Iterative Reweighted Linear Mixed-Model Updates

From a Nonlinear Mixed Model to Repeated Linear Fits

Penalized quasi-likelihood (PQL) replaces a nonlinear mixed-model calculation with a sequence of weighted linear mixed-model fits. At each outer iteration, the current fitted means determine a pseudo-response and working weights; those working data define a linear mixed model that updates the fixed effects, random effects, and, when being estimated, their covariance parameters. The updated quantities then generate the next set of working data.

This is an iterative approximation, not an exact transformation of the original generalized linear mixed model. Its useful computational feature is that each local update uses familiar linear mixed-model machinery. Its important limitation is that the linear approximation changes as the fit changes, and convergence of the updates does not guarantee that the approximation is accurate.

Definition: For observations \(i=1,\ldots,n\), let \(y_i\) have conditional mean \(\mu_i\) and conditional variance \(\phi V(\mu_i)\), given the random-effects vector \(b\). Write \(\eta_i=x_i^\mathsf{T}\beta+z_i^\mathsf{T}b+o_i\), where \(o_i\) is a known offset, and let \(g(\mu_i)=\eta_i\), with inverse link \(h=g^{-1}\). The random effects are modeled as \(b\sim N(0,G(\theta))\). At a current value \(\eta_i^{(r)}\), define \(\mu_i^{(r)}=h(\eta_i^{(r)})\), \(d_i^{(r)}=h'(\eta_i^{(r)})\), and \(z_i^{*(r)}=\eta_i^{(r)}+(y_i-\mu_i^{(r)})/d_i^{(r)}\), \(w_i^{(r)}=(d_i^{(r)})^2/\{\phi V(\mu_i^{(r)})\}\). The PQL update fits these pseudo-responses using weights \(w_i^{(r)}\) and the random-effects covariance structure.

The definitions assume a differentiable inverse link and \(d_i^{(r)}\ne 0\), with a positive conditional variance at the current fitted mean. The displayed weights apply when there are no additional observation-level prior weights; if observation \(i\) has prior weight \(a_i\), its working weight is multiplied by \(a_i\). For binomial data recorded as counts of successes, both the mean and variance must be on the count scale. Equivalently, one can work with proportions using the corresponding trial-count weights.

Why the Working Data Have This Form

The pseudo-response comes from a first-order Taylor expansion of the conditional mean around the current linear predictor. Holding the current value fixed for the expansion,

\mu_i \approx \mu_i^{(r)}+d_i^{(r)}\bigl(\eta_i-\eta_i^{(r)}\bigr).

Rearranging this local relationship suggests a response on the linear-predictor scale:

z_i^{*(r)} =\eta_i^{(r)}+\frac{y_i-\mu_i^{(r)}}{d_i^{(r)}}.

The same linearization gives the working variance. Since the response has conditional variance \(\phi V(\mu_i^{(r)})\), dividing its local error by \(d_i^{(r)}\) gives

\operatorname{Var}\bigl(z_i^{*(r)}\mid b\bigr) \approx \frac{\phi V(\mu_i^{(r)})}{(d_i^{(r)})^2} =\frac{1}{w_i^{(r)}}.

Thus the local approximation treats the pseudo-responses as independent Gaussian observations with means \(x_i^\mathsf{T}\beta+z_i^\mathsf{T}b+o_i\) and variances \(1/w_i^{(r)}\), conditional on \(b\). “Gaussian” here describes the working linear mixed model; it does not assert that the original responses are Gaussian.

Key result: With \(W^{(r)}=\operatorname{diag}(w_1^{(r)},\ldots,w_n^{(r)})\), the conditional penalized quadratic for an outer update is, up to constants, \[ \frac{1}{2}\bigl(z^{*(r)}-X\beta-Zb-o\bigr)^\mathsf{T}W^{(r)} \bigl(z^{*(r)}-X\beta-Zb-o\bigr) +\frac{1}{2}b^\mathsf{T}G(\theta)^{-1}b. \] For fixed \(\theta\), minimizing this expression gives the penalized working fit in \(\beta\) and \(b\). Estimating \(\theta\) requires a linear mixed-model variance-component procedure; the quadratic penalty by itself is not a variance-component estimation criterion.

The penalty is the negative log-density contribution, apart from constants, for a Gaussian random effect with covariance \(G(\theta)\). It discourages random-effect values unsupported by the working data. In a full mixed-model fit, variance parameters are estimated using the chosen linear mixed-model method, such as restricted maximum likelihood or maximum likelihood, rather than by minimizing the conditional quadratic alone.

The Outer Iteration and Its Stopping Quantities

An outer iteration consists of recalculating working data and refitting the resulting linear mixed model. It is important to separate that cycle from the inner calculations used by a linear mixed-model fitting algorithm. The inner fit must itself converge sufficiently for its output to be a meaningful update.

Initialize.
Choose starting values \(\beta^{(0)}\), \(b^{(0)}\), and \(\theta^{(0)}\), and calculate \(\eta_i^{(0)}=x_i^\mathsf{T}\beta^{(0)}+z_i^\mathsf{T}b^{(0)}+o_i\). Starting values must give valid fitted means and finite, positive working weights.

Form the working data.
For each observation, calculate \(\mu_i^{(r)}=h(\eta_i^{(r)})\), its derivative \(d_i^{(r)}\), the pseudo-response \(z_i^{*(r)}\), and the working weight \(w_i^{(r)}\).

Fit the working mixed model.
Fit \(z^{*(r)}=X\beta+Zb+o+\varepsilon\), using residual covariance \((W^{(r)})^{-1}\) and random-effects covariance \(G(\theta)\). Update the fixed effects and random-effects predictions, and estimate \(\theta\) if it is not held fixed.

Recalculate and assess convergence.
Use the updated values to calculate the next linear predictors, fitted means, and working data. Stop only when relevant parameter and fitted-value changes are small and the inner mixed-model fit has converged; otherwise, continue the outer iteration.

No single stopping threshold is universal. A practical check compares successive fixed-effect estimates, random-effect predictions, covariance parameters, linear predictors, and fitted means. For example, a relative change for a parameter vector \(v\) can be summarized as \(\|v^{(r+1)}-v^{(r)}\|_\infty/\max(1,\|v^{(r)}\|_\infty)\). The scaling by at least one avoids an unstable relative change when the current value is near zero. A small change in the working criterion can be monitored too, but it should not be the only check: the criterion is built from iteration-specific working data, not a single exact likelihood that must decrease at every outer update.

Worked Example: Calculating One Binary-Response Working Observation

Suppose one hypothetical Bernoulli observation has \(y_i=1\). At the current iteration, its logit-scale linear predictor is \(\eta_i^{(r)}=0.4000\). For this calculation, the conditional mean is \(p_i=h(\eta_i)\), the inverse link is logistic, \(\phi=1\), and \(V(p_i)=p_i(1-p_i)\). Find its fitted probability, pseudo-response, and weight.

The fitted probability is \(p_i^{(r)}=1/\{1+\exp(-0.4000)\}=0.5987\) (to four significant figures). For the logistic inverse link, \(d_i^{(r)}=p_i^{(r)}(1-p_i^{(r)})\), so \(d_i^{(r)}=(0.5987)(0.4013)\approx0.2403\).

Substitution into the pseudo-response gives \(z_i^{*(r)}=0.4000+(1-0.5987)/0.2403 =0.4000+1.670\approx2.070\). The working weight is \(w_i^{(r)}=(0.2403)^2/\{(0.5987)(0.4013)\}\approx0.2403\). Thus, this observation contributes a pseudo-response of about 2.070 with weight about 0.2403 to this iteration’s linear mixed-model fit. These are working quantities, not a new observed outcome or a final fitted probability.

A Penalized Working Fit in Numbers

The next example isolates the mixed-model calculation after the working data have already been formed. Assume a Gaussian working model with two pseudo-responses, independent residual errors with known weights, and one random coefficient with known variance. These assumptions make it possible to solve the update directly without an inner variance-component estimation step.

Worked Example: Solving a Small Penalized Mixed-Model Update

Suppose the current working pseudo-responses are \(z^*=(1,2)^\mathsf{T}\), with \(W=\operatorname{diag}(2,1)\). Let the fixed-effect design matrix have a column of ones and let the random-effect design be \(Z=(0,1)^\mathsf{T}\). The random coefficient \(b\) has variance \(G=1\), and the offset is zero. Find the penalized working estimates of the fixed intercept \(\beta\) and \(b\).

The quadratic to minimize is \(2(1-\beta)^2+(2-\beta-b)^2+b^2\). Differentiating with respect to \(\beta\) and \(b\), and setting both derivatives to zero, gives the normal equations

\begin{pmatrix}3&1\\1&2\end{pmatrix} \begin{pmatrix}\beta\\b\end{pmatrix} = \begin{pmatrix}4\\2\end{pmatrix}.

The determinant is \(3(2)-1(1)=5\), so \(\beta=\{2(4)-1(2)\}/5=6/5=1.200\), and \(b=\{3(2)-1(4)\}/5=2/5=0.4000\). The fitted working means are therefore \(\beta=1.200\) for the first observation and \(\beta+b=1.600\) for the second. Substituting into the objective gives \(2(1-1.2)^2+(2-1.6)^2+(0.4)^2=0.0800+0.1600+0.1600=0.4000\). The penalty \(b^2=0.1600\) is part of the fit, not an additional observation.

This solves the conditional update for the specified known variance. If the variance were unknown, it would be estimated through the linear mixed-model procedure, and the resulting covariance estimate would feed into the next PQL iteration.

Checking Whether the Iteration Has Settled

A useful convergence assessment asks whether the algorithm has stabilized in quantities that matter for the fitted model, not only whether a software routine returned a convergence message. Consider a hypothetical run with the following changes from outer iteration 6 to iteration 7: the largest absolute fixed-effect change is 0.0006, the largest random-effect change is 0.0010, the largest absolute linear-predictor change is 0.0020, and the largest fitted-probability change is 0.0005. Suppose a variance parameter changes from 0.6400 to 0.6420.

Worked Example: Applying Multiple Stopping Checks

For the variance parameter, the relative change is \(|0.6420-0.6400|/\max(1,|0.6400|)=0.0020\). Suppose the stated stopping rules require fixed-effect and random-effect changes at most 0.0010, predictor changes at most 0.0030, fitted-probability changes at most 0.0010, and variance-parameter relative changes at most 0.0050.

The fixed-effect change, 0.0006, is below its threshold; the random-effect change, 0.0010, meets its threshold. The predictor change, 0.0020, and probability change, 0.0005, are also below their thresholds, as is the variance relative change, 0.0020. Under these specified rules, the outer iteration meets all five checks. This conclusion depends on the chosen thresholds and presumes that the inner mixed-model fit converged; it does not establish that PQL closely approximates the target generalized mixed-model estimates.

Diagnostics, Alternatives, and Limitations

Several warning signs can indicate that an update is unreliable or that the iteration is not behaving well:

Invalid or extreme working quantities. A fitted mean near a boundary can make the derivative small, producing a very large pseudo-response or an extremely small weight. Check for non-finite values and inspect the fitted means and weights.
Oscillation or slow movement. Large alternating changes in coefficients or fitted values suggest that the local linearizations are not settling. A small change in only one summary should not obscure movement in other quantities.
Singular or unstable mixed-model fits. Check whether the inner fitting step reports convergence, whether the covariance estimate is admissible, and whether the fixed-effect design and weighted design provide adequate information.
Convergence under only a weak check. Coefficients can change little while fitted values, variance parameters, or other random-effect components still move materially. Monitor several quantities on meaningful scales.

PQL’s distinctive approximation is to replace each nonlinear response contribution locally with a Gaussian working observation and then use a linear mixed-model fit. Direct likelihood methods instead target an integrated likelihood, using methods such as Laplace approximation or adaptive quadrature. Those methods can be more computationally demanding, but they avoid treating the sequence of working Gaussian models as the target likelihood itself. A generalized linear model without random effects also uses iterative working quantities, but it does not include the mixed-model random-effects covariance and penalty.

The main limitation is therefore conceptual as well as computational: a stable PQL solution is a stable solution to its iterative approximation. It is not, by convergence alone, a guarantee of accurate fixed-effect or random-effect estimates, nor an exact maximum-likelihood fit. Compare methods when the approximation may be consequential, and consider whether the fitted means, weights, and variance estimates are plausible for the application.

Common Mistakes

Using the observed response directly as the working response. The update uses \(z_i^{*(r)}\), which adjusts the current predictor by a scaled residual; it is generally not equal to \(y_i\).
Refitting without recalculating the weights. The weights depend on the current fitted means and must be recalculated at each outer iteration.
Ignoring the random-effect penalty. The working fit includes the covariance-based penalty; an ordinary weighted regression does not perform the same update.
Calling convergence proof of accuracy. Convergence only indicates that the selected iteration has stabilized to its stopping tolerances. It does not remove the approximation made by PQL.

Key takeaway: Each PQL outer step uses the current linear predictor to create a pseudo-response and weight, fits a penalized linear mixed model to those working data, and then recalculates the nonlinear means. Monitor both parameter and fitted-value changes, and keep convergence distinct from approximation accuracy.

Check Your Understanding

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

What information from the current iteration is needed to calculate a working response and its weight?
For the small penalized fit above, what would happen to the normal equations if the random-effect penalty were omitted?
Suppose fixed-effect estimates change by less than a chosen tolerance but the fitted probabilities do not. Is it sufficient to stop based on the fixed effects alone?
Why does a converged PQL iteration not establish that the exact integrated likelihood has been maximized?

Solutions. (1) Calculate the current predictor, its fitted mean, the inverse-link derivative, the conditional variance function, and any dispersion or prior-weight factors. These determine both \(z_i^*\) and \(w_i\). (2) The penalty contributes \(G^{-1}\) to the random-effect block of the normal equations; omitting it removes shrinkage and no longer gives the stated penalized mixed-model update. In this example, the lower-right entry would change from 2 to 1. (3) No. Fitted values are relevant stopping quantities, so a material change in probabilities means the fit has not met all the specified convergence checks. (4) PQL repeatedly fits locally linearized working models; it does not thereby turn their criterion into the exact likelihood integrated over the random effects. Stabilization of the approximation and accuracy relative to the exact target are different questions.