如何使用多项 Logit 模型的标准误差获得平均边际效应 (AME)？

我想要获得具有标准误差的多项 Logit 模型的平均边际效应 (AME)。为此我尝试了不同的方法，但到目前为止还没有达到目标。

最好的尝试

我最好的尝试是使用手动获取 AMEmlogit我在下面展示了这一点。

library(mlogit)
ml.d <- mlogit.data(df1, choice="Y", shape="wide")  # shape data for `mlogit()`
ml.fit <- mlogit(Y ~ 1 | D + x1 + x2, reflevel="1", data=ml.d)  # fit the model

# coefficient names
c.names <- all.vars(ml.fit$call)[2:4]

# get marginal effects
ME.mnl <- sapply(c.names, function(x) 
  stats::effects(ml.fit, covariate=x, data=ml.d), 
  simplify=FALSE) 

# get AMEs
(AME.mnl <- t(sapply(ME.mnl, colMeans)))
#              1            2            3           4          5
# D  -0.03027080 -0.008806072 0.0015410569 0.017186531 0.02034928
# x1 -0.02913234 -0.015749598 0.0130577842 0.013240212 0.01858394
# x2 -0.02724650 -0.005482753 0.0008575982 0.005331181 0.02654047

我知道这些价值观是正确的。但是，我无法通过简单地计算列的标准差来获得正确的标准误差：

# standard errors - WRONG!
(AME.mnl.se <- t(sapply(E.mnl, colSdColMeans)))

(Note: colSdColMeans()提供了色谱柱的 SDhere https://stackoverflow.com/a/17551600/6574038.)

因此这也导致我错误t-values:

# t values - WRONG!
AME.mnl / AME.mnl.se
#             1          2          3         4         5
# D  -0.7110537 -0.1615635 0.04013228 0.4190057 0.8951484
# x1 -0.7170813 -0.2765212 0.33325968 0.3656893 0.8907836
# x2 -0.7084573 -0.1155825 0.02600653 0.1281190 0.8559794

虽然我知道正确的t- 这种情况下的值如下：

# D  -9.26 -1.84  0.31 4.29 8.05   
# x1 -6.66 -2.48  1.60 1.50 3.22  
# x2 -2.95 -0.39  0.06 0.42 3.21 

我了解到应该有一个“增量方法”，但我只找到了一些用于非常特殊情况的代码，其交互位于交叉验证 https://stats.stackexchange.com/a/122537/163114.

失败的尝试

1.) 包装margins似乎无法处理"mlogit"对象：

library(margins)
summary(margins(ml.fit))

2.) 还有另一个 mlogits 包，nnet,

library(nnet) 
ml.fit2 <- multinom(Y ~ D + x1 + x2, data=df1)
summary(ml.fit2)

but margins也无法正确处理这个问题：

> summary(margins(ml.fit2))
 factor     AME SE  z  p lower upper
      D -0.0303 NA NA NA    NA    NA
     x1 -0.0291 NA NA NA    NA    NA
     x2 -0.0272 NA NA NA    NA    NA

3.）还有一个软件包声称可以计算“多项 Logistic 回归模型的平均效应”,

library(DAMisc)
mnlChange2(ml.fit2, varnames="D", data=df1)

但我无法从中得到一滴牛奶，因为该函数什么也没有产生（即使没有该函数的示例）。

现在我们如何获得具有标准误差/多项 Logit 模型的 t 统计量的 AMER?

Data

df1 <- structure(list(Y = c(3, 4, 1, 2, 3, 4, 1, 5, 2, 3, 4, 2, 1, 4, 
1, 5, 3, 3, 3, 5, 5, 4, 3, 5, 4, 2, 5, 4, 3, 2, 5, 3, 2, 5, 5, 
4, 5, 1, 2, 4, 3, 1, 2, 3, 1, 1, 3, 2, 4, 2, 2, 4, 1, 5, 3, 1, 
5, 2, 3, 4, 2, 4, 5, 2, 4, 1, 4, 2, 1, 5, 3, 2, 1, 4, 4, 1, 5, 
1, 1, 1, 4, 5, 5, 3, 2, 3, 3, 2, 4, 4, 5, 3, 5, 1, 2, 5, 5, 1, 
2, 3), D = c(12, 8, 6, 11, 5, 14, 0, 22, 15, 13, 18, 3, 5, 9, 
10, 28, 9, 16, 17, 14, 26, 18, 18, 23, 23, 12, 28, 14, 10, 15, 
26, 9, 2, 30, 18, 24, 27, 7, 6, 25, 13, 8, 4, 16, 1, 4, 5, 18, 
21, 1, 2, 19, 4, 2, 16, 17, 23, 15, 13, 21, 24, 14, 27, 6, 20, 
6, 19, 8, 7, 23, 11, 11, 1, 22, 21, 4, 27, 6, 2, 9, 18, 30, 26, 
22, 10, 1, 4, 7, 26, 15, 26, 18, 30, 1, 11, 29, 25, 3, 19, 15
), x1 = c(13, 12, 4, 3, 16, 16, 15, 13, 1, 15, 10, 16, 1, 17, 
7, 13, 12, 6, 8, 16, 16, 11, 7, 16, 5, 13, 12, 16, 17, 6, 16, 
9, 14, 16, 15, 5, 7, 2, 8, 2, 9, 9, 15, 13, 9, 4, 16, 2, 11, 
13, 11, 6, 4, 3, 7, 4, 12, 2, 16, 14, 3, 13, 10, 11, 10, 4, 11, 
16, 8, 12, 14, 9, 4, 16, 16, 12, 9, 10, 6, 1, 3, 8, 7, 7, 5, 
16, 17, 10, 4, 15, 10, 8, 3, 13, 9, 16, 12, 7, 4, 11), x2 = c(12, 
19, 18, 19, 15, 12, 15, 16, 15, 11, 12, 16, 17, 14, 12, 17, 17, 
16, 12, 20, 11, 11, 15, 14, 18, 10, 14, 13, 10, 14, 18, 18, 18, 
17, 18, 14, 16, 19, 18, 16, 18, 14, 17, 10, 16, 12, 16, 15, 11, 
18, 19, 15, 19, 11, 16, 10, 20, 14, 10, 12, 10, 15, 13, 15, 11, 
20, 11, 12, 16, 16, 11, 15, 11, 11, 10, 10, 16, 11, 20, 17, 20, 
17, 16, 11, 18, 19, 18, 14, 17, 11, 16, 11, 18, 14, 15, 16, 11, 
14, 11, 13)), class = "data.frame", row.names = c(NA, -100L))

我们可以做一些与您链接的答案中所做的非常相似的事情。特别是，首先我们需要一个能够计算给定系数向量的 AME 的函数。为此我们可以定义

AME.fun <- function(betas) {
  tmp <- ml.fit
  tmp$coefficients <- betas
  ME.mnl <- sapply(c.names, function(x) 
    effects(tmp, covariate = x, data = ml.d), simplify = FALSE)
  c(sapply(ME.mnl, colMeans))
}

后半部分是你的，而在前半部分我用了一个技巧来采取同样的方法ml.fit对象并更改其系数。接下来我们找到雅可比

require(numDeriv)
grad <- jacobian(AME.fun, ml.fit$coef)

并应用delta方法。对角线的平方根grad %*% vcov(ml.fit) %*% t(grad)就是我们想要的。因此，

(AME.mnl.se <- matrix(sqrt(diag(grad %*% vcov(ml.fit) %*% t(grad))), nrow = 3, byrow = TRUE))
#             [,1]        [,2]        [,3]        [,4]        [,5]
# [1,] 0.003269320 0.004788536 0.004995723 0.004009762 0.002527462
# [2,] 0.004375795 0.006348496 0.008168883 0.008844684 0.005763966
# [3,] 0.009233616 0.014048212 0.014713090 0.012702188 0.008261734
AME.mnl / AME.mnl.se
#            1          2          3         4        5
# D  -9.259050 -1.8389907 0.30847523 4.2861720 8.051269
# x1 -6.657611 -2.4808393 1.59847852 1.4969683 3.224159
# x2 -2.950794 -0.3902812 0.05828811 0.4197057 3.212458

与Stata的结果一致。