Robust standard errors for a Poisson regression with/without an offset

Question

The following is a question posted to Stack Overflow, but the answer is more to-do with statistical theory than software. I am reposting it here because this is a more appropriate venue, and because similar questions have remained unanswered on this site for 8+ years.

---

Poisson regression can be fit to raw data or can be fit by summarizing data and then using an offset. The model should return the same coefficient estimates (as well as covariance estimates), regardless of the approach (full data or offset).

However it seems the sandwich function from {sandwich} returns different covariance estimates when the model uses an offset. This is particularly troubling because this means the confidence intervals and p values are effected.

Shown below is such an example. The two covariance matrices differ when they are expected to be the same.

Is this intended behaviour?

library("tidyverse")
grouped_data <- tibble::tribble(
  ~treatment,        ~g,    ~y,     ~N,
  "A",      "a", 1338L, 20669L,
  "A",      "b",   36L,  1237L,
  "A",      "c", 2555L, 39438L,
  "A",      "d",  402L,  5713L,
  "B",      "a", 1281L, 19986L,
  "B",      "b",   38L,  1224L,
  "B",      "c", 2495L, 36749L,
  "B",      "d",  382L,  5646L
)

Re-organize data to be as if we hadn't grouped and summarized by uncounting the number of successes (y) and number of failures (N - y)

yes_outcomes <-grouped_data %>% 
  mutate(yy=1) %>% 
  uncount(y)

no_outcomes <-grouped_data %>% 
  mutate(yy=0) %>% 
  uncount(N-y)

# This is equivalent to the data we had before grouping
unit_data <- bind_rows(yes_outcomes, no_outcomes) %>% 
  select(-y) %>% 
  rename(y=yy)

Fit one model on the ungrouped data without offset and one model on the grouped data with offset.

fit <- glm(y ~ treatment + g, data = unit_data, family = poisson)
offset_fit <- glm(y ~ treatment + g, data = grouped_data, family = poisson,
  offset = log(N))

The corresponding sandwich results are:

sandwich::sandwich(fit)
#>               (Intercept)    treatmentB            gb            gc
#> (Intercept)  0.0004689434 -2.213726e-04 -3.599604e-04 -3.621573e-04
#> treatmentB  -0.0002213726  4.385850e-04  5.526098e-06  9.843754e-06
#> gb          -0.0003599604  5.526098e-06  1.346409e-02  3.572512e-04
#> gc          -0.0003621573  9.843754e-06  3.572512e-04  5.422671e-04
#> gd          -0.0003545652 -5.331337e-06  3.572202e-04  3.571526e-04
#>                        gd
#> (Intercept) -3.545652e-04
#> treatmentB  -5.331337e-06
#> gb           3.572202e-04
#> gc           3.571526e-04
#> gd           1.544801e-03
sandwich::sandwich(offset_fit)
#>               (Intercept)    treatmentB            gb            gc
#> (Intercept)  1.672258e-04 -8.616787e-05 -1.238920e-04 -1.263405e-04
#> treatmentB  -8.616787e-05  1.742927e-04 -1.483730e-06  3.456605e-06
#> gb          -1.238920e-04 -1.483730e-06  3.554277e-04  1.246024e-04
#> gc          -1.263405e-04  3.456605e-06  1.246024e-04  2.040329e-04
#> gd          -1.237209e-04 -1.827624e-06  1.246399e-04  1.245998e-04
#>                        gd
#> (Intercept) -1.237209e-04
#> treatmentB  -1.827624e-06
#> gb           1.246399e-04
#> gc           1.245998e-04
#> gd           5.878420e-04

^{Created on 2023-02-19 by the reprex package (v2.0.1)}

Answer 1

How could they possibly be different?

The two sandwich estimators are estimators for different semiparametric models. The individual-observation one is for a model where observations are all independent. The aggregated one is for a model where groups are all independent. Those are different assumptions. Since this is all frequentist, the standard errors describe how $\hat\beta$ would vary in replicates of the data set, and while for any specific data set the $\hat\beta$ is the same in the two models, the set of replicate data sets you need to consider is different.

If the mean model is correctly specified and also the observations really are independent, I think the two sandwich estimators are consistent for the same true variance. But if the mean model isn't correctly specified, it matters whether there's really just one observation per covariate pattern or many observations per covariate pattern. In ANOVA terms, the model with aggregated data has the interaction confounded with the residual error, and the model with individual observations doesn't.

Another way of saying the same thing is to think about bootstraps. The sandwich estimator can be seen as an approximation to the bootstrap. The sandwich estimator in the offset model is approximating a bootstrap that resamples whole groups; the one in the individual model is approximating a bootstrap that resamples individual observations.

We can verify that the grouping in the sandwich estimator is the issue, by using a clustered sandwich estimator on the individual-level data, thus assuming only cluster-level independence in both settings. The results are now the same (they aren't the same as either of the previous estimators, but they could be if I'd worked on choosing the same versions of the estimators: the offset_fit one differs from vcovHC(offset_fit, type="HC0") by a degrees-of-freedom factor of 1.14)

> unit_data$cluster<-with(unit_data, interaction(treatment,g))
> vcovCL(fit, cluster=as.numeric(unit_data$cluster))
            (Intercept) treatmentB        gb        gc        gd
(Intercept)    1.91e-04  -9.85e-05 -1.42e-04 -1.44e-04 -1.41e-04
treatmentB    -9.85e-05   1.99e-04 -1.70e-06  3.95e-06 -2.09e-06
gb            -1.42e-04  -1.70e-06  4.06e-04  1.42e-04  1.42e-04
gc            -1.44e-04   3.95e-06  1.42e-04  2.33e-04  1.42e-04
gd            -1.41e-04  -2.09e-06  1.42e-04  1.42e-04  6.72e-04
> vcovCL(offset_fit, cluster=NULL)
            (Intercept) treatmentB        gb        gc        gd
(Intercept)    1.91e-04  -9.85e-05 -1.42e-04 -1.44e-04 -1.41e-04
treatmentB    -9.85e-05   1.99e-04 -1.70e-06  3.95e-06 -2.09e-06
gb            -1.42e-04  -1.70e-06  4.06e-04  1.42e-04  1.42e-04
gc            -1.44e-04   3.95e-06  1.42e-04  2.33e-04  1.42e-04
gd            -1.41e-04  -2.09e-06  1.42e-04  1.42e-04  6.72e-04

So the issue genuinely is about what units are independent, and it's not a coding error

Who's right?

The individual-observation model is definitely right, because it actually describes how the data were collected. So either they are both right, or the individual-level model is right and the aggregate model is wrong.

I'm inclined to say that the aggregate model is wrong: that the hypothetical replicates it calculates the variance over are not the relevant set of hypothetical replicates.

An improvement?

The usual sandwich estimator estimates the variance of the score function using an empirical variance: it has to, the point is to not assume any distribution for $Y|X$. In this setting, though, we have a distribution for the individual $Y$; they're Bernoulli. It wouldn't be unreasonable to calculate the variance assuming the mean model is correct; that's what we'd routinely do for Poisson regression or logistic regression.

So it would be possible to work out an estimator for the middle of the sandwich (the variability matrix) that's based on the Bernoulli distribution of $Y$, while keeping the usual estimator for the outside of the sandwich (the sensitivity matrix). In maximum likelihood estimation this is kind of pointless, because the sensitivity matrix is just the inverse of the variability matrix and the whole thing collapses back to the model-based estimate. For a binomial model and the working Poisson estimator it could make more sense. Further Research Is Needed.