PQL Under Large Random-Effect Variance

Why Random-Effect Spread Matters

In a binary generalized linear mixed model, a random intercept allows the event probability to vary among clusters. When the random-effect variance is large, that variation can span much of the logistic curve: some clusters have conditional probabilities near zero, others near one, and many may lie in the steep transition between them. Penalized quasi-likelihood (PQL) replaces the nonlinear conditional mean with a local first-order approximation. That is useful when the approximation is local in practice, but a broad random-effect distribution can place substantial probability mass far from its expansion point.

The issue is not simply that the random effect is variable. It is that the logistic mean is curved and bounded, while its first-order Taylor approximation is a straight line that is unbounded. A line can represent a smooth curve near one point; it cannot represent the curve accurately across an arbitrarily wide range of random effects.

Model and First-Order Approximation

Definition: Consider a binary random-intercept logistic model. Conditional on the cluster random intercept $b_j$, observations $Y_{ij}$ are independent Bernoulli variables with probability $p_{ij}(b_j)$, where $$ \operatorname{logit}\{p_{ij}(b_j)\}=x_{ij}^{T}\beta+b_j,\qquad b_j\sim N(0,\tau^2). $$ Here $i$ indexes observations, $j$ indexes clusters, $x_{ij}$ is the fixed-effect covariate vector, $\beta$ is the fixed-effect coefficient vector, and $\tau^2$ is the random-intercept variance. The conditional probability is $$ p_{ij}(b_j)=\frac{\exp(x_{ij}^{T}\beta+b_j)}{1+\exp(x_{ij}^{T}\beta+b_j)}. $$

This model assumes a normal random-intercept distribution with mean zero and finite variance, and conditional independence of the binary outcomes given the random effects and covariates. Its conditional probabilities are not the same as probabilities averaged over the random-effect distribution. For a fixed covariate vector $x$, the latter is

P(Y=1\mid x)=\int \frac{\exp(x^{T}\beta+b)}{1+\exp(x^{T}\beta+b)} \frac{1}{\tau}\phi\!\left(\frac{b}{\tau}\right)\,db,

where $\phi$ is the standard normal density. PQL’s working construction instead uses a local expansion of the conditional mean. Write $\eta=x^{T}\beta+b$, and expand around a current value $\eta_0$. For the logistic function $p(\eta)=\exp(\eta)/\{1+\exp(\eta)\}$,

p(\eta_0+\delta)\approx p(\eta_0)+p'(\eta_0)\delta =p_0+p_0(1-p_0)\delta,

where $p_0=p(\eta_0)$ and $\delta=\eta-\eta_0$. The slope $p_0(1-p_0)$ is positive and at most $1/4$. As established in the earlier material on binary PQL, this linearization leads to a working response and weight for fitting a penalized linear mixed model. The concern here is how well that local line represents the conditional mean across the random-effect values that matter.

Why a Broad Distribution Strains the Line

The derivative changes with $\eta$: it is largest at $p=0.5$, and it approaches zero as $p$ approaches either zero or one. The second derivative is

p''(\eta)=p(\eta)\{1-p(\eta)\}\{1-2p(\eta)\}.

Thus the logistic curve is curved except at its inflection point. Taylor’s theorem makes the cost of ignoring that curvature explicit:

p(\eta_0+\delta) =p_0+p_0(1-p_0)\delta+\frac{1}{2}p''(\xi)\delta^2,

for some $\xi$ between $\eta_0$ and $\eta_0+\delta$. The first-order approximation omits the final term. If the random-effect distribution puts substantial mass at values for which $\delta$ is large, this omitted term need not be small. In particular, a random-effect standard deviation comparable to or larger than a unit on the logit scale can span a substantial part of the curve’s transition from near zero to near one.

There is also a direct warning sign: the linear approximation can predict values below zero or above one, even though a probability must lie between zero and one. Such predictions do not mean the fitted model literally assigns an invalid Bernoulli probability. They show that the local line is a poor representation of the conditional mean at those random-effect values.

Worked Example: A Linearization Across a Wide Range

Suppose, hypothetically, that a binary outcome follows the model above with fixed linear predictor $x^{T}\beta=0.5$, random intercept $b\sim N(0,4)$, and an expansion at $b=0$. The expansion point has $\eta_0=0.5$.

At that point, $p_0=\operatorname{logit}^{-1}(0.5)=0.6225$, and $p_0(1-p_0)=0.6225(0.3775)=0.2350$. The first-order approximation is therefore

p_{\mathrm{lin}}(b)=0.6225+0.2350b.

The following calculations compare it with the exact conditional probability $p(b)=\operatorname{logit}^{-1}(0.5+b)$.

Random intercept $b$	Exact conditional probability	First-order approximation
-4	$\operatorname{logit}^{-1}(-3.5)=0.0293$	$0.6225+0.2350(-4)=-0.3175$
-2	$\operatorname{logit}^{-1}(-1.5)=0.1824$	$0.6225+0.2350(-2)=0.1525$
0	$\operatorname{logit}^{-1}(0.5)=0.6225$	$0.6225$
2	$\operatorname{logit}^{-1}(2.5)=0.9241$	$0.6225+0.2350(2)=1.0925$
4	$\operatorname{logit}^{-1}(4.5)=0.9890$	$0.6225+0.2350(4)=1.5625$

The line matches at the expansion point and is reasonably close at some nearby values. Farther away, it gives impossible probabilities and increasingly misses the saturation of the logistic curve. Since the random-effect standard deviation is $2$, these are not all remote values on the scale of the assumed distribution.

Worked Example: How Much Random-Effect Mass Has an Impossible Linearized Probability?

Continue with the hypothetical model in the previous calculation. The linearized probability is $0.6225+0.2350b$. To find where it is between zero and one, solve both inequalities.

0\leq 0.6225+0.2350b\leq 1 \quad\Longrightarrow\quad -2.649\leq b\leq 1.606.

For example, the lower endpoint is $-0.6225/0.2350=-2.649$, and the upper endpoint is $(1-0.6225)/0.2350=1.606$. Under $b\sim N(0,4)$, standardize by dividing by the standard deviation $2$. The probability that the linearized value is between zero and one is

P(-2.649\leq b\leq 1.606) =\Phi(1.606/2)-\Phi(-2.649/2) \approx \Phi(0.803)-\Phi(-1.325) \approx 0.696.

Thus approximately $1-0.696=0.304$, or 30.4%, of the assumed random-effect distribution lies where this particular linearization predicts a value outside the probability range. This calculation is a diagnostic of the expansion’s reach, not an estimate of model bias or a test of PQL. The expansion point in an actual PQL iteration can differ from zero, but the same check can be made using its current expansion point and slope.

Curvature and Averaging Over Random Effects

The problem is not limited to invalid linearized probabilities. Even where the line remains between zero and one, averaging a nonlinear conditional probability over random effects is generally not equivalent to evaluating the probability at the mean random effect. A short Taylor calculation shows why variance matters. For fixed $x$, write $\eta=x^{T}\beta+b$ and expand around $\eta_0=x^{T}\beta$, so $\delta=b$. Since $E(b)=0$, the linear term averages to zero, but the quadratic term depends on $E(b^2)=\tau^2$:

E\{p(x^{T}\beta+b)\} \approx p(x^{T}\beta)+\frac{1}{2}p''(x^{T}\beta)\tau^2.

This is a local Taylor approximation to the marginal probability, not an exact identity. It indicates that the effect of curvature is on the scale of $\tau^2$; when $\tau^2$ is large, the approximation itself may also be inadequate. The sign of the displayed second-order term depends on the fixed linear predictor, because $p''$ changes sign at zero.

Worked Example: The Omitted Curvature Term Grows With Variance

For a hypothetical covariate pattern with $x^{T}\beta=1$, the conditional probability at $b=0$ is $p(1)=0.7311$. Its derivatives are

p'(1)=0.7311(1-0.7311)=0.1966,\qquad p''(1)=0.7311(1-0.7311)(1-2(0.7311))\approx -0.0909.

With $\tau=0.5$, the second-order term in the Taylor calculation is

\frac{1}{2}p''(1)\tau^2 =\frac{1}{2}(-0.0909)(0.25) =-0.0114.

This term is small compared with $0.7311$, though the size of the term alone does not prove that the entire approximation is accurate. If instead $\tau=2$, the same calculation gives

\frac{1}{2}(-0.0909)(4)=-0.1818.

That much larger curvature contribution signals that a first-order description deserves scrutiny. It should not be reported as a reliable marginal probability: at $\tau=2$, higher-order terms can also matter, so the full nonlinear average is needed for an accurate marginal probability.

Assumptions, Diagnostics, and Alternatives

The argument concerns a binary logit random-intercept model with a normally distributed random effect, conditional Bernoulli outcomes, and a PQL step that linearizes the conditional mean at a current predictor. Large $\tau^2$ is a reason for caution, not a guarantee that every PQL estimate will be poor. The practical concern depends on how much of the random-effect distribution lies far from the expansion point and on how quickly the conditional probabilities change there.

Several checks help reveal when the local approximation may be strained:

Compare the estimated random-effect standard deviation with the logit-scale range over which fitted conditional probabilities change substantially. A spread of several logit units traverses much of the logistic curve.
At representative covariate patterns, evaluate the exact conditional probability and the current first-order approximation over a range of plausible random-effect values, not only at the expansion point.
Check the fraction of the assumed random-effect distribution for which the linearized probability is outside zero and one. A large fraction indicates that the line is being used beyond a plausible local neighborhood.
Examine whether conclusions change materially under a method that retains the nonlinear conditional mean and integrates over the random-effect distribution, such as likelihood-based integration using quadrature.

These checks assess approximation sensitivity; they do not replace an analysis of sampling uncertainty, and a large estimated variance alone does not quantify PQL bias. Conversely, a plausible-looking fit does not establish that the linearization is adequate. PQL can be computationally convenient, but its working linear mixed model is not the original nonlinear likelihood. Likelihood-based integration is an important alternative when the local approximation is not credible, though it too requires appropriate computation and model assumptions.

Key takeaway: Large random-effect variance broadens the range of conditional logits. Because the logistic mean is curved and bounded but its first-order approximation is an unbounded line, PQL can be strained when substantial random-effect probability mass lies far from the expansion point. Evaluate the approximation across that distribution rather than judging it only at its point of contact.

Common Pitfall

Wrong reasoning: “The random intercept has mean zero, so the first-order term averages to zero; therefore PQL’s approximation is adequate even when its variance is large.” The first clause is true for an expansion centered at the mean random effect, but it does not remove curvature. The Taylor remainder contains a term proportional to the squared deviation, and its scale grows with $\tau^2$. In addition, averaging the first-order line does not make it match the nonlinear conditional probabilities away from the expansion point. The correction is to inspect the conditional probability curve across plausible random effects and use a method that retains the nonlinear mean when the linear approximation is inadequate.

Check Your Understanding

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

For a logit model with fixed linear predictor zero, what is the first-order approximation to $p(b)$ at $b=0$, and for which $b$ values does it fall outside zero and one?
Why does a zero-mean random intercept not make the first-order PQL approximation exact after averaging?
At a fixed linear predictor, what quantity controls the scale of the second-order Taylor contribution when $b\sim N(0,\tau^2)$?
Name one diagnostic and one alternative approach for a binary random-intercept model with substantial random-effect spread.

Solutions. (1) Since $p(0)=0.5$ and $p'(0)=0.25$, the line is $0.5+0.25b$. It is below zero for $b<-2$ and above one for $b>2$. (2) The zero mean removes the expected first-order term, but not the curvature remainder, which depends on squared deviations; the nonlinear conditional mean is not generally equal to its tangent line away from the expansion point. (3) The second-order contribution is $\frac{1}{2}p''(x^{T}\beta)\tau^2$, so it scales with $\tau^2$, as well as with the curvature at the expansion point. (4) One can calculate the share of the random-effect distribution where the linearized probability is outside zero and one, or compare conditional probabilities with their linearized values across plausible effects. A likelihood-based method that integrates the nonlinear conditional mean, for example using quadrature, is an alternative.

Random intercept \(b\)	Exact conditional probability	First-order approximation
-4	\(\operatorname{logit}^{-1}(-3.5)=0.0293\)	\(0.6225+0.2350(-4)=-0.3175\)
-2	\(\operatorname{logit}^{-1}(-1.5)=0.1824\)	\(0.6225+0.2350(-2)=0.1525\)
0	\(\operatorname{logit}^{-1}(0.5)=0.6225\)	\(0.6225\)
2	\(\operatorname{logit}^{-1}(2.5)=0.9241\)	\(0.6225+0.2350(2)=1.0925\)
4	\(\operatorname{logit}^{-1}(4.5)=0.9890\)	\(0.6225+0.2350(4)=1.5625\)