Score Equations For Poisson Regression

From the Poisson Likelihood to a Score Vector

In Poisson regression, each observation contributes information about the regression coefficients through its count and its covariate values. The score vector collects the derivatives of the log-likelihood with respect to all coefficients. Setting this vector to zero gives the likelihood equations that a coefficient estimate must satisfy when it is an interior maximum.

Assume \(Y_1,\ldots,Y_n\) are conditionally independent Poisson responses given fixed covariate vectors \(\mathbf{x}_i\), with \(\mathbf{x}_i\in\mathbb{R}^p\). Let \(\boldsymbol{\beta}\in\mathbb{R}^p\) be the coefficient vector, and use the log link: \(\mu_i=E[Y_i\mid\mathbf{x}_i]=\exp(\mathbf{x}_i^\mathsf{T}\boldsymbol{\beta})\). The observed counts \(y_i\) are nonnegative integers. The predictors are treated as fixed while differentiating with respect to \(\boldsymbol{\beta}\).

Definition: The score vector for \(\boldsymbol{\beta}\) is the vector of partial derivatives of the log-likelihood: \(U(\boldsymbol{\beta})=\nabla_{\boldsymbol{\beta}}\ell(\boldsymbol{\beta})\). Its \(j\)th entry is \(\partial\ell(\boldsymbol{\beta})/\partial\beta_j\). The likelihood equations are the equations \(U(\boldsymbol{\beta})=\mathbf{0}\).

For one observation, the Poisson probability mass function is \(P(Y_i=y_i\mid\mathbf{x}_i)=\exp(-\mu_i)\mu_i^{y_i}/y_i!\). Conditional independence makes the full likelihood the product of these contributions. Taking logarithms and substituting the log-link mean gives

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

The factorial term is part of the log-likelihood, but it does not depend on \(\boldsymbol{\beta}\). The exponential term depends on every coefficient through the linear predictor, so its derivative must include the corresponding covariate value.

Deriving the Score

Write \(\beta_j\) for the \(j\)th coefficient and \(x_{ij}\) for the \(j\)th component of \(\mathbf{x}_i\). The chain rule gives \(\partial\exp(\mathbf{x}_i^\mathsf{T}\boldsymbol{\beta})/\partial\beta_j =\exp(\mathbf{x}_i^\mathsf{T}\boldsymbol{\beta})x_{ij} =\mu_i x_{ij}\). Differentiating each term of the log-likelihood therefore yields the following result.

Theorem: Under the independent Poisson regression model with \(\mu_i=\exp(\mathbf{x}_i^\mathsf{T}\boldsymbol{\beta})\), the \(j\)th score component is \(U_j(\boldsymbol{\beta})=\sum_{i=1}^n x_{ij}(y_i-\mu_i)\), for \(j=1,\ldots,p\).

Proof. For each \(j\), differentiate the displayed log-likelihood term by term. The derivative of \(y_i\mathbf{x}_i^\mathsf{T}\boldsymbol{\beta}\) is \(y_i x_{ij}\), the derivative of the negative exponential term is \(-\mu_i x_{ij}\), and the derivative of \(-\log(y_i!)\) is zero. Summing these contributions over \(i\) gives

\frac{\partial\ell(\boldsymbol{\beta})}{\partial\beta_j} =\sum_{i=1}^n(y_i x_{ij}-\mu_i x_{ij}) =\sum_{i=1}^n x_{ij}(y_i-\mu_i). \qquad\square

This component-by-component expression can be assembled into a matrix formula. Let \(X\) be the \(n\times p\) design matrix whose \(i\)th row is \(\mathbf{x}_i^\mathsf{T}\). Let \(\mathbf{y}=(y_1,\ldots,y_n)^\mathsf{T}\), and let \(\boldsymbol{\mu}(\boldsymbol{\beta})\) be the \(n\)-vector with entries \(\mu_i=\exp(\mathbf{x}_i^\mathsf{T}\boldsymbol{\beta})\). The exponential is applied to each component of \(X\boldsymbol{\beta}\).

Score in matrix form: \(U(\boldsymbol{\beta})=X^\mathsf{T}\{\mathbf{y}-\boldsymbol{\mu}(\boldsymbol{\beta})\}\), where \(\boldsymbol{\mu}(\boldsymbol{\beta})=\exp(X\boldsymbol{\beta})\) componentwise.

To verify the matrix expression, the \(j\)th entry of \(X^\mathsf{T}(\mathbf{y}-\boldsymbol{\mu})\) is the inner product of the \(j\)th column of \(X\) with \(\mathbf{y}-\boldsymbol{\mu}\). That inner product is exactly \(\sum_i x_{ij}(y_i-\mu_i)\), the component just derived. Thus the matrix formula is a compact arrangement of the same \(p\) partial derivatives, not a different set of equations.

The Likelihood Equations

Replacing the score vector by zero produces \(p\) equations for the \(p\) coefficients. In component form, each equation uses the corresponding predictor column; in matrix form, all equations are written together.

\sum_{i=1}^n x_{ij}(y_i-\mu_i)=0,\quad j=1,\ldots,p, \qquad\text{or equivalently}\qquad X^\mathsf{T}\{\mathbf{y}-\exp(X\boldsymbol{\beta})\}=\mathbf{0}.

This equivalence follows directly from the component calculation: a vector is zero precisely when each of its entries is zero. The equations are generally nonlinear in \(\boldsymbol{\beta}\), because each fitted mean is an exponential of a linear function of the coefficients. A solution must therefore be understood as a coefficient vector at which the fitted means make all the score components zero.

Result: If \(\boldsymbol{\beta}^\ast\in\mathbb{R}^p\) is a local maximum of the Poisson regression log-likelihood, then it satisfies \(X^\mathsf{T}\{\mathbf{y}-\exp(X\boldsymbol{\beta}^\ast)\}=\mathbf{0}\). A solution of these equations is a stationary point; the equation alone does not establish that it is a maximum.

Justification. The log-likelihood is differentiable for every \(\boldsymbol{\beta}\in\mathbb{R}^p\), because exponentials and linear functions are differentiable and the factorial terms are constant in \(\boldsymbol{\beta}\). The parameter space \(\mathbb{R}^p\) is open. At an interior local maximum of a differentiable function on an open set, every coordinate partial derivative must be zero: fixing all other coordinates gives a differentiable one-dimensional function with a local maximum, whose derivative is zero. Consequently every component of \(U(\boldsymbol{\beta}^\ast)\) is zero, which gives the displayed matrix equation. The reverse implication does not follow from this argument: zero derivative identifies a stationary point, not its type. \(\square\)

Worked Calculations

Worked Example: Solving the Intercept-Only Score Equation

Suppose \(n=4\), the observed counts are \(y=(0,2,3,3)^\mathsf{T}\), and the model has an intercept only. Thus \(x_i=1\) and \(\mu_i=\exp(\beta_0)\) for all four observations. Find the score equation and solve it.

The score for the intercept is \(U_0(\beta_0)=\sum_i(y_i-\mu_i)\). The observed counts sum to \(0+2+3+3=8\), and the four fitted means are each \(\exp(\beta_0)\). Therefore

U_0(\beta_0)=8-4\exp(\beta_0).

Setting the score equal to zero gives \(8-4\exp(\beta_0)=0\), so \(\exp(\beta_0)=2\) and \(\beta_0=\log 2\). Substitution checks the equation: each fitted mean is \(2\), their sum is \(8\), and the score is \(8-8=0\).

Worked Example: Evaluating a Two-Coefficient Score

Consider three observations with an intercept and a numerical predictor. Let \(X=\begin{pmatrix}1&0\\1&1\\1&2\end{pmatrix}\), \(\mathbf{y}=(1,1,5)^\mathsf{T}\), and evaluate the score at \(\boldsymbol{\beta}=(0,\log 2)^\mathsf{T}\).

First calculate the linear predictors: \(X\boldsymbol{\beta}=(0,\log 2,2\log 2)^\mathsf{T}\). Exponentiating componentwise gives the fitted means \(\boldsymbol{\mu}=(1,2,4)^\mathsf{T}\). Hence \(\mathbf{y}-\boldsymbol{\mu}=(0,-1,1)^\mathsf{T}\). The score is

X^\mathsf{T}(\mathbf{y}-\boldsymbol{\mu}) = \begin{pmatrix} 1&1&1\\ 0&1&2 \end{pmatrix} \begin{pmatrix}0\\-1\\1\end{pmatrix} = \begin{pmatrix}0\\1\end{pmatrix}.

Thus the intercept component is \(0\), while the slope component is \(1\). This coefficient vector does not solve the likelihood equations, because the full score vector is not zero. In particular, checking only the intercept equation would not be enough when the model has two coefficients.

Worked Example: Solving the Equations for Two Groups

Four independent counts have predictor values \(x=(0,0,1,1)^\mathsf{T}\) and observed values \(y=(1,2,2,4)^\mathsf{T}\). Use a model with linear predictor \(\beta_0+\beta_1x_i\). Write and solve both score equations.

The design matrix has rows \((1,0),(1,0),(1,1),(1,1)\). The fitted mean is \(\exp(\beta_0)\) for the first two observations and \(\exp(\beta_0+\beta_1)\) for the last two. The intercept equation sums all four observed-minus-fitted contributions:

(1+2)-2\exp(\beta_0)+(2+4)-2\exp(\beta_0+\beta_1)=0.

The slope equation weights each contribution by \(x_i\), so the first two terms are omitted: \(6-2\exp(\beta_0+\beta_1)=0\). Solving this equation gives \(\exp(\beta_0+\beta_1)=3\). Substituting into the intercept equation gives \(3-2\exp(\beta_0)+6-6=0\), and therefore \(\exp(\beta_0)=1.5\). Taking logarithms yields \(\beta_0=\log(1.5)\) and \(\beta_1=\log(3)-\log(1.5)=\log 2\).

Check both equations using these fitted means. The four means are \(1.5,1.5,3,3\). The intercept score is \((1+2+2+4)-(1.5+1.5+3+3)=9-9=0\); the slope score is \((2+4)-(3+3)=0\). Both likelihood equations are satisfied.

Common Mistakes and Practical Use

A frequent differentiation error is to treat \(\mu_i\) as though it were unrelated to the coefficients. In the log-likelihood, \(\mu_i=\exp(\mathbf{x}_i^\mathsf{T}\boldsymbol{\beta})\) changes with \(\beta_j\), and its derivative is \(\mu_i x_{ij}\). Omitting \(x_{ij}\) would incorrectly give the same contribution to every coefficient score, regardless of the design matrix.

Using the wrong sign. Differentiating the negative exponential term gives \(-\mu_i x_{ij}\), so the score contribution is \(x_{ij}(y_i-\mu_i)\), not \(x_{ij}(\mu_i-y_i)\).
Differentiating the factorial term. The quantity \(\log(y_i!)\) depends on the observed count, not on \(\boldsymbol{\beta}\), and contributes zero to the score.
Checking only one equation. With \(p\) coefficients there are \(p\) score components. A candidate must satisfy all of them, not just the intercept equation.
Treating a score root as proof of a maximum. The score equations specify a stationary point. Additional analysis is required to determine whether it is a maximum.

The matrix form is useful because it separates the ingredients of the calculation: the design matrix records the covariates, and the fitted-mean vector is obtained by applying the log-link model at the proposed coefficients. Together they give a concise system for evaluating the score and formulating the likelihood equations. Numerical methods can use this system to search for coefficient values where the score vanishes; the equations themselves state the target, not a guarantee that such a solution exists or is a maximum.

Check Your Understanding

Use the score formulas to answer the following questions.

For one observation with covariate vector \(\mathbf{x}_i\), show the contribution to the \(j\)th score component.
What are the dimensions of \(X\), \(\mathbf{y}-\boldsymbol{\mu}\), and \(X^\mathsf{T}(\mathbf{y}-\boldsymbol{\mu})\) when there are \(n\) observations and \(p\) coefficients?
For an intercept-only model with three observations, write the score equation in terms of the observed counts and \(\beta_0\).
Why does a zero score vector not, by itself, prove that the coefficient vector maximizes the log-likelihood?
In the model with an intercept and predictor \(x_i\), which design-matrix entries weight the slope score contributions?