如何解决三明治 + mlogit:当使用 `vcovHC()` 来计算稳健/集群标准错误时,`ef/X 中的错误:不一致的数组` > 稳健标准误差> 聚类标准误差为什么 vcovHC 对 mlogit 不起作用基本的“稳健”三明治协方差聚类协方差复制Stata结果
在离散选择问题中使用 mlogit()
拟合多项式 Logit (MNL) 后,我正在尝试计算稳健/集群标准误差。不幸的是,我怀疑我遇到了问题,因为我使用的是 long
格式的数据(在我的情况下这是必须的),并且在 #Error in ef/X : non-conformable arrays
之后收到错误 sandwich::vcovHC(,"HC0")
。>
数据
为了说明,请仔细考虑以下数据。它表示来自 5 个个体 (id_ind
) 的数据,这些个体在 3 个备选方案 (altern
) 中进行选择。五个人各选择三遍;因此我们有 15 种选择情况 (id_choice
)。每个备选方案由两个通用属性(x1
和 x2
)表示,并且选项注册在 y
中(1
如果选择,则为 0
否则)。
df <- read.table(header = TRUE,text = "
id_ind id_choice altern x1 x2 y
1 1 1 1 1.586788801 0.11887832 1
2 1 1 2 -0.937965347 1.15742493 0
3 1 1 3 -0.511504401 -1.90667519 0
4 1 2 1 1.079365680 -0.37267925 0
5 1 2 2 -0.009203032 1.65150370 1
6 1 2 3 0.870474033 -0.82558651 0
7 1 3 1 -0.638604013 -0.09459502 0
8 1 3 2 -0.071679538 1.56879334 0
9 1 3 3 0.398263302 1.45735788 1
10 2 4 1 0.291413453 -0.09107974 0
11 2 4 2 1.632831160 0.92925495 0
12 2 4 3 -1.193272276 0.77092623 1
13 2 5 1 1.967624379 -0.16373709 1
14 2 5 2 -0.479859282 -0.67042130 0
15 2 5 3 1.109780885 0.60348187 0
16 2 6 1 -0.025834772 -0.44004183 0
17 2 6 2 -1.255129594 1.10928280 0
18 2 6 3 1.309493274 1.84247199 1
19 3 7 1 1.593558740 -0.08952151 0
20 3 7 2 1.778701074 1.44483791 1
21 3 7 3 0.643191170 -0.24761157 0
22 3 8 1 1.738820924 -0.96793288 0
23 3 8 2 -1.151429915 -0.08581901 0
24 3 8 3 0.606695064 1.06524268 1
25 3 9 1 0.673866953 -0.26136206 0
26 3 9 2 1.176959443 0.85005871 1
27 3 9 3 -1.568225496 -0.40002252 0
28 4 10 1 0.516456176 -1.02081089 1
29 4 10 2 -1.752854918 -1.71728381 0
30 4 10 3 -1.176101700 -1.60213536 0
31 4 11 1 -1.497779616 -1.66301234 0
32 4 11 2 -0.931117325 1.50128532 1
33 4 11 3 -0.455543630 -0.64370825 0
34 4 12 1 0.894843784 -0.69859139 0
35 4 12 2 -0.354902281 1.02834859 0
36 4 12 3 1.283785176 -1.18923098 1
37 5 13 1 -1.293772990 -0.73491317 0
38 5 13 2 0.748091387 0.07453705 1
39 5 13 3 -0.463585127 0.64802031 0
40 5 14 1 -1.946438667 1.35776140 0
41 5 14 2 -0.470448172 -0.61326604 1
42 5 14 3 1.478763383 -0.66490028 0
43 5 15 1 0.588240775 0.84448489 1
44 5 15 2 1.131731049 -1.51323232 0
45 5 15 3 0.212145247 -1.01804594 0
")
问题
因此,我们可以使用 mlogit()
拟合 MNL 并提取它们的稳健方差-协方差,如下所示:
library(mlogit)
library(sandwich)
mo <- mlogit(formula = y ~ x1 + x2|0,method ="nr",data = df,idx = c("id_choice","altern"))
sandwich::vcovHC(mo,"HC0")
#Error in ef/X : non-conformable arrays
正如我们所见,sandwich::vcovHC
产生了一个错误,表示 ef/X
不符合标准。其中 X <- model.matrix(x)
和 ef <- estfun(x,...)
。在查看 the source code on the mirror on GitHub 后,我发现了一个问题,因为数据是长格式的,ef
的维度为 15 x 2
,X
的维度为 {{1 }}。
我的解决方法
鉴于节目必须继续进行,我正在使用从 sandwich
和 I adjusted to accommodate the Stata's output. 中借用的一些函数手动计算稳健和集群标准误差
> 稳健标准误差
这些行的灵感来自 sandwich::meat() 函数。
45 x 2
Stata 等价物
psi<- estfun(mo)
k <- NCOL(psi)
n <- NROW(psi)
rval <- (n/(n-1))* crossprod(as.matrix(psi))
vcov(mo) %*% rval %*% vcov(mo)
# x1 x2
# x1 0.23050261 0.09840356
# x2 0.09840356 0.12765662
> 聚类标准误差
这里,鉴于每个人回答了 3 个问题,个人之间很可能存在某种程度的相关性;因此在这种情况下应该首选集群校正。下面我计算了这种情况下的聚类校正,并展示了与 qui clogit y x1 x2,group(id_choice) r
mat li e(V)
symmetric e(V)[2,2]
y: y:
x1 x2
y:x1 .23050262
y:x2 .09840356 .12765662
的 Stata 输出的等效性。
clogit,cluster()
Stata 等效
id_ind_collapsed <- df$id_ind[!duplicated(mo$model$idx$id_choice,)]
psi_2 <- rowsum(psi,group = id_ind_collapsed )
k_cluster <- NCOL(psi_2)
n_cluster <- NROW(psi_2)
rval_cluster <- (n_cluster/(n_cluster-1))* crossprod(as.matrix(psi_2))
vcov(mo) %*% rval_cluster %*% vcov(mo)
# x1 x2
# x1 0.1766707 0.1007703
# x2 0.1007703 0.1180004
问题:
我想在 qui clogit y x1 x2,group(id_choice) cluster(id_ind)
symmetric e(V)[2,2]
y: y:
x1 x2
y:x1 .17667075
y:x2 .1007703 .11800038
生态系统中容纳我的计算,这意味着不是手动计算矩阵,而是实际使用 sandwich
函数。是否可以使其与此处描述的长格式模型一起使用?例如,直接提供 sandwich
和 meat
对象来执行计算?提前致谢。
PS:我注意到了 there is a dedicated bread
function in sandwich
for mlogit,但我无法发现 bread
的 meat
之类的东西,但无论如何我可能在这里遗漏了一些东西......
解决方法
为什么 vcovHC 对 mlogit 不起作用
HC 协方差估计器的类别只能应用于具有单个线性预测器的模型,其中评分函数又名估计函数是所谓的“工作残差”和回归矩阵的乘积。 Zeileis (2006) 论文(参见公式 7)对此进行了一些详细解释,在包中以 vignette("sandwich-OOP",package = "sandwich")
的形式提供。 ?vcovHC
也指出了这一点,但没有很好地解释。我现在在 http://sandwich.R-Forge.R-project.org/reference/vcovHC.html 的文档中对此进行了改进:
meatHC 函数是用于估计 HC 三明治估计器的肉的真正工作马 - 默认的 vcovHC 方法是调用三明治和面包的包装器。有关更多实现细节,请参阅 Zeileis (2006)。下面和 Zeileis (2004) 描述了线性回归模型的理论背景。其他类型的模型采用类似的公式,前提是它们依赖于单个线性预测器,并且估计函数可以表示为“工作残差”和回归向量的乘积(Zeileis 2006,方程 7)。
这意味着 vcovHC()
不适用于多项 logit 模型,因为它们通常对单独的响应类别使用单独的线性预测变量。同样,不支持两部分或跨栏模型等。
基本的“稳健”三明治协方差
通常,为了计算基本的 Eicker-Huber-White 三明治协方差矩阵估计量,最好的策略是使用 sandwich()
函数而不是 vcovHC()
函数。前者适用于具有 estfun()
和 bread()
方法的任何模型。
对于线性模型,sandwich(...,adjust = FALSE)
(默认)和 sandwich(...,adjust = TRUE)
分别对应于 HC0 和 HC1。在具有 n
观测值和 k
回归系数的模型中,前者使用 1/n
进行标准化,后者使用 1/(n-k)
进行标准化。
然而,Stata 在 logit 模型中除以 1/(n-1)
,参见:
Different Robust Standard Errors of Logit Regression in Stata and R。据我所知,没有明确的理论理由专门使用一种或另一种调整。并且已经在中等规模的样本中,这无论如何都没有区别。
备注: 1/(n-1)
的调整不能直接在 sandwich()
中用作选项。然而,巧合的是,它是 vcovCL()
中的默认值,没有指定 cluster
变量(即,将每个观察值视为一个单独的集群)。因此,如果您想获得与 Stata 完全相同的结果,这是一个方便的“技巧”。
聚类协方差
这可以通过vcovCL(...,cluster = ...)
“照常”计算。对于 mlogit
模型,您只需要考虑 cluster
变量只需要提供一次(而不是在长格式中多次堆叠)。
复制Stata结果
使用您帖子中的数据和模型:
vcovCL(mo)
## x1 x2
## x1 0.23050261 0.09840356
## x2 0.09840356 0.12765662
vcovCL(mo,cluster = df$id_choice[1:15])
## x1 x2
## x1 0.1766707 0.1007703
## x2 0.1007703 0.1180004
版权声明:本文内容由互联网用户自发贡献,该文观点与技术仅代表作者本人。本站仅提供信息存储空间服务,不拥有所有权,不承担相关法律责任。如发现本站有涉嫌侵权/违法违规的内容, 请发送邮件至 dio@foxmail.com 举报,一经查实,本站将立刻删除。