Rare Events And PQL Bias

Why Rare Events Are a Special Problem for PQL

Penalized quasi-likelihood (PQL) turns a generalized linear mixed model into a sequence of weighted linear mixed-model calculations. For binary data, each calculation uses a local linear approximation to the conditional mean. When the fitted event probability is small, the inverse slope of that mean function is large. A single observed event can therefore correspond to a very extreme working response, even though its working weight is small.

That combination matters. The pseudo-response is not a new observed outcome, and its magnitude alone does not determine its influence. PQL uses the response and weight together inside a penalized Gaussian working-model calculation. With rare events, this local Gaussian representation can be a fragile stand-in for the original Bernoulli model. It can contribute to bias in fixed-effect and random-effect estimates, but the direction and size of bias depend on the model and data; rarity alone does not determine them.

Here \(x_{ij}\) is the covariate vector, \(\beta\) is the fixed-effect coefficient vector, and \(\tau^2\) is the random-intercept variance. The probabilities and working quantities are conditional on the current random-effect and fixed-effect values. The assumptions include a correctly specified conditional Bernoulli mean and logit link, and the stated random-effects distribution. PQL adds an approximation: it treats the locally linearized responses as a Gaussian working outcome to estimate the model. It does not calculate the exact marginal Bernoulli likelihood.

Where the Extreme Response Comes From

For the logit link, \(dp_{ij}/d\eta_{ij}=p_{ij}(1-p_{ij})\). A first-order Taylor expansion of the conditional mean around the current linear predictor gives

Solving this local relation for the linear predictor that would correspond to the observed response \(Y_{ij}\) gives the working response. The Bernoulli variance is \(p_{ij}(1-p_{ij})\), so the inverse-variance weight for this linearized response is \(w_{ij}=p_{ij}(1-p_{ij})\). In particular, the response’s departure from the current predictor is

As \(p_{ij}\) approaches zero, a non-event produces a response shift close to \(-1\), while an event produces a shift of \(1/p_{ij}\), which can be enormous. At the same time, the weight \(p_{ij}(1-p_{ij})\) approaches zero. The calculation therefore does not say that a rare event receives an enormous weight. It says that the event creates an extreme pseudo-response with a small weight.

Why the Combination Can Affect Estimation

In a weighted least-squares calculation, a response’s contribution is understood through its weighted residual, not its raw size. For one observation, the weighted squared residual from the current predictor is \(w(z-\eta)^2\). Substituting the binary PQL quantities gives

This quantity is small for a non-event when \(p\) is small, but large for an event. Equivalently, the event has a score contribution \(Y-p\) close to one while the local information \(p(1-p)\) is small. A rare observed event can thus be highly surprising relative to the model’s local information. Its extreme working response should not be read as either an automatically dominant observation or a harmless one: how it affects an estimate depends on the full weighted, penalized fit.

PQL’s bias is connected to using a local Gaussian approximation to represent a nonlinear Bernoulli model. The iteration updates the linearization as fitted values change, but the resulting penalized working-model calculation is still an approximation to the original model. Rare-event observations make the working quantities particularly asymmetric between \(Y=0\) and \(Y=1\). When random effects are estimated, this approximation also affects how the data inform the random-effect distribution. It can distort fixed effects, random-effect predictions, or the estimated random-effect variance. No single bias direction follows from the formulas alone.

A Numerical Check on Event Surprise

The weighted squared residual makes the asymmetry concrete. At \(p=0.02\), the weight is \(0.02(0.98)=0.0196\). If \(Y=1\), then \(z-\eta=1/0.02=50\), so the weighted squared residual is

If \(Y=0\), then \(z-\eta=-1/0.98\approx -1.0204\), giving

This large contrast reflects the surprise of seeing an event under a two-percent conditional probability. It is not, by itself, a bias calculation: bias is the difference between an estimator’s expected value and its target across repeated samples. A single weighted residual cannot establish the sign or amount of that difference.

Diagnostics, Alternatives, and Limits

When event probabilities are low, inspect the quantities PQL actually uses rather than judging the fit only from the scale of its pseudo-responses. Very large \(z_{ij}\) values can be expected under rare events, but they should prompt a check that the fitted probabilities and calculations are valid. Probabilities numerically equal to zero or one make the formulas undefined or numerically unstable. Also examine whether a small number of observed events carries most of the information about particular coefficients or clusters, and whether estimates change materially under reasonable model or estimation choices.

A useful comparison is with likelihood-based estimation that integrates over random effects, using numerical methods such as Laplace approximation or adaptive quadrature. These methods target the marginal likelihood rather than the PQL Gaussian working model, though they too rely on computational approximations and can be difficult to fit in challenging data. Agreement across methods is reassuring, but disagreement is evidence to investigate—not proof that one method is automatically correct. PQL may remain computationally convenient, yet rare-event settings reduce the comfort one should take from its working-Gaussian approximation.

The formulas here concern the conditional probabilities used in an iteration. A low marginal event rate does not, by itself, imply that every conditional probability is low; random effects can make probabilities vary across clusters. Nor do these calculations establish that rare events always produce bias in the same direction, or that a large working response alone invalidates a fit. The size of the approximation error also depends on the covariate pattern, number and arrangement of events, and random-effect structure.

A common mistake is to reason, “The event’s working response is huge, so PQL must give it a huge weight and pull the estimate sharply upward.” The correction is that its weight is \(p(1-p)\), which is small when \(p\) is rare. The event can still be consequential because it is surprising relative to the model’s low local information, but its effect must be assessed within the full penalized fit. Another mistake is to claim that rare events necessarily make every PQL estimate too small; the formulas show an asymmetry and an approximation, not a universal bias direction.