Weighted-Residual Forms Of The Score Equations

From likelihood derivatives to residual constraints

A regression score is a vector of derivatives, one for each coefficient. Its zero is often easier to interpret when written in terms of response residuals \(y_i-\mu_i\). In a generalized linear model (GLM), each residual is multiplied by a factor that depends on the fitted mean, the variance function, the link, and the dispersion. The design matrix then combines these weighted residuals into the coefficient score.

This form makes the estimating equations interpretable: at an interior likelihood stationary point, each design column has a zero weighted sum against the residuals. The resulting constraints are not generally statements that every residual, or even their unweighted total, is zero. Which sums are constrained depends on the design matrix and the weights.

Deriving the weighted-residual score

Suppose \(Y_1,\ldots,Y_n\) are independent observations from an exponential-dispersion family with log-likelihood contributions

\ell_i(\theta_i,\phi) = \frac{y_i\theta_i-b(\theta_i)}{a(\phi)} +c(y_i,\phi).

Here \(\phi\) is treated as fixed while differentiating with respect to the regression coefficients, \(a(\phi)>0\), and the response support does not depend on the parameters. The mean is \(\mu_i=b'(\theta_i)\), and the variance has the form \(\operatorname{Var}(Y_i)=a(\phi)V(\mu_i)\). The GLM specifies \(\eta_i=x_i^\mathsf{T}\beta\) and \(\mu_i=g^{-1}(\eta_i)\), where \(x_i\) is a \(p\)-vector of known covariates and \(g^{-1}\) is differentiable at \(\eta_i\). Assume the parameter value is interior and the derivatives below exist.

Definition: The response residual for observation \(i\) is \(r_i=y_i-\mu_i\). Its score weight is \[ s_i=\frac{d\mu_i/d\eta_i}{a(\phi)V(\mu_i)}. \] With \(X\) the \(n\times p\) design matrix whose \(i\)th row is \(x_i^\mathsf{T}\), and \(D=\operatorname{diag}(s_1,\ldots,s_n)\), the coefficient score is \(U_\beta=X^\mathsf{T}Dr\), where \(r=(r_1,\ldots,r_n)^\mathsf{T}\).

The score weight is not the residual itself: it scales the residual before the design matrix combines contributions to the coefficients. To derive it, first differentiate the \(i\)th log-likelihood contribution with respect to \(\theta_i\):

\frac{\partial\ell_i}{\partial\theta_i} =\frac{y_i-b'(\theta_i)}{a(\phi)} =\frac{y_i-\mu_i}{a(\phi)}.

Since \(d\mu_i/d\theta_i=b''(\theta_i)=V(\mu_i)\), it follows that \(d\theta_i/d\mu_i=1/V(\mu_i)\). The chain rule therefore gives

\frac{\partial\ell_i}{\partial\eta_i} = \frac{\partial\ell_i}{\partial\theta_i} \frac{d\theta_i}{d\mu_i} \frac{d\mu_i}{d\eta_i} = (y_i-\mu_i)\frac{d\mu_i/d\eta_i}{a(\phi)V(\mu_i)}.

Finally, \(\partial\eta_i/\partial\beta=x_i\), so summing the contributions gives the score and its matrix form.

Result: For an interior parameter value under the stated differentiability and exponential-dispersion assumptions, \[ U_\beta(\beta)=\sum_{i=1}^n x_i s_i r_i=X^\mathsf{T}Dr. \] Thus the coefficient score equations are \(X^\mathsf{T}Dr=0\). These equations characterize stationary points with respect to \(\beta\); by themselves, they do not establish that a solution exists or is a maximum.

What the zero equations constrain

Write \(x_{ij}\) for the \(j\)th component of \(x_i\). The \(j\)th score equation is

\sum_{i=1}^n x_{ij}s_i(y_i-\mu_i)=0,\qquad j=1,\ldots,p.

In words, each column of the design matrix has a zero weighted sum against the response residuals. Equivalently, the vector \(Dr\) has zero ordinary dot product with every column of \(X\), and hence with every vector in the column space of \(X\). If all \(s_i>0\), this can also be expressed as ordinary orthogonality between \(D^{1/2}r\) and each column of \(D^{1/2}X\). This latter description uses a positive-weight geometry; the algebraic equation \(X^\mathsf{T}Dr=0\) remains the basic statement.

The result concerns the fitted design columns, not every possible pattern of covariates. With \(p\) coefficients there are \(p\) score equations, while a data set may contain many distinct covariate patterns. A small number of equations can constrain sums across patterns without requiring each pattern's residual sum to vanish.

Result: Suppose observations fall into \(G\) covariate patterns, with common row vector \(z_g^\mathsf{T}\) for pattern \(g\). If the fitted mean and score weight are consequently common within each pattern, write them as \(\mu_g\) and \(s_g\), and let \(R_g=\sum_{i\text{ in pattern }g}(y_i-\mu_g)\). Then \[ U_\beta=\sum_{g=1}^G z_gs_gR_g. \] If the design includes a separate indicator column for every pattern and each \(s_g\ne0\), zero score forces \(R_g=0\) for every pattern. For a design that does not include those separate indicators, the equations constrain only the design-weighted combination of the \(R_g\)'s.

To verify the pattern statement, group the individual score contributions by their common design row, fitted mean, and weight. Within pattern \(g\), each contribution is \(z_gs_g(y_i-\mu_g)\), so their sum is \(z_gs_gR_g\). Summing over the patterns gives the displayed score. When the design contains the pattern indicator for \(g\), its coefficient score is exactly \(s_gR_g\); if this score is zero and \(s_g\ne0\), division by \(s_g\) gives \(R_g=0\). Without those indicator columns, only the stated linear combinations are required to vanish. \(\square\)

Worked examples

Worked Example: A Bernoulli Probit Score at a Fixed Fit

Let \(Y_i\) be independent Bernoulli responses with success probability \(\mu_i\), and use the probit link \(\eta_i=\Phi^{-1}(\mu_i)\), where \(\Phi\) and \(\varphi\) are the standard normal distribution function and density. For a Bernoulli response, \(a(\phi)=1\) and \(V(\mu_i)=\mu_i(1-\mu_i)\). Since \(d\mu_i/d\eta_i=\varphi(\eta_i)\), the score weight is \[ s_i=\frac{\varphi(\eta_i)}{\mu_i(1-\mu_i)}. \]

Consider an intercept-only model at \(\beta=0\) with four responses \(1,0,1,0\). Each \(\eta_i=0\), so each fitted probability is \(\Phi(0)=0.5\), and each residual is \(y_i-0.5\). The four residuals are \(0.5,-0.5,0.5,-0.5\). At zero, \(\varphi(0)=1/\sqrt{2\pi}\), so the common score weight is \[ s_i=\frac{1/\sqrt{2\pi}}{0.5(0.5)} =\frac{4}{\sqrt{2\pi}} \approx1.5958. \]

The intercept score is therefore \[ U_\beta=1.5958(0.5-0.5+0.5-0.5)=0. \] The raw residuals cancel, and the equal weights preserve that cancellation. This calculation checks that \(\beta=0\) satisfies the score equation for these data; it does not, on its own, establish any broader property of the likelihood.

Worked Example: Pattern-Specific Residual Sums Under a Saturated Design

Let \(Y_i\) be normal with mean \(\mu_i\) and known variance \(\sigma^2=2\). Use a log link, so \(\mu_i=\exp(\eta_i)\), with an intercept and an indicator for a second covariate pattern. The two design rows are \((1,0)\) and \((1,1)\). At \(\beta=(\log 2,\log 2)^\mathsf{T}\), the fitted means are \(2\) in the first pattern and \(4\) in the second.

For the normal model, \(a(\phi)=\sigma^2=2\), \(V(\mu_i)=1\), and \(d\mu_i/d\eta_i=\mu_i\). Thus the score weights are \(s_1=2/2=1\) and \(s_2=4/2=2\). Suppose the two observed responses in the first pattern are \(1,3\), and those in the second pattern are \(3,5\). Their residual sums are \[ R_1=(1-2)+(3-2)=0,\qquad R_2=(3-4)+(5-4)=0. \]

The intercept score is \(s_1R_1+s_2R_2=1(0)+2(0)=0\). The indicator-coefficient score is \(s_2R_2=2(0)=0\). Individual residuals are not all zero, but the residuals cancel within each pattern. Because the design has one indicator for each pattern, the score equations require each pattern's residual sum to be zero whenever its score weight is nonzero.

Worked Example: Zero Score Without a Zero Raw Residual Sum

Again use normal responses with known variance \(1\) and a log link, now with an intercept and a numeric covariate \(x\). Take three observations at \(x=0,1,2\), and set \(\beta=(0,\log 2)^\mathsf{T}\). The fitted means are \(1,2,4\), respectively, and the score weights are \(s_i=\mu_i\).

Choose observed responses \(y=(2,1,4.25)^\mathsf{T}\). The residuals are \(r=(1,-1,0.25)^\mathsf{T}\), and their weighted values are \[ (s_1r_1,s_2r_2,s_3r_3)=(1,-2,1). \] The intercept score is their sum, \(1-2+1=0\). The slope score is the covariate-weighted sum, \[ 0(1)+1(-2)+2(1)=0. \] So the coefficient score is zero at this parameter value, even though the raw residual sum is \(1-1+0.25=0.25\), not zero. The equations constrain the residuals after multiplication by their score weights and the design columns; they do not generally constrain the unweighted total.

Reading the equations carefully

The score weights are determined by the link and the mean-variance relationship. In the general formula, the factor \(d\mu_i/d\eta_i\) accounts for how a change in linear predictor changes the mean, while \(V(\mu_i)\) reflects the response variance structure and \(a(\phi)\) supplies the dispersion scale. A change in any of these factors can change which residual combinations must cancel.

These score weights should not be confused with the “working weights” used in iterative fitting algorithms. A commonly used working weight has a squared derivative, \((d\mu_i/d\eta_i)^2/[a(\phi)V(\mu_i)]\), and appears in a different representation involving an adjusted response. The direct score here uses one factor of \(d\mu_i/d\eta_i\) and multiplies the response residual \(y_i-\mu_i\).

A zero score is a local first-order condition for the likelihood with respect to the regression coefficients, provided the parameter is interior and differentiability holds. It does not say that the fitted means equal every observed response. Nor does it alone guarantee a maximum, uniqueness, or even that an interior solution exists. Those questions require additional analysis.

Common mistake

Wrong reasoning: “The intercept score is zero, so the residuals must sum to zero.” This silently treats all score weights as equal.
Correction: The intercept column is a vector of ones, so its score equation is \(\sum_i s_i r_i=0\). The unweighted sum \(\sum_i r_i\) follows only under extra conditions, such as a common nonzero score weight. With other design columns, the score imposes their own weighted sums as well; it does not generally impose one equation per covariate pattern.

Key takeaway: Write the coefficient score as \(X^\mathsf{T}Dr\) before interpreting it. Its zero says that the weighted residual vector has zero dot product with every fitted design column; stronger pattern-by-pattern conclusions require corresponding columns in the design.

Check Your Understanding

Use the weighted-residual form to answer each question.

For a GLM with an intercept, write the intercept score equation in terms of \(s_i\) and \(r_i\). Under what condition does it imply that the raw residuals sum to zero?
Why does a design with one indicator column for each covariate pattern force each pattern's residual sum to zero when the pattern score weights are nonzero?
A model has three covariate patterns but only an intercept and one slope. Why need the score equations not force each pattern's residual sum to vanish separately?
In the normal log-link example with fitted means \(1,2,4\), verify that residuals \(1,-1,0.25\) give zero intercept and slope scores but a nonzero raw residual sum.
Explain why \(X^\mathsf{T}Dr=0\) is a first-order condition and not, by itself, a proof that the likelihood is maximized.