Canonical Links and Simplified IRLS Weights

Why Canonical Links Simplify the Update

In a generalized linear model, an IRLS step is a weighted least-squares calculation built from the current fitted means. Both the working response and the weights depend on how the mean changes with the linear predictor. For a canonical link, that derivative has a particularly simple form: the variance function itself. This causes a cancellation in the coefficient score and makes the working weights proportional to the variance function.

The word proportional matters. The variance function describes how the variance changes with the mean, while the full response variance also includes dispersion. Consequently, depending on the weight convention, an IRLS weight is either the variance function divided by dispersion or just the variance function when a common dispersion factor is omitted.

Exponential-Family Setup and the Canonical Link

Consider independent responses $Y_1,\ldots,Y_n$ with densities or mass functions in exponential-dispersion form $f(y_i;\theta_i,\phi)=\exp\{[y_i\theta_i-b(\theta_i)]/\phi+c(y_i,\phi)\}$. Assume $\phi>0$ is common to the observations and is held fixed when forming the coefficient score. The natural parameter $\theta_i$ lies in the interior of its parameter space, and $b$ is twice differentiable there. The mean and variance are $E[Y_i]=b'(\theta_i)=\mu_i$ and $\operatorname{Var}(Y_i)=\phi b''(\theta_i)=\phi V(\mu_i)$, where $V(\mu_i)=b''(\theta_i)>0$ is the variance function expressed as a function of the mean.

Let $x_i^T$ be row $i$ of the fixed design matrix $X$, and let $\beta$ be the coefficient vector. The linear predictor is $\eta_i=x_i^T\beta$. A canonical link sets the linear predictor equal to the natural parameter: $\eta_i=\theta_i$. Since $d\mu_i/d\theta_i=b''(\theta_i)=V(\mu_i)$, a canonical link satisfies $d\eta_i/d\mu_i=1/V(\mu_i)$.

Definition: A link is canonical for an exponential-family model when it maps the mean to the natural parameter: $g(\mu_i)=\theta_i$, so that $\eta_i=g(\mu_i)=\theta_i$. Where $V(\mu_i)>0$, the inverse derivative is $d\mu_i/d\eta_i=V(\mu_i)$, equivalently $d\eta_i/d\mu_i=1/V(\mu_i)$.

The Coefficient Score Under a Canonical Link

For the stated model, the log-likelihood, up to its dependence on the parameters, is $\ell(\beta)=\sum_i \{[y_i\theta_i-b(\theta_i)]/\phi+c(y_i,\phi)\}$. Differentiating one term with respect to $\theta_i$ gives $(y_i-b'(\theta_i))/\phi=(y_i-\mu_i)/\phi$. Under the canonical link, $d\theta_i/d\beta=x_i$. Thus the coefficient score has no factor involving $V(\mu_i)$.

Theorem: Under the independent exponential-dispersion model above, with $\eta_i=x_i^T\beta=\theta_i$, the coefficient score is $$ U(\beta)=\frac{1}{\phi}X^T(y-\mu), $$ where $y=(y_1,\ldots,y_n)^T$ and $\mu=(\mu_1,\ldots,\mu_n)^T$.

Proof. By the chain rule, the contribution from observation $i$ to the score is $$ \frac{\partial \ell_i}{\partial \beta} =\frac{\partial \ell_i}{\partial \theta_i}\frac{\partial\theta_i}{\partial\beta} =\frac{y_i-\mu_i}{\phi}x_i. $$ Summing over the observations gives $$ U(\beta)=\frac{1}{\phi}\sum_{i=1}^n x_i(y_i-\mu_i) =\frac{1}{\phi}X^T(y-\mu). $$ This differentiation uses only the interior, differentiable parameterization specified above. $\square$

For comparison, with a differentiable noncanonical link the same chain rule leaves the factor $d\theta_i/d\eta_i=(d\mu_i/d\eta_i)/V(\mu_i)$ in the score. Under the canonical link, $d\mu_i/d\eta_i=V(\mu_i)$, so that factor is one. The simplification is therefore a direct consequence of using the natural parameter as the predictor, not a general property of every link.

Worked Example: Bernoulli Response With a Logit Link

Let $Y$ be Bernoulli with mean $\mu$, and use the logit link $\eta=\log\{\mu/(1-\mu)\}$. The Bernoulli natural parameter is this logit, and its variance function is $V(\mu)=\mu(1-\mu)$. For one observation with $y=1$, current mean $\mu=0.2$, and $\phi=1$, calculate the derivative, working response increment, and weight.

First, $V(0.2)=0.2(0.8)=0.16$, so $d\eta/d\mu=1/0.16=6.25$. The current predictor is $\eta=\log(0.2/0.8)=\log(0.25)\approx-1.3863$. The working-response increment is $(y-\mu)d\eta/d\mu=(1-0.2)(6.25)=5$.

Thus the working response is $z\approx-1.3863+5=3.6137$, and the working weight is $W=V(\mu)/\phi=0.16$. To check the cancellation, the score contribution for a scalar predictor is $y-\mu=0.8$. The weighted working residual gives $W(z-\eta)=0.16(5)=0.8$, the same value.

Fisher Information and the Working Weights

The score formula also identifies the information used in Fisher scoring. In the canonical model, differentiating $U(\beta)$ with respect to $\beta$ gives $-X^T\operatorname{diag}(V(\mu_i)/\phi)X$, since $d\mu_i/d\beta=V(\mu_i)x_i$. This Hessian does not depend on the observed responses. Its negative is therefore also the expected Fisher information.

Theorem: Under the same assumptions, the Fisher information for $\beta$ is $$ \mathcal I(\beta)=X^T W X,\qquad W=\operatorname{diag}\left(\frac{V(\mu_1)}{\phi},\ldots,\frac{V(\mu_n)}{\phi}\right). $$ If $X^TWX$ is nonsingular, the Fisher-scoring increment $\delta$, defined by $\mathcal I(\beta)\delta=U(\beta)$, is the solution of the IRLS weighted least-squares normal equations with working response $z_i=\eta_i+(y_i-\mu_i)/V(\mu_i)$ and weights $W_{ii}=V(\mu_i)/\phi$.

Proof. Differentiating the score gives $$ -\frac{\partial U}{\partial\beta^T} =\frac{1}{\phi}\sum_{i=1}^n V(\mu_i)x_ix_i^T =X^TWX. $$ This is the expected information as well, either because the Hessian is already independent of $y$, or by taking its expectation. Define $z_i=\eta_i+(y_i-\mu_i)/V(\mu_i)$. Since $\eta=X\beta$, the weighted least-squares normal equations for the increment $\delta$ are $$ X^TWX\delta=X^TW(z-\eta). $$ For each observation, $$ \frac{V(\mu_i)}{\phi}(z_i-\eta_i) =\frac{V(\mu_i)}{\phi}\frac{y_i-\mu_i}{V(\mu_i)} =\frac{y_i-\mu_i}{\phi}. $$ Therefore $X^TW(z-\eta)=X^T(y-\mu)/\phi=U(\beta)$. The normal equations are exactly the Fisher-scoring equation. Nonsingularity ensures a unique increment. $\square$

This derivation displays the weight simplification: the general working-weight expression is the reciprocal of the response variance times the squared derivative of the link, namely $1/\{\phi V(\mu_i)(d\eta_i/d\mu_i)^2\}$. Substituting the canonical derivative $1/V(\mu_i)$ gives $V(\mu_i)/\phi$. If the common factor $1/\phi$ is omitted from all weights, the working weights are simply $V(\mu_i)$; multiplying every weight by the same positive constant does not change the weighted least-squares coefficient solution.

Worked Example: Poisson Response With a Log Link

For a Poisson response, $\phi=1$, $\mu=e^\eta$, and $V(\mu)=\mu$. With observed count $y=3$ and current mean $\mu=2$, the predictor is $\eta=\log 2\approx0.6931$. The derivative $d\eta/d\mu=1/\mu=0.5$, the working-response increment is $(3-2)(0.5)=0.5$, and hence $z\approx1.1931$. The weight is $W=V(\mu)=2$.

The score contribution for a predictor with design value $x$ is $x(y-\mu)=x$. From the working least-squares expression it is $xW(z-\eta)=x(2)(0.5)=x$. The two calculations agree, and the weight is the current mean, not its reciprocal.

Worked Example: Gamma Response With Its Canonical Link

Use the Gamma exponential-dispersion parameterization with natural parameter $\theta=-1/\mu$, cumulant function $b(\theta)=-\log(-\theta)$, and variance $\phi\mu^2$. Then $b'(\theta)=-1/\theta=\mu$, so $V(\mu)=b''(\theta)=\mu^2$. The canonical link is $\eta=-1/\mu$, which has $d\eta/d\mu=1/\mu^2$.

Take $y=3$, current mean $\mu=2$, and $\phi=0.25$. The current predictor is $\eta=-1/2=-0.5$, and $d\eta/d\mu=1/4=0.25$. The working-response increment is $(3-2)(0.25)=0.25$, so $z=-0.25$. The working weight is $W=\mu^2/\phi=4/0.25=16$. The variance is $\phi\mu^2=1$, consistent with this parameterization. For design value $x$, the score contribution is $x(y-\mu)/\phi=4x$; the weighted working residual is $xW(z-\eta)=x(16)(0.25)=4x$, as required.

What the Simplification Does and Does Not Say

Canonical links simplify both the score and the information, but they do not make the fitted means known or remove the need for iteration. The weights $V(\mu_i)/\phi$ are evaluated at the current fitted means, which generally change after each coefficient update. Nor does “proportional to the variance function” mean “equal to the response variance”: the response variance is $\phi V(\mu_i)$, whereas the canonical IRLS weight is $V(\mu_i)/\phi$ under the full-dispersion convention.

A common mistake is to use the reciprocal variance function as the canonical weight, reasoning that weights should be inversely related to variance. That reasoning overlooks the squared derivative of the link in the working-response scale. The correct calculation is $1/\{\phi V(\mu_i)(d\eta_i/d\mu_i)^2\}$; for a canonical link, $d\eta_i/d\mu_i=1/V(\mu_i)$, so the result is $V(\mu_i)/\phi$, not $1/V(\mu_i)$. Another mistake is to forget dispersion: with a common known $\phi$, dropping its common factor preserves the coefficient solution, but the full Fisher information still contains $1/\phi$.

The simplification is useful as an algebraic check on a score or IRLS derivation. In a canonical model, each observation contributes $(y_i-\mu_i)/\phi$ times its design vector to the score, and the corresponding full working weight is $V(\mu_i)/\phi$. If a derivation produces a different dependence on the mean, recheck the link derivative and whether the variance or its reciprocal was used.

Check Your Understanding

Use the definitions and derivations above to answer these questions.

What condition on the link makes it canonical, and what derivative identity follows when $V(\mu)>0$?
For a canonical model with common dispersion $\phi$, write the coefficient score and explain why the variance function cancels.
If a Poisson model has current mean $\mu=4$, what is its full canonical IRLS weight when $\phi=1$?
For a Gamma model with $V(\mu)=\mu^2$, what is the working weight at mean $\mu=3$ and dispersion $\phi=0.5$?
Why is the canonical working weight not the reciprocal of the variance function?