夹心的错误VCOV矩阵：：vcovHC，R中有观察计数

qyzbxkaa 于 2023-04-27 发布在其他

关注(0)|答案(1)|浏览(109)

我正在尝试估计具有重复观测的加权LM。自然地，可以利用此冗余信息来保存空间和内存。我正在开发具有400k观测的应用程序，其中许多是重复的，并且说LM（以及更复杂的方法）被多次调用，这就是为什么加权OLS等价于具有重复观测的OLS的原因。如果没有这种节省内存的技巧分析不能继续进行，因为这个鲁棒的VCOV被调用了数百次。一个很好的例子是Owen, 2017（“经验可能性的加权自协调优化”）。
模型的两个变体（未加权（m）和加权（m.w））中的点估计值是相同的。但是，sandwich::vcovHC不会将权重视为观测值权重进行处理。
MWE：

rm(list = ls())
x <- c(1, 1, 1, 1, 2:10) # Observation 1 has a count of 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) # Indentical
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) # Very different

这是理论上正确的VCOV计算方法：

# White's asymptotic sandwich formula manually:
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
# Same formula, with observation counts
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

问题可以在sandwich::meatHC函数中看到。以下代码手动计算肉（最简单的情况，对应于type = "HC0"）

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向量？
你真诚的安德烈

r

来源：https://stackoverflow.com/questions/76103997/wrong-vcov-matrix-from-sandwichvcovhc-with-observation-counts-in-r

1条答案

按热度按时间

dgsult0t1#

用例1：无观测加权，仅计数

V.good <- sandwich::vcovHC(m.w, type = "HC0", omega = resid(m.w)^2 * weights(m.w))
round(V.good, 2)
#            (Intercept) x.uniq
#(Intercept)        5.06  -0.53
#x.uniq            -0.53   0.11

这里唯一必要的改变是提供正确的omega向量。在这个应用程序中，它足以提供加权平方残差。我在网上找不到任何类似的例子;也许这可以作为一个例子添加到小插曲，因为Stata用户已经享受他们的[f/a/p]weights太久了，并且具有观察权重的OLS在应用统计学中非常常见。

案例2：计数和观测权重（如IPW）

现在假设有观察权重（抽样权重，逆方差权重，倾向得分等），除此之外，一些观察是重复的，这就是为什么用作lm()参数的最终权重是它们的乘积：

x <- c(1, 1, 1, 1, 2:10) # Observation 1 has a count of 4
y <- c(15, 15, 15, 15, 13, 36, 31, 37, 63, 66, 80, 91, 94)
iv <- 1 / x # If Var(U | X) ~ X, then, these weights are optimal
x.uniq <- x[-(1:3)]
y.uniq <- y[-(1:3)]
iv.uniq <- iv[-(1:3)]
w <- c(4, rep(1, 9))
mw   <- lm(y ~ x, weights = iv)
mwc <- lm(y.uniq ~ x.uniq, weights = iv.uniq * w)
coef(mw) - coef(mwc) # Indentical
round(sandwich::vcovHC(mw, type = "HC0"), 3)
#             (Intercept)      x
# (Intercept)       1.344 -0.189
# x                -0.189  0.151
round(sandwich::vcovHC(mwc, type = "HC0"), 3) # Very different
#             (Intercept) x.uniq
# (Intercept)       2.621 -0.389
# x.uniq           -0.389  0.183

旧的解决方案将不起作用，因为现在重量应该以不同的方式进入面包和肉基质：

round(sandwich::vcovHC(mwc, type = "HC0", omega = resid(mwc)^2 * weights(mwc)), 3)
#             (Intercept) x.uniq
# (Intercept)       3.148 -0.933
# x.uniq           -0.933  0.912

现在，面包矩阵应仅考虑观察计数，而肉矩阵应考虑计数和权重（在本例中为逆方差）：

n <- sum(w) # True number of observations, 13
mm <- model.matrix(mwc)
ef <- mm * resid(mwc) * iv.uniq * sqrt(w) # estfun
me <- crossprod(ef) / n # Meat matrix
max(abs(me - meat(mw))) # 7.1e-15, good
br <- solve(crossprod(mm * sqrt(weights(mwc))) / n) # Bread matrix
max(abs(br - bread(mw))) # 2.2e-16, good
V <- br %*% me %*% br / n
round(V, 3)
#             (Intercept) x.uniq
# (Intercept)       1.344 -0.189
# x.uniq           -0.189  0.151

赞(0）回复(0）举报 2023-04-27

我来回答

夹心的错误VCOV矩阵：：vcovHC，R中有观察计数

1条答案

用例1：无观测加权，仅计数

案例2：计数和观测权重（如IPW）

相关问题

热门标签

最新问答