Score Equations for Inverse-Gaussian Regression

Why the Link Appears in the Score

In inverse-Gaussian regression, the response variance grows with the cube of its mean. Consequently, the coefficient score does not weight each raw residual \(y_i-\mu_i\) equally. Its weight also depends on the fitted mean and on how the chosen link changes the mean when the linear predictor changes. Deriving that score makes the likelihood equations interpretable: at a stationary fit, a particular weighted residual vector is orthogonal to the covariate columns.

The aim here is to derive that condition for a general differentiable link. The focus is the coefficient score and its zero set, not an algorithm for finding its roots or an assessment of whether a root is a maximum.

Model and Score Definition

Suppose \(Y_1,\ldots,Y_n\) are independent inverse-Gaussian responses, with observed values \(y_i>0\), means \(\mu_i>0\), and a common dispersion \(\phi>0\). Use the parameterization in which \(\operatorname{Var}(Y_i)=\phi\mu_i^3\). The covariates are fixed, and the linear predictor and link satisfy \(\eta_i=x_i^\mathsf{T}\beta=g(\mu_i)\), where \(\beta\) is a \(p\)-vector. Assume \(g\) is differentiable on the relevant mean range and \(g'(\mu_i)\ne 0\), so locally \(d\mu_i/d\eta_i=1/g'(\mu_i)\).

Definition: The coefficient score is the gradient of the log-likelihood with respect to \(\beta\), holding \(\phi\) fixed. An interior stationary point for \(\beta\) satisfies the likelihood equations \(U_\beta(\beta)=0\), where \(U_\beta(\beta)=\partial\ell/\partial\beta\). A stationary point need not be a maximum.

For one observation, the inverse-Gaussian density under this parameterization is

f(y_i;\mu_i,\phi) = \frac{1}{\sqrt{2\pi\phi y_i^3}} \exp\left\{-\frac{(y_i-\mu_i)^2}{2\phi\mu_i^2y_i}\right\}, \qquad y_i>0.

Independence gives the log-likelihood by summing the individual log densities. Terms that depend on \(\phi\) or \(y_i\) but not on \(\mu_i\) do not contribute to the derivative with respect to \(\mu_i\).

Deriving the General-Link Coefficient Score

First differentiate the log density with respect to the mean. Writing the exponent's quadratic term as \(y_i^2/\mu_i^2-2y_i/\mu_i+1\) makes the derivative explicit.

Theorem: Under the model and differentiability conditions above, the coefficient score is

U_\beta(\beta) = \sum_{i=1}^n x_i \frac{y_i-\mu_i}{\phi\mu_i^3} \frac{d\mu_i}{d\eta_i} = \frac{1}{\phi}\sum_{i=1}^n x_i (y_i-\mu_i)\mu_i^{-3}\frac{d\mu_i}{d\eta_i}.

Proof. The part of the \(i\)th log density involving \(\mu_i\) is

-\frac{1}{2\phi y_i} \left(\frac{y_i^2}{\mu_i^2}-\frac{2y_i}{\mu_i}+1\right).

Differentiating with respect to \(\mu_i\) gives

\frac{\partial\ell_i}{\partial\mu_i} = -\frac{1}{2\phi y_i} \left(-\frac{2y_i^2}{\mu_i^3}+\frac{2y_i}{\mu_i^2}\right) = \frac{y_i-\mu_i}{\phi\mu_i^3}.

By the chain rule, \(d\mu_i/d\eta_i=1/g'(\mu_i)\), and \(d\eta_i/d\beta=x_i\). Therefore

\frac{\partial\ell_i}{\partial\beta} = \frac{y_i-\mu_i}{\phi\mu_i^3} \frac{d\mu_i}{d\eta_i}x_i.

Summing these contributions over the independent observations proves the stated score formula. \(\square\)

This derivation separates the score into three factors: the raw residual \(y_i-\mu_i\), the inverse-Gaussian mean factor \(\mu_i^{-3}\), and the link factor \(d\mu_i/d\eta_i\). The dispersion supplies the common multiplier \(1/\phi\). The link factor is the derivative of the mean with respect to the predictor, not the derivative of the predictor with respect to the mean.

Likelihood Equations as Weighted Residual Conditions

Define the scalar weighted residual

r_i^{(w)} = (y_i-\mu_i)\mu_i^{-3}\frac{d\mu_i}{d\eta_i}.

The superscript emphasizes that this is a residual multiplied by a model-dependent factor; it is not a standardized residual. Let \(X\) be the \(n\)-by-\(p\) design matrix whose \(i\)th row is \(x_i^\mathsf{T}\), and let \(r^{(w)}\) collect these \(n\) values. The score can then be written as \(U_\beta=X^\mathsf{T}r^{(w)}/\phi\).

Result: For \(\phi>0\), the coefficient likelihood equations are equivalent to

X^\mathsf{T}r^{(w)}=0, \qquad\text{or, equivalently,}\qquad \sum_{i=1}^n x_i (y_i-\mu_i)\mu_i^{-3} \frac{d\mu_i}{d\eta_i}=0.

Justification. The score equation is \(U_\beta=0\). Since \(\phi\) is a positive scalar, multiplying both sides by \(\phi\) does not change its solutions. Substituting the definition of \(r_i^{(w)}\) gives \(X^\mathsf{T}r^{(w)}=0\). In components, for each coefficient \(j=1,\ldots,p\), this says

\sum_{i=1}^n x_{ij}(y_i-\mu_i)\mu_i^{-3} \frac{d\mu_i}{d\eta_i}=0.

Thus each covariate column has zero inner product with the weighted residual vector at an interior stationary fit. This is an algebraic consequence of the score equation; it does not require \(X\) to have full column rank. If columns are linearly dependent, some coefficient directions may not be identifiable, but the displayed score condition remains valid.

Worked Example: Computing a Log-Link Score

Use two observations with covariate rows \(x_1^\mathsf{T}=(1,0)\) and \(x_2^\mathsf{T}=(1,1)\), responses \(y_1=3\) and \(y_2=2\), and the log link \(g(\mu)=\log\mu\). Evaluate the score at \(\beta=(\log 2,\log 2)^\mathsf{T}\).

For the log link, \(\mu_i=\exp(\eta_i)\) and \(d\mu_i/d\eta_i=\mu_i\). The predictors at the stated coefficient values are \(\eta_1=\log 2\) and \(\eta_2=2\log 2\), so the fitted means are \(\mu_1=2\) and \(\mu_2=4\). The scalar factor multiplying \(x_i\), including \(1/\phi\), is

\frac{1}{\phi}(y_i-\mu_i)\mu_i^{-2}.

For observation 1, this factor is \((3-2)/(4\phi)=1/(4\phi)\). For observation 2, it is \((2-4)/(16\phi)=-1/(8\phi)\). Hence

U_\beta = \frac{1}{4\phi} \begin{pmatrix}1\\0\end{pmatrix} - \frac{1}{8\phi} \begin{pmatrix}1\\1\end{pmatrix} = \begin{pmatrix}1/(8\phi)\\-1/(8\phi)\end{pmatrix}.

The score is not zero at this coefficient vector. In particular, the intercept and slope score components are obtained by summing the same scalar factors against the corresponding covariate values.

Worked Example: Reciprocal-Link Weighted Residuals

Use covariate rows \(x_1^\mathsf{T}=(1,0)\) and \(x_2^\mathsf{T}=(1,1)\), responses \(y_1=3\) and \(y_2=2\), and the reciprocal link \(g(\mu)=1/\mu\). Evaluate the score at \(\beta=(1/2,1/2)^\mathsf{T}\).

Here \(\eta_i=1/\mu_i\), so \(\mu_i=1/\eta_i\) and \(d\mu_i/d\eta_i=-1/\eta_i^2=-\mu_i^2\). The two predictors are \(1/2\) and \(1\), giving \(\mu_1=2\) and \(\mu_2=1\). The scalar factor in the score is therefore

\frac{1}{\phi}(y_i-\mu_i)\mu_i^{-3}(-\mu_i^2) = -\frac{y_i-\mu_i}{\phi\mu_i}.

For observation 1, it equals \(-(3-2)/(2\phi)=-1/(2\phi)\). For observation 2, it equals \(-(2-1)/\phi=-1/\phi\). The score is

U_\beta = -\frac{1}{2\phi} \begin{pmatrix}1\\0\end{pmatrix} -\frac{1}{\phi} \begin{pmatrix}1\\1\end{pmatrix} = \begin{pmatrix}-3/(2\phi)\\-1/\phi\end{pmatrix}.

The negative link derivative matters: omitting it would give the wrong sign and the wrong magnitude for both score components.

Worked Example: An Intercept-Only Inverse-Square Link

Consider an intercept-only model with responses \(y_1=1\), \(y_2=2\), and \(y_3=3\), and link \(g(\mu)=1/\mu^2\). The predictor is the intercept \(\beta\), so \(\eta_i=\beta\) and all fitted means are the same value \(\mu\). Evaluate the score equation and its root.

Since \(\eta=\mu^{-2}\), implicit differentiation gives \(d\eta/d\mu=-2\mu^{-3}\), hence \(d\mu/d\eta=-\mu^3/2\). The score for the intercept is

U_\beta = \sum_{i=1}^3 \frac{y_i-\mu}{\phi\mu^3} \left(-\frac{\mu^3}{2}\right) = -\frac{1}{2\phi}\sum_{i=1}^3(y_i-\mu).

Setting this equal to zero yields \(1+2+3-3\mu=0\), so \(\mu=2\) and \(\beta=g(2)=1/4\). At this value the residuals are \(-1,0,1\), and their sum is zero; the score is indeed zero. Here the inverse-square link makes the link and mean factors cancel, leaving an ordinary sum of residuals.

When Does the Ordinary Residual Condition Apply?

The inverse-square example is a special case, not a general rule. In an intercept-only model with a common mean \(\mu\), the score for any differentiable link with nonzero derivative is

U_\beta = \frac{1}{\phi\mu^3} \frac{d\mu}{d\eta} \sum_{i=1}^n(y_i-\mu).

The multiplier outside the residual sum is nonzero under the stated conditions. Thus this intercept-only score equation has the root \(\mu=\bar y\), where \(\bar y=n^{-1}\sum_i y_i\), regardless of the link. This observation concerns the root of this particular score equation; with multiple covariates and varying means, the link-dependent factors generally vary across observations and cannot be removed from the sums.

More generally, the coefficient equations do not assert that raw residuals sum to zero, nor that their unweighted covariate sums vanish. They assert orthogonality for the residuals after multiplication by \(\mu_i^{-3}d\mu_i/d\eta_i\). The factor \(\mu_i^{-3}\) comes from the inverse-Gaussian likelihood, while the derivative factor comes from the regression link.

Common Mistake

A frequent error is to write the score as \(\sum_i x_i(y_i-\mu_i)/(\phi\mu_i^3)\) for every link, treating the mean as though it changed one-for-one with the linear predictor. That expression omits the chain-rule factor \(d\mu_i/d\eta_i\), so it is only appropriate when that derivative is one. Another error is to use \(d\eta_i/d\mu_i\) where the score requires \(d\mu_i/d\eta_i\); these are reciprocals when the link derivative is nonzero.

The safe procedure is to differentiate the log-likelihood with respect to \(\mu_i\) first, then multiply by \(d\mu_i/d\eta_i\), and finally by \(x_i\). Keeping these three steps separate helps identify the source of each weight and prevents accidental sign changes.

Key takeaway: For inverse-Gaussian regression, each observation contributes a coefficient-score factor proportional to \((y_i-\mu_i)\mu_i^{-3}d\mu_i/d\eta_i\). The likelihood equations require this link-weighted residual vector to be orthogonal to every covariate column.

Check Your Understanding

Use the score formula to answer the following questions.

For a general link \(g\), what is \(d\mu_i/d\eta_i\) in terms of \(g'(\mu_i)\), and what condition makes this reciprocal well-defined?
For the log link, simplify the observation-specific factor \((y_i-\mu_i)\mu_i^{-3}d\mu_i/d\eta_i\).
Why does multiplying the coefficient score equation by \(\phi\) leave its solutions unchanged?
Under what model condition do the coefficient equations reduce to an ordinary sum of residuals, and why?
Does \(X^\mathsf{T}r^{(w)}=0\) by itself prove that a stationary point is a likelihood maximum? Explain.