Worked Asymptotics for Poisson Regression

Poisson Regression for Rates

In a Poisson regression for event counts, the amount of time or population at risk can differ from one observation to another. A count of 12 events over 4 person-years and a count of 12 events over 40 person-years do not indicate the same event rate. An offset accounts for this known exposure while the regression models the rate.

Let $Y_1,\ldots,Y_n$ be independent counts, let $t_i>0$ be the known exposure for observation $i$, and let $x_i$ be a fixed $p$-dimensional covariate vector. The log-link Poisson model is

Y_i\sim\operatorname{Poisson}(\mu_i),\qquad \log\mu_i=\log t_i+x_i^\mathsf{T}\beta,\qquad \mu_i=t_i\exp(x_i^\mathsf{T}\beta).

Here $\beta\in\mathbb R^p$ is the coefficient vector. The offset $\log t_i$ enters with coefficient fixed at one; it is not an additional coefficient to estimate. The rate per unit of exposure is $\mu_i/t_i=\exp(x_i^\mathsf{T}\beta)$. Thus, holding other covariates fixed, a one-unit increase in covariate $j$ multiplies the rate by $\exp(\beta_j)$.

Likelihood, Score, and Information

For observed counts $y_1,\ldots,y_n$, independence and the Poisson probability mass function give the log-likelihood

\ell(\beta) =\sum_{i=1}^n\left\{ y_i\bigl(\log t_i+x_i^\mathsf{T}\beta\bigr) -t_i\exp(x_i^\mathsf{T}\beta)-\log(y_i!) \right\}.

Differentiating with respect to the coefficient vector gives the score $U_n(\beta)$. Differentiating again gives the Hessian. The expected information is the negative expected Hessian, and in this model it equals the observed information for every possible data set.

Definition: For the log-link Poisson regression with known offsets $t_i$, the score and information at $\beta$ are $$ U_n(\beta)=\sum_{i=1}^n x_i(y_i-\mu_i),\qquad \mathcal I_n(\beta)=\sum_{i=1}^n \mu_i x_i x_i^\mathsf{T}, \qquad \mu_i=t_i\exp(x_i^\mathsf{T}\beta). $$ The score is a $p$-vector and the information is a $p\times p$ matrix.

To verify these expressions, the derivative of $y_i x_i^\mathsf{T}\beta-\mu_i$ is $x_i(y_i-\mu_i)$, because $\partial\mu_i/\partial\beta=\mu_i x_i$. Its derivative is $-\mu_i x_i x_i^\mathsf{T}$. Therefore the negative Hessian is $\sum_i\mu_i x_i x_i^\mathsf{T}$, which does not depend on the observed counts. Taking its expectation changes nothing, proving the stated information formula.

Worked Example: Two Groups and a Rate-Ratio Coefficient

Suppose independent control and treated group counts are $Y_0=30$ over exposure $t_0=10$, and $Y_1=48$ over exposure $t_1=12$. Use an intercept and a treatment indicator, so that $\beta=(\alpha,\gamma)^\mathsf{T}$, where $\alpha$ is the log control rate and $\gamma$ is the log rate ratio.

The group means are $\mu_0=10e^\alpha$ and $\mu_1=12e^{\alpha+\gamma}$. Apart from terms that do not involve the parameters, the log-likelihood is

\ell(\alpha,\gamma) =78\alpha+48\gamma-10e^\alpha-12e^{\alpha+\gamma}+\text{constant}.

The score equations are $78-\mu_0-\mu_1=0$ and $48-\mu_1=0$. Hence $\hat\mu_1=48$, $\hat\mu_0=30$, $\hat\alpha=\log(30/10)=\log 3$, and $\hat\gamma=\log(48/12)-\hat\alpha=\log(4/3)$. The fitted control and treated rates are 3 and 4 per exposure unit, respectively; the fitted rate ratio is $e^{\hat\gamma}=4/3$.

At a general parameter value the information matrix is

\mathcal I(\alpha,\gamma) = \begin{pmatrix} \mu_0+\mu_1 & \mu_1\\ \mu_1 & \mu_1 \end{pmatrix}.

At the fitted values this is $\begin{pmatrix}78&48\\48&48\end{pmatrix}$. Its determinant is $78(48)-48^2=1440$, so its inverse is

\mathcal I(\hat\beta)^{-1} =\frac{1}{1440} \begin{pmatrix}48&-48\\-48&78\end{pmatrix}.

The information-based variance estimate for the log rate ratio is the bottom-right entry, $78/1440=13/240$, giving standard error $\sqrt{13/240}\approx0.233$. This calculation illustrates the plug-in formula; with only two group counts and two parameters, the numerical approximation should not be mistaken for a large-sample guarantee.

Asymptotic Distribution of the Coefficients

The information matrix describes local curvature of the log-likelihood. To use it asymptotically, suppose the true parameter $\beta_0$ is an interior point of $\mathbb R^p$, the design and offsets are fixed, and the observations are independent Poisson variables with means $\mu_{i0}=t_i\exp(x_i^\mathsf{T}\beta_0)$. Assume the maximum-likelihood estimator $\hat\beta$ exists with probability tending to one and is consistent. Also assume

$\mathcal I_n(\beta_0)/n\to J$, where $J$ is positive definite;
the normalized score satisfies $U_n(\beta_0)/\sqrt n\xrightarrow{d}N(0,J)$;
along the segment from $\beta_0$ to $\hat\beta$, the normalized information converges uniformly in probability to $J$.

These are regular large-sample conditions. For example, bounded covariates and bounded offsets away from zero and infinity, together with suitable convergence of the design averages and the stated consistency and nonsingularity conditions, provide a familiar setting for them. The score covariance is $\mathcal I_n(\beta_0)$, since independent Poisson variables have variance equal to their means.

Theorem: Under these conditions, $$ \sqrt n(\hat\beta-\beta_0)\xrightarrow{d}N(0,J^{-1}). $$ Consequently, the asymptotic covariance of $\hat\beta$ is $J^{-1}/n$, equivalently approximated in practice by $\mathcal I_n(\hat\beta)^{-1}$.

Proof. With probability tending to one, the estimator is an interior likelihood maximizer, so $U_n(\hat\beta)=0$. Apply the integral form of Taylor's theorem to the vector-valued score; this avoids assuming a multivariate mean-value identity with a single intermediate point. Define

A_n=\int_0^1\mathcal I_n\bigl(\beta_0+s(\hat\beta-\beta_0)\bigr)\,ds.

Since the derivative of the score is $-\mathcal I_n(\beta)$, the integral expansion gives

0=U_n(\hat\beta) =U_n(\beta_0)-A_n(\hat\beta-\beta_0).

The assumed uniform convergence along the segment implies $A_n/n\xrightarrow{p}J$. Because $J$ is positive definite, the continuous mapping theorem gives $(A_n/n)^{-1}\xrightarrow{p}J^{-1}$. Rearranging the expansion,

\sqrt n(\hat\beta-\beta_0) =(A_n/n)^{-1}\frac{U_n(\beta_0)}{\sqrt n}.

The score central limit theorem and Slutsky's theorem now yield the stated limiting distribution. $\square$

For a coefficient $\beta_j$, let $e_j$ be the vector with a one in position $j$ and zeros elsewhere. Its asymptotic variance is $e_j^\mathsf{T}J^{-1}e_j/n$; the estimated standard error uses the square root of the corresponding diagonal entry of $\mathcal I_n(\hat\beta)^{-1}$. In particular, a coefficient for a binary treatment is a log rate ratio, so its standard error is on the log-rate-ratio scale.

Worked Example: An Intercept-Only Rate

Suppose $Y=18$ events are observed over $t=6$ exposure units, with model $Y\sim\operatorname{Poisson}(t e^\alpha)$. The log-likelihood, ignoring a parameter-free term, is $\ell(\alpha)=18\alpha-6e^\alpha$. Its score is $18-6e^\alpha$, so the maximum-likelihood estimate solves $6e^{\hat\alpha}=18$. Thus $\hat\alpha=\log 3$, and the estimated rate is $e^{\hat\alpha}=3$ events per exposure unit.

The information at $\alpha$ is $6e^\alpha$, which at the estimate is 18. The estimated variance of $\hat\alpha$ is therefore $1/18$, and its estimated standard error is $1/\sqrt{18}\approx0.236$. To obtain uncertainty for the rate $\lambda=e^\alpha$, use the derivative $d\lambda/d\alpha=e^\alpha$. The delta-method variance estimate is

\widehat{\operatorname{Var}}(\hat\lambda) \approx \hat\lambda^2\,\widehat{\operatorname{Var}}(\hat\alpha) =3^2\left(\frac{1}{18}\right)=\frac12.

Thus the estimated standard error of the rate is $\sqrt{1/2}\approx0.707$. This is an asymptotic, information-based estimate, not an exact finite-sample variance.

Uncertainty in a Fitted Mean

For a new covariate vector $x_*$ and known exposure $t_*>0$, the fitted mean is $m(\beta)=t_*\exp(x_*^\mathsf{T}\beta)$. Its gradient is $m(\beta)x_*$. The delta method therefore gives the following result whenever the coefficient estimator has the asymptotic normal distribution above.

Result: At the true parameter, the asymptotic variance of the fitted mean $m(\hat\beta)$ is $$ \frac{m(\beta_0)^2}{n}\,x_*^\mathsf{T}J^{-1}x_*. $$ The usual plug-in estimate is $$ \widehat{\operatorname{Var}}\{m(\hat\beta)\} \approx \hat m_*^2\,x_*^\mathsf{T}\mathcal I_n(\hat\beta)^{-1}x_*, \qquad \hat m_*=t_*\exp(x_*^\mathsf{T}\hat\beta). $$

Indeed, the gradient of $m$ at $\beta_0$ is $m(\beta_0)x_*$. Applying the delta method to the limiting covariance $J^{-1}/n$ gives the displayed variance. Replacing unknown quantities by their estimates gives the plug-in expression.

Worked Example: Uncertainty for a Treated-Group Mean

Return to the two-group fit, and consider the expected count for treated-group exposure $t_*=5$. The covariate vector is $x_*=(1,1)^\mathsf{T}$. The fitted mean is $\hat m_*=5\exp(\hat\alpha+\hat\gamma)=5(4)=20$.

Using the inverse information matrix calculated above,

x_*^\mathsf{T}\mathcal I(\hat\beta)^{-1}x_* =\frac{48-48-48+78}{1440} =\frac{30}{1440} =\frac{1}{48}.

The estimated variance of the fitted mean is therefore $20^2/48=400/48=25/3\approx8.33$, and its standard error is $\sqrt{25/3}\approx2.89$. This measures uncertainty in the estimated expected count at exposure 5. It does not include the additional Poisson variation of a future observed count, whose conditional variance, if the mean were known, would equal that mean.

Interpretation and a Common Pitfall

Offsets preserve the rate interpretation when exposures differ: the regression models $\log(\mu_i/t_i)=x_i^\mathsf{T}\beta$, not $\log\mu_i=x_i^\mathsf{T}\beta$ unless every exposure is one. Omitting $\log t_i$ can make the model attribute exposure differences to covariates and distort the estimated rate effects.

A second common mistake is to report the Poisson outcome variance as the variance of the fitted mean. The fitted mean is a function of estimated coefficients, so its parameter uncertainty is obtained from the coefficient covariance and the delta method. A future count has additional sampling variation: conditional on the parameters its variance is its mean. These are distinct quantities, and combining them is appropriate only when the target is prediction of a future count rather than uncertainty about the mean itself.

Key takeaway: With a log-link Poisson model and a known exposure offset, the score is $\sum_i x_i(y_i-\mu_i)$ and the information is $\sum_i\mu_i x_i x_i^\mathsf{T}$. Inverting that information gives coefficient uncertainty; applying the delta method gives uncertainty for rates and fitted means.

Check Your Understanding

Work through these questions using the formulas and distinctions developed above.

For one observation with exposure $t_i$ and covariate vector $x_i$, write its mean and explain why $\log t_i$ is an offset rather than a coefficient.
In a treatment-indicator model, what does the exponentiated treatment coefficient represent?
For an intercept-only model with count $y$ and exposure $t$, derive the maximum-likelihood estimate of the rate.
What assumptions allow the inverse information matrix to estimate the asymptotic covariance of the coefficient estimator?
Explain the difference between uncertainty in a fitted mean and the conditional variance of a future Poisson count.