Tutorials › Generalized Linear Mixed Models (GLMM) › Worked Comparison Of PQL Estimates Across Cluster Designs

Penalized Quasi-Likelihood and Its Limitations · Tutorial 239 of 698

Worked Comparison Of PQL Estimates Across Cluster Designs

Learn to compare PQL estimates across cluster designs without mistaking changes in fixed effects or random-intercept variance for effects caused by imbalance alone.

Advanced 12 min read

What You'll Learn

  • Compare balanced and highly imbalanced designs while holding the total number of observations fixed.
  • Calculate a working approximation to the information contributed by clusters of different sizes.
  • Interpret changes in fixed-effect estimates on the log-odds and odds-ratio scales.
  • Assess how event concentration across clusters can affect the estimated random-intercept variance.
  • Distinguish design patterns that suggest weak information from evidence that a PQL fit is unreliable.
  • Compare PQL with likelihood-based integration and with a model that ignores clustering.

Why Compare Cluster Designs?

Two binary random-intercept models can use the same total number of observations and still provide quite different information about cluster effects. If observations are concentrated in a few large clusters, those clusters can estimate their own outcomes well, while many small clusters may provide little information about their individual random intercepts. The balance between clusters also affects how covariate information is distributed. Penalized quasi-likelihood (PQL) uses a local working approximation, so changes in that information pattern can be reflected in its estimates as well as their uncertainty.

The purpose of a design comparison is not to infer that imbalance must move an estimate in one particular direction. Rather, compare the actual cluster sizes, covariate patterns, and event counts, then ask whether the resulting PQL estimates are consistent with those patterns and stable enough to interpret. The earlier material on Unequal Cluster Sizes And Binary PQL and Separating Approximation Bias From Sampling Variability establishes why these issues matter; here the emphasis is on worked comparisons.

Model and assumptions: Let \(Y_{ij}\) be a binary outcome for observation \(i\) in cluster \(j\), with \(Y_{ij}=1\) for an event and \(Y_{ij}=0\) otherwise. A random-intercept logistic GLMM has $$ Y_{ij}\mid b_j \sim \operatorname{Bernoulli}(\mu_{ij}), \qquad \operatorname{logit}(\mu_{ij})=\mathbf{x}_{ij}^{\mathsf T}\boldsymbol{\beta}+b_j, \qquad b_j\sim N(0,\sigma_b^2). $$ Conditional on the random effects and covariates, observations are assumed independent, the conditional mean follows the stated logit model, and the cluster effects are independent draws from the stated normal distribution. PQL estimates the model using a local linearization and working responses and weights; it does not maximize the exact integrated likelihood. Comparisons are most informative when the model specification, outcome definition, covariate coding, and fitting conventions are held fixed.

A fixed-effect coefficient \(\beta_k\) is a conditional log-odds contrast for a one-unit increase in its covariate, holding the other covariates and the random intercept fixed. Its exponentiation, \(\exp(\beta_k)\), is the corresponding conditional odds ratio. The variance \(\sigma_b^2\) describes the spread of cluster intercepts on the log-odds scale; it is not a percentage or a variance in the binary outcome itself.

How Cluster Size Enters the Working Information

For a Bernoulli observation with fitted mean \(\mu_{ij}\), the usual PQL working weight is \(w_{ij}=\mu_{ij}(1-\mu_{ij})\). It is largest at \(\mu_{ij}=0.5\), where it equals 0.25, and becomes smaller as fitted probabilities approach zero or one. For a simple information comparison, suppose a cluster has a common working weight \(w\). Its approximate working information about its random intercept is \(n_jw\), where \(n_j\) is its size.

A useful way to understand the effect of that information is to combine it with the random-intercept distribution. In the working Gaussian approximation, the cluster-level estimate has working sampling variance approximately \(1/(n_jw)\), while the random intercept has variance \(\sigma_b^2\). Their sum gives a marginal variance for the cluster-level signal, and the reciprocal is an approximate marginal information measure.

$$ I_j^{*}\approx\frac{1}{\sigma_b^2+1/(n_jw)}. $$

This is a design heuristic at a specified working iteration, not the exact information for all fixed effects, nor a replacement for fitting the model. It shows that the marginal precision of the noisy cluster-level signal does not grow in direct proportion to cluster size: as \(n_jw\) grows, the term \(1/(n_jw)\) shrinks, but the random-intercept variance remains. This is not information about the cluster's own latent intercept; in the working Gaussian model, sampling precision for that intercept is \(n_jw\), and its posterior precision after combining with the \(N(0,\sigma_b^2)\) distribution is \(1/\sigma_b^2+n_jw\), both of which increase without bound as \(n_jw\) grows. Fixed-effect information also depends on the covariate values and their distribution across clusters.

Worked Example: Equal Total Sample Size, Different Cluster Sizes

Consider two hypothetical designs with 20 clusters and 200 observations each in total. The balanced design has 20 clusters of size 10. The imbalanced design has 10 clusters of size 2 and 10 of size 18. At one representative PQL iteration, suppose the common working weight is \(w=0.16\) and the fitted random-intercept variance is \(\sigma_b^2=0.64\). Compare the approximate cluster-level information \(I_j^*\).

For a size-2 cluster, \(n_jw=2(0.16)=0.32\), so $$ I_j^{*}=\frac{1}{0.64+1/0.32} =\frac{1}{0.64+3.125} =\frac{1}{3.765} \approx 0.2656. $$

For a size-18 cluster, \(n_jw=18(0.16)=2.88\), and $$ I_j^{*}=\frac{1}{0.64+1/2.88} =\frac{1}{0.64+0.3472} =\frac{1}{0.9872} \approx 1.013. $$

For comparison, a size-10 cluster has \(n_jw=1.60\), giving $$ I_j^{*}=\frac{1}{0.64+1/1.60} =\frac{1}{1.265} \approx 0.7905. $$

Summing this heuristic over clusters gives about \(10(0.2656+1.013)=12.79\) for the imbalanced design and \(20(0.7905)=15.81\) for the balanced design. Thus, under these particular working values, the balanced design has greater summed approximate information by this measure, despite having the same total number of observations. This is not a universal ranking: the result depends on the working weights, the variance, and the information quantity being considered. It does illustrate that equal total sample size does not mean equal information about cluster-level variation.

Comparing Fixed Effects and Their Uncertainty

A fixed-effect estimate can change across designs for several reasons. The covariate may be distributed differently across clusters; a small number of large clusters may account for much of the covariate contrast; or the working weights may differ because fitted probabilities differ. When covariates vary both within and between clusters, those patterns can matter even if the model formula is unchanged. Consequently, imbalance alone does not predict whether a coefficient increases or decreases.

Worked Example: Interpreting a Fixed-Effect Difference

Suppose two hypothetical PQL analyses use the same binary outcome, model formula, and total of 200 observations. One has balanced clusters and the other highly imbalanced clusters. Their stipulated fitting summaries for a treatment indicator are:

DesignTreatment coefficientStandard error
Balanced0.600.20
Highly imbalanced0.440.24

The balanced-design conditional odds ratio is \(\exp(0.60)\approx1.822\). The imbalanced-design conditional odds ratio is \(\exp(0.44)\approx1.553\). The coefficient difference is \(0.44-0.60=-0.16\) log-odds units; the ratio of the two odds-ratio estimates is \(\exp(-0.16)\approx0.8521\). Thus the imbalanced-design estimate is about 14.8% lower on the odds-ratio scale, since \(1-0.8521=0.1479\).

The reported standard error increases from 0.20 to 0.24, a 20% increase because \((0.24-0.20)/0.20=0.20\). These calculations describe the stipulated summaries; they do not establish that imbalance caused the coefficient change. To assess that claim, inspect cluster-specific treatment allocation and event patterns, check that the model and coding are identical, and consider whether influential large clusters dominate the contrast. The coefficient is conditional on the random intercept, not a population-averaged odds ratio.

The standard error is also not a diagnostic of approximation bias. A small standard error can accompany a biased PQL estimate, and a large standard error can reflect limited information rather than a faulty algorithm. Comparing PQL with likelihood-based integration can help assess sensitivity to the approximation, while comparing with a GLM that omits the random intercept addresses a different question: that GLM does not account for within-cluster dependence and does not estimate \(\sigma_b^2\).

Event Patterns and the Random-Intercept Variance

Cluster sizes alone do not summarize the information in a Bernoulli outcome. With the same size and the same total number of events, events may be spread evenly across clusters or concentrated in a few. Such patterns can suggest different degrees of between-cluster heterogeneity. They do not determine the fitted variance by themselves: fitted covariates, the assumed distribution, PQL's working approximation, and other data features all contribute.

Worked Example: Events Spread Across Clusters or Concentrated

Suppose two hypothetical data sets each contain eight clusters of size 10 and eight events overall. In the first data set, each cluster has one event. In the second, four clusters have two events each and four have none. Both have \(80\) observations and an overall event proportion of \(8/80=0.10\), but the second distributes events less evenly across clusters.

For illustration, suppose comparable random-intercept PQL fits report \(\hat{\sigma}_b^2=0.18\) when events are spread evenly and \(\hat{\sigma}_b^2=0.70\) when events are concentrated. The difference is \(0.70-0.18=0.52\) log-odds-squared units, and the variance estimate in the concentrated-events fit is \(0.70/0.18\approx3.889\) times as large. The corresponding random-intercept standard deviations are \(\sqrt{0.18}\approx0.4243\) and \(\sqrt{0.70}\approx0.8367\) log-odds units.

These are hypothetical fitting summaries used to illustrate interpretation, not estimates calculated from the event counts alone. The event pattern motivates checking the variance estimate and cluster-level fitted values; it does not prove that the larger estimate is correct. In particular, a concentrated pattern can coexist with sparse information, and PQL may be sensitive in such settings. Review convergence warnings, fitted probabilities near zero or one, and whether the estimated variance changes materially under a likelihood-based fit or a reasonable change in specification.

Derivation: Why the Marginal Signal Precision Levels Off

The approximation above follows directly from combining the working sampling variance and random-effect variance. In the local working model, an estimate of a cluster's intercept based on its observations has variance approximately \(1/(n_jw)\). If the true cluster intercept varies with variance \(\sigma_b^2\), the marginal variance of that noisy cluster-level signal is the sum of these independent variance components:

$$ \operatorname{Var}(\text{cluster-level signal}) \approx \sigma_b^2+\frac{1}{n_jw}. $$

The reciprocal variance is the corresponding marginal Gaussian precision of the noisy cluster-level signal, which gives \(I_j^*\). As \(n_j\) increases without bound while \(w\) stays positive, \(1/(n_jw)\) tends to zero, and \(I_j^*\) tends to \(1/\sigma_b^2\), rather than increasing without limit. This leveling-off is for the marginal precision of the noisy cluster-level signal, not information about one cluster's latent intercept: its sampling precision \(n_jw\) and posterior precision \(1/\sigma_b^2+n_jw\) continue to increase without bound. Additional observations in a large cluster can also contribute to fixed-effect estimation, especially when covariates vary within that cluster.

Diagnostics, Alternatives, and Limitations

A comparison is most useful when it is treated as a structured sensitivity analysis rather than a contest to select the most favorable estimate. Examine the cluster-size distribution, event counts by cluster, and covariate distributions both within and between clusters. Then compare coefficient estimates, standard errors, and the random-intercept variance on their appropriate scales. A change should be described in context and related to the information pattern, not attributed to imbalance without supporting evidence.

  • Check whether large clusters dominate. Compare results with and without especially large or influential clusters as a sensitivity analysis, while being explicit that removing clusters changes the analyzed data.
  • Check sparse and extreme patterns. Clusters with no events or all events can yield weak information about their fitted probabilities; very small working weights can also reduce local information.
  • Check the fit itself. Nonconvergence, extreme fitted values, or large changes under small analytic choices weaken confidence in the approximation. A plausible-looking estimate alone does not establish reliable PQL behavior.
  • Separate distinct uncertainties. Sampling variability concerns which observations were obtained; approximation bias concerns the PQL method's departure from the integrated likelihood. Neither can be read off from the standard error alone.

Likelihood-based integration is a principal alternative because it evaluates the contribution of the random effects through integration rather than relying on PQL's local linearization. Agreement is reassuring but does not guarantee either method is adequate; disagreement signals a need to investigate approximation sensitivity, data sparsity, and model specification. A GLM without random effects may be useful for a different estimand or as a rough contrast, but it does not preserve the assumed cluster structure. PQL comparisons also have limits: changing cluster design changes the data information, and a working-weight calculation cannot reproduce the full fitted model or separate every source of estimation error.

Common mistake: “The imbalanced design has fewer small-cluster events, so its random-intercept variance must be larger, and any fixed-effect shift must be caused by imbalance.” This reasoning overstates what the design pattern establishes. Event concentration can be informative, but the fitted variance also depends on covariates, the random-effects assumption, and the PQL approximation; fixed-effect changes can reflect differing covariate support or influential clusters. The correction is to report the observed pattern, calculate the estimate changes, and treat them as a sensitivity finding unless a controlled comparison supports a causal explanation.

Check Your Understanding

Answer each question using the stated working approximation and distinguish calculations from conclusions that require a full model fit.

  1. For \(n_j=5\), \(w=0.20\), and \(\sigma_b^2=0.50\), calculate \(I_j^*\).
  2. Two fitted treatment coefficients are 0.30 and 0.10. Calculate their conditional odds ratios and the ratio of the second odds ratio to the first.
  3. Why can equal total sample sizes in balanced and imbalanced designs yield different information about cluster effects?
  4. Two data sets have the same overall event proportion, but events are more concentrated in clusters in one. What can and cannot be concluded about the fitted random-intercept variance?

Solutions. (1) \(n_jw=5(0.20)=1\), so \(I_j^*=1/(0.50+1/1)=1/1.50\approx0.6667\). (2) The odds ratios are \(\exp(0.30)\approx1.350\) and \(\exp(0.10)\approx1.105\); their ratio is \(\exp(0.10-0.30)=\exp(-0.20)\approx0.8187\). The second estimate is about 18.1% lower on this odds-ratio comparison. (3) Cluster sizes determine how working information is distributed, and covariate and outcome patterns can differ across clusters even when the total number of observations matches. (4) Concentration is a reason to inspect the variance estimate and PQL sensitivity; it does not determine that estimate or prove that it is larger or more accurate without fitting the model and considering its assumptions.

Clinical Biostats