Working Responses For Binary PQL

Why Construct a Working Response?

A binary GLMM models outcomes on a probability scale but represents covariate and random-effect contributions on a linear-predictor scale. Penalized quasi-likelihood (PQL) connects these scales locally: at a current set of parameter values, it replaces the nonlinear conditional mean by its first-order Taylor approximation. The resulting working response, also called a pseudo-response, is expressed in linear-predictor units and can be used in a linear mixed-model calculation.

The key distinction is between the response actually calculated from one observed binary outcome and its conditional expectation under the model. The realized working response contains sampling noise; it is not simply an estimate equal to the underlying predictor. The earlier material on Taylor linearization of the binary mean supplies the local approximation used here.

Definition: For observation $i$, suppose $Y_i\mid b$ is Bernoulli with conditional mean $\mu_i=h(\eta_i)$, where $\eta_i=x_i^\mathsf{T}\beta+z_i^\mathsf{T}b$. At current values $\eta_i^{(0)}$, set $\mu_i^{(0)}=h(\eta_i^{(0)})$. If $h$ is differentiable and $h'(\eta_i^{(0)})\ne 0$, the binary PQL working response is $$ z_i^*=\eta_i^{(0)}+\frac{Y_i-\mu_i^{(0)}}{h'(\eta_i^{(0)})}. $$

Here $Y_i$ is the observed zero-or-one outcome, while $\eta_i$ and $\mu_i$ denote the predictor and conditional mean at the parameter values being considered. The superscript $(0)$ marks the current linearization point. In a GLMM, the predictor includes the random effects; this construction concerns the conditional mean given $b$, not a marginal mean obtained by integrating over the random-effects distribution.

Deriving the Pseudo-Response

Assume the binary observations are conditionally Bernoulli given $b$, with $P(Y_i=1\mid b)=\mu_i=h(\eta_i)$. At the current predictor, a first-order Taylor expansion gives

h(\eta_i)\approx h(\eta_i^{(0)})+ h'(\eta_i^{(0)})(\eta_i-\eta_i^{(0)}).

The observed outcome differs from its conditional mean by a random residual $Y_i-\mu_i$. Substituting the linearized mean into the decomposition $Y_i=\mu_i+(Y_i-\mu_i)$ gives the local representation

Y_i\approx \mu_i^{(0)}+ h'(\eta_i^{(0)})(\eta_i-\eta_i^{(0)}) +(Y_i-\mu_i).

Rearranging this approximation motivates a response on the predictor scale:

\eta_i^{(0)}+ \frac{Y_i-\mu_i^{(0)}}{h'(\eta_i^{(0)})} \approx \eta_i+\frac{Y_i-\mu_i}{h'(\eta_i^{(0)})}.

The left side is $z_i^*$. The right side shows why it is noisy: it is approximately the predictor plus a scaled conditional residual, not the predictor alone. Taking conditional expectations makes the approximation precise in the relevant sense. Since $E(Y_i\mid b)=h(\eta_i)$,

E(z_i^*\mid b) = \eta_i^{(0)}+ \frac{h(\eta_i)-h(\eta_i^{(0)})}{h'(\eta_i^{(0)})} \approx \eta_i.

The approximation follows from the same Taylor expansion. It is the conditional expectation of the working response that is locally approximated by the predictor; a particular realized $z_i^*$ need not be close to $\eta_i$.

Justification and remainder: If $h$ has a continuous second derivative between $\eta_i^{(0)}$ and $\eta_i$, Taylor's theorem gives $$ h(\eta_i)=h(\eta_i^{(0)})+ h'(\eta_i^{(0)})(\eta_i-\eta_i^{(0)}) +\frac{1}{2}h''(\xi_i)(\eta_i-\eta_i^{(0)})^2 $$ for some $\xi_i$ between those predictors. Therefore $$ E(z_i^*\mid b)-\eta_i = \frac{h''(\xi_i)}{2h'(\eta_i^{(0)})} (\eta_i-\eta_i^{(0)})^2. $$ Thus the conditional-mean discrepancy is second order in the predictor change when the derivative is nonzero and the curvature remains controlled.

The Logit Working Response

For a logit link, $h(\eta)=\exp(\eta)/\{1+\exp(\eta)\}$. Write $p_i^{(0)}=\mu_i^{(0)}$, the current fitted conditional probability. The derivative is $h'(\eta_i^{(0)})=p_i^{(0)}\{1-p_i^{(0)}\}$, so the binary logit working response is

z_i^*=\eta_i^{(0)} +\frac{Y_i-p_i^{(0)}}{p_i^{(0)}\{1-p_i^{(0)}\}}.

It has three inputs: the observed outcome $Y_i$, the current fitted probability $p_i^{(0)}$, and the current predictor $\eta_i^{(0)}$. The correction is positive when $Y_i=1$ exceeds the fitted probability and negative when $Y_i=0$ falls below it. Because the denominator is small when $p_i^{(0)}$ is near zero or one, the correction can be large. A working response can be negative or greater than one: it is measured in linear-predictor units, not probability units.

Worked Example: An Event at a Current Probability of 0.40

Suppose a hypothetical observation has $Y_i=1$ and current fitted probability $p_i^{(0)}=0.40$ under a logit model. Calculate its predictor and working response.

The current predictor is the logit of 0.40: $$ \eta_i^{(0)}=\log\left(\frac{0.40}{0.60}\right) =-0.4055. $$ The derivative at this point is $0.40(1-0.40)=0.24$. Substitution into the logit formula gives $$ z_i^*=-0.4055+\frac{1-0.40}{0.24} =-0.4055+2.500 =2.095. $$

This is a realized working response of about 2.095 predictor units, not a new probability and not a claim that the observation's underlying predictor is 2.095. For comparison, if the same current fit had $Y_i=0$, then $$ z_i^*=-0.4055+\frac{0-0.40}{0.24} =-0.4055-1.667 =-2.072. $$ The change comes from the observed outcome; both calculations use the same current fitted probability and predictor.

Worked Example: A Surprising Response When the Fitted Probability Is High

In a second hypothetical case, let $p_i^{(0)}=0.80$ and observe $Y_i=0$. Calculate the working response and explain its size.

The logit predictor is $$ \eta_i^{(0)}=\log\left(\frac{0.80}{0.20}\right)=1.3863, $$ and the derivative is $0.80(0.20)=0.16$. Therefore, $$ z_i^*=1.3863+\frac{0-0.80}{0.16} =1.3863-5.000 =-3.614. $$

The large negative value results from a zero outcome that is substantially below the fitted probability, divided by a derivative of only 0.16. This calculation is not an arithmetic error, and it does not mean the original outcome was $-3.614$. It is a local, linearized response on the predictor scale. The example also warns that a working response can be far from the predictor when the observed outcome is surprising under the current fit.

Worked Example: Checking the Conditional-Mean Approximation

Suppose the current logit predictor is $\eta_i^{(0)}=0$, so $p_i^{(0)}=0.50$. For illustration, consider a conditional predictor $\eta_i=0.80$, giving $\mu_i=\operatorname{expit}(0.80)\approx0.6900$. Find the conditional expectation of $z_i^*$ and compare it with $\eta_i$.

The derivative at the expansion point is $0.50(1-0.50)=0.25$. Thus $$ E(z_i^*\mid b) =0+\frac{0.6900-0.5000}{0.25} =0.7600. $$ The predictor is $0.8000$, so the conditional expectation of the pseudo-response is about $0.0400$ lower. It is close, but not equal: the change in predictor is 0.80, and the expansion is local rather than exact over that distance.

If the observed outcome in this setup is $Y_i=1$, its realized working response is instead $$ z_i^*=0+\frac{1-0.50}{0.25}=2.000. $$ That value is not the conditional expectation 0.7600. The expectation averages over the binary outcome under probability 0.6900; a single realization retains outcome noise.

Assumptions, Diagnostics, and Failure Modes

The construction requires a correctly specified conditional mean link $h(\eta_i)$, a differentiable link function at the current predictor, and a nonzero derivative there. For the logit case, finite current predictors imply $0<p_i^{(0)}<1$ and hence a positive derivative. The response is defined observation by observation, conditional on the random effects; in the usual GLMM formulation, outcomes are conditionally independent given those effects and covariates. Conditional independence is part of the model for the joint outcomes, although the algebra defining each response uses its own conditional mean.

Useful checks focus on whether the local approximation is credible and numerically stable:

Inspect extreme current probabilities. When $p_i^{(0)}$ is near zero or one, $p_i^{(0)}\{1-p_i^{(0)}\}$ is small. Even a modest difference $Y_i-p_i^{(0)}$ can then produce an extreme response.
Compare the current and updated predictor scales. Large predictor changes move the calculation away from its expansion point, where curvature can make the conditional-mean approximation poor.
Check whether extreme responses dominate the calculation. A few large values can be a warning of near-boundary fitted probabilities or outcomes poorly represented by the current fit; they are not, by themselves, proof of a coding error.
Verify the link and derivative together. Substituting a probability-scale derivative from one link while using a predictor from another produces an invalid response.

For logit, the conditional variance of the realized working response follows directly from the Bernoulli variance: $$ \operatorname{Var}(z_i^*\mid b) = \frac{\operatorname{Var}(Y_i\mid b)} {\{h'(\eta_i^{(0)})\}^2} = \frac{\mu_i(1-\mu_i)} {\{p_i^{(0)}(1-p_i^{(0)})\}^2}. $$ At the expansion point, where $\mu_i=p_i^{(0)}$, this reduces to $1/\{p_i^{(0)}(1-p_i^{(0)})\}$. This makes the noise particularly clear: near probabilities zero and one, the working response can be highly variable as well as sensitive to the observed outcome.

How This Approximation Fits Among Approaches

PQL uses the working response as one ingredient in an iterative approximation that treats the nonlinear conditional mean locally as linear. The response itself does not provide a complete fitting procedure: it does not specify how observations should contribute relative to one another, nor does it by itself determine fixed effects or random effects. Those require the rest of the penalized working mixed-model calculation.

Likelihood-based integration instead targets the likelihood after accounting for the random-effects distribution, using methods such as Laplace approximation or adaptive quadrature. These methods do not replace each binary outcome with this Taylor-derived pseudo-response. They can be computationally more demanding, but avoid relying on this particular response-scale linearization as their central approximation. Standard GLM iterative fitting is another nearby idea, but it does not by itself account for GLMM random effects. The choice among approaches therefore involves both the nonlinear mean and the random-effects structure, not just a preference for one response formula.

A common mistake is to write $Y_i-\mu_i^{(0)}\approx h'(\eta_i^{(0)})(\eta_i-\eta_i^{(0)})$, then rearrange as though the realized outcome residual were only a change in the predictor. That drops the random discrepancy between $Y_i$ and its conditional mean. The correction is to distinguish $Y_i-\mu_i$, the conditional sampling residual, from $\mu_i-\mu_i^{(0)}$, the mean change approximated by Taylor expansion. Consequently, $z_i^*$ is a noisy working response, while its conditional expectation is locally approximated by $\eta_i$.

The main limitation is local validity. A small derivative amplifies the outcome residual, and a large move from the expansion point can make curvature important. The formula is therefore best understood as a computational approximation at current values, not as a transformation that converts binary data into exact Gaussian observations.

Key takeaway: For a binary logit model, $z_i^*=\eta_i^{(0)}+\{Y_i-p_i^{(0)}\}/[p_i^{(0)}(1-p_i^{(0)})]$. It depends on the outcome, fitted probability, and current predictor; its conditional expectation, not each realized value, is locally approximated by the predictor.

Check Your Understanding

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

For a logit model with $p_i^{(0)}=0.25$, calculate $\eta_i^{(0)}$, the derivative, and $z_i^*$ when $Y_i=1$.
Explain why a realized working response is not generally equal to the conditional predictor.
At $\eta_i^{(0)}=0$, suppose the conditional mean at a target predictor is 0.60. Calculate $E(z_i^*\mid b)$ and compare it with the target predictor $\operatorname{logit}(0.60)$.
Name two conditions that make the working-response construction valid, and one warning sign that its local approximation may be poor.

Solutions. (1) $\eta_i^{(0)}=\log(0.25/0.75)=-1.0986$, and the derivative is $0.25(0.75)=0.1875$. Thus $z_i^*=-1.0986+(1-0.25)/0.1875=-1.0986+4=- -$ Wait, the arithmetic is $-1.0986+4=2.9014$. The working response is about 2.901 predictor units. (2) The response includes the scaled random residual $Y_i-\mu_i$; only its conditional expectation is locally approximated by the predictor. (3) The derivative at zero is $0.25$, so $E(z_i^*\mid b)=0+(0.60-0.50)/0.25=0.4000$. The target is $\log(0.60/0.40)=0.4055$, a close but not identical value. (4) The mean link must be differentiable and its derivative nonzero at the expansion point; the conditional binary mean must be modeled as $h(\eta_i)$. A current probability near zero or one, or a large predictor change from the expansion point, warns of instability or substantial linearization error.