Common Interpretive Errors With Binary PQL

What PQL Output Does—and Does Not—Say

For a binary generalized linear mixed model, a fitted coefficient can be easy to misstate. A positive PQL slope does not, by itself, give the population-averaged change in outcome probability. A random-intercept variance estimate close to zero does not prove that clusters are homogeneous. And an algorithm that reports convergence has only passed a numerical stopping check; that is not a certificate that its estimates are accurate.

These distinctions matter because PQL (penalized quasi-likelihood) estimates a mixed model through an approximation. Its output is often useful, but interpretation must match the model scale, the approximation, and the information in the data. This tutorial concentrates on three errors that can turn a plausible-looking binary PQL summary into an unjustified conclusion.

Model and assumptions: For individuals \(i\) in clusters \(j\), suppose \(Y_{ij}\mid b_j\) is Bernoulli with probability \(p_{ij}\), and \[ \operatorname{logit}(p_{ij})=\beta_0+\beta_1x_{ij}+b_j,\qquad b_j\sim N(0,\sigma_b^2). \] The random intercepts are independent across clusters, and outcomes are conditionally independent given the random effects and covariates. The covariate values are treated as given. PQL approximates estimation for this conditional model; it does not change the distinction between conditional and marginal quantities.

Here \(\beta_1\) describes a one-unit covariate change in the conditional log odds, holding the cluster random effect fixed. Thus \(\exp(\beta_1)\) is a conditional, or cluster-specific, odds ratio under this logit random-intercept model. A population-averaged probability instead averages \(p_{ij}\) over the distribution of \(b_j\). Because the inverse-logit function is nonlinear, that average generally is not the inverse-logit evaluated at the average random effect.

Error 1: Calling a Conditional Coefficient Population-Averaged

The key distinction is between evaluating the response at a specified random effect and averaging the response over random effects. For covariate value \(x\), the model-implied marginal probability is

P(Y=1\mid x)=\int \operatorname{logit}^{-1}(\beta_0+\beta_1x+b)\,f_b(b)\,db,

where \(f_b\) is the normal density with mean zero and variance \(\sigma_b^2\). In contrast, the conditional probability at \(b=0\) is \(\operatorname{logit}^{-1}(\beta_0+\beta_1x)\). These quantities need not agree. Consequently, \(\exp(\beta_1)\) is not generally the population-averaged odds ratio, and \(\beta_1\) is not generally a marginal log-odds slope.

Worked Example: From a Conditional Odds Ratio to Marginal Effects

Suppose a hypothetical binary random-intercept model has \(\hat\beta_0=-0.5185\), \(\hat\beta_1=1.000\), and \(\hat\sigma_b^2=1.000\). The PQL fit reports a conditional odds ratio of \(\exp(1)=2.718\). Approximate the marginal probabilities at \(x=0\) and \(x=1\), then calculate a marginal risk difference and odds ratio.

For a five-point normal quadrature approximation with nodes \(0,\ \pm1.355626,\ \pm2.856970\) and respective weights \(0.5333333,\ 0.2220759\) for each of the first pair, and \(0.0112574\) for each of the second pair, the marginal probability at predictor \(\eta\) is approximated by

\sum_k w_k\operatorname{logit}^{-1}(\eta+b_k).

At \(x=1\), \(\eta=-0.5185+1.000=0.4815\). Substitution gives

0.5333333\operatorname{logit}^{-1}(0.4815) +0.2220759\{\operatorname{logit}^{-1}(1.837126)+\operatorname{logit}^{-1}(-0.874126)\} +0.0112574\{\operatorname{logit}^{-1}(3.338470)+\operatorname{logit}^{-1}(-2.375470)\} \approx 0.5984.

At \(x=0\), the predictor is \(\eta=-0.5185\), so the same calculation gives

0.5333333\operatorname{logit}^{-1}(-0.5185) +0.2220759\{\operatorname{logit}^{-1}(0.837126)+\operatorname{logit}^{-1}(-1.874126)\} +0.0112574\{\operatorname{logit}^{-1}(2.338470)+\operatorname{logit}^{-1}(-3.375470)\} \approx 0.3942.

The marginal risk difference for \(x=1\) versus \(x=0\) is therefore \(0.5984-0.3942=0.2042\), or about 20.42 percentage points. The corresponding marginal odds ratio is \[ \frac{0.5984/(1-0.5984)}{0.3942/(1-0.3942)} =\frac{0.5984/0.4016}{0.3942/0.6058} \approx 2.29. \] Thus the conditional odds ratio is 2.718, while the marginal odds ratio is about 2.29. They answer different questions. The calculations assume the fitted normal random-intercept distribution is an appropriate basis for averaging; they do not remove possible PQL approximation error.

A tempting shortcut is to use the symmetry of a mean-zero random-effect distribution to infer complementary marginal probabilities. Symmetry and \(\operatorname{logit}^{-1}(-z)=1-\operatorname{logit}^{-1}(z)\) imply that complement only when the two fixed predictors are negatives of one another. In the example, \(0.4815\) is not the negative of \(-0.5185\), so symmetry does not justify replacing one marginal probability by one minus the other.

Interpretive rule: State whether an effect is conditional or marginal. For a logit random-intercept model, exponentiating a conditional slope gives a conditional odds ratio. A marginal odds ratio or risk difference requires marginal probabilities, obtained by integrating over the random-effect distribution.

Error 2: Treating a Near-Zero Variance as Proof of No Heterogeneity

The random-intercept variance \(\sigma_b^2\) describes between-cluster variation on the latent log-odds scale under the stated model. A very small fitted value may reflect little detectable variation, but it may also arise from limited information, sparse binary outcomes, too few informative clusters, model misspecification, or the behavior of the PQL approximation. The estimate is constrained to be nonnegative, so a value at or near zero is also a boundary estimate rather than an ordinary estimate well inside the parameter space.

Worked Example: Reading a Small Variance Estimate

Suppose a hypothetical binary PQL analysis of 18 clinics reports \(\hat\sigma_b^2=0.0025\) on the log-odds scale. Convert this to a random-intercept standard deviation and state what can be concluded.

The estimated standard deviation is \[ \hat\sigma_b=\sqrt{0.0025}=0.0500 \] log-odds units. Under the fitted model, this is a small estimated spread in clinic intercepts. It does not establish that the true variance is exactly zero: the data contain only 18 clusters, and binary outcomes may provide limited information about their variation. A PQL variance estimate can also be affected by approximation bias.

The appropriate conclusion is therefore descriptive and qualified: this fit estimates a small between-clinic random-intercept variance, but the estimate alone does not rule out meaningful heterogeneity. To assess how stable that conclusion is, one would examine the data's cluster-level information, relevant model diagnostics, uncertainty or boundary behavior, and results from an appropriate alternative estimation method. A point estimate by itself cannot supply those checks.

A variance component is not the same as a test result. In particular, “the estimate is nearly zero” and “there is no heterogeneity” are not equivalent statements. If the fitted value is exactly zero, that says the optimization or estimation procedure landed at the boundary for this dataset and model. It does not prove the population variance is zero. Inference about a variance component must account for the boundary and the method used; a routine interpretation of a symmetric interval or an unexamined standard error may be unsuitable.

Error 3: Equating Convergence With Accuracy

PQL is iterative. Convergence means that the fitting procedure met its numerical stopping rule—for example, successive parameter updates or an objective change became sufficiently small. It does not establish that the approximation is close to the exact integrated likelihood, that the assumed random-effects distribution is correct, or that the model is well supported by the observed binary data. A fit can converge precisely to a biased or otherwise misleading answer.

Worked Example: A Converged Fit Can Still Be Approximate

Consider hypothetical data for which a PQL fit converges and reports \(\hat\beta_1=0.850\). The covariate is measured in one-unit increments. An analyst says: “The algorithm converged, so the population-averaged odds are multiplied by \(\exp(0.850)\) per unit.” Identify and correct both errors.

First, the coefficient is on the conditional log-odds scale. Its conditional odds ratio is \[ \exp(0.850)=2.340. \] Thus, under the fitted model and for the same random-effect value, a one-unit increase in the covariate multiplies the conditional odds by about 2.34. It does not directly give the population-averaged odds ratio; that requires averaging predicted probabilities over the random-effect distribution and then forming the marginal odds ratio.

Second, convergence only confirms a numerical stopping condition. For illustration, suppose a separate likelihood-based integration fit to the same hypothetical data gives \(\hat\beta_1=0.620\). Its conditional odds ratio is \[ \exp(0.620)=1.859. \] The difference between 2.340 and 1.859 does not alone prove which method is correct, but it demonstrates that numerical convergence of PQL cannot establish agreement with another estimation approach. The discrepancy calls for investigation of the approximation and data support; it is not resolved by the convergence flag.

Diagnostics, Alternatives, and Limits

Interpretation should be accompanied by checks aimed at the particular failure mode, rather than by treating a software status message as a general quality assessment. Useful checks include:

Check the estimand. Label reported coefficients as conditional effects, and state the link and random-effects structure. If the scientific question concerns population-averaged risk or odds, calculate marginal predictions by integrating over the fitted random-effect distribution.
Inspect the variance estimate in context. Record whether it is near the boundary, and examine how many clusters and within-cluster events and non-events inform it. Sparse outcomes or a small number of clusters make a precise heterogeneity conclusion especially difficult.
Inspect the fitting output and data support. Review warnings, extreme fitted probabilities, very large coefficients, sensitivity to starting values, and whether some clusters contain only zeros or only ones. Such features can destabilize estimation or make an approximation less reliable even when a stopping rule is met.
Compare suitable estimation approaches. Likelihood-based integration, including numerical quadrature when appropriate, targets the integrated likelihood rather than PQL's approximation. Generalized estimating equations target marginal mean relationships under their own assumptions. These approaches do not answer identical questions automatically, so compare like estimands and compatible models.

The comparison is diagnostic, not a contest in which one method is always best. Likelihood-based integration can be computationally demanding or unstable in difficult data configurations; PQL can be attractive computationally but may incur appreciable bias for binary outcomes. Generalized estimating equations address marginal means rather than estimating the same random-effects distribution. Model assumptions, data structure, and the desired estimand determine which results are relevant.

A further limitation is that integrating PQL estimates yields model-based marginal predictions, not guaranteed accurate population effects. If the fixed effects or variance estimate are biased, or the normal random-intercept model is unsuitable, the marginal quantities inherit those limitations. Reporting should keep the distinction clear: an estimate is conditional or marginal, model-based or directly targeted, and approximate or obtained from a different fitting criterion.

A Short Justification for the Conditional–Marginal Distinction

The distinction follows directly from nonlinearity. Let \(g(z)=\operatorname{logit}^{-1}(z)\). The conditional probability is \(g(\beta_0+\beta_1x+b)\), whereas the marginal probability is its expectation over \(b\). In general, averaging a nonlinear function is not the same as applying that function to an average:

E_b\{g(\beta_0+\beta_1x+b)\} \ne g\{\beta_0+\beta_1x+E_b(b)\} = g(\beta_0+\beta_1x).

The inequality is not universal for every parameter choice, but equality should not be assumed. Even when the marginal probabilities are calculated at two covariate values, their odds ratio need not equal \(\exp(\beta_1)\). This is why the conditional coefficient cannot simply be relabelled as a population-averaged effect.

Key takeaway: In binary PQL, interpret a fixed slope as a conditional effect under the fitted GLMM. Treat a near-zero variance as an estimate—not proof—of absent heterogeneity, and treat convergence as a numerical property—not proof—of statistical accuracy.

Check Your Understanding

Answer each question, then compare your reasoning with the solutions.

A fitted logit random-intercept model gives \(\hat\beta_1=0.40\). What does \(\exp(0.40)\) describe, and what does it not automatically describe?
A PQL fit reports \(\hat\sigma_b^2=0\). Give one conclusion that is justified and one that is not.
Why does a converged PQL algorithm not establish that its fixed-effect estimate is accurate?
If marginal probabilities at two covariate values are 0.30 and 0.45, calculate the risk difference. Is that calculation the conditional odds ratio?

Solutions. (1) \(\exp(0.40)=1.492\) is the conditional odds ratio for a one-unit increase, holding the random effect fixed. It is not automatically the population-averaged odds ratio. (2) The estimate is at the zero boundary for this fitted model and dataset; it does not prove that the population has no heterogeneity. (3) Convergence means the numerical stopping rule was met, not that the PQL approximation, model assumptions, or data information guarantee accuracy. (4) The risk difference is \(0.45-0.30=0.15\), or 15 percentage points. It is a marginal probability contrast, not a conditional odds ratio.