I think the default vcovHC in R's sandwich package does not handle offsets in Poisson models. We see this because robust (heteroscedasticity consistent) standard errors, using the sandwich package, are vastly different. I believe that the formulation of the bread and meat of the regression model are not correctly specified in this case. The two different variance-covariance matrices (model based versus robust) are ludicrously different.
When we generate data from the following relative risk model, the Poisson model for rare outcomes should be correct still.
set.seed(1)
library(sandwich)
x <- 0:1
n <- rpois(2, exp(10 + .4*x))
y <- rbinom(2, n, exp(-4 + .2*x))
model <- glm(y ~ x + offset(log(n)), family=poisson)
vcov(model)
Gives the following information:
> vcov(model)
(Intercept) x
(Intercept) 0.002415459 -0.002415459
x -0.002415459 0.003779715
> vcovHC(model, type='HC0')
(Intercept) x
(Intercept) 4.713025e-31 -4.713025e-31
x -4.713025e-31 4.953579e-31
vcovHC(model, type='HC0')
However, expanding the array to produce equivalent regression without the use of offsets produces different robust error estimates. By my understanding of offsets, all aspects of these two analyses should be the same:
xlong <- rep.int(x, times=n)
ylong <- rep.int(c(1,0,1,0), times=c(y[1], n[1]-y[1], y[2], n[2]-y[2]))
modellong <- glm(ylong ~ xlong, family=poisson)
and doing the same
> vcov(modellong)
(Intercept) xlong
(Intercept) 0.002415402 -0.002415402
xlong -0.002415402 0.003779643
> vcovHC(modellong, type='HC0')
(Intercept) xlong
(Intercept) 0.002369866 -0.002369866
xlong -0.002369866 0.003703911
