Why Binary PQL Shrinks Random-Effect Variance Estimates

Why the Variance Can Be Too Small

Penalized quasi-likelihood (PQL) turns a non-Gaussian mixed model into a sequence of weighted linear mixed-model calculations. For binary responses, that sequence uses a local approximation to the mean and treats the resulting pseudo-data as Gaussian. This is computationally useful, but it can make evidence for between-cluster variation look weaker than it is. The effect is particularly important when clusters contain few responses or their fitted probabilities give little working information.

The claim is a tendency, not a rule that every PQL fit must underestimate the variance. It is also important to distinguish two kinds of shrinkage. A cluster’s estimated random effect is pulled toward zero by its Gaussian prior in any mixed model. The concern here is different: the estimated population variance of the random intercept, often written \(\tau^2\), may itself be too small because PQL’s approximate working model does not represent the binary likelihood’s information and curvature adequately.

Where the Attenuation Enters

At a current value of the linear predictor \(\eta_{ij}=x_{ij}^{T}\beta+b_i\), write \(\mu_{ij}=\operatorname{logit}^{-1}(\eta_{ij})\). Binary PQL forms the working response and weight

The weighted working model treats \(z_{ij}\) as a noisy observation of the linear predictor, with working variance \(1/w_{ij}\). Thus a small weight means large working noise. The weight is at most \(1/4\), attained at \(\mu_{ij}=1/2\), and it gets smaller as the fitted probability approaches either 0 or 1.

The same weight is the local curvature of the Bernoulli log likelihood with respect to the linear predictor. For one response, the log likelihood, omitting a constant, is \(\ell(\eta)=Y\eta-\log(1+\exp(\eta))\). Its first derivative is \(Y-\mu\), and its second derivative is \(-\mu(1-\mu)=-w\). Consequently, around a current predictor, the data provide approximately \(w\) units of curvature about that predictor. For a cluster, the local curvature contributed by its responses is \(W_i=\sum_j w_{ij}\).

To obtain the factor, expand at \(b_i=0\), so evaluate \(\mu_{ij}\) and \(w_{ij}\) at \(\eta_{ij}=x_{ij}^{T}\beta\). Let \(r_i=\sum_j(Y_{ij}-\mu_{ij})\) be the cluster’s score for its intercept at this expansion point. The quadratic approximation to the log likelihood has the form \(r_i b_i-\tfrac12W_i b_i^2\), apart from terms not involving \(b_i\). Adding the Gaussian log-prior contribution \(-b_i^2/(2\tau^2)\) gives a quadratic whose maximum solves

When \(W_i\) is small, the cluster contributes little curvature relative to the prior precision \(1/\tau^2\), so its estimated intercept is strongly pulled toward zero. Sparse clusters often have small \(W_i\) simply because they contain few responses; probabilities close to zero or one can make each response’s weight small as well. This reduces the separation between cluster-to-cluster signal and working noise in the Gaussian calculations used by PQL.

Why Shrinkage Can Affect the Variance Estimate

PQL does not estimate \(\tau^2\) by simply averaging the squared conditional modes. It estimates variance components through a working Gaussian mixed-model criterion, accounting for the working residual variance and the random-effects distribution. Nevertheless, the working-data representation helps explain the tendency. If the cluster-specific pseudo-data are noisy relative to their between-cluster spread, the Gaussian criterion has limited evidence that the spread is due to random intercepts rather than working error. The variance estimate can then be pulled toward the boundary value zero.

A simple calculation makes this distinction visible. Suppose, just for illustration, that each cluster is summarized by a working mean \(Z_i\) satisfying the approximate model \(Z_i=b_i+e_i\), with independent \(e_i\) of variance \(v\). Then \(\operatorname{Var}(Z_i)=\tau^2+v\). A moment estimate would subtract the known working-noise variance \(v\) from the observed variance of the \(Z_i\), truncating at zero because a variance cannot be negative. This is a teaching approximation, not the exact PQL variance-component estimator.

There is a second source of concern: linearization is local, whereas random effects can move the conditional probabilities well away from the expansion point. For example, at \(\eta=0\), the logistic mean is \(0.5\) and its derivative is \(0.25\). The tangent approximation is \(0.5+0.25b\). At \(b=-2\) and \(b=2\), it gives \(0\) and \(1\), while the exact logistic probabilities are approximately \(0.1192\) and \(0.8808\). Beyond those values the tangent can even leave the probability range. The linearized Gaussian calculation therefore need not preserve how the binary likelihood accumulates information across the random-effects distribution. This mismatch can flatten or distort the evidence for a nonzero \(\tau^2\).

What to Check in Practice

A small PQL variance estimate should be treated as an approximation result to investigate, not automatically as proof that clusters are homogeneous. Useful checks include the following:

• Inspect the boundary and convergence information. A variance estimate at or near zero may indicate weak information in the working model. It is not by itself proof that the true variance is zero.
• Examine cluster sizes and fitted weights. Few observations per cluster, or many fitted probabilities near zero or one, produce small \(W_i\) and noisier working responses.
• Compare with a likelihood method that integrates over random effects. Laplace approximation and adaptive quadrature evaluate an approximation to the integrated likelihood rather than relying on Gaussian pseudo-responses as the final data representation. Agreement is reassuring; disagreement signals sensitivity to approximation choice. Neither method is automatically reliable if its numerical diagnostics are poor.
• Check sensitivity to the working fit. Changes in starting values or iterations that materially change the variance estimate merit attention, as do estimates that remain at a boundary despite apparent cluster-level pattern.

This comparison does not mean that PQL is always unusable or that integration-based estimates are always correct. PQL can be convenient and can perform adequately when its approximation is suitable. Its key limitation here is structural: a local linearization followed by a Gaussian variance calculation may not retain the binary model’s full information about heterogeneity. With limited within-cluster data, there may also be little information to estimate \(\tau^2\) precisely even under a more accurate likelihood method.

Common Mistake

A second overstatement is to claim that binary PQL always underestimates \(\tau^2\). The direction and size of approximation error depend on the design, the fitted probabilities, the amount of information, and the fitting details. “Can shrink” or “often underestimates in sparse settings” is more defensible than a universal claim.