拟合多分类逻辑回归后，用R的margins复现Stata操作遇报错求助

问题背景 #

使用nnet包拟合多分类逻辑回归模型后，尝试通过R的margins命令复现Stata中设定预测变量特定取值以获取因变量obesity_E预测值的操作，但ChatGPT生成的代码运行时出现报错：

Error in find_terms_in_model.default(model, variables = variables): Some values in 'variables' are not in the model terms.

相关代码 #

1. R中多分类逻辑回归拟合代码 #

library(nnet)
model <- multinom(obesity_E ~ age_100 + I(age_100^2) + obesity + age_100:obesity +
                  covid:race + age_100:race + education + education:race +
                  rabplace_5, data = female_98_2020, maxit=1000)

2. ChatGPT生成的报错代码 #

margins_model <- margins(model, variables = "obesity:race:rabplace_5",
                         at = list(age_100 = seq(0, 25, 1)), atMethod = "mean",
                         method = "probs", force = TRUE, noSe = TRUE,
                         save = "tran_point_F", replace = TRUE)

3. 目标复现的Stata代码 #

margins  , at (age_100=(0 (1) 25) obesity=(1 (1) 4) race=(0 (1) 3) rabplace_5=(1 (1) 5)) atmeans force nose saving(tran_point_F, replace ) 

报错原因 #

ChatGPT指定的variables = "obesity:race:rabplace_5"是一个三变量交互项，但该交互项并未出现在拟合的回归模型中。原模型仅包含二元交互项age_100:obesity、covid:race、age_100:race、education:race，因此margins包无法识别该变量，触发报错。

修正后的R代码 #

要完全复现Stata的margins操作，需对齐Stata的逻辑：指定目标变量的取值序列，将其余变量设为均值，同时匹配参数对应关系。修正后的代码如下：

# 加载margins包
library(margins)

# 生成符合要求的边际预测
margins_model <- margins(
  model,
  at = list(
    age_100 = seq(0, 25, 1),  # 对应Stata的age_100=(0 (1) 25)
    obesity = 1:4,            # 对应Stata的obesity=(1 (1) 4)
    race = 0:3,               # 对应Stata的race=(0 (1) 3)
    rabplace_5 = 1:5          # 对应Stata的rabplace_5=(1 (1) 5)
  ),
  atmethod = "mean",  # 对应Stata的atmeans，将未指定变量设为均值
  method = "probs",   # 获取因变量各分类的预测概率
  force = TRUE,       # 对应Stata的force，强制计算边际效应
  noSE = TRUE         # 对应Stata的nose，不计算标准误
)

# 将结果转为数据框并保存（替代Stata的save参数）
tran_point_F <- as.data.frame(margins_model)
write.csv(tran_point_F, "tran_point_F.csv", row.names = FALSE)

关键注意事项 #

1. 变量匹配：必须确保variables参数指定的变量/交互项存在于模型公式中，否则会触发识别错误；
2. 参数对应：R的margins包与Stata的margins命令参数并非完全同名，例如atmethod对应Stata的atmeans，noSE对应nose；
3. 结果保存：R中没有直接对应Statasaving的参数，需将margins对象转为数据框后导出。

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