在R中复现Stata逻辑回归结果时的meno变量OR值差异问题

在R中复现Stata的基础逻辑回归模型时，除meno变量外，其余变量的优势比（OR）及系数完全一致，但meno的OR值在R中显著偏高：Stata中1.meno的OR为1.680891，R中factor(meno)1的OR为2.83。

所用代码 #

Stata代码 #

logistic sfcancer2yr c.age i.bminewcat i.n_rels i.race_eth i.meno i.bdnewcat i.biopsycoll i.birads i.bdnewcat#i.meno

R代码 #

trainmodel <- glm(sfcancer2yr ~ age + factor(bminewcat) + n_rels + factor(race_eth) + factor(meno) + bdnewcat + biopsycoll + birads + bdnewcat:meno, family = "binomial"(link = "logit"), data = traindata)

coef.orig <- coef(trainmodel)

exp(trainmodel$coefficients[-1])

已排查操作 #

• 将meno在R中设为因子或非因子，问题依旧
• 确认两个软件导入的数据集完全一致
• 移除交互项后，两个软件结果完全一致，推测差异来自交互项的处理方式

差异原因与解决方法 #

核心原因 #

差异源于交互项的变量类型与参考组定义不一致：

1. Stata中i.bdnewcat会自动将变量转为分类因子，并以第一个水平为参考组；而你的R代码中bdnewcat未用factor()包裹，默认被当作连续变量处理，导致交互项bdnewcat:meno的计算逻辑和Stata完全不同。
2. 即使bdnewcat是分类变量，若R与Stata的参考组设置不一致（比如Stata用0作参考，R默认用1作参考），会导致meno的主效应解释变化——有交互项时，主效应代表的是参考组下的效应，而非全局平均效应，这直接影响系数及OR值。

修正步骤 #

1. 将bdnewcat转为因子，确保和Stata的分类逻辑一致，修改R代码：

trainmodel <- glm(sfcancer2yr ~ age + factor(bminewcat) + n_rels + factor(race_eth) + factor(meno) + factor(bdnewcat) + biopsycoll + birads + factor(bdnewcat):factor(meno), family = binomial(link = "logit"), data = traindata)

2. 对齐参考组：将R中meno和bdnewcat的参考组调整为与Stata一致（以下示例以参考组为0为例，需根据你的实际数据调整）：

traindata$meno <- relevel(factor(traindata$meno), ref = "0")
traindata$bdnewcat <- relevel(factor(traindata$bdnewcat), ref = "0")

3. 重新计算OR值，此时factor(meno)1的OR应与Stata结果一致。

示例数据提供建议 #

无需提供全量7万条数据，只需生成最小复现数据集：

• 提取包含模型所有变量的100-200条随机样本
• 确保样本覆盖meno和bdnewcat的所有水平组合
• 导出为CSV格式即可，便于他人快速复现问题

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