如何用emmeans输出多项Logistic回归系数？OR结果为何与multinom有差异？

问题背景 #

我用以下代码拟合多项Logistic回归模型：

ml <- foreign::read.dta("https://stats.idre.ucla.edu/stat/data/hsbdemo.dta")
model.fit <- multinom(ses~female*prog+read, data=ml)

提取模型系数后计算比值比（OR），得到exp(1.271990)=3.567946和exp(2.092619)=8.106117。

但用emmeans包运行以下代码：

EMM <- emmeans(model.fit, revpairwise~ses|female*prog)
REMM <- regrid(EMM, transform="logit")
contrast(REMM, interaction="revpairwise", type="response", by="female")

得到的OR结果（如4.926、15.314）和multinom系数指数化后的结果不一致，请问原因是什么？怎么让emmeans输出和multinom系数指数化后一致的结果？

不一致的原因 #

• 模型系数的OR是单条件比值比：multinom输出的系数，是ses各水平相对于默认参考组（通常是ses的第一个水平），在固定其他协变量（read取均值、female和prog取特定组合）下的对数比值，指数化后得到的是单一协变量组合下ses水平间的条件OR。
• 当前emmeans代码计算的是交互对比OR：你写的代码先按female*prog分组计算ses的反向两两对比，再通过regrid和交互对比，本质是在计算不同female*prog组之间的ses对比差异，并非直接对应multinom系数所代表的单一组合下的OR，逻辑上和模型系数的定义完全不同。

解决方法 #

要让emmeans输出与multinom系数指数化一致的OR，需直接计算特定协变量组合下ses各水平相对于参考组的条件OR：

1. 先通过summary(model.fit)确认ses的参考水平（比如默认是low）。
2. 针对对应协变量组合计算OR：

# 指定female、prog的特定水平，read取均值（与模型系数计算的默认设定一致）
emm <- emmeans(model.fit, ~ ses, at = list(female = 0, prog = "general", read = mean(ml$read)))
# 计算相对于参考组的对比，指数化得到OR
contrast(emm, method = "pairwise", ref = 1, type = "response")

如果需要输出所有female*prog组合下的对应OR，可按分组计算：

# 按female*prog分组，read取均值
emm <- emmeans(model.fit, ~ ses | female*prog, at = list(read = mean(ml$read)))
# 每组内计算ses相对于参考组的OR
contrast(emm, method = "pairwise", ref = 1, type = "response")

这样得到的结果就会和multinom系数指数化后的OR完全匹配，因为两者都是基于相同的协变量固定值，计算ses水平相对于参考组的条件比值比。

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