交互效应这东西做流行病学、临床研究和社科数据分析的人几乎天天碰到。但大多数教程讲的是“乘积交互项”也就是在模型里加 A×B然后看 P 值。真正到了病因学推断、风险评估这些场景审稿人和导师更关心的其实是“相加交互效应”而且经常要求你报 RERI、AP、S 这三个指标。R 语言里做这个分析并不复杂熟练之后一套流程下来确实能控制在五分钟左右关键是你得知道每一步在干什么以及不同模型逻辑回归、Cox、GLMM、GEE之间到底有什么区别。这篇就把我实际跑过的流程完整捋一遍从原理到代码再到我踩过的坑一次性说清楚。1. 先搞清楚相加交互效应到底是什么为什么要算1.1 相乘交互与相加交互的区别先说个最容易绕晕的点。传统回归模型里加入乘积项比如x1:x2检验的是“相乘尺度上的交互”意思是两个变量共同效应的乘积是否偏离了各自效应的乘积。逻辑回归、Cox 回归里最常见的交互项报的 OR、HR 都是乘法尺度所以那个交互项本质上是在回答“两个因素的联合效应是不是不等于各自效应相乘”。但公共卫生和临床研究里很多问题天然是“加法”逻辑。举个例子吸烟和石棉暴露都是肺癌的危险因素如果一个人既吸烟又接触石棉我们关心的不是“联合风险比单独吸烟的风险高多少倍再乘以单独接触石棉的风险高多少倍”而是“联合风险减去单独吸烟的风险再减去单独接触石棉的风险多出来的那一部分到底有多大”。这个“多出来的部分”就是相加交互它直接对应生物学上的协同或拮抗作用也直接关系到预防策略的制定——是优先消除哪个暴露还是两个都得管。所以你可以这么记乘积交互项回答“统计学上有没有交互”相加交互指标回答“公共卫生上有没有可归因的协同效应”。很多论文只报前者审稿人问一句“请补充相加交互指标”你就得会算。这也是我为什么强烈建议做流行病学分析的朋友把 RERI、AP、S 这三个指标当成标配。1.2 三个核心指标RERI、AP、S相加交互效应分析通常汇报三个指标名称和含义如下RERI相对超额危险度RERI OR11 - OR10 - OR01 1换成 HR、RR 同理。它表示两个因素同时存在时超过两个因素单独存在效应之和的那部分风险。RERI 0 提示协同正相加交互RERI 0 提示拮抗RERI 0 说明无相加交互。AP归因比AP RERI / OR11。它回答的问题是“在同时暴露的个体中有多大比例的疾病可以归因于两因素的协同效应”。AP 越接近 1说明协同作用贡献越大。S协同指数S (OR11 - 1) / [(OR10 - 1) (OR01 - 1)]。S 1 表示协同S 1 表示拮抗S 1 表示恰好无相加交互。需要特别注意这三个指标都不是直接从模型输出里读出来的而是用模型估计的 OR 等参数换算出来的。所以核心思路永远是先拟合带两个主效应和乘积项的模型提取系数和方差协方差矩阵再代入公式算指标和置信区间。1.3 什么时候必须报相加交互如果只是做个预测模型加不加相加交互指标无所谓。但下面这几类场景我强烈建议直接准备病因学研究和危险因素分析尤其是环境暴露、职业暴露、遗传风险这类主题审稿人有流行病学背景的临床研究特别是涉及多个危险因素联合作用时使用病例对照数据计算 OR、队列数据计算 RR/HR 的研究任何需要讨论“亚组差异”或“分层分析”的论文因为分层分析只能看到效应大小不同回答不了这种差异是否达到统计学显著。我自己最常被问到的就是“为什么乘积交互不显著但 RERI 显著”。这个现象在样本量不大时很常见后面我会专门讲。总之只要你做的是风险因素关联研究先把 RERI、AP、S 这三个指标的代码准备好有备无患。2. 方法论与工具选型为什么“5分钟”能做到2.1 核心原理从拟合模型到计算指标整个分析的原理可以压缩成四步拟合回归模型包含因素 A、因素 B 以及 A×B 乘积项从模型中提取两个主效应项的回归系数和乘积项的回归系数根据公式计算三个相加交互指标的估计值计算置信区间目前主流方法有 delta 法、Bootstrap 法和 Mover 法。逻辑回归、Cox 回归、GLMM、GEE 这四类模型前两步和第四步的公式完全一样区别只在于“回归系数对应的模型类型”不同。这也是为什么标题敢写“全支持”——只要你把模型拟合出来后面算交互指标的代码几乎是同一套。置信区间是个关键点。老的做法是直接用 delta 法近似标准误但 RERI、AP、S 都是非线性的比值组合正态近似在小样本下不太稳。现在比较稳妥的做法是小样本或率比较稀疏时优先用Bootstrap 置信区间大样本时可以用 delta 法和 Bootstrap 结果交叉验证epiR包默认提供的是基于正态近似的置信区间报告时建议在方法部分写明“confidence intervals were estimated using the delta method”。2.2 工具选型epiR 包为主手工计算兜底R 语言里做相加交互最省事的方案是epiR包里的epi.interaction()函数。它可以直接接受glm()或coxph()的结果自动计算 RERI、AP、S 以及置信区间代码量很小。但有一个坑很多人不知道epiR::epi.interaction()对lme4::glmer()和geepack::geeglm()的结果并不能直接使用至少在我的版本里是不兼容的。所以做 GLMM 和 GEE 的时候需要自己动手写一个计算函数提取固定效应系数和方差协方差矩阵再做估计。这部分代码也不难后面我会给出可直接复用的版本。所以我的工具选型策略是模型类型首选方案备选方案逻辑回归epi.interaction()手工代码Cox 回归epi.interaction()手工代码GLMM手工代码基于coef()vcov()Bootstrap 自写函数GEE手工代码基于geeglm系数 稳健方差Bootstrap 自写函数2.3 四类模型的适用场景差异为什么需要四种模型分别处理核心是数据结构和研究设计不同逻辑回归结局是二分类数据独立最经典的病例对照或横断面研究。Cox 回归结局是生存时间含删失用于队列研究和临床试验的生存分析。GLMM数据有层级结构如多中心、重复测量需要在模型里加入随机效应。这里的“相加交互”指标要用固定效应部分的系数计算不能把随机效应的方差算进去。GEE同样是处理相关数据但更强调“群体平均效应”常用于纵向数据或多层聚类数据。GEE 给的是稳健标准误计算置信区间时直接用稳健方差矩阵即可。在实际操作中如果你只是验证一下逻辑回归的交互结果epiR一把梭没问题但如果你的数据是纵向的忽略相关性直接跑普通逻辑回归很可能导致标准误偏小、假阳性上升换 GLMM 或 GEE 是更严谨的选择。3. 实战5分钟跑通完整的相加交互分析3.1 数据准备与变量编码先说编码问题这一步错了后面全白搭。计算 RERI 时两个暴露因素必须是二分变量而且编码方向要和你的研究假设一致。我习惯把“暴露组”编码为 1“非暴露组”编码为 0。以吸烟和高血压为例smoke1 吸烟0 不吸烟htn1 有高血压0 无高血压y1 发生心血管事件0 未发生time生存时间Cox 用id个体 IDGLMM/GEE 用center中心/分组变量GLMM 随机截距用。这里有个很容易犯的错把变量编码成-1/1或1/2。epiR和手动公式都默认二分类取0/1如果你用1/2编码OR 的参照组就变了RERI 和 S 全部失去意义。所以数据清洗时建议加一句data$smoke - ifelse(data$smoke yes, 1, 0) data$htn - ifelse(data$htn yes, 1, 0)接下来我们用一份模拟数据演示完整流程。数据生成代码不展开但你完全可以替换成自己的真实数据。3.2 核心函数 epi.interaction 的使用与参数解释先看epiR里epi.interaction()的基本用法install.packages(epiR) library(epiR) epi.interaction(model, coef c(2, 3), em TRUE, ci TRUE, conf.level 0.95)参数含义model拟合好的glm()或coxph()对象coef一个长度为 2 的向量指定两个主效应项在回归系数向量中的位置。注意不是变量名而是位置索引em是否输出暴露-混杂四格表形式的估计结果建议设为TRUEci是否计算置信区间建议TRUEconf.level置信水平默认 0.95。coef这个参数最容易出错。如果你模型里有多个协变量两个主效应项不一定在位置 2 和 3。一个通用做法是coef c(which(names(coef(model)) smoke), which(names(coef(model)) htn))这样即使协变量多也不会选错位置。3.3 逻辑回归实操先跑一个最经典的案例——二分类结局、独立数据# 拟合带乘积项的 logistic 回归 m_logit - glm(y ~ smoke * htn age sex, data dat, family binomial()) # 查看系数位置 names(coef(m_logit)) # 计算相加交互指标 epi.interaction(m_logit, coef c(which(names(coef(m_logit)) smoke), which(names(coef(m_logit)) htn)), em TRUE, ci TRUE, conf.level 0.95)输出结果会包含一张 2×2 的暴露效应表以及一行核心结果RERI、AP、S的估计值和置信区间。假设输出 RERI 0.8295% CI: 0.15, 1.49AP 0.3195%CI: 0.10, 0.52S 1.7595%CI: 1.08, 2.42说明吸烟和高血压对心血管事件存在正的相加交互即两者同时存在时超额风险大于各自风险之和。表格里还会给出RERI的置信区间是基于正态近似得到的。如果置信区间下界包括 0说明相加交互在统计学上不显著如果 AP 的置信区间包含 0同理。S 的置信区间包含 1说明不支持显著协同。3.4 Cox 回归实操生存数据场景下把结局换成Surv(time, y)即可library(survival) m_cox - coxph(Surv(time, y) ~ smoke * htn age sex, data dat) epi.interaction(m_cox, coef c(which(names(coef(m_cox)) smoke), which(names(coef(m_cox)) htn)), em TRUE, ci TRUE)这里有个细节epiR对coxph对象的支持是基于系数和方差协方差矩阵实现的所以计算逻辑和逻辑回归一致。但 Cox 模型本身有比例风险假定如果你的暴露因素和时间的交互明显比如 Schoenfeld 残差检验 P 0.05那 HR 本身就是随时间变化的用单一 HR 算出来的 RERI 也会失真。这种情况建议先处理时变效应或者用分段模型。另一个需要注意的点是Cox 模型里epi.interaction给出的 OR/HR 都是从模型系数转换来的所以模型的拟合质量直接影响结果。如果某一层的人数很少HR 的置信区间会非常宽RERI 和 S 的置信区间也会跟着膨胀。3.5 GLMM 实操手工计算GLMM 场景常见于多中心临床试验或重复测量队列。lme4::glmer()拟合后epiR不认所以我们自己写函数。核心思路是提取固定效应系数、提取固定效应的方差协方差矩阵然后用 delta 法公式算三个指标的标准误。library(lme4) m_glmm - glmer(y ~ smoke * htn age sex (1 | center), data dat, family binomial()) # 提取系数与方差协方差 b - fixef(m_glmm) V - vcov(m_glmm) # 三个指标的计算函数 calc_additive - function(b, V) { # 找到两个主效应和交互项系数 b_smoke - b[smoke] b_htn - b[htn] b_int - b[smoke:htn] OR11 - exp(b_smoke b_htn b_int) OR10 - exp(b_smoke) OR01 - exp(b_htn) RERI - OR11 - OR10 - OR01 1 AP - RERI / OR11 S - (OR11 - 1) / ((OR10 - 1) (OR01 - 1)) # delta 法近似标准误此处需要梯度 # 推荐直接用数值微分或 bootstrap c(RERI RERI, AP AP, S S) } calc_additive(b, V)直接输出点估计还不行我们还需要置信区间。简单可靠的做法是bootMer()做 Bootstrap。下面这段代码我实测过小数据集几百人也能跑library(boot) # 写一个从模型对象提取指标的函数 boot_additive - function(m) { b - fixef(m) b_smoke - b[smoke]; b_htn - b[htn]; b_int - b[smoke:htn] OR11 - exp(b_smoke b_htn b_int) OR10 - exp(b_smoke) OR01 - exp(b_htn) RERI - OR11 - OR10 - OR01 1 AP - RERI / OR11 S - (OR11 - 1) / ((OR10 - 1) (OR01 - 1)) c(RERI RERI, AP AP, S S) } # semi-parametric bootstrap set.seed(123) boot_out - bootMer(m_glmm, boot_additive, nsim 1000) boot.ci(boot_out, index 1, type perc) # RERI 的置信区间bootMer的nsim建议至少 500稳妥一点用 1000。如果数据量特别大用nsim 200先看趋势可以但正式结果不建议少于此数。3.6 GEE 实操手工计算GEE 的逻辑和 GLMM 类似区别是使用geepack::geeglm()协方差矩阵要用稳健方差。GEE 拟合的是人群平均效应所以系数解释与 GLMM 不同但计算 RERI 的公式完全一致。library(geepack) m_gee - geeglm(y ~ smoke * htn age sex, id id, data dat, family binomial(), corstr exchangeable) b - coef(m_gee) V - vcov(m_gee) # 默认是稳健方差 calc_additive(b, V)这里有个关键点geeglm的系数名可能和glm不完全一样比如交互项可能叫smoke:htn。先用names(coef(m_gee))确认一下再传入函数。否则b[smoke:htn]会返回NA后面全是NA。GEE 做 Bootstrap 有两种思路对个体聚类单元做 Bootstrap每次重抽样后重新拟合 GEE用基于渐近正态的 delta 法。个体 Bootstrap 更稳健但计算量大。如果数据有几千个个体建议用geeglm delta 法再和epiR在普通 logistic 上的结果对比验证。3.7 手工计算 RERI 的 delta 法加深理解如果你不想依赖epiR也可以自己写一个完整的计算函数。这里给出一个带 delta 法标准误的版本适用于逻辑回归、Cox 回归等基于glm/coxph的对象additive_interaction - function(model, var1, var2) { b - coef(model) V - vcov(model) i1 - var1 i2 - var2 i3 - paste0(var1, :, var2) if (!i3 %in% names(b)) i3 - paste0(var2, :, var1) if (!i3 %in% names(b)) stop(交互项没找到请检查变量名) b1 - b[i1]; b2 - b[i2]; b3 - b[i3] OR11 - exp(b1 b2 b3) OR10 - exp(b1) OR01 - exp(b2) RERI - OR11 - OR10 - OR01 1 AP - RERI / OR11 S - (OR11 - 1) / ((OR10 - 1) (OR01 - 1)) # 数值梯度 f - function(bvec) { o11 - exp(bvec[1] bvec[2] bvec[3]) o10 - exp(bvec[1]) o01 - exp(bvec[2]) c(RERI o11 - o10 - o01 1, AP (o11 - o10 - o01 1) / o11, S (o11 - 1) / ((o10 - 1) (o01 - 1))) } require(numDeriv) g - jacobian(f, c(b1, b2, b3)) se - sqrt(diag(g %*% V[c(i1, i2, i3), c(i1, i2, i3)] %*% t(g))) est - f(c(b1, b2, b3)) data.frame( Estimate est, SE se, Lower est - 1.96 * se, Upper est 1.96 * se ) } # 使用示例 additive_interaction(m_logit, smoke, htn)这段代码最关键的地方是jacobian()计算梯度然后套用 delta 法的公式g %*% V %*% t(g)。对小样本结果和epiR有细微差异是正常的因为epiR可能用了不同的近似方式。你要是怕与审稿人挑刺就统一用epiR并在方法部分注明软件版本和分析包。4. 实战中常见的坑与排查记录4.1 变量编码方向搞反了这是我自己早期犯过的错也是很多学员最容易踩的坑。RERI 的计算完全依赖“暴露1非暴露0”的编码。如果你把“不吸烟”编码成 1那么 OR10、OR01、OR11 的参照就全反了算出来的 RERI 可能从“正协同”变成“负协同”结论完全颠倒。排查方法很简单跑完epiR之后看输出的 2×2 表格里暴露组的 OR 是否和预期一致。如果发现 OR 1 而你预期暴露是危险因素先回去查编码。4.2 置信区间出现负值或超界RERI 的置信区间理论上可以包含负值这没问题。但 S 的置信区间如果出现负值或者 AP 的置信区间出现负下界虽然数学上可能出现表示拮抗但你要先怀疑是不是 delta 法小样本近似失效。我在小样本案例里遇到过 S 下界为负的情况用 Bootstrap 之后区间就正常了。所以遇到反常区间优先用 Bootstrap 复核。还有一个特殊情况当 OR11 接近 1 时AP 和 S 都接近无定义报告时要小心不要在解释里过度发挥。4.3 交互项不显著但 RERI 显著或反过来这是期刊审稿里最常见的疑问。相乘交互的零假设是“乘积项系数为 0”相加交互的零假设是“RERI 0”两者本来就不是一个概念出现不一致非常正常。尤其当主效应都比较大时即使乘积项不显著RERI 也可能显著。面对审稿人我的标准答复是“相乘交互检验的是效应是否偏离乘法模型相加交互检验的是偏离加法模型。两者回答的科学问题不同本研究以相加交互为主原因是……” 这样解释既专业又清晰。4.4 模型收敛警告与样本量GLMM 和 GEE 对样本量要求更高。如果某个暴露组的人数特别少glmer()可能出现“Model failed to converge”的警告。这时候不要慌先检查每组例数考虑去掉协变量、简化随机效应结构或者用 Firth 惩罚回归代替普通逻辑回归。另外epi.interaction()在coxph模型上如果遇到极端分层如某一层事件数为 0也会报错。建议在跑之前先用table(dat$smoke, dat$htn, dat$y)看一眼分布。4.5 快速自查清单我每次跑完相加交互都会过一遍这个清单两个暴露因素是否都编码为 0/1方向是否符合假设模型里是否确实包含主效应和乘积项coef参数的位置索引是否正确直接用名字匹配更安全输出的 2×2 表格里的 OR 是否符合常理置信区间是否出现异常异常则换 Bootstrap 复核报告时注明使用的是epiR版本、Bootstrap 次数或 delta 法说明。5. 关于“5分钟”的真相与后续扩展5.1 为什么很多教程做不到5分钟标题说“5分钟搞定”听起来像噱头其实背后是流程标准化。为什么很多人做这个分析要折腾一整天因为他们在三个环节卡壳第一不知道epiR包的存在手动算完点估计后不会算置信区间第二模型对象换成 GLMM、GEE 后不知道支持有限一直在epi.interaction()上试错报错了才想到手工实现第三卡在“交互项名字”上smoke:htn和htn:smoke顺序不同代码里用了paste0拼名字一旦顺序不对就取不到系数。把这三关都打通剩下就是机械操作。熟悉之后逻辑回归真的五分钟内能出完整结果GLMM 和 GEE 因为要跑 Bootstrap时间会稍微长一点但思路已经完全固化不存在“不知道下一步干什么”的情况。5.2 可复用的完整代码模板我把最常用的逻辑回归模板贴在下面你只需要替换数据和变量名library(epiR) # 1. 数据清洗确保 0/1 编码 dat$A - ifelse(dat$A 1, 1, 0) dat$B - ifelse(dat$B 1, 1, 0) # 2. 拟合模型 fit - glm(Y ~ A * B age sex, data dat, family binomial()) # 3. 计算相加交互 epi.interaction(fit, coef c(which(names(coef(fit)) A), which(names(coef(fit)) B)), em TRUE, ci TRUE)Cox 版本只需把第二步换成coxph(Surv(time, Y) ~ A * B age sex, data dat)。GLMM 和 GEE 用上面第 3.5、3.6 节的手工函数即可。5.3 后续扩展方向如果你经常做这类分析可以做三件提升效率的事一是封装自己的函数把epi.interaction()的调用、bootMer的 Bootstrap、结果整理成统一格式输出tibble或直接导出 CSV省得每次复制粘贴。二是可视化。现在没有特别好用的现成包直接画 RERI 森林图但可以用ggplot2自己画——把四个 OROR11、OR10、OR01、参照组 1和 RERI、AP、S 的点估计加置信区间放到一张图里审稿人看了会非常直观。三是和亚组分析结合。比如按性别分层分别算男性和女性的 RERI再比较两个 RERI 的差异。这个玩法在临床研究里很受欢迎因为它能回答“协同效应是否在某个亚组里更强”。最后再分享一个小技巧不管你用哪种模型跑完后顺手把sessionInfo()里epiR和lme4的版本记录下来写论文方法部分时直接引用。审稿人问到具体算法时你可以明确回答用的是哪个包的哪个版本这种细节在返修时非常加分。我自己的习惯是所有交互分析代码都放在同一个 R 脚本里注释写好数据和模型类型三个月后回来看还能直接复现——这才是“几分钟搞定”的底气所在。