A Hand-Worked First Poisson IRLS Update

From Starting Coefficients to a First Update

An IRLS update turns the current Poisson fit into a weighted least-squares calculation. For a small dataset, it is possible to show every step: evaluate the fitted means at the starting coefficients, form the working responses and weights, and solve the resulting equations for the next coefficients. The calculation is useful both for understanding the algorithm and for checking software output.

This tutorial considers independent counts $Y_i$ with $Y_i\sim\operatorname{Poisson}(\mu_i)$, and a log-linear model $\log(\mu_i)=\eta_i=x_i^\mathsf{T}\beta$. There is no offset in the examples. The design vector $x_i$ includes an intercept, and $\beta$ denotes the current coefficient vector unless a superscript identifies an iteration. The examples use the Poisson dispersion value $1$.

Definition: At a current coefficient vector $\beta^{(0)}$, the fitted mean is $\mu_i^{(0)}=\exp(x_i^\mathsf{T}\beta^{(0)})$. For a Poisson log-link model, the IRLS working response and weight at observation $i$ are $$ z_i^{(0)}=x_i^\mathsf{T}\beta^{(0)}+\frac{y_i-\mu_i^{(0)}}{\mu_i^{(0)}}, \qquad w_i^{(0)}=\mu_i^{(0)}. $$ The first updated coefficient vector $\beta^{(1)}$ is the weighted least-squares fit of $z^{(0)}$ on the design matrix $X$, using weights $w_i^{(0)}$.

The superscript $(0)$ marks quantities calculated at the starting coefficients; it does not indicate that the observed count is zero. The working response can be fractional or negative, even though the observed count $y_i$ is a nonnegative integer. The weight is computed from the current fitted mean, not from the observed count.

The Hand Calculation

For $n$ observations and $p$ coefficients, put the design vectors into the rows of the $n\times p$ matrix $X$, let $z^{(0)}$ be the column vector of working responses, and let $W^{(0)}$ be the diagonal matrix of working weights. If $X$ has full column rank, then $X^\mathsf{T}W^{(0)}X$ is invertible because every Poisson fitted mean is positive. The update solves

X^\mathsf{T}W^{(0)}X\beta^{(1)} = X^\mathsf{T}W^{(0)}z^{(0)}.

Thus the arithmetic has three stages. First calculate $\eta_i^{(0)}=x_i^\mathsf{T}\beta^{(0)}$ and $\mu_i^{(0)}=\exp(\eta_i^{(0)})$. Next calculate each $z_i^{(0)}$ and $w_i^{(0)}$. Finally solve the weighted normal equations. In a small example, writing those equations explicitly makes sign errors and incorrect coefficients easier to catch.

Worked Example: An Intercept-Only Update

Suppose the observed counts are $y=(0,2,4)^\mathsf{T}$, and the model has an intercept only: $\eta_i=\beta_0$. Start at $\beta_0^{(0)}=0$. Calculate the first update.

For every observation, $\eta_i^{(0)}=0$, so $\mu_i^{(0)}=\exp(0)=1$, and $w_i^{(0)}=1$. The working responses are $$ z_1^{(0)}=0+\frac{0-1}{1}=-1,\qquad z_2^{(0)}=0+\frac{2-1}{1}=1,\qquad z_3^{(0)}=0+\frac{4-1}{1}=3. $$

With an intercept-only design, the weighted least-squares estimate is the weighted average of the working responses. Here all three weights are one, so $$ \beta_0^{(1)} = \frac{(-1)+1+3}{1+1+1} = \frac{3}{3} = 1. $$

The updated linear predictor is $1$ for each observation, and the updated fitted mean is $\exp(1)$ for each. That new mean is calculated after the coefficient update; it is not used to construct the working responses for this first step.

Worked Example: An Intercept and a Slope

Consider three observations with covariate values $x=(-1,0,1)$, counts $y=(0,1,2)$, and model $\eta_i=\beta_0+\beta_1x_i$. Start from $\beta^{(0)}=(0,0)^\mathsf{T}$.

Each starting linear predictor is zero, so each fitted mean and weight is one. The working responses are $z_i^{(0)}=0+(y_i-1)/1$, giving $z^{(0)}=(-1,0,1)^\mathsf{T}$. The design matrix and weighted normal equations are $$ X= \begin{pmatrix} 1&-1\\ 1&0\\ 1&1 \end{pmatrix}, \qquad X^\mathsf{T}X= \begin{pmatrix} 3&0\\ 0&2 \end{pmatrix}, \qquad X^\mathsf{T}z^{(0)} = \begin{pmatrix} 0\\ 2 \end{pmatrix}. $$

Consequently, $3\beta_0^{(1)}=0$ and $2\beta_1^{(1)}=2$, so $\beta^{(1)}=(0,1)^\mathsf{T}$. Direct substitution checks the result: the fitted line to the working responses takes values $-1,0,1$ at $x=-1,0,1$, exactly matching those responses. The updated Poisson means are $\exp(-1),1,\exp(1)$.

Why This Is Fisher Scoring

Theorem: Suppose the observations are independent Poisson variables with means $\mu_i=\exp(x_i^\mathsf{T}\beta)$, and $X$ has full column rank. At any finite current coefficient vector $\beta$, the Poisson log-link IRLS update is the Fisher scoring update $$ \beta_{\mathrm{new}} = \beta+(X^\mathsf{T}WX)^{-1}X^\mathsf{T}(y-\mu), \qquad W=\operatorname{diag}(\mu_1,\ldots,\mu_n). $$

To verify this identity, write $\eta=X\beta$, $\mu_i=\exp(\eta_i)$, and $z_i=\eta_i+(y_i-\mu_i)/\mu_i$. Multiplying each working response by its weight gives $$ w_i z_i = \mu_i\eta_i+(y_i-\mu_i). $$ Since $\eta=X\beta$, collecting these identities across observations yields $$ X^\mathsf{T}Wz = X^\mathsf{T}WX\beta+X^\mathsf{T}(y-\mu). $$ The IRLS normal equations give $X^\mathsf{T}WX\beta_{\mathrm{new}}=X^\mathsf{T}Wz$. Substituting the preceding identity and multiplying by $(X^\mathsf{T}WX)^{-1}$ gives the asserted update.

For completeness, the log-likelihood, ignoring terms that do not depend on $\beta$, is $\ell(\beta)=\sum_i\{y_i x_i^\mathsf{T}\beta-\exp(x_i^\mathsf{T}\beta)\}$. Differentiating gives the score $U(\beta)=X^\mathsf{T}(y-\mu)$. Its negative expected derivative, the expected information, is $X^\mathsf{T}WX$: the derivative of the score is $-X^\mathsf{T}WX$, which does not depend on the observed counts. Thus the expression just derived is the Fisher scoring step. The full-rank condition ensures this information matrix is invertible, since for any nonzero vector $a$, $$ a^\mathsf{T}X^\mathsf{T}WXa = \sum_{i=1}^n \mu_i(x_i^\mathsf{T}a)^2>0. $$ At least one $x_i^\mathsf{T}a$ is nonzero by full column rank, and all $\mu_i>0$. $\square$

This identity provides an independent arithmetic check: the change in coefficients can be obtained by solving the information equations with right-hand side $X^\mathsf{T}(y-\mu)$. It does not mean that the working response is an additional observation; it is the value used to construct this particular weighted least-squares update.

Worked Example: Unequal Starting Weights

Take $x=(0,1,2)$, $y=(1,2,5)$, and $\eta_i=\beta_0+\beta_1x_i$. Start at $\beta^{(0)}=(\log 2,0)^\mathsf{T}$. Calculate the first updated coefficients.

The starting linear predictor is $\log 2$ at each observation, so each fitted mean and weight is $2$. The working responses are $$ z_1^{(0)}=\log 2+\frac{1-2}{2}=\log 2-\frac12,\quad z_2^{(0)}=\log 2,\quad z_3^{(0)}=\log 2+\frac{5-2}{2}=\log 2+\frac32. $$

The design matrix has rows $(1,0),(1,1),(1,2)$. Since all weights are two, $$ X^\mathsf{T}WX= \begin{pmatrix} 6&6\\ 6&10 \end{pmatrix}, \qquad X^\mathsf{T}Wz^{(0)} = \begin{pmatrix} 6\log 2+2\\ 6\log 2+6 \end{pmatrix}. $$

The normal equations are $6\beta_0^{(1)}+6\beta_1^{(1)}=6\log 2+2$ and $6\beta_0^{(1)}+10\beta_1^{(1)}=6\log 2+6$. Subtracting the first equation from the second gives $4\beta_1^{(1)}=4$, hence $\beta_1^{(1)}=1$. Substitution into the first gives $6\beta_0^{(1)}=6\log 2-4$, so $$ \beta^{(1)} = \begin{pmatrix} \log 2-\frac23\\ 1 \end{pmatrix} \approx \begin{pmatrix} 0.0265\\ 1 \end{pmatrix}. $$

As a separate check, the starting residual vector $y-\mu$ is $(-1,0,3)^\mathsf{T}$, so the score is $X^\mathsf{T}(y-\mu)=(2,6)^\mathsf{T}$. The coefficient increment $(-2/3,1)^\mathsf{T}$, multiplied by the information matrix above, gives $(2,6)^\mathsf{T}$. This agrees with the Fisher scoring equation.

What the Weighted Fit Guarantees

Result: If $\beta^{(1)}$ solves the IRLS weighted normal equations, then the working-response residuals $r=z^{(0)}-X\beta^{(1)}$ satisfy $X^\mathsf{T}W^{(0)}r=0$. Thus they are orthogonal to every column of the design matrix under the weighted inner product.

Indeed, the normal equations state $X^\mathsf{T}W^{(0)}X\beta^{(1)}=X^\mathsf{T}W^{(0)}z^{(0)}$. Moving the left side to the right and factoring gives $$ X^\mathsf{T}W^{(0)}\{z^{(0)}-X\beta^{(1)}\}=0. $$ The expression in braces is $r$, proving the result. $\square$

This is a check on the weighted least-squares solution, not a claim that every working response equals its fitted value. With more observations than coefficients, individual residuals will generally remain. It is their weighted sums against the design columns that equal zero.

Common Mistakes and a Reliable Check

A frequent error is to use $y_i$ as the working response, or to use $y_i-\mu_i$ without dividing by $\mu_i$ and adding the current linear predictor. For the Poisson log link, the working response is specifically $\eta_i+(y_i-\mu_i)/\mu_i$. Another error is to calculate the weights from the counts; they are the current fitted means $\mu_i$. Either error changes the weighted regression and therefore changes the update.

A different source of error is to solve the normal equations incorrectly. For example, after calculating a slope, substitute it into an original equation and solve carefully for the intercept; do not divide a coefficient term or constant by the wrong quantity. In the unequal-weight example, the first equation after substitution is $6\beta_0^{(1)}+6=6\log 2+2$, which gives $\beta_0^{(1)}=\log 2-2/3$, not an intercept obtained by treating $6\beta_0^{(1)}$ as $3\beta_0^{(1)}$.

For a hand calculation, check the fitted means and weights first, then verify both sides of the weighted normal equations. When practical, also compare the coefficient increment with the score and information calculation. Agreement between these routes is a strong check that the first update has been computed consistently.

Check Your Understanding

Try these questions using the Poisson log-link formulas from this tutorial.

For one observation with starting linear predictor $\eta^{(0)}=0$ and count $y=3$, calculate its starting mean, weight, and working response.
In an intercept-only model, explain why the IRLS update is a weighted average of the working responses.
For a current fitted mean $\mu=4$, linear predictor $\eta=\log 4$, and observed count $y=2$, calculate the working response and weight.
State the score and expected information for a Poisson log-link model in matrix notation, and explain how they check an IRLS coefficient increment.
What weighted residual equation must hold after the least-squares update, and which matrix defines the weighting?