Worked Example Of PQL With Rare Binary Outcomes

Why Rare Events Need Careful Reading

In a binary generalized linear mixed model, a rare event can produce a surprisingly large PQL pseudo-response. That number is not an observed outcome or a predicted probability. It is a local working quantity formed from the current model fit. To understand its role, it must be read together with its working weight and the cluster-level penalty.

This tutorial calculates those quantities for hypothetical clustered data. The aim is to distinguish three things that are easy to conflate: the size of a pseudo-response, the weight assigned to it in the working mixed model, and the effect of the observation on an update. All calculations refer to a specified current iteration; they are not presented as final fitted estimates.

Model and Working Quantities

Assume independent binary outcomes conditional on a normally distributed random intercept. For individual $i$ in cluster $j$, let $Y_{ij}$ be 1 for an event and 0 otherwise. At the current PQL iteration, the conditional mean is $\mu_{ij}$, and the linear predictor is $\eta_{ij}=\mathbf{x}_{ij}^{\mathsf T}\boldsymbol{\beta}+b_j$.

Definition: Under the binary logit random-intercept model, $Y_{ij}\mid b_j$ is Bernoulli with mean $\mu_{ij}$, $\operatorname{logit}(\mu_{ij})=\eta_{ij}=\mathbf{x}_{ij}^{\mathsf T}\boldsymbol{\beta}+b_j$, and $b_j\sim N(0,\tau^2)$. Conditional on the random effects, the outcomes are assumed independent. At a current iterate, the PQL working weight and pseudo-response are $$ w_{ij}=\mu_{ij}(1-\mu_{ij}), \qquad z_{ij}=\eta_{ij}+\frac{Y_{ij}-\mu_{ij}}{w_{ij}}. $$ The working linear mixed model treats $z_{ij}$ locally as having mean $\eta_{ij}$ and residual variance $1/w_{ij}$, while retaining the random-effect distribution and its penalty.

These quantities use the current $\eta_{ij}$, not a predictor recomputed from the pseudo-response. For the logit link, $d\mu/d\eta=\mu(1-\mu)$. A first-order Taylor expansion of the mean around the current predictor gives $$ \mu(\eta+\Delta)\approx\mu(\eta)+w\Delta. $$ Setting this local approximation equal to the observed $Y$ gives $\Delta\approx(Y-\mu)/w$, which yields $z=\eta+\Delta$. The Bernoulli variance is $\mu(1-\mu)=w$; dividing this variance by the squared derivative $w^2$ gives the working residual variance $1/w$. This is the local linearization underlying the PQL working quantities.

The normal random-intercept assumption is used for the Gaussian random-effect penalty, and conditional independence gives the observation-level Bernoulli variances used here. The linearization is local: it is most credible when the mean response changes nearly linearly over the relevant predictor range. With rare events, the derivative $w$ can be small, making the correction in $z$ large. That is a property of the approximation, not evidence that the observed outcome itself is large.

Worked Example: One Rare Event at a Given Predictor

Suppose a hypothetical individual has observed outcome $Y=1$ and current linear predictor $\eta=-2.3918$. Calculate the probability, weight, and pseudo-response using the same predictor throughout.

Mean. The inverse-logit gives $$ \mu=\frac{1}{1+\exp(-\eta)} =\frac{1}{1+\exp(2.3918)} \approx 0.08380. $$

Weight. $$ w=\mu(1-\mu)=0.08380(0.91620)\approx 0.07678. $$

Pseudo-response. For an event, $Y-\mu=1-\mu$, so $(Y-\mu)/w=1/\mu$. Therefore $$ z=-2.3918+\frac{1}{0.08380}\approx -2.3918+11.9332=9.5414. $$ The pseudo-response is about 9.541, although the observed outcome is 1 and the current fitted probability is about 0.08380. It is a working response on the linear-predictor scale, not a probability.

Why the Rare-Event Pseudo-Response Can Be Extreme

The event and non-event formulas make the mechanism explicit. Since $w=\mu(1-\mu)$, $$ z= \begin{cases} \eta+\dfrac{1}{\mu}, & Y=1,\\[6pt] \eta-\dfrac{1}{1-\mu}, & Y=0. \end{cases} $$ When $\mu$ is close to zero, the event correction $1/\mu$ is large and positive. The non-event correction is close to $-1$, not equally extreme. Conversely, when $\mu$ is close to one, a non-event can have a large negative correction. Thus it is an outcome that is surprising relative to the current fitted mean, combined with a small derivative, that produces a large pseudo-response.

A large $z$ alone does not show that the observation has a large working weight. In fact, for a rare event the weight $w=\mu(1-\mu)$ is small. In a weighted least-squares representation, influence also depends on the design and the leverage of the observation, and in a mixed model it depends on the cluster structure and random-effect penalty. The pseudo-response and weight must be inspected together.

Worked Example: Rare Events at Two Fitted Probabilities

Consider two hypothetical individuals, each with $Y=1$. For the first, let $\mu=0.01$; for the second, let $\mu=0.05$. Calculate their current predictors, weights, and pseudo-responses, then compare the raw response correction with its weighted contribution.

For $\mu=0.01$, $\eta=\log\{0.01/0.99\}=-4.5951$, and $w=0.01(0.99)=0.009900$. Thus $$ z=-4.5951+\frac{1-0.01}{0.009900} =-4.5951+100 =95.4049. $$ The pseudo-response is very large, while the working weight is only 0.009900.

For $\mu=0.05$, $\eta=\log\{0.05/0.95\}=-2.9444$, and $w=0.05(0.95)=0.04750$. Therefore $$ z=-2.9444+\frac{1-0.05}{0.04750} =-2.9444+20 =17.0556. $$ This pseudo-response is also large, but smaller than 95.4049, and its weight is larger.

The weighted correction from the current predictor is $w(z-\eta)=Y-\mu$. It is $0.009900(100)=0.9900$ for the first individual and $0.04750(20)=0.9500$ for the second. The first raw pseudo-response is much larger, but its weighted correction is not correspondingly larger. This identity is useful for detecting arithmetic errors and for avoiding the mistaken conclusion that the largest $z$ automatically dominates the update.

From Individual Quantities to a Cluster Update

For a random-intercept model, consider one PQL working-model update with the fixed-effect predictor held at a value $\mathbf{x}_{ij}^{\mathsf T}\boldsymbol{\beta}$. The penalized weighted least-squares criterion for a single cluster’s intercept contains $$ \sum_i w_{ij}\bigl(z_{ij}-\mathbf{x}_{ij}^{\mathsf T}\boldsymbol{\beta}-b_j\bigr)^2+\frac{b_j^2}{\tau^2}. $$ Differentiating with respect to $b_j$, setting the derivative to zero, and solving gives $$ b_j^{\mathrm{new}}= \frac{\sum_i w_{ij}\bigl(z_{ij}-\mathbf{x}_{ij}^{\mathsf T}\boldsymbol{\beta}\bigr)} {\sum_i w_{ij}+1/\tau^2}. $$ The denominator includes the Gaussian random-effect penalty. This conditional update illustrates why a cluster’s event count and its size matter: individual weighted residuals aggregate, while the penalty pulls the update toward zero. In a full iteration, fixed effects and variance components may also be updated; the calculation below holds them fixed to isolate this one step.

Worked Example: A Three-Person Cluster With One Event

Consider a hypothetical cluster of three individuals with no varying covariates. At the current iterate, set the fixed intercept to $\beta=-\log(19)=-2.9444$, the current random intercept to zero, and $\tau^2=0.25$. The outcomes are $1,0,0$. Thus all three current predictors are $-2.9444$, all fitted probabilities are 0.05, and all weights are 0.04750. The Bernoulli and conditional-independence assumptions hold by construction for this illustrative model.

The event has $z_1=-2.9444+20=17.0556$. Each non-event has $$ z_0=-2.9444-\frac{1}{0.95}=-3.9971. $$ For a cluster intercept update with $\beta$ held fixed, each weighted residual in the numerator is $w(z-\beta)=Y-\mu$. The event contributes $1-0.05=0.9500$, while each non-event contributes $0-0.05=-0.05000$. The numerator is therefore $$ 0.9500-0.05000-0.05000=0.8500. $$ The denominator is $$ 3(0.04750)+\frac{1}{0.25}=0.1425+4=4.1425. $$ Consequently, $$ b^{\mathrm{new}}=\frac{0.8500}{4.1425}\approx 0.2052. $$ The event generated the largest pseudo-response, but the cluster update is not 17.0556. It reflects all three weighted residuals and the random-intercept penalty. This is one conditional working-model step, not a final estimate of the cluster effect or a complete PQL fit.

Reading the Calculations and Diagnosing Problems

A useful check is to recompute each quantity from the same $\eta$: first $\mu=\{1+\exp(-\eta)\}^{-1}$, then $w=\mu(1-\mu)$, and finally $z=\eta+(Y-\mu)/w$. Rounding $\mu$ too early can noticeably alter $z$ when $w$ is small. Keep extra digits during the calculations and round the reported values only at the end.

The identity $w(z-\eta)=Y-\mu$ is an independent arithmetic check. For $Y=1$, its value is $1-\mu$; for $Y=0$, it is $-\mu$. If the displayed weight agrees with the probability but this identity fails materially, the pseudo-response or predictor has likely been calculated inconsistently.

The working weights are not probabilities, and a small weight does not mean an observation is irrelevant to every part of the fit. The weighted residual, the covariate pattern, and the cluster’s other observations all matter. A practical diagnostic calculation therefore records $Y$, $\eta$, $\mu$, $w$, and $z$ together, rather than reporting an extreme pseudo-response in isolation. If fitted probabilities are very close to zero or one, also recognize that the local linearization can be strained: a numerically correct working response does not guarantee a reliable approximation to the nonlinear model.

Comparison and Limitations

PQL replaces the nonlinear conditional mean locally with a Gaussian working mixed model. Its advantage for this calculation is that weighted mixed-model machinery can be used at each iteration. Its limitation is equally relevant: the pseudo-data depend on current fitted values, and extreme corrections can arise when the derivative is small. This tutorial’s arithmetic identifies those corrections; it does not establish that the resulting PQL estimate is unbiased or that the approximation is adequate for a particular dataset.

A direct likelihood-based GLMM analysis is a principal alternative because it targets the model likelihood rather than relying on this same sequence of first-order pseudo-responses. The two approaches can differ in sparse rare-event settings. Choosing between them requires attention to the data structure, the fitting method, and the scientific goal; a large pseudo-response by itself is not a decision rule. PQL’s working weights also should not be confused with inverse-probability weights or with final precision estimates.

Common Pitfalls

Treating $z$ as an outcome. A pseudo-response such as 95.4049 is a local linearization quantity, not an observed number of events or a probability.
Calling the largest pseudo-response the most influential observation. Its weight may be small, and influence also depends on design, cluster structure, and the penalty.
Mixing values from different calculations. Compute $\mu$, $w$, and $z$ from the same predictor and unrounded probability. Check $w(z-\eta)=Y-\mu$.
Interpreting a working update as a final fit. A conditional random-intercept update with other parameters fixed is only one step of the iterative procedure.
Assuming correct arithmetic validates the approximation. The calculations can be internally consistent even when the local linearization is inadequate for the data.

Key takeaway: Rare binary events can create extreme PQL pseudo-responses because the logit mean has a small derivative at a low fitted probability. The corresponding working weight is also small. Interpret both quantities, check their arithmetic together, and account for the cluster-level penalty before drawing conclusions about an update.

Check Your Understanding

Answer each question using the definitions and calculations above. Solutions follow the questions.

At $\mu=0.02$, calculate the logit predictor and weight. For $Y=1$, calculate the pseudo-response to three decimal places.
At $\mu=0.02$ and $Y=1$, what is $w(z-\eta)$? Explain what it represents.
At $\eta=-\log(19)$, calculate $z$ for $Y=0$, using $\mu=0.05$.
Why does a very large event pseudo-response not, by itself, establish that the observation has a very large working weight?
For the three-person cluster example, what would the random-intercept update numerator be if all three outcomes were zero, with all current predictors and weights unchanged?

Solutions. (1) $\eta=\log(0.02/0.98)=-3.8918$, and $w=0.02(0.98)=0.01960$. For an event, $z=\eta+1/\mu=-3.8918+50=46.108$. (2) $w(z-\eta)=Y-\mu=1-0.02=0.9800$. It is the weighted correction from the current predictor, not the pseudo-response itself. (3) $z=\eta-1/(1-\mu)=-2.9444-1/0.95=-3.9971$. (4) The event correction can be large because it divides by a small $w$, while that same $w$ is the observation’s working weight; actual influence also depends on design and the mixed-model structure. (5) Each of the three non-events contributes $Y-\mu=-0.05$, so the numerator is $-0.1500$.