Tutorials › Generalized Linear Mixed Models (GLMM) › Worked Example Of Binary Random-Intercept PQL

Home  /  Tutorials  /  Penalized Quasi-Likelihood and Its Limitations

Penalized Quasi-Likelihood and Its Limitations · Tutorial 227 of 698

Worked Example Of Binary Random-Intercept PQL

Follow one binary random-intercept PQL step from the current linear predictors through the working data to updated mixed-model estimates.

Advanced 12 min read

What You'll Learn

  • Derive the binary-logit PQL pseudo-response by linearizing the mean, not the observed binary outcome.
  • Calculate fitted probabilities, working weights and pseudo-responses from current linear predictors.
  • Solve a small penalized working mixed-model problem for a fixed intercept and cluster random intercepts.
  • Interpret the updated working estimates without treating them as final logistic-model estimates.
  • Recognize numerical and modeling conditions that can make a PQL step unstable or misleading.

What One Binary PQL Iteration Computes

Penalized quasi-likelihood (PQL) turns a nonlinear mixed-model fitting step into a working linear mixed-model calculation. For binary data, that calculation is especially concrete: current fitted probabilities determine a working weight and pseudo-response for every observation, and the working mixed model updates the fixed effects and random effects. This tutorial carries out one such step for a small random-intercept dataset.

The result is an iteration update, not a final fitted logistic model. In particular, the pseudo-responses are constructed from the current parameter values, and the updated working estimates may change the probabilities used in the next iteration.

Model and assumptions: For observation \(j\) in cluster \(i\), let \(Y_{ij}\) be conditionally Bernoulli given the cluster random effect, with \(P(Y_{ij}=1\mid b_i)=p_{ij}\) and \(\operatorname{logit}(p_{ij})=\eta_{ij}=x_{ij}^{T}\beta+b_i\). The random intercepts are independent with \(b_i\sim N(0,\tau^2)\). Outcomes are conditionally independent given the random effects, and the Bernoulli variance is \(p_{ij}(1-p_{ij})\). The working step below uses current values of \(\beta\), \(b_i\) and \(\tau^2\); it treats the resulting working responses as a Gaussian mixed model for the update.

Constructing the Working Responses

Write the current linear predictor and probability as \(\eta_{ij}^{(r)}\) and \(p_{ij}^{(r)}\). For the logit link, \(\eta=\operatorname{logit}(p)\), and the logistic mean as a function of the linear predictor is \(p(\eta)=1/(1+\exp(-\eta))\). Its derivative is \(p(\eta)\{1-p(\eta)\}\). Linearizing this mean function at the current linear predictor gives the local approximation needed for a working response.

$$ p(\eta)\approx p_{ij}^{(r)} +p_{ij}^{(r)}\{1-p_{ij}^{(r)}\} \bigl(\eta-\eta_{ij}^{(r)}\bigr). $$

The observed binary response \(Y_{ij}\) is equated to this local mean approximation plus working error. Solving the resulting equation for the response on the linear-predictor scale gives the pseudo-response \(z_{ij}^{(r)}\). This derivation does not apply the logit function to \(Y_{ij}\): for binary observations, that would involve \(\operatorname{logit}(0)\) or \(\operatorname{logit}(1)\), which are infinite.

$$ z_{ij}^{(r)} =\eta_{ij}^{(r)} +\frac{Y_{ij}-p_{ij}^{(r)}} {p_{ij}^{(r)}\{1-p_{ij}^{(r)}\}}, \qquad w_{ij}^{(r)}=p_{ij}^{(r)}\{1-p_{ij}^{(r)}\}. $$

The working response has approximate conditional mean \(\eta_{ij}\) near the current value. Its working residual variance is approximately \(1/w_{ij}^{(r)}\): the Bernoulli variance at the current probability is \(w_{ij}^{(r)}\), and dividing the response residual by \(w_{ij}^{(r)}\) gives variance \(w_{ij}^{(r)}/\{w_{ij}^{(r)}\}^2=1/w_{ij}^{(r)}\). Thus observations with larger weights have greater precision in the working model.

Working-model update: Let \(z^{(r)}\) collect the pseudo-responses, \(W^{(r)}\) be diagonal with entries \(w_{ij}^{(r)}\), and \(Z\) be the random-intercept design matrix. For fixed current \(\tau^2\), the working fixed effects and random effects minimize \((z^{(r)}-X\beta-Zb)^T W^{(r)}(z^{(r)}-X\beta-Zb)+b^T G^{-1}b\), where \(G=\tau^2 I\). Equivalently, the working residual covariance is \(W^{-1}\) and the random-effect covariance is \(G\). A separate variance-component update, if performed by the fitting algorithm, is not implied by this minimization alone.

Worked Example: Updating Three Clusters

Worked Example: Updating Three Clusters

Consider three hypothetical clusters, each with two binary observations. The outcomes are \((1,1)\), \((1,0)\), and \((0,0)\), respectively. Fit an intercept-only fixed effect and a random intercept for each cluster. For this single iteration, take the current values to be \(\beta^{(r)}=0\), \(b_i^{(r)}=0\) for all clusters, and \(\tau^2=0.25\). These current values give \(\eta_{ij}^{(r)}=0\) and \(p_{ij}^{(r)}=0.5\) for every observation.

Working responses and weights. For every observation, \(w_{ij}^{(r)}=0.5(1-0.5)=0.25\). For an observed 1, \(z=0+(1-0.5)/0.25=2\). For an observed 0, \(z=0+(0-0.5)/0.25=-2\). The working-response pairs are therefore \((2,2)\), \((2,-2)\), and \((-2,-2)\). The working residual variance for each observation is \(1/0.25=4\).

Update the fixed intercept. In working-model notation, the residual covariance within a two-observation cluster is \(V_i=4I_2+0.25\mathbf{1}\mathbf{1}^T\), where \(\mathbf{1}\) is a two-vector of ones. Since \(V_i\mathbf{1}=4.5\mathbf{1}\), it follows that \(V_i^{-1}\mathbf{1}=\mathbf{1}/4.5\). The generalized least-squares intercept is

$$ \hat\beta =\frac{\sum_i \mathbf{1}^T V_i^{-1}z_i} {\sum_i \mathbf{1}^T V_i^{-1}\mathbf{1}} =\frac{(4+0-4)/4.5}{3(2/4.5)} =0. $$

Update the cluster intercepts. Given this fixed-effect estimate, the working-model random-effect estimate for cluster \(i\) is \(\hat b_i=\tau^2\mathbf{1}^TV_i^{-1}(z_i-\hat\beta\mathbf{1})\). For the first cluster, the pseudo-response sum is 4, so \(\hat b_1=0.25(4/4.5)=0.2222\). For the second cluster, the sum is 0 and \(\hat b_2=0\). For the third cluster, the sum is \(-4\), giving \(\hat b_3=0.25(-4/4.5)=-0.2222\). The updated working linear predictors are consequently \(0.2222\), \(0\), and \(-0.2222\), one for each cluster.

As a check on the interpretation, converting these updated linear predictors back to logistic probabilities gives approximately \(0.5553\), \(0.5000\), and \(0.4447\), respectively. These are probabilities implied by the updated working predictors; the example has completed one step, not established convergence or an exact likelihood-based fit. The variance \(\tau^2=0.25\) was held fixed for this calculation.

Why the Random Effects Are Shrunk

The working-model update balances fit to the pseudo-responses against the Gaussian random-effect penalty. For a cluster with two observations and equal working variance 4, the variance of its working-response mean is \(4/2=2\). If the fixed intercept is \(\beta\), its random-intercept estimate is therefore the working residual mean multiplied by the usual variance ratio:

$$ \hat b_i =\frac{\tau^2}{\tau^2+2}(\bar z_i-\beta). $$

Here \(\bar z_i\) is the cluster mean of the pseudo-responses. The multiplier is \(0.25/(0.25+2)=1/9\), so a cluster's working residual mean is not assigned wholly to its random effect. The random-effect penalty pulls the estimate toward zero. In the three-cluster example, the first cluster's residual mean is 2 and its estimate is \(2/9\), matching the matrix calculation.

Worked Example: Computing Pseudo-Responses at a Nonzero Probability

Suppose a current fitted probability is \(p^{(r)}=0.6\), so \(\eta^{(r)}=\log(0.6/0.4)=0.4054651\). For an observation with outcome \(Y=1\), the weight is \(0.6(0.4)=0.24\), and \(z^{(r)}=0.4054651+(1-0.6)/0.24=2.0721318\), or \(2.0721\) to four decimal places.

For an observation with outcome \(Y=0\) at the same current probability, the weight remains \(0.24\), but the response residual changes: \(z^{(r)}=0.4054651+(0-0.6)/0.24=-2.0945349\), or \(-2.0945\) to four decimal places. The two observations have the same working precision but different pseudo-responses. Neither pseudo-response is a possible binary outcome; each is a local linearization device.

Worked Example: A Conditional Random-Intercept Update

Use the same current probability \(0.6\) for two observations in one hypothetical cluster, with outcomes \((1,0)\). Let the fixed intercept be held at \(\beta=0.4054651\) for this conditional random-effect calculation, and let \(\tau^2=0.25\). From the preceding calculations, the pseudo-responses are \(2.0721318\) and \(-2.0945349\); their mean is \(-0.0112016\). The working-response mean has residual variance \(1/(2\times0.24)=2.0833333\).

The conditional random-intercept update is the working residual mean multiplied by the random-effect variance ratio: \(\hat b=\{0.25/(0.25+2.0833333)\}(-0.0112016-0.4054651)\). The ratio is \(3/28\), and the residual mean is \(-0.4166667\), so \(\hat b=(3/28)(-0.4166667)=-0.0446429\). This is a conditional update with \(\beta\) fixed, not a joint fixed- and random-effect update. It illustrates how both the working uncertainty and the random-effect penalty determine the amount of cluster adjustment.

Assumptions, Diagnostics, and Alternatives

The arithmetic of an iteration is straightforward only when the current probabilities are strictly between zero and one, so that the weights are positive and the pseudo-responses are finite. The working mixed-model calculation also depends on the stated design and covariance structure: outcomes are conditionally independent Bernoulli observations, the random intercepts are independent Gaussian effects, and the working residual variances are \(1/w_{ij}^{(r)}\). A calculation that uses different weights, a different random-effect structure, or an updated variance component must state that choice explicitly.

After computing a step, check that every weight equals \(p(1-p)\), every pseudo-response uses its own observed outcome, and the random-effect update balances working fit with the penalty. Very small weights produce large working residual variances and can make pseudo-responses highly influential. Large changes in predictors or random-effect estimates across iterations, or estimates that are sensitive to starting values, are reasons to examine the fit rather than assume the algorithm has produced a reliable approximation. These checks do not remove the known approximation limitations of binary PQL.

PQL is useful because it converts each local update into a familiar weighted mixed-model problem. Its main limitation is that the working linearization approximates the nonlinear binary model; it is not an exact integration over random effects. Likelihood-based integration is a distinct alternative that targets the marginal likelihood rather than repeatedly treating local pseudo-responses as Gaussian. The earlier material on PQL bias explains why a numerically completed iteration or convergence does not by itself establish that the resulting estimates have negligible approximation bias.

Common Pitfall

A tempting but invalid derivation is to Taylor-expand \(\operatorname{logit}(Y)\) around the current probability. For Bernoulli data, \(Y\) is 0 or 1, and the logit is infinite at both endpoints. The valid construction linearizes \(p(\eta)\), the mean as a function of the linear predictor, then solves the local working equation for \(z\). A second mistake is to report a working pseudo-response as a fitted probability or as a transformed observed outcome. It is neither: it is a temporary response on the linear-predictor scale, with precision \(w\), used to calculate an update.

Key takeaway: In a binary random-intercept PQL iteration, calculate \(p^{(r)}\), \(w^{(r)}=p^{(r)}(1-p^{(r)})\), and \(z^{(r)}=\eta^{(r)}+(Y-p^{(r)})/w^{(r)}\) observation by observation. Then solve the penalized working mixed model, keeping clear which parameter values were held fixed and that the update remains an approximation.

Check Your Understanding

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

  1. At current probability \(p^{(r)}=0.5\) and current predictor \(\eta^{(r)}=0\), calculate the working weight and pseudo-response when \(Y=0\).
  2. Why is a Taylor expansion of \(\operatorname{logit}(Y)\) not a valid derivation of the binary PQL pseudo-response?
  3. With two observations per cluster, working residual variance 4, and random-intercept variance \(0.25\), what fraction of a cluster's working residual mean is assigned to its random intercept?
  4. In the three-cluster example, why was the estimated fixed intercept zero even though two cluster random effects were nonzero?

Solutions. (1) \(w=0.5(1-0.5)=0.25\), and \(z=0+(0-0.5)/0.25=-2\). (2) A Bernoulli outcome is 0 or 1, where the logit is infinite; PQL instead linearizes the mean function \(p(\eta)\). (3) The cluster-mean working residual variance is \(4/2=2\), so the fraction is \(0.25/(0.25+2)=1/9\). (4) The three cluster pseudo-response means are \(2,0,-2\), with equal working-model precision; their grand mean is zero, so the generalized least-squares intercept is zero. The positive and negative cluster deviations are then partly represented by shrunken random intercepts.

Previous

Next

Back to Tutorials
Clinical Biostats