R中sandwich::vcovHC处理观测计数时生成错误VCOV矩阵问题

我在处理含大量重复观测的线性模型时，为节省内存将观测计数作为权重（加权OLS等价于重复观测的OLS），数据量达40万条且模型需多次调用。未加权模型m与加权模型m.w的点估计完全一致，但sandwich::vcovHC生成的方差-协方差（VCOV）矩阵差异显著，显然未正确将权重视为观测重复计数。

最小可复现示例 #

rm(list = ls())
x <- c(1, 1, 1, 1, 2:10) # 观测1重复4次
y <- c(15, 15, 15, 15, 13, 36, 31, 37, 63, 66, 80, 91, 94)
x.uniq <- x[-(1:3)]
y.uniq <- y[-(1:3)]
w <- c(4, rep(1, 9))
m <- lm(y ~ x)
m.w <- lm(y.uniq ~ x.uniq, weights = w)

coef(m) - coef(m.w) # 点估计完全一致

round(sandwich::vcovHC(m, type = "HC0"), 2)
#            (Intercept)     x
#(Intercept)        5.06 -0.53
#x                 -0.53  0.11

round(sandwich::vcovHC(m.w, type = "HC0"), 2) # 结果差异极大

理论正确的VCOV计算 #

手动实现White渐近三明治公式，无论是否按重复计数加权，结果均与未加权模型的vcovHC输出一致：

# 未加权版本
br <- solve(crossprod(model.matrix(m)))
me <- crossprod(model.matrix(m) * resid(m))
V <- br %*% me %*% br
round(V, 2)
#            (Intercept) x.uniq
#(Intercept)        5.06  -0.53
#x.uniq            -0.53   0.11

# 加权版本（正确处理重复计数）
br.w <- solve(crossprod(model.matrix(m.w) * sqrt(weights(m.w))))
me.w <- crossprod(model.matrix(m.w) * resid(m.w) * sqrt(weights(m.w)))
V.w <- br.w %*% me.w %*% br.w
round(V.w, 2)
#            (Intercept) x.uniq
#(Intercept)        5.06  -0.53
#x.uniq            -0.53   0.11

问题根源 #

通过手动计算meat部分验证，问题出在sandwich::meatHC函数的处理逻辑：

X <- model.matrix(m)
X.w <- model.matrix(m.w)
res <- rowMeans(estfun(m)/X)
res.w <- rowMeans(estfun(m.w)/X.w)
omega <- res^2
omega.w <- res.w^2

meat <- crossprod(sqrt(omega) * X) / nrow(X)
meat.w <- crossprod(sqrt(omega.w) * X.w) / nrow(X.w)

all.equal(meat, meatHC(m, type = "HC0")) # TRUE
all.equal(meat.w, meatHC(m.w, type = "HC0")) # TRUE

疑问 #

是否遗漏了传递给estfun的参数？或者需要自定义omega向量来修正加权模型的meat计算？

内容的提问来源于stack exchange，提问作者André V. Kostyrka