使用R的nnet包计算多项Logistic回归AME遇阻，求替代方案

针对你使用nnet::multinom拟合多项Logistic回归后，无法用marginaleffects::avg_slopes计算AME的问题，以下是几个可行的替代方案：

方案1：使用margins包直接计算 #

margins包对nnet::multinom模型有较好的支持，只需在拟合模型时指定Hess = TRUE以确保能提取方差协方差矩阵：

# 安装并加载依赖包
install.packages("margins")
library(margins)
library(nnet)

# 拟合多项Logistic回归模型（注意修正原代码的变量名语法，如vote_choice而非vote choice）
model <- multinom(vote_choice ~ tchange + controls, data = your_data, Hess = TRUE)

# 计算tchange的平均边际效应
ame_results <- margins(model, variables = "tchange")

# 查看带标准误、p值的汇总结果
summary(ame_results)

方案2：手动计算AME（含bootstrap标准误） #

如果偏好手动实现逻辑，可以通过对比变量变化前后的预测概率均值来计算AME，还能通过bootstrap估计标准误：

基础计算 #

library(nnet)
library(dplyr)

# 拟合模型（必须加Hess=TRUE）
model <- multinom(vote_choice ~ tchange + controls, data = your_data, Hess = TRUE)

# 生成tchange增加1单位的数据集副本
data_shifted <- your_data %>% mutate(tchange = tchange + 1)

# 计算原数据集和移位后数据集的类别预测概率
pred_orig <- predict(model, newdata = your_data, type = "probs")
pred_shifted <- predict(model, newdata = data_shifted, type = "probs")

# 计算每个类别的平均边际效应
ame_manual <- colMeans(pred_shifted - pred_orig)
names(ame_manual) <- colnames(pred_orig)
ame_manual

Bootstrap估计标准误 #

library(boot)

# 定义bootstrap统计量函数
calc_ame <- function(data, idx) {
  sample_data <- data[idx, ]
  mod <- multinom(vote_choice ~ tchange + controls, data = sample_data, Hess = TRUE, trace = FALSE)
  pred_o <- predict(mod, newdata = sample_data, type = "probs")
  pred_s <- predict(mod, newdata = mutate(sample_data, tchange = tchange + 1), type = "probs")
  colMeans(pred_s - pred_o)
}

# 运行bootstrap（R为重复次数，可根据需求调整）
set.seed(123)
boot_output <- boot(data = your_data, statistic = calc_ame, R = 1000)

# 查看bootstrap结果（包含均值、标准误、置信区间）
print(boot_output)

方案3：改用支持marginaleffects的MNL模型包 #

如果可以更换模型拟合包，mlogit或brmultinom等包的模型能直接被marginaleffects::avg_slopes支持，以mlogit为例：

install.packages("mlogit")
library(mlogit)
library(marginaleffects)

# 将数据转换为mlogit要求的格式
mlogit_data <- mlogit.data(your_data, choice = "vote_choice", shape = "wide")

# 拟合多项Logistic回归
mlogit_model <- mlogit(vote_choice ~ 1 | tchange + controls, data = mlogit_data)

# 直接计算AME
avg_slopes(mlogit_model, variables = "tchange")

内容的提问来源于stack exchange，提问作者Gina