Why Cluster Size Matters for PQL
In a binary generalized linear mixed model, observations within a cluster share a latent random effect. When a cluster contains only a few Bernoulli observations, its data may provide little information about that effect. The random-effect distribution then has an important role in determining fitted values, while the cluster’s likelihood may be too broad for a first-order local approximation to represent it well. This is one reason first-order penalized quasi-likelihood (PQL) can be biased with small cluster sizes.
The issue is not simply that small clusters produce noisy estimates. PQL replaces a nonlinear binary response model with a locally linear working model and estimates its random effects under a Gaussian penalty. If a cluster’s data do not sharply locate its random effect, the approximation is being used over a relatively broad range of plausible values. The resulting bias depends on the model and data; it does not have a universal direction for every coefficient.
Model and Information About a Cluster Effect
Consider independent clusters indexed by \(i=1,\ldots,G\), with observations \(j=1,\ldots,m_i\). Let \(Y_{ij}\) be a binary outcome and \(x_{ij}\) a fixed-effect design vector. For this tutorial, use a logistic random-intercept model:
Assume outcomes are conditionally independent within each cluster given \(b_i\), clusters are independent, and \(\tau^2\) is the random-intercept variance. The logistic mean function is \(g(u)=1/(1+\exp(-u))\). PQL uses a first-order Taylor expansion of this mean around a current linear predictor, then fits a penalized working mixed model. Its local nature is central here: the working response and weights describe the model near the expansion point, not the entire nonlinear response curve.
The Gaussian random-effect density contributes a quadratic penalty \(b_i^2/(2\tau^2)\), with curvature \(1/\tau^2\). Thus, near a current value, the curvature of the penalized log-likelihood in \(b_i\) is approximately \(\sum_jp_{ij}(1-p_{ij})+1/\tau^2\). The prior penalty stabilizes the estimate, but it does not create information in that cluster’s observed outcomes. With few observations, the likelihood component can remain weak relative to the penalty.
A Local Curvature Calculation
The information calculation gives a useful, though local, measure of how much each cluster can say about its random intercept. It also shows why cluster size matters even when the model and the average response probability are held fixed.
Worked Example: Comparing Information in Two Cluster Sizes
Suppose the current conditional probability is \(p_{ij}=0.5\) for every observation and the random-intercept variance is \(\tau^2=1\). Compare clusters with \(m_i=2\) and \(m_i=8\).
For each observation, the conditional information is \(p(1-p)=0.5(0.5)=0.25\). Therefore the two-observation cluster has likelihood information \(2(0.25)=0.50\), while the eight-observation cluster has \(8(0.25)=2.00\). The Gaussian penalty has curvature \(1/\tau^2=1\) in both cases.
The approximate total penalized curvatures are \(0.50+1=1.50\) and \(2.00+1=3.00\). Under a local quadratic approximation, the corresponding standard deviations for the random intercept are approximately \(\sqrt{1/1.50}=0.8165\) and \(\sqrt{1/3.00}=0.5774\). The two-observation cluster therefore leaves substantially more local uncertainty about its intercept.
These are local curvature-based approximations, not exact posterior standard deviations. They show the information imbalance: the prior penalty is as strong for both clusters, but the smaller cluster contributes less likelihood curvature to counter it.
Why Weak Information Challenges First-Order PQL
A first-order Taylor approximation is accurate near its expansion point when the relevant range of the predictor is narrow and the mean function is nearly linear over that range. With a small cluster, the likelihood may not concentrate the random intercept near one value. A broader range of plausible intercepts makes it more likely that the local tangent is a poor summary of the nonlinear logistic mean. The working responses and weights can then misrepresent how the cluster’s binary outcomes inform both \(b_i\) and the population parameters.
The Taylor remainder makes this limitation explicit. For a linear predictor \(\eta\) and a change \(b\),
The first-order approximation omits the curvature term and higher-order remainder. The error is not guaranteed to have one sign across predictor values. Its practical importance depends on the range of random effects and fixed predictors that matter in the fitted model, as well as on how much information each cluster supplies.
Worked Example: Checking the Tangent Approximation
At \(\eta=0\), the logistic mean is \(g(0)=0.5\) and its derivative is \(g'(0)=0.25\). The first-order approximation at this point is \(g(0+b)\approx 0.5+0.25b\). Compare it with the exact mean for two possible random-intercept values.
For \(b=1\), the tangent gives \(0.5+0.25(1)=0.75\), while the exact mean is \(g(1)=1/(1+\exp(-1))=0.7311\). The approximation is high by about \(0.0189\).
For \(b=2\), the tangent gives \(0.5+0.25(2)=1.000\), while the exact mean is \(g(2)=1/(1+\exp(-2))=0.8808\). The approximation is high by about \(0.1192\). It even reaches the boundary of the probability scale, although the logistic mean remains below one.
This calculation does not predict the bias of a complete PQL fit. It demonstrates that a local linearization can deteriorate as the predictor moves away from its expansion point. If a small cluster leaves the random intercept weakly determined, the model has less cluster-specific information to keep the relevant range close to that point.
There is a second linked mechanism. The penalized estimate of a cluster effect is pulled toward zero by the Gaussian penalty. With sparse cluster data, that penalty can have substantial influence relative to the cluster’s likelihood. PQL’s working approximation and the penalized estimate are then combined in estimating \(\beta\) and \(\tau^2\). Distortion in the representation of cluster effects can consequently affect population-level estimates as well. The direction and size of that effect depend on the covariates, response pattern, and variance structure; small cluster size alone does not determine them.
How Shrinkage Appears in a Small Cluster
The next calculation isolates the penalty’s contribution to a cluster-specific mode. It illustrates why a very small number of binary observations may not overcome the pull toward the population mean, even for a response pattern that points strongly in one direction.
Worked Example: A Cluster With Two Successes
Suppose a hypothetical cluster has two successes out of two observations, the fixed linear predictor is zero, and \(\tau^2=1\). The conditional log-likelihood for its intercept \(b\), apart from constants, is \(2\log\{g(b)\}\). Including the Gaussian penalty gives \(\ell_{\mathrm{pen}}(b)=2\log\{g(b)\}-b^2/2\).
Differentiating gives \(\ell_{\mathrm{pen}}'(b)=2\{1-g(b)\}-b\). The penalized mode solves \(b=2\{1-g(b)\}=2/(1+\exp(b))\). Substitution gives a root of approximately \(b=0.675\): \(2/(1+\exp(0.675))\approx0.675\). By contrast, without the Gaussian penalty, the log-likelihood increases toward its maximum as \(b\) tends to positive infinity. The penalty makes the estimate finite and pulls it toward zero.
This is a conditional-mode calculation, not a full PQL fixed-effect or variance estimate. It shows concretely that two responses provide limited leverage against the penalty. If the working response and weights are also based on a local approximation, the uncertainty and nonlinearity in this tiny cluster are not fully represented by treating its fitted mode as if it were precisely known.
Practical Assessment and Alternatives
Small cluster sizes are a reason to scrutinize PQL, not an automatic proof that every estimate is materially biased. Assess the fitted model by considering whether the approximation is plausible for the observed design and by comparing it with an approach that better accounts for the nonlinear random-effect contribution.
- Inspect cluster-level information. Note how many Bernoulli observations each cluster contributes and whether the fitted probabilities imply small conditional weights \(p_{ij}(1-p_{ij})\). These weights are largest at \(p=0.5\) and approach zero near zero or one.
- Check the working fit. Review convergence, fitted probabilities, random-effect modes, and whether working responses or updates become extreme. A converged iteration indicates that the algorithm stopped; it does not establish that the first-order approximation is adequate.
- Compare estimates across methods. Likelihood-based integration, including Laplace approximation or adaptive quadrature, incorporates the nonlinear conditional likelihood through an approximation to the random-effect integral. Meaningful changes in fixed-effect or variance estimates relative to PQL indicate sensitivity to its working approximation; they do not by themselves prove which estimate is closest to the truth.
- Separate sampling uncertainty from approximation concerns. Few observations per cluster can make cluster-specific information weak, while approximation error is a distinct issue. A narrow standard error from a working fit should not be taken as evidence that approximation error is negligible.
| Approach | How it handles the binary mean | Relevance when clusters are small |
|---|---|---|
| First-order PQL | Uses a local Taylor linearization and a penalized working mixed model | Computationally convenient, but vulnerable when local linearity is poor or cluster effects are weakly informed |
| Likelihood-based integration | Accounts for the nonlinear conditional likelihood by approximating integration over random effects | Useful as a sensitivity comparison; accuracy and computational cost depend on the integration method and model |
These alternatives have limitations. Integration methods are still numerical approximations and can require more computation. PQL may nevertheless be useful for exploratory or computationally demanding analyses, provided its estimates are not treated as automatically equivalent to likelihood-based estimates. No single cluster-size threshold guarantees negligible PQL bias across models.
Common Mistake
Check Your Understanding
Answer each question, then compare with the brief solutions below.
- For \(m_i=4\) observations with \(p_{ij}=0.5\), what is the conditional information about the random intercept? If \(\tau^2=1\), what is the local penalized curvature?
- Why can a small cluster make a first-order Taylor approximation less dependable, even if the PQL algorithm converges?
- At \(\eta=0\), calculate the first-order logistic approximation for \(b=1\), and compare it with the exact logistic mean to three significant figures.
- Does the Gaussian penalty create additional data information about \(b_i\)? Explain the difference between penalty curvature and likelihood information.
- Name one useful comparison for evaluating sensitivity to first-order PQL, and state one limitation of that comparison.
Solutions. (1) The information is \(4(0.5)(0.5)=1.00\); adding penalty curvature \(1/\tau^2=1\) gives local penalized curvature \(2.00\). (2) Few observations may leave a broad range of plausible random intercepts, over which the nonlinear logistic mean can differ from its tangent; convergence only concerns the iteration’s stopping rule. (3) The tangent gives \(0.5+0.25(1)=0.750\); the exact mean is \(g(1)=0.731\), so the tangent is high by \(0.0189\). (4) No. The penalty contributes curvature \(1/\tau^2\) because of the assumed random-effect distribution, whereas likelihood information comes from observed outcomes and equals \(\sum_jp_{ij}(1-p_{ij})\). (5) Compare with likelihood-based integration such as adaptive quadrature. It can be computationally demanding and remains subject to numerical approximation.