The Penalized Working Mixed-Model Criterion

Why the Working Criterion Has a Penalty

A binary PQL step replaces a nonlinear conditional mean model, locally, with a weighted linear mixed-model problem. The weights describe how much information each observation contributes at the current fit. The random-effect penalty represents the Gaussian distribution assumed for the random effects. Both parts are essential: omitting the weights discards the local information structure, while omitting the penalty treats random effects as unrestricted coefficients rather than as draws from a specified distribution.

Let $y_i\in\{0,1\}$ be conditionally independent binary responses, with conditional mean $\mu_i$ and linear predictor $\eta_i=\mathbf{x}_i^\top\boldsymbol{\beta}+\mathbf{z}_i^\top\mathbf{b}$. Here $\boldsymbol{\beta}$ is a vector of fixed effects, $\mathbf{b}$ is a vector of random effects, and $\mathbf{X}$ and $\mathbf{Z}$ are their design matrices. Assume a logit link, so $\operatorname{logit}(\mu_i)=\eta_i$, and assume $\mathbf{b}\sim N(\mathbf{0},\mathbf{G})$, with $\mathbf{G}$ positive definite. In a PQL working step, the current values of the parameters and random effects determine a working response and weight for each observation.

Definition: For a current linear predictor $\eta_i^{(0)}$, fitted mean $\mu_i^{(0)}$, and binary logit model, the working response and working weight are $$ z_i=\eta_i^{(0)}+\frac{y_i-\mu_i^{(0)}}{\mu_i^{(0)}(1-\mu_i^{(0)})}, \qquad w_i=\mu_i^{(0)}(1-\mu_i^{(0)}). $$ Let $\mathbf{z}$ be the vector of working responses and let $\mathbf{W}=\operatorname{diag}(w_1,\ldots,w_n)$. The penalized working mixed-model criterion is $$ Q(\boldsymbol{\beta},\mathbf{b}) =\frac{1}{2}\left[ (\mathbf{z}-\mathbf{X}\boldsymbol{\beta}-\mathbf{Z}\mathbf{b})^\top \mathbf{W} (\mathbf{z}-\mathbf{X}\boldsymbol{\beta}-\mathbf{Z}\mathbf{b}) +\mathbf{b}^\top\mathbf{G}^{-1}\mathbf{b} \right]. $$

The factor $1/2$ is conventional and does not change the minimizer. The first term measures weighted discrepancy between the working responses and their linear predictions. The second term is the negative log-density contribution of the Gaussian random effects, up to constants and the same factor of $1/2$. In this step, $\mathbf{W}$ and $\mathbf{G}$ are treated as fixed. The criterion is minimized jointly over $\boldsymbol{\beta}$ and $\mathbf{b}$.

From the Binary Model to the Working Responses

The working quantities are a local approximation, not newly observed outcomes. To see their origin, write $\mu_i=\mu(\eta_i)$ and expand the conditional log likelihood around a current predictor $\eta_i^{(0)}$. For a Bernoulli observation, the score with respect to $\eta_i$ at the current value is $(y_i-\mu_i^{(0)})\mu_i'^{(0)}/[\mu_i^{(0)}(1-\mu_i^{(0)})]$. The Fisher information is $(\mu_i'^{(0)})^2/[\mu_i^{(0)}(1-\mu_i^{(0)})]$. A Fisher-scoring step therefore has information weight $w_i=(\mu_i'^{(0)})^2/[\mu_i^{(0)}(1-\mu_i^{(0)})]$ and pseudo-response $z_i=\eta_i^{(0)}+(y_i-\mu_i^{(0)})/\mu_i'^{(0)}$. For the logit link, $\mu_i'^{(0)}=\mu_i^{(0)}(1-\mu_i^{(0)})$, giving the definition above.

This quadratic approximation says that, near the current predictor, the likelihood contributes approximately a weighted squared distance between $z_i$ and $\eta_i$. Adding the Gaussian random-effect contribution gives $Q$. Thus minimizing $Q$ is the mixed-model working problem for this local approximation; it is not an exact minimization of the original binary likelihood.

Worked Example: Computing One Binary Working Response and Weight

Suppose a hypothetical observation has $y_i=1$, and the current logit predictor is $\eta_i^{(0)}=0$. Find its fitted mean, working response, and weight.

The inverse logit at zero gives $\mu_i^{(0)}=1/(1+e^0)=0.5$. Therefore the logit working weight is $w_i=0.5(1-0.5)=0.25$.

The working response is $$ z_i=0+\frac{1-0.5}{0.5(1-0.5)} =\frac{0.5}{0.25}=2. $$

This observation contributes $0.25(z_i-\eta_i)^2/2$ to the unpenalized part of $Q$. The value $z_i=2$ is a local pseudo-response, not a possible binary outcome. The weight is positive because the fitted mean is strictly between zero and one.

Deriving the Coupled Normal Equations

Set $\mathbf{r}=\mathbf{z}-\mathbf{X}\boldsymbol{\beta}-\mathbf{Z}\mathbf{b}$. Since $\mathbf{W}$ is symmetric, differentiating the criterion with respect to the fixed effects gives $\partial Q/\partial\boldsymbol{\beta}=-\mathbf{X}^\top\mathbf{W}\mathbf{r}$. Differentiating with respect to the random effects gives $\partial Q/\partial\mathbf{b}=-\mathbf{Z}^\top\mathbf{W}\mathbf{r}+\mathbf{G}^{-1}\mathbf{b}$. At a joint minimum, both gradients are zero. Substituting the definition of $\mathbf{r}$ yields the two equations:

\begin{aligned} \mathbf{X}^\top\mathbf{W}\mathbf{X}\boldsymbol{\beta} +\mathbf{X}^\top\mathbf{W}\mathbf{Z}\mathbf{b} &=\mathbf{X}^\top\mathbf{W}\mathbf{z},\\ \mathbf{Z}^\top\mathbf{W}\mathbf{X}\boldsymbol{\beta} +\left(\mathbf{Z}^\top\mathbf{W}\mathbf{Z}+\mathbf{G}^{-1}\right)\mathbf{b} &=\mathbf{Z}^\top\mathbf{W}\mathbf{z}. \end{aligned}

The first equation balances the fixed-effect fit against the contribution of the random effects. The second balances the random-effect fit against both the fixed-effect contribution and the Gaussian penalty. Together they form a block system:

\begin{bmatrix} \mathbf{X}^\top\mathbf{W}\mathbf{X} & \mathbf{X}^\top\mathbf{W}\mathbf{Z}\\ \mathbf{Z}^\top\mathbf{W}\mathbf{X} & \mathbf{Z}^\top\mathbf{W}\mathbf{Z}+\mathbf{G}^{-1} \end{bmatrix} \begin{bmatrix} \boldsymbol{\beta}\\ \mathbf{b} \end{bmatrix} = \begin{bmatrix} \mathbf{X}^\top\mathbf{W}\mathbf{z}\\ \mathbf{Z}^\top\mathbf{W}\mathbf{z} \end{bmatrix}.

The cross-product terms show why the two sets of coefficients are not generally estimated separately. A random effect changes the fitted linear predictor and therefore affects the fixed-effect equation; the reverse is also true. If $\mathbf{X}$ has full column rank, all weights are positive, and $\mathbf{G}$ is positive definite, the criterion is strictly convex in $(\boldsymbol{\beta},\mathbf{b})$ and has a unique minimizer. Indeed, its quadratic curvature in a direction $(\mathbf{u},\mathbf{v})$ is proportional to $(\mathbf{X}\mathbf{u}+\mathbf{Z}\mathbf{v})^\top\mathbf{W}(\mathbf{X}\mathbf{u}+\mathbf{Z}\mathbf{v})+\mathbf{v}^\top\mathbf{G}^{-1}\mathbf{v}$. This can be zero only if $\mathbf{v}=\mathbf{0}$ and $\mathbf{X}\mathbf{u}=\mathbf{0}$, which, under full column rank, implies $\mathbf{u}=\mathbf{0}$.

Worked Example: Solving a Small Penalized Working System

Consider three hypothetical binary observations with current predictors all zero and outcomes $1,0,0$. Use a fixed intercept and one random-effect coefficient with design values $0,1,2$, so $\mathbf{X}=(1,1,1)^\top$ and $\mathbf{Z}=(0,1,2)^\top$. Let the random-effect variance be $1$, so $\mathbf{G}^{-1}=1$.

At predictor zero, each fitted mean is $0.5$ and each weight is $0.25$. The working responses are $2,-2,-2$, since $z_i=0+(y_i-0.5)/0.25$. Consequently, $\mathbf{X}^\top\mathbf{W}\mathbf{X}=0.75$, $\mathbf{X}^\top\mathbf{W}\mathbf{Z}=0.75$, and $\mathbf{Z}^\top\mathbf{W}\mathbf{Z}=0.25(0^2+1^2+2^2)=1.25$. The right-hand sides are $\mathbf{X}^\top\mathbf{W}\mathbf{z}=0.25(2-2-2)=-0.5$ and $\mathbf{Z}^\top\mathbf{W}\mathbf{z}=0.25(0\cdot2+1\cdot(-2)+2\cdot(-2))=-1.5$.

Adding the penalty $1$ to the random-effect diagonal gives $$ \begin{aligned} 0.75\beta+0.75b&=-0.5,\\ 0.75\beta+2.25b&=-1.5. \end{aligned} $$

Subtracting the first equation from the second gives $1.5b=-1$, hence $b=-2/3$. Substitution into the first gives $0.75\beta-0.5=-0.5$, so $\beta=0$. The fitted working predictors are $0,-2/3,-4/3$. These are the minimizers for this working step, not final estimates from a single binary likelihood fit.

What the Penalty Changes

The term $\mathbf{b}^\top\mathbf{G}^{-1}\mathbf{b}$ discourages random-effect values that are large relative to their assumed covariance. Its strength depends on $\mathbf{G}$: smaller random-effect variances produce stronger penalties in the corresponding directions. For a scalar random effect with variance $\tau^2$, the penalty is $b^2/\tau^2$. It is not a penalty on the fixed effects, and it is not an arbitrary numerical adjustment: it comes from the assumed Gaussian distribution.

Worked Example: Comparing Penalized and Unpenalized Fits

Use the working data and design matrices from the preceding example, but temporarily remove the random-effect penalty. The unpenalized equations have lower-right entry $1.25$, not $2.25$: $0.75\beta+0.75b=-0.5$ and $0.75\beta+1.25b=-1.5$.

Subtracting the first equation from the second gives $0.5b=-1$, so $b=-2$. Substitution gives $0.75\beta-1.5=-0.5$, hence $\beta=4/3$. With the Gaussian penalty included, the fit was instead $\beta=0$ and $b=-2/3$. The difference illustrates how the penalty reallocates the fit between fixed and random effects while shrinking the random-effect coefficient toward zero.

For the penalized solution, the weighted residual sum of squares is $0.25[2^2+(-2+2/3)^2+(-2+4/3)^2]=14/9$, and the penalty is $(2/3)^2=4/9$. Thus the criterion is $Q=\tfrac12(14/9+4/9)=1$. The unpenalized fit has a smaller weighted residual sum of squares, but that alone does not make it the appropriate mixed-model fit: it ignores the assumed distribution of the random effect.

Assumptions, Diagnostics, and Limitations

The criterion relies on a well-defined binary conditional model, conditional independence of observations given the random effects, and the specified Gaussian random-effect covariance. At the working step, it also relies on positive finite weights and on treating the current weights and covariance as fixed. In the logit case, finite predictors imply $0<\mu_i^{(0)}<1$ and therefore $w_i>0$. Numerically extreme predictors can nevertheless make weights very small, so calculations may be poorly conditioned even when the model is formally defined.

Useful checks include confirming that every working weight is positive, that the dimensions of $\mathbf{X}$, $\mathbf{Z}$, and $\mathbf{G}$ agree, and that the block system is solved stably. Large working responses can occur when a fitted mean is close to zero or one and the observed outcome points in the opposite direction. Such values can make a local quadratic approximation unreliable. Near-collinearity among design columns or very small weights can also make the fixed-effect information matrix difficult to solve. These are warnings to inspect the fit and approximation, not reasons to alter the criterion silently.

The working criterion differs from both unpenalized weighted least squares and the marginal likelihood. Unpenalized weighted least squares omits the Gaussian contribution and does not implement the assumed random-effect distribution. Marginal likelihood instead integrates over the random effects; minimizing $Q$ for a fixed set of working quantities does not perform that integration. PQL uses this local penalized problem as an approximation within a broader estimation procedure, so the criterion by itself does not establish that the approximation is accurate.

Common mistake: Treating $\mathbf{b}$ as an unrestricted nuisance coefficient and solving the equations without $\mathbf{G}^{-1}$ gives a different, unpenalized problem. In the mixed-model criterion, the random effects are constrained by their Gaussian distribution; the penalty must appear in the random-effect equation, added to $\mathbf{Z}^\top\mathbf{W}\mathbf{Z}$, and not to the fixed-effect block.

Check Your Understanding

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

For a binary logit observation with current fitted mean $0.25$, current predictor $\eta^{(0)}=\operatorname{logit}(0.25)$, and observed outcome $y=0$, find its working weight and working response.
Why does the random-effect normal equation contain $\mathbf{G}^{-1}\mathbf{b}$, while the fixed-effect equation has no corresponding penalty term?
In a scalar random-effect model, what happens to the penalty $b^2/\tau^2$ as $\tau^2$ becomes smaller, for a fixed nonzero $b$?
Name two conditions or numerical features that can make the working system unreliable or difficult to solve.

Solutions. (1) The weight is $0.25(0.75)=0.1875$. Since $y=0$, the working response is $\eta^{(0)}+(0-0.25)/0.1875=\operatorname{logit}(0.25)-4/3$, approximately $-2.432$. (2) The Gaussian random-effect density contributes the quadratic term $\mathbf{b}^\top\mathbf{G}^{-1}\mathbf{b}/2$; no prior penalty for $\boldsymbol{\beta}$ was assumed. (3) The penalty increases, because dividing a fixed positive $b^2$ by a smaller positive variance gives a larger value. (4) Examples include extremely small working weights, near-collinear design columns, an incorrectly specified or non-positive-definite covariance matrix, and predictors so extreme that the local quadratic approximation is suspect.