Poisson Regression Score Equations

Poisson Regression With a Log Link

In Poisson regression, a count observation has a mean that depends on predictors, while its conditional distribution remains Poisson. With the log link, the linear predictor is the logarithm of the mean. This structure makes the likelihood score especially interpretable: at a finite solution, the observed counts and fitted means have equal totals when each is weighted by any column of the design matrix.

The equations are useful both for understanding fitted models and for checking calculations. They do not generally require each fitted mean to equal its observed count. Instead, the design matrix specifies the particular totals that must agree.

Definition: Let \(Y_1,\ldots,Y_n\) be independent observations with \(Y_i\sim\operatorname{Poisson}(\mu_i)\), where \(\mu_i>0\). Let \(x_i=(x_{i1},\ldots,x_{ip})^\mathsf{T}\) be a fixed design vector and \(\beta=(\beta_1,\ldots,\beta_p)^\mathsf{T}\) the regression coefficient vector. Under the log link, \(\log(\mu_i)=x_i^\mathsf{T}\beta\), so \(\mu_i=\exp(x_i^\mathsf{T}\beta)\). The design matrix \(X\) has rows \(x_i^\mathsf{T}\).

This tutorial considers the model without an offset. Write \(\eta_i=x_i^\mathsf{T}\beta\), so that \(\mu_i=\exp(\eta_i)\). Because the exponential is positive, every finite coefficient vector gives a positive fitted mean.

Likelihood and Score

For a Poisson observation, \(P(Y_i=y_i)=\exp(-\mu_i)\mu_i^{y_i}/y_i!\) for a nonnegative integer \(y_i\). Independence lets us multiply these probabilities across observations. Taking logarithms gives the log-likelihood as a function of \(\beta\):

L(\beta)=\prod_{i=1}^n \frac{\exp(-\mu_i)\mu_i^{y_i}}{y_i!}, \qquad \ell(\beta)=\sum_{i=1}^n\left[y_i x_i^\mathsf{T}\beta-\exp(x_i^\mathsf{T}\beta)-\log(y_i!)\right].

The factorial terms matter in the likelihood itself, but do not depend on \(\beta\). To differentiate, use \(\partial\mu_i/\partial\beta_j=\mu_i x_{ij}\). Thus the derivative of the \(i\)th parameter-dependent contribution is \(x_{ij}(y_i-\mu_i)\).

Theorem: For the independent Poisson log-link model, the score vector and observed Hessian are

U(\beta)=\frac{\partial\ell(\beta)}{\partial\beta} =X^\mathsf{T}(y-\mu), \qquad H(\beta)=\frac{\partial^2\ell(\beta)}{\partial\beta\,\partial\beta^\mathsf{T}} =-X^\mathsf{T}\operatorname{diag}(\mu_1,\ldots,\mu_n)X.

Here \(y=(y_1,\ldots,y_n)^\mathsf{T}\) and \(\mu=(\mu_1,\ldots,\mu_n)^\mathsf{T}\), with each \(\mu_i=\exp(x_i^\mathsf{T}\beta)\).

Proof. For coordinate \(j\), differentiating each summand of \(\ell\) gives \(\partial\ell/\partial\beta_j=\sum_i x_{ij}(y_i-\mu_i)\). Collecting these \(p\) coordinate derivatives gives \(U(\beta)=X^\mathsf{T}(y-\mu)\). Differentiating coordinate \(j\) of the score with respect to coordinate \(k\) gives \(\partial^2\ell/(\partial\beta_j\partial\beta_k)=-\sum_i x_{ij}\mu_i x_{ik}\). These are precisely the entries of \(-X^\mathsf{T}\operatorname{diag}(\mu)X\). \(\square\)

The Hessian records the curvature of the log-likelihood. For any vector \(a\in\mathbb{R}^p\),

a^\mathsf{T}H(\beta)a =-\sum_{i=1}^n\mu_i(x_i^\mathsf{T}a)^2\leq 0.

Every \(\mu_i\) is positive, so the inequality holds for every coefficient vector and every design. In particular, the log-likelihood is concave. If \(X\) has full column rank, then \(Xa\ne 0\) for every nonzero \(a\); at least one squared term is positive, making the Hessian negative definite. This is strict concavity in the full-rank case. The statement concerns curvature; by itself it does not guarantee that a finite maximizing coefficient vector exists.

What the Score Equations Match

For a finite interior maximum, differentiability implies that the score is zero. More generally, the score equations are the equations \(U(\beta)=0\). Reading them one coordinate at a time gives a simple interpretation.

Theorem: The \(j\)th Poisson regression score equation is

\sum_{i=1}^n x_{ij}y_i=\sum_{i=1}^n x_{ij}\mu_i, \qquad j=1,\ldots,p.

Thus the observed counts and fitted means have equal totals after weighting by the \(j\)th design-matrix column.

Proof. The \(j\)th score coordinate is \(U_j(\beta)=\sum_i x_{ij}(y_i-\mu_i)\). Setting it equal to zero and moving the fitted-mean sum to the other side gives the displayed equation. This holds separately for every column \(j\), including columns whose entries are not indicators. \(\square\)

An intercept column consists of ones. Its equation therefore says that the sum of the observed counts equals the sum of the fitted means. A group-indicator column says that the observed and fitted totals agree within that group. A continuous predictor column instead matches a predictor-weighted total; it does not say that the unweighted totals agree within each predictor value.

Worked Examples

Worked Example: The Intercept-Only Poisson Model

Suppose three independent counts are \(y_1=1\), \(y_2=2\), and \(y_3=5\), and fit an intercept-only model. Each design value is \(x_i=1\), so \(\mu_i=\exp(\beta)\) for all three observations.

The log-likelihood is \(\ell(\beta)=8\beta-3\exp(\beta)-\log(1!)-\log(2!)-\log(5!)\), because \(1+2+5=8\). Its score is \(U(\beta)=8-3\exp(\beta)\). Setting the score equal to zero gives \(\exp(\hat\beta)=8/3\), hence \(\hat\beta=\log(8/3)\).

The fitted mean for each observation is \(8/3\), so the fitted total is \(3(8/3)=8\), exactly the observed total. The Hessian is \(H(\beta)=-3\exp(\beta)\), which at \(\hat\beta\) equals \(-8\). It is negative, as expected for this one-parameter, full-rank model.

Worked Example: Matching Totals for Two Groups

Suppose four counts are \(2,4,1,5\), with the first two observations in group A and the last two in group B. Use two indicator columns and no separate intercept: \(x_1=x_2=(1,0)^\mathsf{T}\) and \(x_3=x_4=(0,1)^\mathsf{T}\). Write the coefficients as \((\beta_A,\beta_B)^\mathsf{T}\).

The fitted means are \(\mu_1=\mu_2=\exp(\beta_A)\) and \(\mu_3=\mu_4=\exp(\beta_B)\). The score for group A is \(U_A=(2+4)-2\exp(\beta_A)=6-2\exp(\beta_A)\); the score for group B is \(U_B=(1+5)-2\exp(\beta_B)=6-2\exp(\beta_B)\).

Setting both scores to zero gives \(\exp(\hat\beta_A)=3\) and \(\exp(\hat\beta_B)=3\), so \(\hat\beta_A=\hat\beta_B=\log 3\). The observed group totals are \(2+4=6\) and \(1+5=6\). The fitted totals are \(3+3=6\) in each group. The individual fitted values need not equal the observations: in group A, for example, both fitted means are 3, whereas the observed counts are 2 and 4.

Worked Example: A Continuous Predictor Column

Consider three observations with predictor values \(x_i=0,1,2\) and counts \(y_i=1,2,4\). Fit an intercept and a slope, so \(x_i^\mathsf{T}=(1,x_i)\). At \(\beta=(0,\log 2)^\mathsf{T}\), the linear predictors are \(0,\log 2,2\log 2\), and the fitted means are \(1,2,4\).

The intercept-column equation compares unweighted totals: the observed total is \(1+2+4=7\), and the fitted total is \(1+2+4=7\). For the slope column, the observed weighted total is \(0(1)+1(2)+2(4)=10\), and the fitted weighted total is \(0(1)+1(2)+2(4)=10\). Thus both score coordinates are zero at this coefficient vector.

The Hessian there is \[ -\begin{pmatrix} \sum_i\mu_i & \sum_i x_i\mu_i\\ \sum_i x_i\mu_i & \sum_i x_i^2\mu_i \end{pmatrix} = -\begin{pmatrix} 7&10\\ 10&18 \end{pmatrix}. \] The determinant of the matrix inside the minus sign is \(7(18)-10^2=26>0\), and its first diagonal entry is positive, so that matrix is positive definite. Consequently the displayed Hessian is negative definite, consistent with the full-rank design having strict concavity.

Common Pitfall and Practical Interpretation

A frequent mistake is to read \(X^\mathsf{T}(y-\mu)=0\) as \(y_i=\mu_i\) for every observation. That conclusion does not follow: the score gives \(p\) weighted-sum equations, not \(n\) separate equations. In the group example, the fitted means are both 3 in group A, even though its observed counts are 2 and 4. Their group total matches, which is what the indicator-column equation requires.

The exact constraints depend on the columns included in \(X\). With an intercept, the total fitted count equals the total observed count. With group indicators, the totals match for the represented groups. With a continuous predictor, its column imposes equality of predictor-weighted totals. A coefficient's score equation is therefore a direct link between the model's design and the data summaries that the fit must reproduce.

These identities are useful for checking a derivation or a numerical fit: compute \(X^\mathsf{T}(y-\mu)\) and verify that each component is zero, up to numerical tolerance. They also explain why model coding matters. Changing the design columns changes which weighted totals are constrained, even when the underlying counts have not changed.

Check Your Understanding

Use the score and Hessian formulas to answer these questions.

In an intercept-only Poisson model, what equality does the score equation impose on the observed and fitted counts?
For a design column of group indicators, interpret the corresponding score equation.
Why does a continuous predictor's score equation generally not require each fitted mean to equal its observed count?
State the observed Hessian for the Poisson log-link model and explain why it is negative semidefinite.
In the two-group example, calculate the fitted total in each group at the score solution and compare it with the observed total.