泊松回归含offset时sandwich函数协方差矩阵不一致是否为预期行为？

问题背景 #

泊松回归既可以基于原始数据拟合，也可以先汇总数据后使用offset参数拟合。理论上两种方式应返回完全相同的系数估计及协方差估计，但实际使用R语言{sandwich}包的sandwich()函数时，发现带offset的分组模型与无offset的原始数据模型返回的协方差估计不一致，这会直接影响置信区间与p值。

复现代码 #

准备分组数据 #

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
)

还原为原始未分组数据 #

将分组数据中的成功次数（y）和失败次数（N - y）展开，得到等价于分组前的原始数据：

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

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

unit_data <- bind_rows(yes_outcomes, no_outcomes) %>% 
  select(-y) %>% 
  rename(y=yy)

拟合两个泊松回归模型 #

分别在原始数据和分组数据上拟合模型：

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

查看sandwich函数输出结果 #

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

疑问 #

请问这是预期行为吗？

内容的提问来源于stack exchange，提问作者Demetri Pananos