Why Add a Second-Order Correction?
Penalized quasi-likelihood (PQL) replaces a nonlinear mixed model with a sequence of approximate weighted linear mixed-model calculations. In binary models, the response mean changes nonlinearly with the linear predictor. A first-order local approximation can therefore miss how variation in the random effects bends the mean curve. Second-order PQL corrections are intended to account for some of that curvature and reduce a particular approximation error.
The word some matters. A second-order correction is an additional local approximation, not an exact integration over random effects and not a general guarantee of low bias. Its performance still depends on the amount and pattern of information in the binary data, the random-effect variation, and how well the local expansion represents the model over the relevant range.
Model and Curvature Correction
For binary observations \(Y_{ij}\) in cluster \(i\), consider a logistic random-effects model. Conditional on the random effects, the observations are Bernoulli and independent, with \(\eta_{ij}=x_{ij}^{\mathsf T}\beta+z_{ij}^{\mathsf T}b_i\) and \(P(Y_{ij}=1\mid b_i)=h(\eta_{ij})\), where \(h(u)=1/(1+\exp(-u))\). Assume \(b_i\sim N(0,G)\), independently across clusters, and that the design vectors and model structure are specified correctly. The conditional independence assumption is given the cluster random effects; marginally, observations in a cluster can be dependent.
The first-order term in the Taylor expansion is \(h'(\eta)U\), whose expectation is zero because \(E(U)=0\). The second-order term is \(\frac{1}{2}h''(\eta)U^2\), whose expectation is \(\frac{1}{2}h''(\eta)q\). For the logistic function, \[ h'(\eta)=h(\eta)\{1-h(\eta)\},\qquad h''(\eta)=h(\eta)\{1-h(\eta)\}\{1-2h(\eta)\}. \] Thus the adjustment is negative when \(h(\eta)>1/2\), positive when \(h(\eta)<1/2\), and zero at \(h(\eta)=1/2\). It reflects the local concavity or convexity of the logistic curve, not a universal shift in one direction.
Expand the inverse-link mean around the reference predictor: \(h(\eta+U)\approx h(\eta)+h'(\eta)U+\frac{1}{2}h''(\eta)U^2\).
Use \(E(U)=0\) and \(E(U^2)=q\), so the linear term disappears and the quadratic term contributes \(\frac{1}{2}h''(\eta)q\).
The expansion omits higher-order terms. It approximates a feature of the nonlinear mean; it does not evaluate the exact random-effects integral or recover the entire joint distribution.
The justification is Taylor's theorem: for some value \(\xi\) between \(\eta\) and \(\eta+U\), the remainder after the quadratic terms is \(h'''(\xi)U^3/6\). Taking expectations does not generally make that remainder vanish, because \(\xi\) depends on \(U\). Even when a centered normal \(U\) has \(E(U^3)=0\), replacing the varying \(h'''(\xi)\) by a constant is another approximation; fourth- and higher-order contributions can remain important. The second-order term is therefore useful when variation around \(\eta\) is modest and the curve is locally well represented by its quadratic expansion.
Worked Example: Computing a Logistic Curvature Adjustment
Suppose a hypothetical binary-model calculation has reference predictor \(\eta=1.00\) and random-effect contribution \(U\) with variance \(q=0.800\). Calculate the second-order adjusted mean.
First, \(h(1)=1/(1+\exp(-1))=0.7311\). The logistic second derivative is \(h''(1)=0.7311(1-0.7311)(1-2(0.7311))=-0.09086\).
The curvature adjustment is \(\frac{1}{2}h''(1)q=0.5(-0.09086)(0.800)=-0.03634\). Therefore the second-order approximation is \(0.7311-0.03634=0.6947\).
The negative adjustment is consistent with the logistic curve's concavity at a predictor of 1.00. This is an approximate response probability, not a new estimate of the conditional probability at \(b=0\), which remains \(h(1)=0.7311\). The distinction between conditional and random-effects-averaged quantities must be kept clear.
What the Correction Does Not Repair
PQL already uses an approximation to handle the nonlinear response. In a typical first-order construction, the response is represented by a working response and variance obtained from a local linearization, and the random effects are handled through a penalized linear mixed-model calculation. The second-order correction adds curvature information to that approximation. It does not turn the working response into an exact Bernoulli likelihood, nor does it make the random-effects distribution known from the data.
For binary outcomes, important challenges can remain even after adding the curvature term. Sparse clusters, few clusters, rare events, large random-effect variation, or random slopes can leave fixed effects and variance components weakly informed or sensitive to modeling choices. A correction aimed at one local approximation error cannot create information absent from the data, and it need not cancel the several errors that jointly affect a PQL fit.
Worked Example: A Large Variance Can Make the Correction Implausible
Consider a hypothetical calculation with \(\eta=0.500\) and \(q=30.0\), in squared log-odds units. Here \(h(0.500)=0.6225\), and \(h''(0.500)=0.6225(1-0.6225)(1-2(0.6225))=-0.05756\).
The second-order term is \(0.5(-0.05756)(30.0)=-0.8633\). The resulting approximation is \(0.6225-0.8633=-0.2409\), which is outside the allowable range for a probability.
This is not a possible fitted probability; it is evidence that the local expansion is being used beyond a credible range. The large variance allows \(U\) to move far from the expansion point, where the local curvature at \(\eta=0.500\) no longer describes the logistic curve adequately. Truncating the result to zero would hide the failure rather than validate the correction.
A small correction is not, by itself, evidence that the overall approximation is accurate. For instance, at \(\eta=0\), the logistic second derivative is zero, so the displayed second-order adjustment to the marginal mean is zero for any \(q\). Yet random effects can still induce substantial dependence among observations in the same cluster. A mean correction and a joint-model correction are not interchangeable.
Worked Example: A Corrected Mean Does Not Capture Cluster Dependence
In a hypothetical logistic random-intercept model, let the reference predictor be \(\eta=0\) and let the random intercept \(U\) have variance \(q=0.800\). For one observation, symmetry gives \(E\{h(U)\}=0.500\): since \(h(u)+h(-u)=1\), averaging over a symmetric distribution yields one half. The second-order mean adjustment also equals zero because \(h''(0)=0\).
Now take two conditionally independent observations in the same cluster, each with conditional success probability \(h(U)\). Their marginal covariance is \(\operatorname{Cov}(Y_{i1},Y_{i2})=\operatorname{Var}\{h(U)\}\), which is positive when \(U\) varies. To see its approximate size for modest variation, linearize \(h(U)\) at zero: \(h(U)\approx 0.5+0.25U\). Then \(\operatorname{Var}\{h(U)\}\approx (0.25)^2(0.800)=0.0500\).
The covariance value is itself only a first-order approximation, but it demonstrates the point: the marginal mean can be exactly 0.500 while within-cluster dependence remains. A second-order adjustment to the mean alone does not reproduce that dependence or ensure accurate variance-component estimates.
Assessing Use and Comparing Approaches
When examining a second-order PQL result, treat the correction as a diagnostic and approximation component, not as a certificate of reliability. Check whether its assumptions and numerical behavior are plausible in the fitted setting. Warning signs include a large estimated random-effect variance, many fitted probabilities near zero or one, very sparse event counts, small or highly unequal cluster sizes, convergence difficulties, and large changes in estimates when the approximation or model specification changes. These features are especially concerning in light of the established PQL difficulties with binary outcomes.
A practical comparison is about what each method approximates and what it can establish:
| Approach | What it adds or uses | Important limitation |
|---|---|---|
| First-order PQL | Local linearization and a penalized working-model fit | Can miss curvature in binary response means |
| Second-order PQL correction | A local quadratic curvature term, such as \(\frac{1}{2}h''(\eta)q\) | Still a truncated local approximation; does not guarantee reliable fixed effects or variance components |
| Likelihood-based integration | Targets the likelihood after accounting for random effects through numerical or analytic integration approximations | Requires its own computational and model checks; it is not automatically reliable under every data design |
Likelihood-based integration is a principal alternative when PQL warnings are serious, and the next topic treats its comparison with PQL for binary models. This does not mean that changing methods removes weak information or model misspecification. A second-order correction can be useful when the local approximation is reasonable, but it should not be mistaken for a general bias correction, an exact marginalization, or a substitute for checking the fitted model against the data structure.
Common Mistake
Wrong reasoning: “The second-order term accounts for nonlinearity, so the corrected PQL estimates must be unbiased or close to the likelihood estimates.” This conclusion does not follow. The term accounts for a particular local curvature contribution under its expansion assumptions; it leaves higher-order terms, the working-model approximation, and the information limitations of the binary data.
Correction: Describe exactly what was corrected, inspect whether the Taylor expansion is plausible over the random-effect variation, and interpret the fitted estimates with the remaining approximation uncertainty in mind. A correction that improves one component of an approximation is not proof that the full model fit is accurate.
Check Your Understanding
Answer each question, then compare with the brief solutions below.
- For a logistic mean, what determines the sign of the second-order adjustment?
- At \(\eta=0\), why can the second-order marginal-mean adjustment be zero even when random effects induce within-cluster dependence?
- For \(\eta=0.500\) and \(q=0.100\), calculate the second-order adjusted mean using \(h(0.500)=0.6225\) and \(h''(0.500)=-0.05756\).
- What does an adjusted probability below zero indicate in this approximation?
- Name two reasons a second-order correction does not guarantee reliable binary PQL estimates.
Solutions. (1) Its sign is determined by \(h''(\eta)\), since \(q\) is nonnegative: it is negative above probability 0.5 and positive below it. (2) At zero, \(h''(0)=0\), but the shared random effect still varies the conditional probability across clusters, creating positive covariance between conditionally independent outcomes. (3) The adjustment is \(0.5(-0.05756)(0.100)=-0.002878\), so the adjusted mean is \(0.6225-0.002878=0.6196\). (4) It signals that the local Taylor approximation is being used outside a plausible range; it is not a valid probability. (5) Higher-order terms remain, and sparse binary data or large random-effect variation can undermine estimation even if the local curvature term is included.