PQL Bias For Random Slopes In Binary Models

Why Random-Slope PQL Needs Particular Care

A random slope allows the association between a covariate and a binary response to differ across clusters. Estimating that variation is challenging because binary observations carry limited information, and PQL replaces the nonlinear conditional mean with a local linear approximation. When the covariate has little support within clusters, or when random effects move fitted probabilities into regions of strong curvature, the approximation can poorly represent the information available about slope variation. The resulting variance and covariance estimates can be biased, unstable, or weakly identified.

These are risks, not automatic outcomes. Within-cluster covariate variation often gives direct information about differences in slopes, but it is not the only possible source: between-cluster differences in constant covariate values can also inform the random-effect variance function. The amount of information depends on the covariate support, replication, and model assumptions.

Model and PQL Approximation

For cluster $i$ and observation $j$, let $Y_{ij}$ be binary and, conditional on the cluster random effects, assume independent Bernoulli outcomes. In a logit random-intercept and random-slope model, write

Y_{ij}\mid b_i \sim \operatorname{Bernoulli}(\mu_{ij}), \qquad \operatorname{logit}(\mu_{ij})=\beta_0+\beta_1 x_{ij}+b_{0i}+b_{1i}x_{ij}, \qquad \begin{pmatrix}b_{0i}\\b_{1i}\end{pmatrix}\sim N\!\left( \begin{pmatrix}0\\0\end{pmatrix},G\right), \quad G=\begin{pmatrix}\tau_{00}&\tau_{01}\\\tau_{01}&\tau_{11}\end{pmatrix}.

The clusters are assumed independent, the random effects have mean zero and common covariance matrix $G$, and the conditional response model and logit link are correctly specified. The covariance matrix must be positive semidefinite: $\tau_{00}\geq 0$, $\tau_{11}\geq 0$, and $\tau_{00}\tau_{11}-\tau_{01}^{2}\geq 0$. Conditional independence does not mean observations are marginally independent; integrating over shared random effects induces within-cluster dependence.

Definition: At covariate value $x$, the random contribution to the linear predictor is $b_{0i}+xb_{1i}$. Its variance is $$ q(x)=\operatorname{Var}(b_{0i}+xb_{1i})=\tau_{00}+2x\tau_{01}+x^2\tau_{11}. $$ The random-slope variance is $\tau_{11}$; $\tau_{01}$ is the random-intercept/random-slope covariance. Neither parameter alone describes all between-cluster variation across the observed covariate range.

For PQL, the binary mean is locally linearized around a current linear predictor $\eta_{ij}$. With $\mu_{ij}=\operatorname{logit}^{-1}(\eta_{ij})$, the working response and weight are

Y^*_{ij}=\eta_{ij}+\frac{Y_{ij}-\mu_{ij}}{\mu_{ij}(1-\mu_{ij})}, \qquad w_{ij}=\mu_{ij}(1-\mu_{ij}).

The iterative working model treats $Y^*_{ij}$ as approximately Gaussian with working residual variance $1/w_{ij}$, while retaining the mixed-model random effects. Its variance-component estimates therefore use a locally approximated working model, not the exact marginal likelihood of the binary observations. This distinction matters especially for a random slope: the data must inform how cluster-specific linear predictors change over $x$, while the approximation and its weights change with the fitted probabilities.

Where Sparse Covariate Support Enters

If $x$ varies within a cluster, observations at different covariate values can help separate a cluster-specific shift in baseline from a cluster-specific change in slope. With little within-cluster variation, that direct contrast is weak. However, variation among clusters whose $x$ values are constant within cluster can also provide information: clusters at different $x$ values experience random predictor variation with variance $q(x)$. Since $q(x)$ is quadratic, three distinct covariate values can, in principle, provide three constraints on its three covariance parameters, if there is adequate replication and the assumed model is appropriate. Two values alone do not generally determine all three parameters from $q(x)$.

Worked Example: Three Covariate Values Can Identify the Variance Function

Suppose a hypothetical study has many replicated clusters at each of three covariate values, with $x$ constant within each cluster. For illustration, suppose the variance of the random predictor contribution is estimated as $q(0)=0.36$, $q(1)=0.85$, and $q(2)=1.84$. Treat these values as exact for this algebraic demonstration.

Write $a=\tau_{00}$, $c=\tau_{01}$, and $d=\tau_{11}$. The variance-function equations are

q(0)=a=0.36,\qquad q(1)=a+2c+d=0.85,\qquad q(2)=a+4c+4d=1.84.

Subtracting the first equation from the second gives $2c+d=0.49$. Subtracting it from the third gives $4c+4d=1.48$, or $c+d=0.37$. Subtracting $2(c+d)=0.74$ from $2c+d=0.49$ gives $-d=-0.25$, hence $d=0.25$. This is consistent with a nonnegative variance. The arithmetic is: from $q(2)-a=1.84-0.36=1.48$, we have $4c+4d=1.48$, hence $c+d=0.37$. From $q(1)-a=0.49$, $2c+d=0.49$. Subtracting $c+d=0.37$ from $2c+d=0.49$ gives $c=0.12$; then $d=0.25$.

Thus $\tau_{00}=0.36$, $\tau_{01}=0.12$, and $\tau_{11}=0.25$. The matrix is admissible because its determinant is $0.36(0.25)-0.12^2=0.0756>0$. This demonstrates possible identification from between-cluster support alone; it does not guarantee precise PQL estimates in a finite sample.

By contrast, two support values provide only two equations for three unknowns. For example, $q(0)=0.36$ and $q(1)=0.85$ are compatible both with $(\tau_{00},\tau_{01},\tau_{11})=(0.36,0.12,0.25)$ and with $(0.36,0.15,0.19)$. Both give $q(1)=0.85$; the second matrix is valid because $0.36(0.19)-0.15^2=0.0459>0$. Additional information in the design or model may help, but those two values of $q$ alone do not distinguish the parameter triples.

Curvature and the Direction of Approximation Error

For the logistic mean, $\mu(\eta)=1/(1+e^{-\eta})$, the first derivative is $\mu'(\eta)=\mu(1-\mu)$, and the second derivative is $\mu''(\eta)=\mu(1-\mu)(1-2\mu)$. A first-order Taylor expansion around $\eta$ gives $\mu(\eta+\delta)\approx\mu(\eta)+\mu'(\eta)\delta$. The neglected second-order term is approximately $\frac12\mu''(\eta)\delta^2$. Its sign changes around $\mu=0.5$, and its size can matter when random slopes produce substantial changes in the linear predictor over the observed range.

This local expansion explains a mechanism, not a universal formula for the bias in $\hat{\tau}_{11}$. PQL fits the working response and weights iteratively; it does not simply replace the true random-slope variance by a fixed shrunken value. In binary models, PQL often underestimates random-effect variation, but the magnitude and even the direction of bias for the slope variance or intercept-slope covariance depend on the design, parameter values, and fitted probabilities.

Worked Example: A Tangent Approximation Misses a Curved Mean

Consider a hypothetical observation with current probability $\mu=0.90$, so $\eta=\log(0.90/0.10)=\log(9)\approx2.197$. Suppose a random-slope change contributes $\delta=1$ to the linear predictor. The first-order tangent approximation predicts

\mu(\eta)+\mu'(\eta)\delta=0.90+(0.90)(0.10)(1)=0.990.

The actual probability after the shift is $\operatorname{logit}^{-1}(3.197)\approx0.9607$. The tangent prediction is therefore about $0.0293$ higher. The second derivative at the expansion point is $0.90(0.10)(1-1.80)=-0.072$, so the second-order contribution is $\frac12(-0.072)(1)^2=-0.036$, consistent in sign and approximate scale with the observed discrepancy.

The working weight at the starting point is $0.90(0.10)=0.09$, compared with $0.25$ at probability $0.50$. Thus observations near the probability extremes have less local working weight, while curvature can make a large random-slope shift poorly represented by a tangent. Together these features can weaken or distort information about slope heterogeneity. They do not, by themselves, determine the direction or size of the resulting variance bias.

Centering, Covariance, and Interpretation

The intercept-slope covariance depends on where $x=0$ is placed. This is not a defect: the random intercept is the cluster deviation at the chosen zero point. Centering can make that point meaningful and can reduce avoidable dependence between estimated components, but it does not create covariate information or remove PQL’s nonlinear approximation.

Worked Example: Re-centering Changes the Covariance, Not the Slope Variance

Use the hypothetical estimates $\tau_{00}=0.36$, $\tau_{01}=0.12$, and $\tau_{11}=0.25$ from the variance-function example. Define $x^*=x-1$, so that $x=x^*+1$. The random predictor is

$$ b_0+b_1x=(b_0+b_1)+b_1x^*. $$

The centered random intercept is $b_0^*=b_0+b_1$. Therefore,

\operatorname{Var}(b_0^*)=0.36+2(0.12)+0.25=0.85,\qquad \operatorname{Cov}(b_0^*,b_1)=0.12+0.25=0.37,\qquad \operatorname{Var}(b_1)=0.25.

The intercept variance and covariance change because the intercept is now defined at $x=1$, while the slope variance remains $0.25$. The variance at $x=1$ is $0.85$ under either coding. Consequently, a reported covariance is interpretable only alongside the covariate’s origin and scale.

Diagnostics, Alternatives, and Limits

Warning signs include a slope-variance estimate at or near zero despite substantial cluster-to-cluster pattern differences, a covariance estimate near the positive-semidefinite boundary, large uncertainty, or marked changes under reasonable covariate scaling and centering. Check the number and spread of distinct $x$ values within and between clusters, the amount of replication at those values, and whether fitted probabilities cluster near zero or one. Compare the fitted mean curve with observed patterns across the covariate range, while remembering that a visual discrepancy can also reflect sampling variability.

A random-intercept-only model is simpler but rules out cluster-specific slopes; it is not an unbiased substitute when slope variation is real. A standard GLM omits the cluster dependence altogether. Likelihood-based integration methods, including Laplace approximation and adaptive quadrature, retain the nonlinear conditional response model more directly and provide useful comparisons with PQL. They are not guaranteed to solve weak identification or small-sample uncertainty, but disagreement with PQL is a reason to investigate approximation sensitivity rather than to select the more favorable estimate.

The main limitation is that sparse support and nonlinear curvature can make the variance components hard to recover even when the model is correctly specified. If all clusters share nearly the same $x$, the data contain little information about how random predictor variation changes with $x$; if support is limited to two values, $q(x)$ alone cannot separate all three covariance parameters. Conversely, absent within-cluster variation does not prove non-identifiability when there are three or more replicated support values. Conclusions must account for the actual design and the approximation method.

Common Mistake

Wrong reasoning: “There is no within-cluster change in $x$, so random-slope variation cannot be estimated.” This overlooks information from clusters at different constant covariate values. Their variation can inform $q(x)$, and three distinct values can in principle constrain its three parameters under suitable replication and assumptions. The correction is to assess all covariate support and its replication, while recognizing that two support values generally leave the variance function underdetermined and that PQL may still be biased or imprecise.

Key takeaway: Random-slope PQL bias is shaped jointly by covariate support, the random-effect covariance structure, and the curvature and weighting of the binary mean. Treat a small or unstable slope-variance estimate as a design- and approximation-dependent result, not as automatic evidence that slope heterogeneity is absent.

Check Your Understanding

Answer each question, then compare with the brief solutions.

In the random-intercept/random-slope model, write $q(x)$ and identify which parameter is the random-slope variance.
Why can clusters with constant but different $x$ values still inform random-effect variation? What is the limitation if there are only two distinct values?
At $\mu=0.80$, calculate the logistic working weight and the second derivative of the mean with respect to $\eta$.
If $x^*=x-2$, state whether the random-slope variance changes and give the centered intercept in terms of the original random effects.

Solutions. (1) $q(x)=\tau_{00}+2x\tau_{01}+x^2\tau_{11}$; the slope variance is $\tau_{11}$. (2) The variance of the random predictor contribution differs with $x$, following $q(x)$. Three distinct support values can constrain its three parameters; two values alone generally supply only two constraints. (3) The weight is $0.80(0.20)=0.16$, and the second derivative is $0.16(1-1.60)=-0.096$. (4) The slope variance is unchanged; since $x=x^*+2$, the centered intercept is $b_0^*=b_0+2b_1$.