问题描述
在离散选择问题中使用 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 2Stata 等价物
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