Poisson Log-Link Updates in IRLS

From Poisson Means to Working Quantities

In a Poisson regression with a log link, the linear predictor describes the logarithm of the fitted mean. IRLS uses that fitted mean to construct a working response and a weight at each iteration. The response is a temporary quantity for a weighted regression calculation; it is not a transformation of the observed count.

Assume $Y_1,\ldots,Y_n$ are independent, with $Y_i\sim\operatorname{Poisson}(\mu_i)$, where $\mu_i>0$. Let $x_i$ be the $p$-vector of fixed covariates for observation $i$, and write the model as $\log(\mu_i)=\eta_i=x_i^\mathsf{T}\beta$. There are no offsets or additional prior weights in these formulas. At an iteration with finite coefficient vector $\beta^{(t)}$, define $\eta_i^{(t)}=x_i^\mathsf{T}\beta^{(t)}$ and $\mu_i^{(t)}=\exp(\eta_i^{(t)})$. In particular, each fitted mean is strictly positive, even if the observed count is zero.

Definition: For the Poisson log-link model, the IRLS working response and working weight at iteration $t$ are $$ z_i^{(t)}=\eta_i^{(t)}+\frac{Y_i-\mu_i^{(t)}}{\mu_i^{(t)}}, \qquad w_i^{(t)}=\mu_i^{(t)}. $$ The working response is defined for every nonnegative count, including $Y_i=0$.

The response formula comes from locally linearizing the mean as a function of the linear predictor. At the current value, the derivative is $d\mu_i/d\eta_i=\mu_i^{(t)}$. Thus, for a nearby predictor $\eta_i$, the first-order approximation is $\mu_i\approx\mu_i^{(t)}+\mu_i^{(t)}(\eta_i-\eta_i^{(t)})$. Setting this approximation equal to the observed response and solving for the predictor gives $z_i^{(t)}$. This construction does not take $\log(Y_i)$; in particular, it remains defined when $Y_i=0$.

Why the Weight Is the Fitted Mean

The Poisson distribution has $E[Y_i]=\mu_i$ and $\operatorname{Var}(Y_i)=\mu_i$. Since the local slope of the mean against the predictor is $\mu_i^{(t)}$, a response discrepancy of a given size corresponds locally to a predictor discrepancy divided by $\mu_i^{(t)}$. The resulting working variance is approximately $\mu_i^{(t)}/(\mu_i^{(t)})^2=1/\mu_i^{(t)}$. Its inverse, used as the weighted least-squares weight, is $\mu_i^{(t)}$.

The same weight follows directly from expected information. For one observation, the Poisson log-likelihood, omitting a term that does not depend on $\beta$, is $Y_i\eta_i-\exp(\eta_i)$. Differentiating with respect to $\beta$ gives score $x_i(Y_i-\mu_i)$. The negative expected second derivative is $\mu_i x_i x_i^\mathsf{T}$, so the observation’s expected information is its fitted mean times $x_i x_i^\mathsf{T}$.

Worked Example: A Zero Count

Suppose the observed count is $Y_i=0$, and the current fitted mean is $\mu_i^{(t)}=0.8$. Then $\eta_i^{(t)}=\log(0.8)\approx-0.22314$. The working response is

z_i^{(t)} =\log(0.8)+\frac{0-0.8}{0.8} \approx -0.22314-1 =-1.22314.

The working weight is $w_i^{(t)}=0.8$. Both quantities are well-defined. The response is negative here, which causes no problem: it is a value on the linear-predictor scale for the weighted regression, not a possible observed count.

Worked Example: A Positive Count and Its Weight

Let $Y_i=4$ and let the current fitted mean be $\mu_i^{(t)}=2$. Then $\eta_i^{(t)}=\log 2$, and

z_i^{(t)}=\log 2+\frac{4-2}{2}=\log 2+1\approx1.69315, \qquad w_i^{(t)}=2.

The observation’s weighted response contribution is $w_i^{(t)}z_i^{(t)}=2\log 2+2\approx3.38629$. Its score contribution, for a scalar intercept-only predictor, is $Y_i-\mu_i^{(t)}=2$; its information contribution is $2$. These quantities are connected in the regression equations below.

The Weighted Regression Update

Let $X$ be the $n\times p$ design matrix whose $i$th row is $x_i^\mathsf{T}$. At iteration $t$, let $W^{(t)}$ be diagonal with entries $\mu_i^{(t)}$, and let $z^{(t)}$ be the vector of working responses. The weighted least-squares step minimizes the following quadratic function over a coefficient vector $b$:

Q^{(t)}(b)=\sum_{i=1}^n \mu_i^{(t)} \left(z_i^{(t)}-x_i^\mathsf{T}b\right)^2 =\left(z^{(t)}-Xb\right)^\mathsf{T}W^{(t)} \left(z^{(t)}-Xb\right).

If $X^\mathsf{T}W^{(t)}X$ is nonsingular, this minimizer is unique and satisfies the weighted normal equations. Nonsingularity holds, for example, when $X$ has full column rank, because every diagonal entry of $W^{(t)}$ is positive.

Theorem: Under the independent Poisson log-link model, if $X^\mathsf{T}W^{(t)}X$ is nonsingular, the weighted least-squares solution $b$ is exactly the Fisher scoring update $\beta^{(t+1)}$.

Proof. Differentiating $Q^{(t)}(b)$ with respect to $b$ gives the normal equations $X^\mathsf{T}W^{(t)}Xb=X^\mathsf{T}W^{(t)}z^{(t)}$. For the Poisson model, the score and expected information at $\beta^{(t)}$ are

U^{(t)}=\sum_{i=1}^n x_i(Y_i-\mu_i^{(t)}), \qquad I^{(t)}=\sum_{i=1}^n \mu_i^{(t)}x_ix_i^\mathsf{T} =X^\mathsf{T}W^{(t)}X.

The right-hand side of the weighted normal equations can be expanded using the definition of $z_i^{(t)}$:

X^\mathsf{T}W^{(t)}z^{(t)} =\sum_{i=1}^n x_i\mu_i^{(t)} \left(\eta_i^{(t)}+\frac{Y_i-\mu_i^{(t)}}{\mu_i^{(t)}}\right) =I^{(t)}\beta^{(t)}+U^{(t)}.

Therefore the normal equations give $b=\beta^{(t)}+(I^{(t)})^{-1}U^{(t)}$, which is the Fisher scoring update. Uniqueness follows from the assumed nonsingularity. $\square$

This identity explains why an IRLS step can be implemented as weighted regression: the weights and response are chosen so that its normal equations reproduce the Fisher scoring equation. It does not assert that the working responses are newly observed data or that the quadratic function is the original Poisson log-likelihood.

What Each Observation Contributes

Each row contributes separately to both sides of the weighted normal equations. For observation $i$, its contribution to the matrix on the left is $\mu_i^{(t)}x_ix_i^\mathsf{T}$, and its contribution to the vector on the right is $x_i\mu_i^{(t)}z_i^{(t)}$. Substituting the working response gives

x_i\mu_i^{(t)}z_i^{(t)} =x_i\left(\mu_i^{(t)}\eta_i^{(t)}+Y_i-\mu_i^{(t)}\right).

Thus the fitted mean determines the observation’s weight in the regression and also appears in its weighted-response contribution. The observed count enters the latter through the residual $Y_i-\mu_i^{(t)}$. The whole system is obtained by adding these observation-specific contributions.

Worked Example: Two Observations and Their Matrix Contributions

Consider two observations with covariate vectors $x_1=(1,0)^\mathsf{T}$ and $x_2=(1,1)^\mathsf{T}$. Suppose the current coefficients are $\beta^{(t)}=(0,\log 3)^\mathsf{T}$, so the current predictors are $0$ and $\log 3$, and the fitted means are $1$ and $3$. Let the observed counts be $Y_1=2$ and $Y_2=1$.

For observation 1, $z_1^{(t)}=0+(2-1)/1=1$ and $w_1^{(t)}=1$. Its matrix contribution is $1\begin{pmatrix}1&0\\0&0\end{pmatrix}$, and its right-hand-side contribution is $(1,0)^\mathsf{T}$.

For observation 2, $z_2^{(t)}=\log 3+(1-3)/3=\log 3-2/3\approx0.43195$, and $w_2^{(t)}=3$. Its matrix contribution is $3\begin{pmatrix}1&1\\1&1\end{pmatrix}$. Its right-hand-side contribution is $3(\log 3-2/3)(1,1)^\mathsf{T}\approx(1.29584,1.29584)^\mathsf{T}$.

Adding the contributions gives

X^\mathsf{T}W^{(t)}X= \begin{pmatrix}4&3\\3&3\end{pmatrix}, \qquad X^\mathsf{T}W^{(t)}z^{(t)} \approx \begin{pmatrix}2.29584\\1.29584\end{pmatrix}.

The matrix has determinant $4\cdot3-3\cdot3=3$, so it is nonsingular and the weighted regression has a unique solution. This calculation shows explicitly how the fitted means, covariate vectors, and working responses enter the update system.

Common Mistake and Practical Interpretation

A tempting but incorrect derivation begins by writing $z_i^{(t)}$ as a Taylor approximation to $\log(Y_i)$ for every observation. That expression is undefined when $Y_i=0$. The correction is to define the working response directly as $\eta_i^{(t)}+(Y_i-\mu_i^{(t)})/\mu_i^{(t)}$, which is valid for all nonnegative counts. For positive counts, a Taylor approximation involving $\log(Y_i)$ may offer intuition, but it is not needed to define or derive the IRLS working response.

Another potential confusion is to treat $w_i^{(t)}=\mu_i^{(t)}$ as a fixed property of the observation. It is recalculated at every iteration because the fitted mean changes with the current coefficient vector. A larger current fitted mean therefore gives that row a larger weight in the current weighted regression, while the next iteration may assign it a different weight.

Key takeaway: For Poisson regression with a log link, each iteration uses $z_i^{(t)}=\eta_i^{(t)}+(Y_i-\mu_i^{(t)})/\mu_i^{(t)}$ and $w_i^{(t)}=\mu_i^{(t)}$. These choices make the weighted normal equations equal the Fisher scoring equation, with each observation contributing $\mu_i^{(t)}x_ix_i^\mathsf{T}$ and $x_i\mu_i^{(t)}z_i^{(t)}$.

Check Your Understanding

Use the formulas and reasoning developed here to answer the following.

For an observed count of zero and a current fitted mean of $1.5$, calculate the working response and weight.
Why does deriving the working response from a local linearization of the mean avoid the zero-count problem?
For one observation with covariate vector $x_i$, state its contribution to each side of the weighted normal equations.
Explain why the IRLS weight changes from one iteration to the next.
What condition ensures that the weighted least-squares update is unique, and why are the diagonal weights positive in this model?