Deriving the Poisson Regression Score

Why Derive the Score?

Poisson regression specifies how a conditional mean changes with covariates. Its maximum-likelihood estimates are defined by the likelihood, but the score equations reveal what the fitted model must match in the data. The Hessian, in turn, describes the likelihood’s curvature and helps establish whether a solution is a maximum and whether it is unique.

The derivation below is for the ordinary Poisson likelihood under a correctly specified conditional distribution. It is useful to keep that scope clear: the likelihood calculations rely on independent Poisson outcomes, not merely on a chosen conditional mean. The focus here is the likelihood, its derivatives, and the characterization of its maximum.

Model and assumptions: For observations $i=1,\ldots,n$, let $Y_i$ be a count and let $\mathbf{x}_i$ be a fixed $p$-dimensional covariate vector, usually including an intercept. Conditional on the covariates, assume the $Y_i$ are independent and $$ Y_i\mid \mathbf{x}_i \sim \operatorname{Poisson}(\mu_i), \qquad \mu_i=\exp(\mathbf{x}_i^\top\boldsymbol{\beta}), $$ where $\boldsymbol{\beta}$ is a $p$-dimensional coefficient vector. Thus $P(Y_i=y_i\mid\mathbf{x}_i)=\exp(-\mu_i)\mu_i^{y_i}/y_i!$ for nonnegative integer $y_i$.

Log-Likelihood, Score, and Hessian

Write the observed counts as $y_1,\ldots,y_n$. Conditional independence means that the joint likelihood is the product of the individual Poisson probabilities. Taking its logarithm turns that product into a sum. Since $\log(\mu_i)=\mathbf{x}_i^\top\boldsymbol{\beta}$, the log-likelihood is

\ell(\boldsymbol{\beta}) =\sum_{i=1}^n\left[ y_i\mathbf{x}_i^\top\boldsymbol{\beta} -\exp(\mathbf{x}_i^\top\boldsymbol{\beta}) -\log(y_i!) \right].

The factorial term depends on the observed data but not on $\boldsymbol{\beta}$. It must be retained in the full log-likelihood, although it contributes zero to derivatives with respect to the coefficients.

Definition: The score vector is the vector of first derivatives of the log-likelihood with respect to $\boldsymbol{\beta}$. The Hessian is the matrix of second derivatives. For the Poisson regression model above, $$ \mathbf{U}(\boldsymbol{\beta}) =\frac{\partial\ell(\boldsymbol{\beta})}{\partial\boldsymbol{\beta}} =\sum_{i=1}^n\mathbf{x}_i(y_i-\mu_i), \qquad \mathbf{H}(\boldsymbol{\beta}) =\frac{\partial^2\ell(\boldsymbol{\beta})} {\partial\boldsymbol{\beta}\,\partial\boldsymbol{\beta}^\top} =-\sum_{i=1}^n\mu_i\mathbf{x}_i\mathbf{x}_i^\top. $$

To see the first derivative, the contribution of observation $i$ to the log-likelihood is $y_i\mathbf{x}_i^\top\boldsymbol{\beta}-\mu_i-\log(y_i!)$. The derivative of the linear term is $y_i\mathbf{x}_i$, while the derivative of $\mu_i=\exp(\mathbf{x}_i^\top\boldsymbol{\beta})$ is $\mu_i\mathbf{x}_i$. Subtracting gives $\mathbf{x}_i(y_i-\mu_i)$. Differentiating this score contribution again gives $-\mu_i\mathbf{x}_i\mathbf{x}_i^\top$, and summing over observations yields the Hessian.

The maximum-likelihood score equations set each component of $\mathbf{U}(\boldsymbol{\beta})$ to zero. Equivalently, the fitted conditional means must satisfy

\sum_{i=1}^n\mathbf{x}_i y_i = \sum_{i=1}^n\mathbf{x}_i\mu_i.

These are covariate-weighted matching conditions. When the model has an intercept, its score equation says that the sum of observed counts equals the sum of fitted means. Other coefficient equations match the corresponding covariate-weighted totals. They do not generally require every fitted mean to equal its observed count.

Worked Example: An Intercept-Only Poisson Model

Suppose three hypothetical independent observations have counts $y_1=0$, $y_2=1$, and $y_3=5$. Use an intercept-only model, so $\mathbf{x}_i=1$ and $\mu_i=\exp(\beta)$ for every observation.

The log-likelihood and score are $$ \ell(\beta)=6\beta-3\exp(\beta)-\log(0!)-\log(1!)-\log(5!), \qquad U(\beta)=6-3\exp(\beta). $$ Setting the score equal to zero gives $3\exp(\hat\beta)=6$, so $\exp(\hat\beta)=2$ and $\hat\beta=\log(2)$. The fitted conditional mean is $2$ counts per observation, equal to the observed average, $6/3=2$.

At this estimate, the score is $6-3(2)=0$, and the Hessian is $H(\hat\beta)=-3(2)=-6$. The negative second derivative confirms downward curvature at the solution. In this intercept-only model the log-likelihood is strictly concave, so this stationary point is its unique maximum.

Curvature and the Maximum-Likelihood Characterization

The Hessian has a useful quadratic-form interpretation. For any $p$-dimensional vector $\mathbf{a}$,

\mathbf{a}^\top\mathbf{H}(\boldsymbol{\beta})\mathbf{a} =-\sum_{i=1}^n\mu_i(\mathbf{a}^\top\mathbf{x}_i)^2 \leq 0.

Every $\mu_i$ is positive, so the log-likelihood is concave in $\boldsymbol{\beta}$. If the covariate vectors span all $p$ dimensions (the design matrix has full column rank), then for every nonzero $\mathbf{a}$, at least one $\mathbf{a}^\top\mathbf{x}_i$ is nonzero. The quadratic form is then strictly negative, making the Hessian negative definite and the log-likelihood strictly concave.

Consequently, if a finite parameter vector solves the score equations, it is a global maximum; under full column rank it is the unique finite maximum. The qualification “if a finite solution exists” matters. For example, if all observed counts are zero in an intercept-only model, then $\ell(\beta)=-n\exp(\beta)$ up to a constant. Its supremum is approached as $\beta$ tends to negative infinity, and there is no finite score-zero estimate.

Worked Example: Binary Treatment and the Hessian

Consider four hypothetical independent observations with an intercept and a treatment indicator $D_i$. The untreated counts are $1$ and $3$; the treated counts are $2$ and $6$. The model is $\mu_i=\exp(\beta_0+\beta_1D_i)$. The untreated sample mean is $4/2=2$, and the treated sample mean is $8/2=4$.

The score equations require the fitted means for each group to equal its group average. Thus $\exp(\hat\beta_0)=2$, and $\exp(\hat\beta_0+\hat\beta_1)=4$. It follows that $\hat\beta_0=\log(2)$ and $\hat\beta_1=\log(4)-\log(2)=\log(2)$. At these estimates the fitted means are $2$ for each untreated observation and $4$ for each treated observation.

For the intercept score, the sum of residuals is $(1-2)+(3-2)+(2-4)+(6-4)=0$. For the treatment score, only treated observations contribute, and $(2-4)+(6-4)=0$. Hence the score vector is zero.

The untreated observations have covariate vector $(1,0)^\top$, and their fitted means sum to $2+2=4$. The treated observations have vector $(1,1)^\top$, and their fitted means sum to $4+4=8$. Therefore $$ \mathbf{H}(\hat{\boldsymbol{\beta}}) =-\left[ 4\begin{pmatrix}1&0\\0&0\end{pmatrix} +8\begin{pmatrix}1&1\\1&1\end{pmatrix} \right] =-\begin{pmatrix}12&8\\8&8\end{pmatrix}. $$ The positive weights $4$ and $8$ are the total fitted means in the untreated and treated groups, respectively. The design has full column rank because the two groups have different covariate vectors; the Hessian is negative definite, so this finite score-zero estimate is the unique maximum.

Using the Derivatives in Computation

When the score equations do not have a closed-form solution, an iterative method can use the score and Hessian to update the coefficient vector. A Newton step for maximizing the log-likelihood is

\boldsymbol{\beta}_{\mathrm{new}} =\boldsymbol{\beta}_{\mathrm{old}} -\mathbf{H}(\boldsymbol{\beta}_{\mathrm{old}})^{-1} \mathbf{U}(\boldsymbol{\beta}_{\mathrm{old}}).

This is a computational update, not a guarantee that one step reaches the maximum. It requires an invertible Hessian at the current value; practical fitting procedures may also use safeguards such as step-size control.

Worked Example: One Newton Step

Take three hypothetical observations with covariates $x_i=0,1,2$, counts $y_i=0,1,3$, and model $\mu_i=\exp(\beta_0+\beta_1x_i)$. Start at $\boldsymbol{\beta}_{\mathrm{old}}=(0,0)^\top$, so each fitted mean is $1$. The residuals $y_i-\mu_i$ are $-1,0,2$.

The score is $$ \mathbf{U}(\boldsymbol{\beta}_{\mathrm{old}}) =\begin{pmatrix}1&0\\1&1\\1&2\end{pmatrix}^{\!\top} \begin{pmatrix}-1\\0\\2\end{pmatrix} =\begin{pmatrix}1\\4\end{pmatrix}. $$ Since each $\mu_i=1$, the negative Hessian is $$ -\mathbf{H}(\boldsymbol{\beta}_{\mathrm{old}}) =\sum_{i=1}^3\mathbf{x}_i\mathbf{x}_i^\top =\begin{pmatrix}3&3\\3&5\end{pmatrix}. $$ Its determinant is $3(5)-3(3)=6$, so its inverse is $\frac{1}{6}\begin{pmatrix}5&-3\\-3&3\end{pmatrix}$. Substituting into the Newton update gives $$ \boldsymbol{\beta}_{\mathrm{new}} =\begin{pmatrix}0\\0\end{pmatrix} +\frac{1}{6}\begin{pmatrix}5&-3\\-3&3\end{pmatrix} \begin{pmatrix}1\\4\end{pmatrix} =\begin{pmatrix}-7/6\\3/2\end{pmatrix} \approx\begin{pmatrix}-1.167\\1.500\end{pmatrix}. $$ This is one update from the starting value, not a claim that the score equations are now solved. A fitting algorithm would evaluate the score and Hessian again and continue until an appropriate convergence criterion is met.

Common Pitfall and Practical Meaning

A common mistake is to read the score equations as saying that every residual $y_i-\mu_i$ must be zero. The equations instead say that the residuals sum to zero after weighting by each covariate vector: $\sum_i\mathbf{x}_i(y_i-\mu_i)=\mathbf{0}$. For an intercept this implies equality of the observed and fitted totals, but individual observations can have positive or negative residuals. A second pitfall is to treat a zero score as sufficient without checking whether it corresponds to a maximum or whether a finite solution exists. Concavity and the Hessian provide the relevant curvature check; full rank gives strict concavity and uniqueness when a finite solution is attained.

The score and Hessian are also the building blocks of numerical estimation and likelihood-based uncertainty calculations. Their formulas make explicit how each observation contributes: the score contribution is its covariate vector times its count residual, while the Hessian contribution is a negative weighted outer product of its covariates. Thus observations with larger fitted means contribute more curvature, with weights $\mu_i$.

Key takeaway: Under independent Poisson sampling with $\log(\mu_i)=\mathbf{x}_i^\top\boldsymbol{\beta}$, the score is $\sum_i\mathbf{x}_i(y_i-\mu_i)$, and the Hessian is $-\sum_i\mu_i\mathbf{x}_i\mathbf{x}_i^\top$. A finite score-zero solution is a global maximum; with full column rank it is the unique finite maximum.

Check Your Understanding

Answer each question using the model and derivative results above.

Why does the term $\log(y_i!)$ appear in the log-likelihood but not in the score?
In an intercept-only model, what condition does the score equation impose on the sum of fitted means?
Explain why the Hessian is negative semidefinite for every coefficient vector.
What additional condition makes the log-likelihood strictly concave, and what does that imply for a finite score-zero solution?
In the binary-treatment example, why does the treatment score involve only the treated observations?