Fixed-Effect Bias Induced By Binary PQL

Why fixed effects can be distorted

In a binary mixed model, fixed effects describe conditional relationships: they relate predictors to the outcome while holding a cluster’s random effect fixed. Penalized quasi-likelihood (PQL) estimates these effects by repeatedly replacing the nonlinear binary model with a weighted linear mixed model. That replacement is computationally useful, but it is approximate. In particular, a good-looking working linear fit does not guarantee that its fixed-effect estimates closely match those from the binary mixed model’s likelihood.

The distortion is not a single constant shift. A slope can be attenuated because averaging over random intercepts flattens the relationship between a predictor and the marginal probability. The intercept can shift in either direction, depending on the baseline probability at the covariate value used to define it. These are related effects, but they are not interchangeable, and neither direction is guaranteed for every data set or PQL fit.

Definition: Consider independent clusters $j=1,\ldots,J$, with binary outcomes $Y_{ij}$ and covariates $x_{ij}$. A binary random-intercept logistic model assumes $$ Y_{ij}\mid b_j \ \text{are conditionally independent},\qquad Y_{ij}\mid b_j\sim\operatorname{Bernoulli}(\mu_{ij}),\qquad \operatorname{logit}(\mu_{ij})=\beta_0+\beta_1x_{ij}+b_j,\qquad b_j\sim N(0,\tau^2), $$ with random effects independent across clusters and independent of the covariates. The fixed effects $\beta_0,\beta_1$ describe conditional log odds: at a fixed value of $b_j$, a one-unit increase in $x$ changes the log odds by $\beta_1$. PQL fixed-effect bias here means systematic discrepancy caused by the PQL approximation, not ordinary sampling variation.

What the working response and weights do

At a current parameter value, write $\eta_{ij}=\beta_0+\beta_1x_{ij}+b_j$, $p_{ij}=\operatorname{logit}^{-1}(\eta_{ij})$, and $v_{ij}=p_{ij}(1-p_{ij})$. The binary-logit working response and its weight are

z_{ij}=\eta_{ij}+\frac{Y_{ij}-p_{ij}}{v_{ij}}, \qquad w_{ij}=v_{ij}.

The response is a first-order linearization: the residual $Y_{ij}-p_{ij}$ is divided by the derivative of the mean with respect to the linear predictor. The weight is the conditional Bernoulli variance divided by the squared derivative, which for the logit link simplifies to $v_{ij}$. PQL uses these quantities to fit a penalized weighted linear mixed model, updating the means, weights, and random effects iteratively.

At the true current predictor, the working response has conditional mean $\eta_{ij}$, since $E(z_{ij}\mid b_j,x_{ij})=\eta_{ij}+E(Y_{ij}-p_{ij}\mid b_j,x_{ij})/v_{ij}=\eta_{ij}$. That identity matters: it would be inaccurate to say that the working response is automatically biased at each observation just because the outcome is binary. But its conditional variance is $1/v_{ij}$, and treating the transformed response as Gaussian is an approximation. Further, the weights and predictors are estimated iteratively, and the random-effect distribution and nonlinear likelihood are not represented exactly by one weighted linear fit. The mismatch can therefore affect the fitted fixed effects even though the conditional mean identity holds at a fixed, correct predictor.

Worked Example: Calculating a Binary Working Response

Suppose a hypothetical observation has current linear predictor $\eta=-1.3863$, so $p=0.2000$, and its observed outcome is $Y=1$. Then $v=p(1-p)=0.2000(0.8000)=0.1600$. Its working response and weight are

z=-1.3863+\frac{1-0.2000}{0.1600} =-1.3863+5.000 =3.6137,\qquad w=0.1600.

The large working response is not a predicted probability; it is a linearized response on the logit scale. If instead $Y=0$, then $z=-1.3863+(0-0.2000)/0.1600=-2.6363$, with the same weight. Under the assumed conditional model, the expected working response is $0.2000(3.6137)+0.8000(-2.6363)=-1.3863$, the current predictor. The calculation illustrates both the mean identity and how a binary outcome can yield a highly variable transformed response when $v$ is small.

Why slope and intercept effects differ

A useful way to see the distinction is to compare the conditional model with its marginal mean, averaged over the random intercept. Let $f(\eta)=\operatorname{logit}^{-1}(\eta)$, and let $\eta=\beta_0+\beta_1x$. For small $\tau^2$, a Taylor expansion around $b=0$ gives a local approximation to the marginal probability $m(x)=E_b\{f(\eta+b)\}$. Since $E(b)=0$, $E(b^2)=\tau^2$, and the odd central moments of a normal random variable are zero, the second-order terms give

m(x)\approx f(\eta)+\frac{\tau^2}{2}f''(\eta) =f(\eta)+\frac{\tau^2}{2}f(\eta)\{1-f(\eta)\}\{1-2f(\eta)\}.

To translate this probability shift to the logit scale, use the first-order expansion of $\operatorname{logit}(u)$ around $u=f(\eta)$. Its derivative is $1/[f(\eta)\{1-f(\eta)\}]$. Thus

\operatorname{logit}\{m(x)\} \approx \eta+\frac{\tau^2}{2}\{1-2f(\eta)\}.

This is a local approximation, not an exact identity or a formula for every PQL estimate. Its value is that it exposes two distinct mechanisms. At a fixed reference value $x$, the approximate marginal-logit shift is $\frac{\tau^2}{2}(1-2p)$, where $p=f(\beta_0+\beta_1x)$. It is positive when $p<0.5$, negative when $p>0.5$, and approximately zero when $p=0.5$. The intercept shift therefore depends on where $x=0$ lies and on the baseline probability there.

Differentiating the same approximation with respect to $x$ gives an approximate marginal-logit slope:

\frac{d}{dx}\operatorname{logit}\{m(x)\} \approx \beta_1\left[1-\tau^2p(1-p)\right].

For a positive $\beta_1$ and sufficiently small $\tau^2$, this local slope is smaller than the conditional slope. The marginal relationship is flatter because mixing over random intercepts smooths the probability curve. This helps explain why binary PQL can show fixed-slope attenuation, but it does not prove that a particular PQL estimate must be attenuated by this amount. A PQL estimate depends on the design, fitted weights, cluster information, and variance-component estimate as well as on this marginal relationship.

Worked Example: A Local Slope Attenuation Calculation

Consider a hypothetical model with $\beta_1=1.00$, $\tau^2=0.200$, and baseline probability $p=0.200$ at the covariate value of interest. The second-order marginal-logit slope approximation is

1.00\left[1-(0.200)(0.200)(0.800)\right] =1.00(0.968) =0.968.

The local marginal-logit slope is therefore approximately $0.968$ per one-unit increase in $x$, compared with the conditional slope of $1.00$. This is a reduction of about $0.032$ log-odds units per covariate unit in this approximation. It illustrates a possible direction and scale of slope attenuation; it is not a computed PQL estimate, and higher-order terms in $\tau^2$ have been omitted.

Worked Example: The Intercept Shift Changes With Baseline Risk

Suppose $\tau^2=0.360$, and compare two hypothetical reference points. At a reference point where $p=0.200$, the approximate intercept-scale shift is

\frac{0.360}{2}(1-2(0.200)) =0.180(0.600) =0.108.

At a reference point where $p=0.800$, it is

\frac{0.360}{2}(1-2(0.800)) =0.180(-0.600) =-0.108.

Thus the same random-intercept variance gives opposite approximate shifts at these two baseline probabilities. The intercept is the log odds at $x=0$, so its interpretation depends on the covariate’s reference point. Centering $x$ at a different value changes which baseline the intercept represents. It is therefore not sound to claim that binary PQL always biases the intercept upward or always biases it toward zero.

Interpreting discrepancies and diagnosing trouble

A fitted PQL coefficient is an approximation to a conditional-model coefficient, not automatically a marginal effect. A comparison is meaningful only when both methods use the same outcome, covariates, random-effects structure, and interpretation. If a PQL slope is smaller than a likelihood-based estimate, attenuation is plausible, but the difference could also reflect sampling variability or a poorly determined fit. Likewise, an intercept difference should be interpreted at the stated covariate reference value rather than as a universal direction of bias.

Useful checks focus on whether the approximation is consequential and stable:

Compare PQL fixed effects with estimates from a likelihood-based method using Laplace approximation or adaptive quadrature, when feasible. Large changes are a warning that the working-model approximation matters.
Check whether the fitted random-intercept variance is sensitive to the estimation method. Underestimating that variance can alter how much heterogeneity the fixed effects must accommodate, though it does not by itself determine the direction of either fixed-effect bias.
Inspect working weights and fitted probabilities. Very small $p_{ij}(1-p_{ij})$ values give low weights and can produce extreme working responses, making the local Gaussian approximation fragile.
Check convergence and sensitivity to reasonable starting values. An apparently stable iteration is not evidence that the approximation has negligible bias; it only indicates that the algorithm has settled on its approximate solution.

These checks do not establish that one method is correct in every finite sample. Likelihood-based integration is generally preferable when its computational cost is manageable because it targets the specified likelihood rather than the PQL working approximation. PQL can be attractive for speed and for obtaining initial estimates, but agreement with a likelihood-based fit should be assessed when fixed-effect accuracy is important. A marginal model is another possible analysis, but it targets a population-averaged relationship rather than the conditional coefficients in the random-intercept model; its coefficients should not be treated as a direct alternative estimate of the same parameter without clarifying the estimand.

Common mistake and limitations

Common pitfall: “The working response has conditional mean equal to the linear predictor, so PQL cannot bias the fixed effects.” The equality holds at a fixed, correct predictor. It does not make the Gaussian working fit an exact likelihood analysis, nor does it remove the effects of estimated weights, iterative linearization, or approximated random-effect variation. The correction is to distinguish the working-response identity from the accuracy of the complete PQL estimator.

The local expansion above assumes a normal random intercept and a variance small enough for a second-order approximation to be informative. It explains why slope attenuation and baseline-dependent intercept shifts are plausible; it is not a universal bias correction. In a finite sample, the sign and size of PQL’s fixed-effect discrepancy also depend on the covariate pattern, cluster information, and estimation variability. The expansion should therefore support interpretation and diagnosis, not replace comparison with an appropriate likelihood-based analysis.

Check Your Understanding

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

At a fixed correct predictor, what is $E(z_{ij}\mid b_j,x_{ij})$ for the binary-logit working response, and why does this not establish that PQL is unbiased?
For $\tau^2=0.100$ and $p=0.500$, calculate the second-order intercept-scale shift and the approximate marginal-logit slope factor relative to $\beta_1$.
For $\tau^2=0.400$ and baseline $p=0.100$, calculate the approximate intercept-scale shift. State its direction.
Why is “PQL always biases the intercept downward” not a defensible general statement?

Solutions. (1) The conditional mean is $\eta_{ij}$, because $E(Y_{ij}-p_{ij}\mid b_j,x_{ij})=0$. This is an identity for the working response at a fixed correct predictor, not a proof that the iteratively fitted weighted Gaussian model reproduces the binary mixed-model likelihood. (2) The shift is $(0.100/2)(1-2(0.500))=0$. The slope factor is $1-(0.100)(0.500)(0.500)=0.975$, so the approximate slope is $0.975\beta_1$. (3) The shift is $(0.400/2)(1-2(0.100))=0.200(0.800)=0.160$, upward on the log-odds scale at this reference point. (4) The local shift changes sign with baseline probability: it is positive below $0.5$ and negative above $0.5$. In addition, the intercept refers to $x=0$, so its baseline depends on covariate coding. These approximations do not guarantee the exact direction of a PQL discrepancy in any particular sample.