Binary Logit PQL Equations In A Random-Intercept Model

From a Binary Logit Model to a Working Mixed Model

For clustered binary outcomes, the logit model is nonlinear in its mean: changing a cluster intercept changes every observation's probability through the inverse-logit function. Penalized quasi-likelihood (PQL) handles this by repeatedly constructing a local, weighted working response and fitting a penalized linear mixed model to it. This tutorial makes the equations explicit for one Gaussian random intercept per cluster.

Model and assumptions: For clusters \(i=1,\ldots,m\) and observations \(j=1,\ldots,n_i\), let \(y_{ij}\in\{0,1\}\). Conditional on the cluster random intercepts, outcomes are independent and \(y_{ij}\mid b_i\sim\operatorname{Bernoulli}(\mu_{ij})\), where \(\operatorname{logit}(\mu_{ij})=x_{ij}^{T}\beta+b_i\). The covariate vector \(x_{ij}\) is observed, \(\beta\) is the fixed-effect vector, and the cluster intercepts are independent with \(b_i\sim N(0,\tau^2)\), \(\tau^2>0\), and independent of the covariates. Thus observations in a cluster may be dependent marginally, through their shared \(b_i\).

Write \(\eta_{ij}=x_{ij}^{T}\beta+b_i\) and \(\mu_{ij}=\exp(\eta_{ij})/[1+\exp(\eta_{ij})]\). In a PQL iteration, use current values \(\eta_{ij}^{(0)}\), \(\mu_{ij}^{(0)}\), and \(\tau^2\). The superscript \(0\) marks the values at which the local working quantities are calculated; it does not mean the model is restricted to a single iteration.

Definition: The binary logit PQL pseudo-response and working weight at an observation are \(z_{ij}=\eta_{ij}^{(0)}+(y_{ij}-\mu_{ij}^{(0)})/[\mu_{ij}^{(0)}(1-\mu_{ij}^{(0)})]\) and \(w_{ij}=\mu_{ij}^{(0)}(1-\mu_{ij}^{(0)})\). The working linear mixed model treats \(z_{ij}\) as having mean \(x_{ij}^{T}\beta+b_i\) and working variance \(1/w_{ij}\). This is a local working approximation, not a claim that the original Bernoulli response is normally distributed.

These quantities follow from a first-order Taylor expansion of the link around the current mean. For the logit link \(g(\mu)=\log[\mu/(1-\mu)]\), the derivative is \(g'(\mu)=1/[\mu(1-\mu)]\). If the current predictor is \(\eta^{(0)}=g(\mu^{(0)})\), then \(g(y)\) is not a useful quantity for binary \(y\) (it is infinite at 0 and 1). Instead, the standard working-response construction linearizes the mean as a function of the predictor. Since \(d\mu/d\eta=\mu(1-\mu)\), a local response-scale displacement \(y-\mu^{(0)}\) corresponds to predictor-scale displacement \((y-\mu^{(0)})/[\mu^{(0)}(1-\mu^{(0)})]\). Adding this displacement to \(\eta^{(0)}\) gives \(z\); the Bernoulli variance \(\mu^{(0)}(1-\mu^{(0)})\), divided by the squared derivative \( [d\mu/d\eta]^2\), gives working variance \(1/w\).

Worked Example: One Observation's Pseudo-Response and Weight

Suppose a hypothetical observation has outcome \(y=1\), and the current linear predictor is \(\eta^{(0)}=0.4\). The fitted probability is \(\mu^{(0)}=1/(1+\exp(-0.4))\approx0.5987\). Using the displayed rounded value,

w=0.5987(1-0.5987)=0.2403,\qquad z=0.4+\frac{1-0.5987}{0.2403} =0.4+\frac{0.4013}{0.2403} \approx2.070.

The corresponding working variance, using the same displayed rounded weight, is \(1/0.2403\approx4.161\). The observation therefore contributes a working response of about \(2.070\) with weight \(0.2403\); it does not contribute an ordinary response of 2.070 to a Bernoulli likelihood.

For comparison, at the same current predictor but with \(y=0\), using those displayed rounded values gives \(z=0.4+(0-0.5987)/0.2403\approx-2.091\), with the same weight \(0.2403\). The outcome changes the pseudo-response, whereas the current fitted probability determines the weight.

The Penalized Equations for Cluster Intercepts

Let \(X\) be the fixed-effect design matrix, with one row \(x_{ij}^{T}\) per observation, and let \(Z\) be the cluster-incidence matrix: row \((i,j)\) has a 1 in column \(i\) and zeros elsewhere. Collect the pseudo-responses into \(z\), set \(W=\operatorname{diag}(w_{ij})\), and write \(b=(b_1,\ldots,b_m)^{T}\). The working linear predictor is \(X\beta+Zb\).

For fixed current weights and a specified \(\tau^2\), the PQL working step minimizes the weighted residual sum of squares plus the Gaussian random-effect penalty. The factor \(1/2\) does not change the minimizer.

Q(\beta,b)= (z-X\beta-Zb)^{T}W(z-X\beta-Zb) +\frac{1}{\tau^2}b^{T}b.

The penalty is the negative-log-density contribution for independent \(N(0,\tau^2)\) intercepts, up to constants and a factor of 2. Differentiating \(Q\) with respect to \(\beta\) and \(b\), setting each gradient equal to zero, and rearranging gives the penalized normal equations:

\begin{pmatrix} X^{T}WX & X^{T}WZ\\ Z^{T}WX & Z^{T}WZ+\tau^{-2}I_m \end{pmatrix} \begin{pmatrix}\widehat\beta\\ \widehat b\end{pmatrix} = \begin{pmatrix}X^{T}Wz\\ Z^{T}Wz\end{pmatrix}.

In detail, the fixed-effect gradient is \(-2X^{T}W(z-X\beta-Zb)\), while the random-effect gradient is \(-2Z^{T}W(z-X\beta-Zb)+2b/\tau^2\). Their zero conditions yield the two block rows above. The term \(\tau^{-2}I_m\) is essential: it shrinks each cluster intercept toward zero, with stronger shrinkage when the assumed random-intercept variance is smaller.

Because the random effect is a random intercept, the equations can also be written cluster by cluster. The fixed-effect equation is \(\sum_i\sum_j w_{ij}x_{ij}(z_{ij}-x_{ij}^{T}\beta-b_i)=0\). The equation for cluster \(i\) is

\sum_{j=1}^{n_i}w_{ij} \left(z_{ij}-x_{ij}^{T}\beta-b_i\right)-\frac{b_i}{\tau^2}=0, \qquad \widehat b_i= \frac{\sum_j w_{ij}(z_{ij}-x_{ij}^{T}\beta)} {\sum_j w_{ij}+\tau^{-2}}.

The second expression is conditional on the current or fitted \(\beta\), \(w_{ij}\), and \(\tau^2\). Its denominator combines the cluster's total working information, \(\sum_jw_{ij}\), with the penalty precision, \(\tau^{-2}\). It is not a standalone estimate of \(\tau^2\).

Worked Example: Solving the Joint Equations for Two Clusters

Consider hypothetical binary data with one observation in each of two clusters. At the current iteration, take \(\beta=0\), \(b_1=b_2=0\), and \(\tau^2=1\). Both current predictors are 0, so both fitted probabilities are 0.5 and both weights are \(0.5(1-0.5)=0.25\). Let \(y_1=1\) and \(y_2=0\). Their pseudo-responses are \(z_1=0+(1-0.5)/0.25=2\) and \(z_2=0+(0-0.5)/0.25=-2\).

The working objective for a common fixed intercept \(\beta\) and two cluster intercepts is \(Q=0.25(2-\beta-b_1)^2+0.25(-2-\beta-b_2)^2+b_1^2+b_2^2\). Differentiating with respect to each parameter gives

2\beta+b_1+b_2=0,\qquad 0.25(2-\beta-b_1)-b_1=0,\qquad 0.25(-2-\beta-b_2)-b_2=0.

The second equation gives \(b_1=0.4-0.2\beta\), and the third gives \(b_2=-0.4-0.2\beta\). Substituting their sum, \(b_1+b_2=-0.4\beta\), into the first equation gives \(2\beta-0.4\beta=0\), hence \(\widehat\beta=0\). Therefore \(\widehat b_1=0.400\) and \(\widehat b_2=-0.400\) for this working step. The opposite outcomes are balanced by opposite penalized cluster intercepts; the fixed intercept remains zero.

Reading the Cluster-Level Update

The cluster-specific equation makes the shrinkage mechanism visible without requiring a full matrix solve. If \(\beta\) is held fixed for a working update, each cluster's intercept is a weighted residual average, shrunk toward zero by the penalty precision. The amount of shrinkage depends both on \(\tau^2\) and on the working weights and observation count in that cluster.

Worked Example: Updating One Random Intercept

Suppose a hypothetical cluster contains two Bernoulli observations, both equal to 1. At the start of this step, let their current predictors both be 0, hold the fixed intercept at \(\beta=0\), and set \(\tau^2=1\). Each fitted probability is 0.5, each weight is 0.25, and each pseudo-response is 2. The cluster equation gives

\widehat b= \frac{0.25(2-0)+0.25(2-0)}{0.25+0.25+1} =\frac{1}{1.5} =0.667.

This is the minimizer for the cluster intercept while \(\beta\), the working quantities, and \(\tau^2\) are held fixed. To check it directly, the derivative condition is \(0.25(2-\widehat b)+0.25(2-\widehat b)-\widehat b=0\), which becomes \(1-1.5\widehat b=0\). In the working model, the fitted linear predictor for each observation is therefore about 0.667. This calculation is not a new Bernoulli probability calculation at the original iteration; probabilities and working quantities are refreshed in subsequent PQL iterations.

Iteration, Diagnostics, and Alternatives

A PQL iteration alternates between calculating the current working responses and weights and solving the resulting penalized working mixed model. In practice, that mixed-model fit may update the fixed effects, random-effect predictions, and variance component. The displayed normal equations specify the coefficient update conditional on the weights and \(\tau^2\); variance-component estimation requires the mixed-model fitting procedure and its chosen likelihood convention, such as ML or REML. After an update, recompute \(\eta_{ij}\), \(\mu_{ij}\), \(z_{ij}\), and \(w_{ij}\), then continue until the specified convergence criteria are met.

Initialize.
Choose starting values for \(\beta\), \(b_i\), and \(\tau^2\), and calculate each current predictor and fitted probability.

Construct the working data.
For every observation, calculate \(z_{ij}\) and \(w_{ij}\) using the binary logit formulas. Check that \(0<\mu_{ij}<1\) and that the calculations use the same iteration's current values.

Fit the penalized working model.
Solve the block equations for fixed effects and cluster intercepts, while the mixed-model fitting procedure handles the variance component.

Recompute and assess.
Refresh probabilities and working quantities, then assess changes in parameters and the fitting criterion. Stop only when the chosen convergence tolerances are satisfied.

Several checks help locate a faulty or unstable calculation. First, verify \(w_{ij}=\mu_{ij}(1-\mu_{ij})\), which is positive and at most 0.25; a larger value signals an arithmetic or coding error. Verify that the pseudo-response uses the current predictor and the matching current weight. Compare the cluster equation's residual contribution with its penalty term: at a solution, their difference is zero. Finally, check that fitted parameters and the objective stabilize rather than oscillate or diverge. Very small weights can produce large working variances and extreme pseudo-responses, so inspect them rather than silently treating every working observation as equally informative.

The calculation has clear boundaries. A working normal approximation is used to turn the nonlinear binary model into a sequence of linear mixed-model problems; PQL is therefore not the same as maximizing the exact integrated Bernoulli likelihood. Its result can depend on the quality of the approximation and on whether the iterations converge. The equations above do not by themselves establish that a particular fit is adequate, nor do they supply a universal variance-component update. They also depend on the stated random-intercept structure; random slopes or a different link require corresponding changes to the model and working quantities.

An ordinary logistic regression is a simpler alternative when the data structure does not call for a cluster random effect, but it does not represent the shared cluster intercept in this model. Exact marginal-likelihood methods, including numerical integration approaches, target the likelihood after integrating over the random effects rather than using the PQL working response. These alternatives can be more computationally demanding. The choice is a modeling and computational decision: the PQL equations are useful to understand and implement the working approximation, but they should not be mistaken for exact likelihood equations.

Common Mistakes

Using the response-scale variance as the working weight. The Bernoulli variance is \(\mu(1-\mu)\), and in this logit working model it is the weight; the corresponding working variance is its reciprocal.
Forgetting the random-effect penalty. Omitting \(b_i/\tau^2\) from the cluster equation eliminates shrinkage and no longer gives the Gaussian random-intercept working fit.
Updating the pseudo-response but not the weight. Both are functions of the current fitted probability and must be recomputed together at each iteration.
Calling a working response an observed outcome. The pseudo-response can be outside \([0,1]\); it is an algebraic device for a local weighted fit, not a predicted probability or binary measurement.
Treating the conditional equations as a variance estimate. They give fixed-effect and random-intercept equations conditional on \(\tau^2\); estimating \(\tau^2\) is a separate part of the working mixed-model fit.

Key takeaway: In binary logit PQL, each current probability determines a working weight \(w_{ij}=\mu_{ij}(1-\mu_{ij})\) and pseudo-response \(z_{ij}=\eta_{ij}^{(0)}+(y_{ij}-\mu_{ij})/w_{ij}\). The weighted residual equations are coupled across clusters through the fixed effects, while each Gaussian random intercept receives the precision penalty \(1/\tau^2\).

Check Your Understanding

Answer each question, then compare with the brief solutions below.

If a current fitted probability is 0.2, what is the binary logit PQL working weight and working variance?
For \(y=0\) and current predictor \(\eta^{(0)}=0\), calculate the pseudo-response and weight.
What term in the cluster equation represents Gaussian shrinkage, and how does it change if \(\tau^2\) decreases?
Why can the pseudo-response be negative without contradicting the binary nature of the outcome?
Do the penalized normal equations alone estimate \(\tau^2\)? Explain.

Solutions. (1) \(w=0.2(0.8)=0.16\), so the working variance is \(1/0.16=6.25\). (2) At \(\eta^{(0)}=0\), \(\mu^{(0)}=0.5\) and \(w=0.25\), giving \(z=0+(0-0.5)/0.25=-2\). (3) The term is \(b_i/\tau^2\); decreasing \(\tau^2\) increases the penalty precision and, other quantities fixed, strengthens shrinkage toward zero. (4) The pseudo-response is a local working quantity on the predictor scale, not an observed Bernoulli value. (5) No. They give coefficient equations conditional on the working weights and \(\tau^2\); a mixed-model variance-component procedure must estimate \(\tau^2\).