如何在 mlogit 中获得分类变量的边际效应？
How to get marginal effects for categorical variables in mlogit?
我想计算 "mlogit" 对象的边际效应，其中解释变量是分类变量（因子）。对于数值数据 effects() 会抛出一些东西，而对于分类数据则不会。
为简单起见，我在下面展示了一个双变量示例。
数值变量
# with mlogit
library(mlogit)
ml.dat <- mlogit.data(df3, choice="y", shape="wide")
fit.mnl <- mlogit(y ~ 1 | x, data=ml.dat)
head(effects(fit.mnl, covariate="x", data=ml.dat))
# FALSE TRUE
# 1 -0.01534581 0.01534581
# 2 -0.01534581 0.01534581
# 3 -0.20629452 0.20629452
# 4 -0.06903946 0.06903946
# 5 -0.24174312 0.24174312
# 6 -0.39306240 0.39306240
# with glm
fit.glm <- glm(y ~ x, df3, family = binomial)
head(effects(fit.glm))
# (Intercept) x
# -0.2992979 -4.8449254 2.3394989 0.2020127 0.4616640 1.0499595
因子变量
# transform to factor
df3F <- within(df3, x <- factor(x))
class(df3F$x) == "factor"
# [1] TRUE
虽然 glm() 仍然抛出一些东西，
# with glm
fit.glmF <- glm(y ~ x, df3F, family = binomial)
head(effects(fit.glmF))
# (Intercept) x2 x3 x4 x5 x6
# 0.115076511 -0.002568206 -0.002568206 -0.003145397 -0.003631992 -0.006290794
mlogit() 方法
# with mlogit
ml.datF <- mlogit.data(df3F, choice="y", shape="wide")
fit.mnlF <- mlogit(y ~ 1 | x, data=ml.datF)
head(effects(fit.mnlF, covariate="x", data=ml.datF))
抛出此 错误:
Error in `contrasts<-`(`*tmp*`, value = contr.funs[1 + isOF[nn]]) :
contrasts can be applied only to factors with 2 or more levels
In addition: Warning message:
In Ops.factor(data[, covariate], eps) :
Error in `contrasts<-`(`*tmp*`, value = contr.funs[1 + isOF[nn]]) :
contrasts can be applied only to factors with 2 or more levels
我该如何解决这个问题？
我已经尝试过用 this answer 来操纵 effects.mlogit() 但它并没有帮助解决我的问题。
注意：这个问题与有关，我想将其应用于分类解释变量。
编辑
（在将给定解决方案应用于与上面链接的问题相关的潜在问题时演示问题。请参阅评论。）
# new example ----
library(mlogit)
ml.d <- mlogit.data(df1, choice="y", shape="wide")
ml.fit <- mlogit(y ~ 1 | factor(x), reflevel="1", data=ml.d)
AME.fun2 <- function(betas) {
aux <- model.matrix(y ~ x, df1)[, -1]
ml.datF <- mlogit.data(data.frame(y=df1$y, aux),
choice="y", shape="wide")
frml <- mFormula(formula(paste("y ~ 1 |", paste(colnames(aux),
collapse=" + "))))
fit.mnlF <- mlogit(frml, data=ml.datF)
fit.mnlF$coefficients <- betas # probably?
colMeans(effects(fit.mnlF, covariate="x2", data=ml.datF)) # first co-factor?
}
(AME.mnl <- AME.fun2(ml.fit$coefficients))
require(numDeriv)
grad <- jacobian(AME.fun2, ml.fit$coef)
(AME.mnl.se <- matrix(sqrt(diag(grad %*% vcov(ml.fit) %*% t(grad))),
nrow=3, byrow=TRUE))
AME.mnl / AME.mnl.se
# doesn't work yet though...
# probably "true" values, obtained from Stata:
# # ame
# 1 2 3 4 5
# 1. NA NA NA NA NA
# 2. -0.400 0.121 0.0971 0.113 0.0686
# 3. -0.500 -0.179 0.0390 0.166 0.474
#
# # z-values
# 1 2 3 4 5
# 1. NA NA NA NA NA
# 2. -3.86 1.25 1.08 1.36 0.99
# 3. -5.29 -2.47 0.37 1.49 4.06
数据
df3 <- structure(list(x = c(11, 11, 7, 10, 9, 8, 9, 6, 9, 9, 8, 9, 11,
7, 8, 11, 12, 5, 8, 8, 11, 6, 13, 12, 5, 8, 7, 11, 8, 10, 9,
10, 7, 9, 2, 10, 3, 6, 11, 9, 7, 8, 4, 12, 8, 12, 11, 9, 12,
9, 7, 7, 7, 10, 4, 10, 9, 6, 7, 8, 9, 13, 10, 8, 10, 6, 7, 10,
9, 6, 4, 6, 6, 8, 6, 9, 3, 7, 8, 2, 8, 6, 7, 9, 10, 8, 6, 5,
5, 7, 9, 1, 6, 11, 11, 9, 7, 8, 9, 9), y = c(TRUE, TRUE, TRUE,
TRUE, TRUE, TRUE, TRUE, FALSE, FALSE, TRUE, FALSE, TRUE, TRUE,
TRUE, FALSE, TRUE, TRUE, FALSE, TRUE, TRUE, TRUE, FALSE, TRUE,
TRUE, FALSE, TRUE, FALSE, TRUE, FALSE, TRUE, TRUE, TRUE, FALSE,
TRUE, FALSE, TRUE, FALSE, FALSE, TRUE, TRUE, TRUE, FALSE, FALSE,
TRUE, FALSE, TRUE, TRUE, FALSE, TRUE, TRUE, TRUE, FALSE, FALSE,
TRUE, FALSE, TRUE, TRUE, FALSE, FALSE, TRUE, TRUE, TRUE, TRUE,
FALSE, TRUE, FALSE, FALSE, TRUE, FALSE, FALSE, FALSE, FALSE,
FALSE, FALSE, FALSE, TRUE, FALSE, FALSE, FALSE, FALSE, TRUE,
FALSE, FALSE, TRUE, TRUE, FALSE, FALSE, FALSE, FALSE, FALSE,
TRUE, FALSE, FALSE, TRUE, TRUE, TRUE, FALSE, FALSE, TRUE, FALSE
)), class = "data.frame", row.names = c(NA, -100L))
> summary(df3)
x y
Min. : 1.00 Mode :logical
1st Qu.: 7.00 FALSE:48
Median : 8.00 TRUE :52
Mean : 8.08
3rd Qu.:10.00
Max. :13.00
df1 <- structure(list(y = c(5, 4, 2, 2, 2, 3, 5, 4, 1, 1, 2, 4, 1, 4,
5, 5, 2, 3, 3, 5, 5, 3, 2, 4, 5, 1, 3, 3, 4, 3, 5, 2, 4, 4, 5,
5, 5, 2, 1, 5, 1, 3, 1, 4, 1, 2, 2, 4, 3, 1, 4, 3, 1, 1, 5, 2,
5, 4, 2, 2, 4, 2, 3, 5, 4, 1, 2, 2, 3, 5, 2, 5, 3, 3, 3, 1, 3,
1, 1, 4, 3, 4, 5, 2, 1, 1, 3, 1, 5, 4, 4, 2, 5, 3, 4, 4, 3, 1,
5, 2), x = structure(c(2L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 2L, 2L,
2L, 1L, 1L, 2L, 2L, 3L, 2L, 2L, 2L, 2L, 3L, 2L, 2L, 3L, 3L, 2L,
3L, 2L, 2L, 2L, 3L, 2L, 1L, 3L, 2L, 3L, 3L, 1L, 1L, 3L, 2L, 2L,
1L, 2L, 1L, 1L, 1L, 2L, 2L, 1L, 1L, 2L, 1L, 1L, 2L, 2L, 3L, 2L,
2L, 2L, 3L, 2L, 3L, 1L, 2L, 1L, 2L, 2L, 1L, 3L, 2L, 2L, 1L, 2L,
2L, 1L, 3L, 1L, 1L, 2L, 2L, 3L, 3L, 2L, 2L, 1L, 1L, 1L, 3L, 2L,
3L, 2L, 3L, 1L, 2L, 3L, 3L, 1L, 2L, 2L), .Label = c("1", "2",
"3"), class = "factor")), row.names = c(NA, -100L), class = "data.frame")
预计 effects 不适用于因子，否则输出将包含另一个维度，使结果有些复杂，这是非常合理的，就像我下面的解决方案一样，人们可能只希望对某个因素水平而不是所有水平产生影响。此外，正如我在下面解释的那样，分类变量的边际效应并不是唯一定义的，因此这将是 effects.
的一个额外并发症。
一个自然的解决方法是手动将因子变量转换为一系列虚拟变量，如
aux <- model.matrix(y ~ x, df3F)[, -1]
head(aux)
# x2 x3 x4 x5 x6 x7 x8 x9 x10 x11 x12 x13
# 1 0 0 0 0 0 0 0 0 0 1 0 0
# 2 0 0 0 0 0 0 0 0 0 1 0 0
# 3 0 0 0 0 0 1 0 0 0 0 0 0
# 4 0 0 0 0 0 0 0 0 1 0 0 0
# 5 0 0 0 0 0 0 0 1 0 0 0 0
# 6 0 0 0 0 0 0 1 0 0 0 0 0
这样数据就是
ml.datF <- mlogit.data(data.frame(y = df3F$y, aux), choice = "y", shape = "wide")
我们还需要手动构造公式
frml <- mFormula(formula(paste("y ~ 1 |", paste(colnames(aux), collapse = " + "))))
到目前为止一切顺利。现在如果我们 运行
fit.mnlF <- mlogit(frml, data = ml.datF)
head(effects(fit.mnlF, covariate = "x2", data = ml.datF))
# FALSE TRUE
# 1 -1.618544e-15 0.000000e+00
# 2 -1.618544e-15 0.000000e+00
# 3 -7.220891e-08 7.221446e-08
# 4 -1.618544e-15 0.000000e+00
# 5 -5.881129e-08 5.880851e-08
# 6 -8.293366e-08 8.293366e-08
那么结果不正确。 effects 在这里所做的是将 x2 视为一个 连续 变量，并计算这些案例的通常边际效应。即，如果对应于 x2 的系数是 b2 而我们的模型是 f(x,b2)，effects 计算 f 关于 b2 的导数并在每个观测向量 x我。这是错误的，因为 x2 只取值 0 和 1，而不是 0 或 1 左右的值，这是取导数的假设（极限的概念）！例如，考虑您的其他数据集 df1。在那种情况下，我们错误地得到
colMeans(effects(fit.mnlF, covariate = "x2", data = ml.datF))
# 1 2 3 4 5
# -0.25258378 0.07364406 0.05336283 0.07893391 0.04664298
这里有另一种方法（使用导数近似）得到这个不正确的结果：
temp <- ml.datF
temp$x2 <- temp$x2 + 0.0001
colMeans(predict(fit.mnlF, newdata = temp, type = "probabilities") -
predict(fit.mnlF, newdata = ml.datF, type = "probabilities")) / 0.0001
# 1 2 3 4 5
# -0.25257597 0.07364089 0.05336032 0.07893273 0.04664202
我没有使用 effects，而是通过两次使用 predict 手动计算了错误的边际效应：结果是平均值（{x2new = x2old + 0.0001 的拟合概率} - {拟合概率x2new = x2old}) / 0.0001。也就是说，我们通过将 x2 向上移动 0.0001（从 0 到 0.0001 或从 1 到 0.0001）来查看预测概率的变化。这两个都没有意义。当然，我们不应该期望 effects 有任何其他内容，因为数据中的 x2 是数字。
那么接下来的问题就是如何计算正确的（平均）边际效应。正如我所说，分类变量的边际效应并不是唯一定义的。假设 x_i 是个人 i 是否有工作，y_i 是他们是否有车。那么，那么至少有以下六点需要考虑。
- 从 x_i=0 到 x_i=1 时对 y_i = 1 的概率的影响。
- 从x_i=0 到x_i（观测值）时。
- 从 x_i 到 1.
现在，当我们对平均边际效应感兴趣时，我们可能只想对 1-3 的变化产生影响的那些个体进行平均。也就是说，
- 如果观测值不是 1，则从 x_i=0 到 x_i=1。
- 如果观测值不为 0，则从 x_i=0 到 x_i。
- 如果观测值不是 1，则从 x_i 到 1。
根据您的结果，Stata 使用选项 5，因此我将重现相同的结果，但实施任何其他选项都很简单，我建议考虑哪些选项对您的特定应用程序感兴趣。
AME.fun2 <- function(betas) {
aux <- model.matrix(y ~ x, df1)[, -1]
ml.datF <- mlogit.data(data.frame(y = df1$y, aux), choice="y", shape="wide")
frml <- mFormula(formula(paste("y ~ 1 |", paste(colnames(aux), collapse=" + "))))
fit.mnlF <- mlogit(frml, data = ml.datF)
fit.mnlF$coefficients <- betas
aux <- ml.datF # Auxiliary dataset
aux$x3 <- 0 # Going from 0 to the observed x_i
idx <- unique(aux[aux$x3 != ml.datF$x3, "chid"]) # Where does it make a change?
actual <- predict(fit.mnlF, newdata = ml.datF)
counterfactual <- predict(fit.mnlF, newdata = aux)
colMeans(actual[idx, ] - counterfactual[idx, ])
}
(AME.mnl <- AME.fun2(ml.fit$coefficients))
# 1 2 3 4 5
# -0.50000000 -0.17857142 0.03896104 0.16558441 0.47402597
require(numDeriv)
grad <- jacobian(AME.fun2, ml.fit$coef)
AME.mnl.se <- matrix(sqrt(diag(grad %*% vcov(ml.fit) %*% t(grad))), nrow = 1, byrow = TRUE)
AME.mnl / AME.mnl.se
# [,1] [,2] [,3] [,4] [,5]
# [1,] -5.291503 -2.467176 0.36922 1.485058 4.058994
我想计算 "mlogit" 对象的边际效应，其中解释变量是分类变量（因子）。对于数值数据 effects() 会抛出一些东西，而对于分类数据则不会。
为简单起见，我在下面展示了一个双变量示例。
数值变量
# with mlogit
library(mlogit)
ml.dat <- mlogit.data(df3, choice="y", shape="wide")
fit.mnl <- mlogit(y ~ 1 | x, data=ml.dat)
head(effects(fit.mnl, covariate="x", data=ml.dat))
# FALSE TRUE
# 1 -0.01534581 0.01534581
# 2 -0.01534581 0.01534581
# 3 -0.20629452 0.20629452
# 4 -0.06903946 0.06903946
# 5 -0.24174312 0.24174312
# 6 -0.39306240 0.39306240
# with glm
fit.glm <- glm(y ~ x, df3, family = binomial)
head(effects(fit.glm))
# (Intercept) x
# -0.2992979 -4.8449254 2.3394989 0.2020127 0.4616640 1.0499595
因子变量
# transform to factor
df3F <- within(df3, x <- factor(x))
class(df3F$x) == "factor"
# [1] TRUE
虽然 glm() 仍然抛出一些东西，
# with glm
fit.glmF <- glm(y ~ x, df3F, family = binomial)
head(effects(fit.glmF))
# (Intercept) x2 x3 x4 x5 x6
# 0.115076511 -0.002568206 -0.002568206 -0.003145397 -0.003631992 -0.006290794
mlogit() 方法
# with mlogit
ml.datF <- mlogit.data(df3F, choice="y", shape="wide")
fit.mnlF <- mlogit(y ~ 1 | x, data=ml.datF)
head(effects(fit.mnlF, covariate="x", data=ml.datF))
抛出此 错误:
Error in `contrasts<-`(`*tmp*`, value = contr.funs[1 + isOF[nn]]) :
contrasts can be applied only to factors with 2 or more levels
In addition: Warning message:
In Ops.factor(data[, covariate], eps) :
Error in `contrasts<-`(`*tmp*`, value = contr.funs[1 + isOF[nn]]) :
contrasts can be applied only to factors with 2 or more levels
我该如何解决这个问题？
我已经尝试过用 this answer 来操纵 effects.mlogit() 但它并没有帮助解决我的问题。
注意：这个问题与
编辑
（在将给定解决方案应用于与上面链接的问题相关的潜在问题时演示问题。请参阅评论。）
# new example ----
library(mlogit)
ml.d <- mlogit.data(df1, choice="y", shape="wide")
ml.fit <- mlogit(y ~ 1 | factor(x), reflevel="1", data=ml.d)
AME.fun2 <- function(betas) {
aux <- model.matrix(y ~ x, df1)[, -1]
ml.datF <- mlogit.data(data.frame(y=df1$y, aux),
choice="y", shape="wide")
frml <- mFormula(formula(paste("y ~ 1 |", paste(colnames(aux),
collapse=" + "))))
fit.mnlF <- mlogit(frml, data=ml.datF)
fit.mnlF$coefficients <- betas # probably?
colMeans(effects(fit.mnlF, covariate="x2", data=ml.datF)) # first co-factor?
}
(AME.mnl <- AME.fun2(ml.fit$coefficients))
require(numDeriv)
grad <- jacobian(AME.fun2, ml.fit$coef)
(AME.mnl.se <- matrix(sqrt(diag(grad %*% vcov(ml.fit) %*% t(grad))),
nrow=3, byrow=TRUE))
AME.mnl / AME.mnl.se
# doesn't work yet though...
# probably "true" values, obtained from Stata:
# # ame
# 1 2 3 4 5
# 1. NA NA NA NA NA
# 2. -0.400 0.121 0.0971 0.113 0.0686
# 3. -0.500 -0.179 0.0390 0.166 0.474
#
# # z-values
# 1 2 3 4 5
# 1. NA NA NA NA NA
# 2. -3.86 1.25 1.08 1.36 0.99
# 3. -5.29 -2.47 0.37 1.49 4.06
数据
df3 <- structure(list(x = c(11, 11, 7, 10, 9, 8, 9, 6, 9, 9, 8, 9, 11,
7, 8, 11, 12, 5, 8, 8, 11, 6, 13, 12, 5, 8, 7, 11, 8, 10, 9,
10, 7, 9, 2, 10, 3, 6, 11, 9, 7, 8, 4, 12, 8, 12, 11, 9, 12,
9, 7, 7, 7, 10, 4, 10, 9, 6, 7, 8, 9, 13, 10, 8, 10, 6, 7, 10,
9, 6, 4, 6, 6, 8, 6, 9, 3, 7, 8, 2, 8, 6, 7, 9, 10, 8, 6, 5,
5, 7, 9, 1, 6, 11, 11, 9, 7, 8, 9, 9), y = c(TRUE, TRUE, TRUE,
TRUE, TRUE, TRUE, TRUE, FALSE, FALSE, TRUE, FALSE, TRUE, TRUE,
TRUE, FALSE, TRUE, TRUE, FALSE, TRUE, TRUE, TRUE, FALSE, TRUE,
TRUE, FALSE, TRUE, FALSE, TRUE, FALSE, TRUE, TRUE, TRUE, FALSE,
TRUE, FALSE, TRUE, FALSE, FALSE, TRUE, TRUE, TRUE, FALSE, FALSE,
TRUE, FALSE, TRUE, TRUE, FALSE, TRUE, TRUE, TRUE, FALSE, FALSE,
TRUE, FALSE, TRUE, TRUE, FALSE, FALSE, TRUE, TRUE, TRUE, TRUE,
FALSE, TRUE, FALSE, FALSE, TRUE, FALSE, FALSE, FALSE, FALSE,
FALSE, FALSE, FALSE, TRUE, FALSE, FALSE, FALSE, FALSE, TRUE,
FALSE, FALSE, TRUE, TRUE, FALSE, FALSE, FALSE, FALSE, FALSE,
TRUE, FALSE, FALSE, TRUE, TRUE, TRUE, FALSE, FALSE, TRUE, FALSE
)), class = "data.frame", row.names = c(NA, -100L))
> summary(df3)
x y
Min. : 1.00 Mode :logical
1st Qu.: 7.00 FALSE:48
Median : 8.00 TRUE :52
Mean : 8.08
3rd Qu.:10.00
Max. :13.00
df1 <- structure(list(y = c(5, 4, 2, 2, 2, 3, 5, 4, 1, 1, 2, 4, 1, 4,
5, 5, 2, 3, 3, 5, 5, 3, 2, 4, 5, 1, 3, 3, 4, 3, 5, 2, 4, 4, 5,
5, 5, 2, 1, 5, 1, 3, 1, 4, 1, 2, 2, 4, 3, 1, 4, 3, 1, 1, 5, 2,
5, 4, 2, 2, 4, 2, 3, 5, 4, 1, 2, 2, 3, 5, 2, 5, 3, 3, 3, 1, 3,
1, 1, 4, 3, 4, 5, 2, 1, 1, 3, 1, 5, 4, 4, 2, 5, 3, 4, 4, 3, 1,
5, 2), x = structure(c(2L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 2L, 2L,
2L, 1L, 1L, 2L, 2L, 3L, 2L, 2L, 2L, 2L, 3L, 2L, 2L, 3L, 3L, 2L,
3L, 2L, 2L, 2L, 3L, 2L, 1L, 3L, 2L, 3L, 3L, 1L, 1L, 3L, 2L, 2L,
1L, 2L, 1L, 1L, 1L, 2L, 2L, 1L, 1L, 2L, 1L, 1L, 2L, 2L, 3L, 2L,
2L, 2L, 3L, 2L, 3L, 1L, 2L, 1L, 2L, 2L, 1L, 3L, 2L, 2L, 1L, 2L,
2L, 1L, 3L, 1L, 1L, 2L, 2L, 3L, 3L, 2L, 2L, 1L, 1L, 1L, 3L, 2L,
3L, 2L, 3L, 1L, 2L, 3L, 3L, 1L, 2L, 2L), .Label = c("1", "2",
"3"), class = "factor")), row.names = c(NA, -100L), class = "data.frame")
预计 effects 不适用于因子，否则输出将包含另一个维度，使结果有些复杂，这是非常合理的，就像我下面的解决方案一样，人们可能只希望对某个因素水平而不是所有水平产生影响。此外，正如我在下面解释的那样，分类变量的边际效应并不是唯一定义的，因此这将是 effects.
一个自然的解决方法是手动将因子变量转换为一系列虚拟变量，如
aux <- model.matrix(y ~ x, df3F)[, -1]
head(aux)
# x2 x3 x4 x5 x6 x7 x8 x9 x10 x11 x12 x13
# 1 0 0 0 0 0 0 0 0 0 1 0 0
# 2 0 0 0 0 0 0 0 0 0 1 0 0
# 3 0 0 0 0 0 1 0 0 0 0 0 0
# 4 0 0 0 0 0 0 0 0 1 0 0 0
# 5 0 0 0 0 0 0 0 1 0 0 0 0
# 6 0 0 0 0 0 0 1 0 0 0 0 0
这样数据就是
ml.datF <- mlogit.data(data.frame(y = df3F$y, aux), choice = "y", shape = "wide")
我们还需要手动构造公式
frml <- mFormula(formula(paste("y ~ 1 |", paste(colnames(aux), collapse = " + "))))
到目前为止一切顺利。现在如果我们 运行
fit.mnlF <- mlogit(frml, data = ml.datF)
head(effects(fit.mnlF, covariate = "x2", data = ml.datF))
# FALSE TRUE
# 1 -1.618544e-15 0.000000e+00
# 2 -1.618544e-15 0.000000e+00
# 3 -7.220891e-08 7.221446e-08
# 4 -1.618544e-15 0.000000e+00
# 5 -5.881129e-08 5.880851e-08
# 6 -8.293366e-08 8.293366e-08
那么结果不正确。 effects 在这里所做的是将 x2 视为一个 连续 变量，并计算这些案例的通常边际效应。即，如果对应于 x2 的系数是 b2 而我们的模型是 f(x,b2)，effects 计算 f 关于 b2 的导数并在每个观测向量 x我。这是错误的，因为 x2 只取值 0 和 1，而不是 0 或 1 左右的值，这是取导数的假设（极限的概念）！例如，考虑您的其他数据集 df1。在那种情况下，我们错误地得到
colMeans(effects(fit.mnlF, covariate = "x2", data = ml.datF))
# 1 2 3 4 5
# -0.25258378 0.07364406 0.05336283 0.07893391 0.04664298
这里有另一种方法（使用导数近似）得到这个不正确的结果：
temp <- ml.datF
temp$x2 <- temp$x2 + 0.0001
colMeans(predict(fit.mnlF, newdata = temp, type = "probabilities") -
predict(fit.mnlF, newdata = ml.datF, type = "probabilities")) / 0.0001
# 1 2 3 4 5
# -0.25257597 0.07364089 0.05336032 0.07893273 0.04664202
我没有使用 effects，而是通过两次使用 predict 手动计算了错误的边际效应：结果是平均值（{x2new = x2old + 0.0001 的拟合概率} - {拟合概率x2new = x2old}) / 0.0001。也就是说，我们通过将 x2 向上移动 0.0001（从 0 到 0.0001 或从 1 到 0.0001）来查看预测概率的变化。这两个都没有意义。当然，我们不应该期望 effects 有任何其他内容，因为数据中的 x2 是数字。
那么接下来的问题就是如何计算正确的（平均）边际效应。正如我所说，分类变量的边际效应并不是唯一定义的。假设 x_i 是个人 i 是否有工作，y_i 是他们是否有车。那么，那么至少有以下六点需要考虑。
- 从 x_i=0 到 x_i=1 时对 y_i = 1 的概率的影响。
- 从x_i=0 到x_i（观测值）时。
- 从 x_i 到 1.
现在，当我们对平均边际效应感兴趣时，我们可能只想对 1-3 的变化产生影响的那些个体进行平均。也就是说，
- 如果观测值不是 1，则从 x_i=0 到 x_i=1。
- 如果观测值不为 0，则从 x_i=0 到 x_i。
- 如果观测值不是 1，则从 x_i 到 1。
根据您的结果，Stata 使用选项 5，因此我将重现相同的结果，但实施任何其他选项都很简单，我建议考虑哪些选项对您的特定应用程序感兴趣。
AME.fun2 <- function(betas) {
aux <- model.matrix(y ~ x, df1)[, -1]
ml.datF <- mlogit.data(data.frame(y = df1$y, aux), choice="y", shape="wide")
frml <- mFormula(formula(paste("y ~ 1 |", paste(colnames(aux), collapse=" + "))))
fit.mnlF <- mlogit(frml, data = ml.datF)
fit.mnlF$coefficients <- betas
aux <- ml.datF # Auxiliary dataset
aux$x3 <- 0 # Going from 0 to the observed x_i
idx <- unique(aux[aux$x3 != ml.datF$x3, "chid"]) # Where does it make a change?
actual <- predict(fit.mnlF, newdata = ml.datF)
counterfactual <- predict(fit.mnlF, newdata = aux)
colMeans(actual[idx, ] - counterfactual[idx, ])
}
(AME.mnl <- AME.fun2(ml.fit$coefficients))
# 1 2 3 4 5
# -0.50000000 -0.17857142 0.03896104 0.16558441 0.47402597
require(numDeriv)
grad <- jacobian(AME.fun2, ml.fit$coef)
AME.mnl.se <- matrix(sqrt(diag(grad %*% vcov(ml.fit) %*% t(grad))), nrow = 1, byrow = TRUE)
AME.mnl / AME.mnl.se
# [,1] [,2] [,3] [,4] [,5]
# [1,] -5.291503 -2.467176 0.36922 1.485058 4.058994