Sandwich Inference for Poisson Regression

Why Sandwich Inference Is Needed

Poisson regression specifies a log-linked conditional mean, $\mu_i=\exp(x_i'\beta)$. Poisson quasi-maximum likelihood (QMLE) estimates $\beta$ using that mean relationship. Its conventional covariance formula also uses the Poisson variance restriction $\operatorname{Var}(Y_i\mid X_i)=\mu_i$. If the conditional variance differs from the mean, that conventional formula can give misleading standard errors even when the conditional mean is correctly specified.

The heteroskedasticity-robust sandwich covariance replaces the Poisson variance assumption with an estimate based on the observed score contributions. Its two main components are called the bread and the meat. The bread describes how quickly the estimating equation changes with the coefficients; the meat describes the variation in the estimating equation across observations.

Definition: Let $Y_i$ be a count outcome and $x_i$ a $k$-vector of regressors, including an intercept if appropriate. Poisson QMLE specifies $\mu_i(\beta)=\exp(x_i'\beta)$ and solves $$ \sum_{i=1}^n x_i\{Y_i-\mu_i(\hat\beta)\}=0. $$ The score contribution for observation $i$ is $s_i(\beta)=x_i\{Y_i-\mu_i(\beta)\}$. The empirical bread and meat are $$ \widehat A=\sum_{i=1}^n \hat\mu_i x_i x_i', \qquad \widehat B=\sum_{i=1}^n \hat u_i^2 x_i x_i', \qquad \hat u_i=Y_i-\hat\mu_i. $$ The heteroskedasticity-robust sandwich covariance estimator is $\widehat V=\widehat A^{-1}\widehat B\widehat A^{-1}$.

Assumptions and Interpretation

For the usual robust inference, observations are independent, the regressors provide enough variation to identify the coefficients, and the relevant moments are finite. For a conditional-mean interpretation of the coefficients, the mean specification must also be correct: $E(Y_i\mid x_i)=\exp(x_i'\beta_0)$. The conditional variance may vary with $x_i$ and need not equal the conditional mean. In that setting, Poisson QMLE can consistently estimate the conditional-mean coefficients while the sandwich estimates their sampling covariance.

Under independent, identically distributed sampling, the sandwich can also support inference about the QMLE pseudo-true parameter when the mean model is misspecified, subject to the usual regularity conditions. That parameter solves the population estimating equation $E[x_i\{Y_i-\exp(x_i'\beta)\}]=0$. It need not describe the true conditional mean. Thus robust standard errors address uncertainty in the estimator; they do not repair a misspecified mean relationship.

The bread in the formula comes from the derivative of the score. Since $$ \frac{\partial s_i(\beta)}{\partial\beta'}=-\mu_i(\beta)x_i x_i', $$ the negative derivative of the total score is $\sum_i\mu_i x_i x_i'$. The meat uses the outer products of the score contributions: $s_i(\hat\beta)s_i(\hat\beta)'=\hat u_i^2x_i x_i'$. If the conditional variance really is Poisson, the expected meat and bread coincide; otherwise, they need not.

A short expansion explains why the two pieces appear on either side of the meat. Write $A=E[\mu_i x_i x_i']$ and $B=E[s_i(\beta_0)s_i(\beta_0)']$, assuming a correctly specified conditional mean. A first-order expansion of the estimating equation around $\beta_0$ gives

0 \approx \frac{1}{\sqrt n}\sum_{i=1}^n s_i(\beta_0) -A\sqrt n(\hat\beta-\beta_0), \qquad \sqrt n(\hat\beta-\beta_0) \approx A^{-1}\frac{1}{\sqrt n}\sum_{i=1}^n s_i(\beta_0).

The score sum has asymptotic covariance $B$, so the coefficient estimator has asymptotic covariance $A^{-1}BA^{-1}$ on the $\sqrt n$ scale. Replacing the population quantities with their sample estimates yields the stated covariance for $\hat\beta$. The two bread factors are why the expression is called a sandwich.

Worked Calculation: An Intercept-Only Model

Worked Example: Robust Inference for a Mean Count

Suppose four independent, hypothetical observations have counts $0,1,2,3$. Fit an intercept-only Poisson QMLE, so $x_i=1$ and $\mu_i=\exp(\beta)$. Compute the robust standard error for $\hat\beta$ and test $H_0:\beta=0$ against a two-sided alternative.

The intercept score equation is $\sum_i(Y_i-\exp(\beta))=0$, giving $\hat\mu=\bar Y=(0+1+2+3)/4=1.5$, and $\hat\beta=\log(1.5)=0.4055$. The fitted residuals are $-1.5,-0.5,0.5,1.5$.

The scalar bread and meat are $$ \widehat A=\sum_i\hat\mu=4(1.5)=6, \qquad \widehat B=\sum_i\hat u_i^2=2.25+0.25+0.25+2.25=5. $$ Therefore $\widehat V=5/6^2=0.1389$, and $SE(\hat\beta)=\sqrt{0.1389}=0.3727$.

The robust Wald statistic is $z=0.4055/0.3727=1.088$, with a two-sided large-sample $p$-value of approximately $0.277$. At the 5% level, this test does not reject $H_0$. This tiny sample illustrates the arithmetic, not a setting in which large-sample calibration should be expected to perform reliably.

Robust Tests and Confidence Intervals

For a single coefficient $\beta_j$, the robust standard error is the square root of the $j$th diagonal element of $\widehat V$. To test $H_0:\beta_j=b_0$, calculate $z=(\hat\beta_j-b_0)/SE(\hat\beta_j)$ and compare it with a standard normal reference distribution in a large sample. A two-sided $100(1-\alpha)\%$ Wald confidence interval is

\hat\beta_j \ \pm\ z_{1-\alpha/2}\,SE(\hat\beta_j).

This interval is for the coefficient on the log-mean scale. Exponentiating both endpoints gives an interval for $\exp(\beta_j)$, a conditional mean ratio for a one-unit increase in $x_j$ when the other regressors are held fixed. It does not turn a confidence interval into a statement that there is a particular probability the fixed coefficient lies inside it.

For $q$ linear restrictions $H_0:R\beta=r$, the robust Wald statistic is

W=(R\hat\beta-r)' \{R\widehat V R'\}^{-1} (R\hat\beta-r).

Under the null and regularity conditions, $W$ is approximately chi-square distributed with $q$ degrees of freedom. A robust covariance matrix changes the estimated uncertainty used in these tests and intervals; it does not change the fitted coefficients.

Worked Calculation: A Robust Slope Test

Worked Example: Comparing Two Conditional Mean Counts

Consider four hypothetical independent observations, with a binary regressor $z_i$: the two observations with $z_i=0$ have counts $0,2$, and the two with $z_i=1$ have counts $1,3$. Fit $\mu_i=\exp(\beta_0+\beta_1z_i)$. Compute the robust covariance matrix and test whether the conditional mean differs between the two groups.

The score equations match each group's fitted mean to its sample mean. Thus $\hat\mu=1$ for $z=0$, $\hat\mu=2$ for $z=1$, $\hat\beta_0=\log(1)=0$, and $\hat\beta_1=\log(2)=0.6931$. With $x_i=(1,z_i)'$, the bread is $$ \widehat A= \begin{pmatrix}2&0\\0&0\end{pmatrix} + \begin{pmatrix}4&4\\4&4\end{pmatrix} = \begin{pmatrix}6&4\\4&4\end{pmatrix}. $$

For $z=0$, the two residuals are $-1,1$; for $z=1$, they are also $-1,1$. Their squared residuals are all 1. The meat is therefore $$ \widehat B= \begin{pmatrix}2&0\\0&0\end{pmatrix} + \begin{pmatrix}2&2\\2&2\end{pmatrix} = \begin{pmatrix}4&2\\2&2\end{pmatrix}. $$

Since $\widehat A^{-1}=\begin{pmatrix}0.5&-0.5\\-0.5&0.75\end{pmatrix}$, multiplication gives $$ \widehat V=\widehat A^{-1}\widehat B\widehat A^{-1} = \begin{pmatrix}0.5&-0.5\\-0.5&0.625\end{pmatrix}. $$ The robust standard error for the slope is $\sqrt{0.625}=0.7906$. The test statistic for $H_0:\beta_1=0$ is $z=0.6931/0.7906=0.8767$, with a two-sided large-sample $p$-value of approximately $0.381$. The 95% robust interval is $0.6931\pm1.96(0.7906)=[-0.8565,2.243]$. On exponentiating, the estimated mean ratio is $2.000$, with an interval of approximately $[0.425,9.42]$. The data provide little precision about the group difference.

Worked Calculation: A Joint Wald Test

Worked Example: Testing Both Regression Coefficients

Using the same hypothetical group data and fitted model, test the joint null $H_0:\beta_0=0,\ \beta_1=0$. This asks whether both the baseline log conditional mean and the group effect are zero. Use the robust covariance matrix calculated above.

Here $\hat\beta=(0,\log 2)'$, $R$ is the $2$-by-$2$ identity matrix, and $r=(0,0)'$. The inverse robust covariance is $$ \widehat V^{-1}= \begin{pmatrix}10&8\\8&8\end{pmatrix}. $$ Thus the Wald statistic is $$ W= \begin{pmatrix}0&\log 2\end{pmatrix} \begin{pmatrix}10&8\\8&8\end{pmatrix} \begin{pmatrix}0\\\log 2\end{pmatrix} =8(\log 2)^2 =3.844. $$

The reference distribution is chi-square with $q=2$ degrees of freedom. The corresponding upper-tail $p$-value is approximately $0.146$. At the 5% level, the joint null is not rejected. This is an asymptotic Wald calculation; with only four observations, its numerical result should be treated as an illustration rather than dependable small-sample evidence.

Limits and Common Pitfalls

The sandwich is robust to an unrestricted conditional variance for independent observations when the mean model and other required conditions support the target of inference. It is not automatically robust to dependence between observations. If outcomes are correlated within, for example, individuals observed repeatedly or hospitals, an observation-level meat treats correlated score contributions as independent and can understate or otherwise misstate uncertainty. Such settings require a covariance estimator that accounts for the dependence structure.

Using the conventional covariance by default. The inverse bread alone relies on the Poisson variance restriction. When that restriction is not warranted, use the sandwich for robust inference.
Thinking robust standard errors repair the model. They do not correct a wrong conditional mean, omitted regressors, or an inappropriate functional form.
Ignoring the sampling structure. The observation-level sandwich does not account for correlated outcomes; independence must be justified or the covariance estimator must match the dependence.
Reading a non-rejection as proof of no effect. A large $p$-value means the data do not provide strong evidence against the stated null under the test assumptions; it does not establish that the null is true.

Some software applies a finite-sample multiplier, such as $n/(n-k)$, to the meat. Such adjustments can affect small-sample standard errors, but they are conventions rather than substitutes for the assumptions behind the covariance estimator. Always identify whether reported robust standard errors use such a correction.

Key takeaway: Poisson QMLE's sandwich covariance combines a bread based on fitted means with a meat based on observed score variation. It permits robust standard errors, tests, and intervals without imposing equality of the conditional variance and mean, but it does not fix mean-model misspecification or unaccounted dependence.

Check Your Understanding

Try these questions without looking back.

What does the bread measure, and how does it differ from the meat in Poisson QMLE?
Which variance restriction does the sandwich avoid imposing?
If the robust standard error for a coefficient is $0.20$ and its estimate is $0.50$, calculate the Wald $z$-statistic for testing that the coefficient is zero.
Does a heteroskedasticity-robust covariance estimator automatically account for correlation among repeated observations from the same person? Explain.
What does exponentiating a confidence interval for a log-linked slope allow you to interpret?