Why PQL Warnings Need Careful Interpretation
Penalized quasi-likelihood (PQL) fits a generalized linear mixed model by repeatedly approximating its non-Gaussian response model with a weighted linear mixed model. For binary outcomes, that approximation can be sensitive to sparse information and to fitted probabilities near zero or one. A fitting routine may report convergence even when its estimates deserve scrutiny: convergence describes the behavior of an algorithm, not necessarily the accuracy of the approximation.
The earlier material on binary PQL bias establishes that PQL can underestimate random-effect variation and can distort fixed effects in some designs. Here the task is diagnostic: identify features of a fitted model that should prompt investigation, explain what they do and do not imply, and use more than one indicator before drawing a conclusion.
Model and Working Quantities
Consider a binary response \(Y_{ij}\) for observation \(i\) in cluster \(j\), with \(Y_{ij}=1\) for an event and \(Y_{ij}=0\) otherwise. A logistic random-intercept model is
The standard working assumptions are that responses are conditionally independent given the random effects and covariates, clusters are independent, and the random intercepts follow the stated normal distribution. PQL additionally relies on iterative local linearization of the mean-variance relationship. Its working response and weights are calculated at the current fitted linear predictors; they are not new observed outcomes or ordinary observation-level precision estimates.
Why Extreme Working Quantities Arise
These formulas follow from a first-order Taylor approximation. Write \(\eta_0\) for the current linear predictor and \(\mu_0\) for its fitted mean. Locally,
For the logistic link, \(d\mu/d\eta=\mu(1-\mu)\). Substituting the observed response \(Y\) for the locally approximated mean and solving for the linear predictor gives
This is the working response. The corresponding weighted least-squares weight is the squared derivative divided by the Bernoulli variance: \(\{\mu(1-\mu)\}^2/\{\mu(1-\mu)\}=\mu(1-\mu)\). Thus, when a fitted probability is very close to 0 or 1, the denominator in the working-response correction is small. A response that disagrees with that fitted probability can produce an exceptionally large correction, while receiving a small working weight. Such a pair of values is a warning to inspect the fit and its data support, not proof by itself that the entire model is unusable.
Worked Example: Calculating an Extreme Working Response
Suppose a hypothetical PQL iteration assigns an observation \(Y=1\) a fitted probability \(\mu=0.02\). Calculate its linear predictor, working response, and weight.
The linear predictor is \(\eta=\log\{0.02/(1-0.02)\}=\log(0.02041)\approx -3.892\). The Bernoulli variance factor is \(0.02(0.98)=0.0196\). Therefore,
The working response, \(46.11\), is far from the observed binary value because the current model assigned a very small probability to an event that occurred. Its weight is also small. The combination can make the local linear approximation consequential for the iteration even though this observation is downweighted. The value \(46.11\) is not a predicted probability; probabilities must lie between 0 and 1.
Small Random-Effect Variance Estimates
A PQL estimate of \(\tau^2\) very close to zero is a potential warning because binary PQL can shrink estimated random-effect variation downward. But the same estimate can also arise when the underlying variation is genuinely small or when the data contain little information about it. The estimate alone cannot distinguish these explanations. In particular, do not interpret a small estimate as proof that clusters are homogeneous or that a random intercept is unnecessary.
Check the estimate alongside the number of clusters, cluster sizes, event counts within clusters, and the distribution of fitted probabilities. Also inspect whether the estimate is at the boundary, whether reported standard errors or convergence diagnostics are reliable near that boundary, and whether small changes to the fitting setup alter the result. A variance estimate is on a squared scale; taking its square root gives the random-effect standard deviation on the linear-predictor scale.
Worked Example: Interpreting a Very Small Variance Estimate
In a hypothetical data set, 24 clinics contribute four binary observations each. There are 14 events among the 96 observations. A logistic random-intercept PQL fit reports \(\hat{\tau}^2=0.014\). Interpret the scale and identify what can be concluded.
The estimated random-intercept standard deviation is \(\sqrt{0.014}\approx 0.118\) on the log-odds scale. A one-standard-deviation difference in random intercept corresponds to an odds multiplier of \(\exp(0.118)\approx 1.125\), conditional on the fixed effects. This is a modest estimated amount of between-clinic variation in this fitted model.
It is not justified to conclude that the true variance is zero: the estimate is positive, and PQL may understate variation. Nor does the total of 14 events alone show whether the information is adequate; their allocation across clinics matters. The small cluster size and the pattern of within-clinic events are reasons to inspect fitted probabilities, estimation stability, and an appropriate alternative fit before making a substantive claim about clinic variation.
Sensitivity to Starting Values
PQL is iterative. An implementation begins from initial fixed-effect and variance-component values, updates working quantities, and repeats until its stopping rule is met. Strong dependence of the final estimates on reasonable starting values is a warning about numerical or approximation stability. It can arise from weak information, boundary behavior, extreme working quantities, or a difficult iteration path. It does not, by itself, show which fit is correct.
A useful check holds the data, model specification, software settings, and convergence tolerances constant while varying plausible starting values. Compare not only whether the routine declares convergence, but also the final variance estimate, fixed-effect estimates, fitted probabilities, and any warnings. If estimates differ materially, verify coding and scaling, inspect the data and working quantities, and seek an independent fit or expert review. “Converged” means only that the stated numerical stopping rule was met.
Worked Example: Comparing Fits From Different Starts
For the same hypothetical binary random-intercept model, three reasonable starting values for \(\tau^2\) produce the following PQL results. Assume the software reports convergence for each run under the same tolerance.
| Starting \(\tau^2\) | Final \(\hat{\tau}^2\) | Final fixed-effect estimate \(\hat{\beta}_1\) |
|---|---|---|
| 0.01 | 0.006 | −0.42 |
| 0.25 | 0.071 | −0.57 |
| 1.00 | 0.240 | −0.69 |
The largest and smallest final variance estimates differ by \(0.240-0.006=0.234\); their ratio is \(0.240/0.006=40\). The fixed-effect estimates also span \(−0.69-(−0.42)=−0.27\). Those differences are substantial relative to the estimates themselves, so the runs do not support a stable single result. Convergence flags do not resolve that disagreement. The appropriate conclusion is that these fitted results are starting-value sensitive and need investigation, not that the largest or smallest estimate should automatically be selected.
Putting the Warning Signs Together
No universal cutoff defines an “extreme” working response or a “too small” variance estimate. Interpret each value relative to the model, the scale, the data, and the other diagnostics. A practical review can proceed in this order:
Confirm the binary outcome coding, covariates, cluster identifiers, and stated random-effects structure. Check cluster sizes and event counts; sparse clusters can contribute weak information.
Find probabilities near 0 or 1, then examine the corresponding working responses and weights. Look for observations whose outcomes strongly disagree with the fitted probabilities.
Compare reasonable starting values and check convergence messages, variance estimates, fixed effects, and fitted probabilities. A single convergence label is not a sufficient diagnostic.
Describe the output as a PQL estimate and state relevant warnings. If results are consequential or unstable, compare with a suitable alternative fitting approach rather than treating one PQL run as definitive.
This diagnostic use of alternatives is deliberately limited: likelihood-based integration methods, including Laplace approximation and adaptive Gaussian quadrature, can provide a comparison, but they have their own computational and modeling considerations. A standard GLM that omits the random effects answers a different modeling question, so it is not a direct validation of the PQL mixed-model fit. Nor does agreement between methods prove that all assumptions are correct.
Common Pitfall and Limitations
A common mistake is to argue, “The algorithm converged, so the small variance estimate is trustworthy.” The correction is that numerical convergence and approximation quality are separate questions. A second mistake is to treat every large working response as an erroneous observation. The correction is to trace it to its fitted probability and outcome: the quantity may be a predictable consequence of logistic linearization, while still flagging a region where the approximation deserves attention.
These indicators are not formal tests and do not supply a probability that PQL is correct or incorrect. Small variance estimates can be valid, extreme working quantities can occur in otherwise useful fits, and starting-value sensitivity can reflect numerical issues rather than approximation bias alone. Conversely, the absence of conspicuous warnings does not guarantee negligible PQL bias. The diagnostics guide scrutiny; they do not replace design knowledge, model checking, or appropriately cautious interpretation.
Check Your Understanding
Answer each question, then compare your reasoning with the solution immediately below it.
- A binary observation has \(Y=0\) and current fitted probability \(\mu=0.98\). Calculate its working response and weight to three significant figures.
Solution: \(\eta=\log(0.98/0.02)=\log(49)\approx3.892\). Thus \(Y^*=3.892+(0-0.98)/\{0.98(0.02)\}=3.892-50=-46.1\), and \(w=0.98(0.02)=0.0196\). The large negative working response accompanies a small weight. - Why does a PQL variance estimate of \(0.0025\) not establish that the true random-intercept variance is zero?
Solution: The estimate is positive, and even a near-zero PQL estimate may reflect downward approximation bias or weak information as well as genuinely small variation. It is a fitted estimate, not proof of a zero population variance. - Three fits meet their convergence criteria but produce \(\hat{\tau}^2=0.01,\ 0.02,\ 0.19\). What is the diagnostic conclusion?
Solution: The estimates differ materially: the largest is \(0.19/0.01=19\) times the smallest. This is evidence of starting-value sensitivity under the stated fitting procedure and warrants investigation. It does not identify which estimate is correct. - Does an extreme working response of \(46.1\) mean the fitted event probability is 46.1?
Solution: No. A working response is on the linear-predictor scale and can fall outside the probability range. In the worked calculation, the fitted probability was \(0.02\); \(46.1\) was the linearized working response, not a probability. - Name two data features to inspect when a binary PQL fit reports a very small random-effect variance.
Solution: Inspect the number of clusters and the event counts within clusters; also examine cluster sizes and fitted probabilities. These features help assess how much information the data provide about between-cluster variation and whether extreme working quantities are present.