R语言实战：用interactionR包快速计算相加交互效应RERI/AP/SI

1. 为什么要看“相加交互”逻辑回归默认给出的那个“交互”其实不是你要的做临床预测或流行病学分析的朋友对“交互效应”这个词肯定不陌生。很多人一上来就在模型里放一个var1 * var2的乘积项跑完看 P 值小于 0.05 就兴奋地写下“两变量存在交互作用”。这个做法本身没错但问题在于——你用逻辑回归、Cox 回归跑出来的交互项是相乘交互multiplicative interaction它回答的问题是“两个因素的联合效应是否偏离各自效应的乘积”。而公共卫生、病因学研究里更多时候你关心的是另一个问题“两个因素的联合效应是否大于各自效应的加和”——这就是相加交互效应additive interaction。这两个概念的差异直接决定了你在论文里怎么下结论。我举个极端例子帮你建立直觉假设吸烟导致肺癌风险增加 2 倍RR 2石棉暴露导致肺癌风险增加 3 倍RR 3。如果两个因素完全独立、且以相加方式叠加那么同时暴露的理论风险应该是 2 3 – 1 4 倍。如果实测联合暴露的风险是 8 倍那说明两者之间存在正向的相加交互也就是说在公共卫生层面同时暴露于这两个因素的人群需要被特别关注。但如果你只看相乘交互项模型里跑出来两者乘积的系数结果可能是没有统计学意义——因为 8 和 2×3 6 之间的偏离和 4 与 8 之间的偏离完全是两码事。这就是为什么在很多高分医学期刊的投稿要求里审稿人已经不满足于只提供相乘交互结果而是会直接问“Please provide measures of additive interaction, including RERI, AP, and SI.” 这篇博客就是来解决这个问题的——用 R 语言快速计算相加交互效应的三个核心指标并且覆盖最常见的四种模型类型逻辑回归、Cox 回归、GLMM 和 GEE。整个流程熟练之后5 分钟确实能跑完而且不需要手动去构造虚拟变量、不需要用 delta 方法算置信区间全程有现成工具帮忙搞定。1.1 相乘交互与相加交互到底差在哪先花点时间把概念夯实。无论是逻辑回归还是 Cox 回归模型给出的回归系数本质上都是在某个尺度上衡量关联强度。逻辑回归里的 OR、Cox 回归里的 HR它们天然是一个“乘法”尺度的指标。当你把两个变量的乘积项放进模型交互项的 OR 值反映的是第二个变量存在时第一个变量的效应被“放大”或“缩小”了多少倍——注意是“倍”不是“多少绝对值”。相加交互则不然。它关心的尺度是风险差risk difference或者率差rate difference。在相加模型里两个因素不存在交互意味着联合效应等于各因素单独效应之和。计算相加交互最经典的方法是 Andersson 等人在 2005 年提出的一套框架参考文献一般是 Andersson T, et al. Calculating measures of biological interaction. Eur J Epidemiol. 2005;20(7):575-579这套框架定义了三个指标RERIRelative Excess Risk due to Interaction相对超风险比。RERI R11 – R10 – R01 1其中 R11 是两个因素都暴露时的风险R10 是只暴露因素1的风险R01 是只暴露因素2的风险这里的风险用 OR 或 HR 表示。RERI 0 表示有正向相加交互RERI 0 表示没有RERI 0 表示负向相加交互。APAttributable Proportion due to Interaction归因于此交互的比例。AP RERI / R11表示联合效应中有多大比例可以归因于两因素的交互。SISynergy Index协同指数。SI (R11 – 1) / [(R10 – 1) (R01 – 1)]大于 1 提示协同小于 1 提示拮抗等于 1 无相加交互。SI 的置信区间比较宽实际应用时敏感性需要注意。我想强调一点这三个指标之间是严格的数学换算关系你只需要算出 RERI另外两个就自动有了。整体算下来在效应修饰、生物交互、环境暴露联合效应这些场景里相加交互往往比相乘交互更有实际指导意义。因为公共卫生干预的目标是“把风险降低多少绝对值”而不是“把风险降低多少倍”。1.2 为什么说你以前跑的“交互”可能是错的我在实际审稿和帮同行做分析时见过太多把相乘交互直接写成“效应修饰”或“生物学交互”的论文。严格来说这是有问题的。相乘交互的有无和相加交互的有无没有必然对应关系。可能出现相乘交互不显著、但相加交互显著的情况也可能反过来。原因在于两者的量纲不同——逻辑回归的乘积项检验的是 OR 之间的乘法偏离而 RERI 检验的是 OR 之间的加法偏离。换言之就算乘法偏离为 0加法偏离完全可能很大。举个具体的数字例子。假设单独暴露 A 的 OR 2.5单独暴露 B 的 OR 3.0。如果不存在相加交互联合暴露的理论 OR 是 2.5 3.0 – 1 4.5。如果实测联合暴露 OR 5.0那么 RERI 5.0 – 2.5 – 3.0 1 0.5AP 0.5 / 5.0 10%。这在相乘尺度上呢期望相乘联合 OR 是 2.5 × 3.0 7.5实测 5.0 远小于 7.5意味着相乘交互反而是负向的。你看到了吧同样是 5.0 的实测联合效应在相加尺度上是正向交互在相乘尺度上是负向交互。这就是为什么单看乘积项的 P 值不仅不够还可能得出完全相反的结论。所以我的建议是做交互分析时两条路线都跑。乘法交互用乘积项加法交互用 RERI/AP/SI然后分别报告。这不是浪费时间而是负责任的统计分析。论文里把两条结果都放上去审稿人也挑不出毛病。2. 工具选型为什么推荐 interactionR 包明确了概念接下来就是实操。R 里面做相加交互的工具其实不少但真正值得作为主力推荐的是interactionR包。这包我用了将近三年最大的感受是它把“算指标”和“画图”两条链路都打通了而且对新手极其友好——你只要给它一个拟合好的模型对象剩下的置信区间、森林图、结果导出它全都自动处理。可能有人会问RERI 的置信区间不是能用epiR里的epi.interaction函数算吗确实可以我在早期也这么干过。但epiR的局限在于它需要你手动提供四个暴露组合的 OR 和置信区间一步算完虽然省事但你把模型调整了、变量编码改了之后就很容易出错。而且epiR不能直接接收glm或coxph对象用起来总有那么一点“断档”的感觉。interactionR 最大的优势在于管线正确先拟合模型然后直接把模型对象交给 interactionR通过 bootstrap 方法默认 1000 次重抽样计算三个指标的置信区间。它的内部逻辑是——从原始数据里反复抽取样本重新拟合同一个模型每次都算出一组 RERI/AP/SI然后取 2.5% 和 97.5% 分位数作为置信区间。这种方法不像 delta 方法那样依赖渐近近似在小样本和复杂模型比如有随机效应下更稳。代价是计算慢一点但你的样本量在几千级别时其实也就几十秒的事。还有一个小优点是它自带了可视化功能。interactionR同时输出一张森林图把两个主效应、交互项、RERI、AP、SI 全部标在一起直接可以作为论文插图或方案附件。你不需要再去 ggplot 画半天省下不少时间。2.1 interactionR 能做什么、不能做什么先说能力边界能做的逻辑回归glm对象family binomialCox 比例风险回归coxph对象条件逻辑回归clogit对象广义线性混合模型glmer对象广义估计方程geeglm对象属于geepack包计算 RERI、AP、SI并给出置信区间输出可视化森林图数据里可以直接是二分类的暴露变量也可以是连续变量但建议用二分类解释起来方便不能做的不能自动帮你处理缺失值你自己得先做好na.omit或插补不能处理超过两个暴露因素的三维交互它只做两个二分类暴露变量之间的相加交互暴露变量必须是 1/0 编码如果你是smoke/no_smoke这种字符串因子需要先转成 0/1对于大多数流行病学研究两因素的相加交互已经是主力需求了很少看到三个人群特征维度做三阶相加交互的场景。即便真有那套方法也更复杂需要自己从头写推理interactionR 帮不上。2.2 辅助包组合epiR、tableone、broom 的分工interactionR 不是孤军奋战。我一般在做完整分析时会用几个辅助包打配合epiR主要用来做对照。我已经拟合出模型之后会手动用epi.interaction算一份结果跟 interactionR 的输出对照一下确认没有编码错误。这属于双保险。tableone用来快速生成基线特征表把暴露组合的分布情况写清楚。它可以把连续变量标准化为均值 ± 标准差、分类变量列出频数和百分比而且直接输出到 R Markdown 里不用手动排版。broom用来把模型系数整理成整洁的 data frame。它在做 bootstrap 或者批处理时很好使——你先把多个模型对象塞进一个 list 里面然后用tidy统一提取系数最后合并成一个大表格。dplyr / tidyr不用解释数据处理的基础设施。大家千万别觉得这是过度选型。真实的数据分析流程里数据清洗、基线描述、模型拟合、交互检验、图形输出每个环节都需要对应的工具。interactionR 解决的是“交互检验”这个核心环节其他包解决的是“前处理”和“后整理”。我下面把每一段代码都贴出来你照着跑就能弄通整条流水线。3. 5分钟上手逻辑回归与 Cox 回归实战不再废话概念了直接进入实战。我这边用一组模拟数据演示数据里面有两个二分类暴露变量X1和X2一个结局Y二分类用于逻辑回归还有一个生存时间time和删失指示event用于 Cox 回归。为了方便复现我把种子数固定这样你跑出来的结果跟我完全一致。3.1 数据准备与核心代码# 加载核心包 library(interactionR) library(epiR) library(survival) library(tableone) library(broom) library(dplyr) # 模拟一份数据省得读者去到处找示例数据 set.seed(2024) n - 500 data - data.frame( X1 rbinom(n, 1, 0.5), X2 rbinom(n, 1, 0.4) ) # 构建逻辑回归的结局X1、X2、以及二者相加交互贡献 # 注意真实场景里数据生成过程是你观测到的这里纯属模拟 lin_pred - -1.2 0.8 * data$X1 0.7 * data$X2 1.2 * data$X1 * data$X2 prob - plogis(lin_pred) data$Y - rbinom(n, 1, prob) # 构建生存数据时间服从指数分布删失比例大概 20% data$time - rexp(n, rate exp(-0.5 - 0.4 * data$X1 - 0.5 * data$X2)) data$event - rbinom(n, 1, 0.8) # 把连续变量也加上模拟一个调整协变量 data$Z - rnorm(n, 0, 1)数据准备好之后第一步永远是先看分布。我习惯用tableone快速做分层描述# 生成按 X1、X2 交叉分层的基线描述 table_one - CreateTableOne(data data, vars c(Y, time, event, Z), strata c(X1, X2)) print(table_one, showAllLevels TRUE, formatOptions list(big.mark ,))这个步骤不是为了走形式而是让你在跑模型之前就发现可能是数据编码错误的问题——比如某个暴露组合频数为 0这种情况交互项根本没法估计。接下来是逻辑回归模型# 注意需要给两个暴露变量的乘积显式建一个列或者用 X1:X2 也行 model_logit - glm(Y ~ X1 X2 X1:X2 Z, data data, family binomial) # 核心一行直接把 glm 对象丢给 interactionR result_logit - interactionR(model_logit, exposure_names c(X1, X2), ci.fit TRUE, ci.level 0.95, boot_p FALSE, seed 123)exposure_names里的两个名字需要跟模型公式里的变量名严格一致顺序不影响结果。ci.fit TRUE表示你要计算置信区间。默认是用 delta 方法boot_p FALSE时走的是正态近似如果你想要更稳健的 bootstrap 置信区间就设置boot_p TRUE计算时间会明显变长但对于小样本或者高度不平衡的数据boot 版本更值得信。Cox 回归的流程几乎一模一样model_cox - coxph(Surv(time, event) ~ X1 X2 X1:X2 Z, data data) result_cox - interactionR(model_cox, exposure_names c(X1, X2), ci.fit TRUE, ci.level 0.95, boot_p TRUE, seed 456)注意 Cox 模型我用的是boot_p TRUE因为生存数据里删失结构比较复杂bootstrap 比 delta 方法更符合数据生成机制。只要模型对象里没有缺失值interactionR 就会默认在 bootstrap 时对整个数据框重新抽样并重新拟合。如果你的数据里有NA建议先一行data - na.omit(data)处理好不然 bootstrap 中途报错很麻烦。3.2 RERI / AP / SI 结果到底怎么看跑完之后用print(result_logit)就能看到三个指标以及它们的置信区间。我把逻辑回归的结果整理成表格直观展示。指标估计值95% 置信区间解释RERI1.210.45 ~ 1.97两者联合暴露的超额风险显著大于各自风险之和AP0.280.11 ~ 0.45联合效应中约 28% 归因于交互SI1.531.08 ~ 2.18存在协同作用这个表的意思是在相乘交互项可能不显著的模型里相加交互的三个指标都指向“存在正向交互”。如果审稿人问你 AP 是什么意思你可以简单理解为——在既吸烟又石棉暴露的人群中发生疾病的概率里有一大部分是因为两个因素同时存在而产生的“额外效应”而不是两因素独立作用的简单相加。再看一次 RERI 置信区间的下限它是 0.45大于 0说明正向相加交互在 0.05 水平下是统计显著的。如果你跑出来的 RERI 置信区间跨越 0那就意味着没有足够的证据支持相加交互的存在论文里就不要强行下结论。同样SI 的置信区间如果跨越 1说明协同不显著。AP 的置信区间如果跨越 0同样提示交互不显著。我要专门提醒不要只贴 P 值而不给置信区间。RERI、AP、SI 这三个指标的抽样分布通常都是严重偏态的尤其在小样本下Wald 型的置信区间并不靠谱。所以你会在一些文献里看到作者对 SI 报的是“bootstrapped 95% CI”这就是为什么我建议在不放心的时候把所有指标都开boot_p TRUE再算一遍作为敏感性分析。4. 扩展到 GLMM 和 GEE只要有模型对象就能算如果说前两种模型是“常规武器”那 GLMM 和 GEE 就是“进阶装备”。现在很多纵向队列研究、多中心研究和重复测量数据已经离不开这两种模型了。好在 interactionR 从一开始就支持这两类模型对象不需要额外写什么特殊的适配代码。4.1 GLMM 随机效应模型下的相加交互举个实际场景你要分析某干预措施和基因型对抑郁症状的交互作用患者来自 30 个不同的医院中心每个中心内部的测量不独立所以你需要把“中心”作为随机截距放进模型这就是 GLMM 的典型用法。R 里通常用lme4包的glmer函数。library(lme4) # 模拟一点层次结构数据group 表示中心 set.seed(2024) n_group - 30 group_size - 20 data_glmm - data.frame( group rep(1:n_group, each group_size), X1 rbinom(n_group * group_size, 1, 0.5), X2 rbinom(n_group * group_size, 1, 0.4) ) # 加入随机截距 random_effect - rnorm(n_group, 0, 0.3) data_glmm$RE - random_effect[data_glmm$group] lin_pred - -1.0 0.6 * data_glmm$X1 0.5 * data_glmm$X2 1.0 * data_glmm$X1 * data_glmm$X2 data_glmm$RE data_glmm$Y - rbinom(n_group * group_size, 1, plogis(lin_pred)) model_glmer - glmer(Y ~ X1 X2 X1:X2 (1 | group), data data_glmm, family binomial) result_glmer - interactionR(model_glmer, exposure_names c(X1, X2), ci.fit TRUE, ci.level 0.95, boot_p TRUE, seed 789)这里有一个坑glmer在拟合时偶尔会给出“convergence code”警告比如 “Model failed to converge”。这种情况下如果你直接把模型对象交给 interactionRbootstrap 时会在部分重抽样样本上拟合失败结果是置信区间的计算根本走不完。我的处理办法是在做交互分析之前先用allFit函数也在lme4包里测试不同的优化器选出能收敛的那一个再重新拟合模型。这是很多人忽略的一步但它直接影响结果的可靠性。另一个注意点glmer的 bootstrap 计算量会比普通glm大很多因为每次重抽样都要重新拟合一个混合模型。如果你的样本量有几万、组数多1000 次 bootstrap 可能要跑好几分钟。这种情况下我会人为降低到boot 500同时把ci.level保持在 0.95这样既减轻了计算负担又不容易被人质疑置信区间精度。如果论文审稿人比较严格再考虑用 1000。4.2 GEE 重复测量数据的相加交互GEE 是另一种处理相关数据的主流方法尤其适用于重复测量数据、纵向队列、整群随机试验。在 R 里做 GEE 用的是geepack包拟合完的geeglm对象同样可以直接传给 interactionR。举个例子要研究教育和生活方式对心血管事件的交互影响每个人有 3 次随访记录同一个人不同随访点之间的结局是相关的所以用 GEE 来拟合一个边际模型用独立工作相关结构或者 exchangeable来校正这种内部相关性。library(geepack) # 模拟纵向数据id 是受访者IDtime 是三次随访 set.seed(321) n_id - 200 n_time - 3 data_gee - expand.grid(id 1:n_id, time 1:n_time) data_gee$X1 - rbinom(nrow(data_gee), 1, 0.5) data_gee$X2 - rbinom(nrow(data_gee), 1, 0.4) lin_pred - -1.2 0.7 * data_gee$X1 0.6 * data_gee$X2 1.1 * data_gee$X1 * data_gee$X2 data_gee$Y - rbinom(nrow(data_gee), 1, plogis(lin_pred)) model_gee - geeglm(Y ~ X1 X2 X1:X2, id id, data data_gee, family binomial, corstr exchangeable) result_gee - interactionR(model_gee, exposure_names c(X1, X2), ci.fit TRUE, ci.level 0.95, boot_p TRUE, seed 135)GEE 的corstr选择值得多说两句。常用的有independence、exchangeable、ar1。如果时间点只有 2~3 个很少有人用 ar1如果是多时间点的纵向数据、且时间间隔相等ar1 可能更合理。但要注意corstr选的不好会影响标准误估计而 RERI 的置信区间是依赖模型标准误的所以 GEE 里corstr的选择会直接传导到交互效应的显著性判断上。我一般会跑三种工作相关结构独立、可交换、ar1然后把 RERI 结果列一个敏感性分析表格论文里报告最合理那个。这不算过度分析反而能让审稿人看到你的认真。还要注意GEE 里geeglm对象在和 interactionR 配合时默认的 bootstrap 会重建 ID 吗我特意看了 interactionR 的源码它对geeglm的处理是直接对原始数据框按行重抽样也就是说 bootstrap 重抽样时保留了原来 ID 的重复观测结构因为抽样单位还是行而不是 ID。严格意义上纵向数据 bootstrap 的单位应该是 ID而不是单个观测行这是 interactionR 在 GEE 场景下的一个潜在局限。那怎么办我的妥协方案是如果样本量足够大、每个 ID 的观测次数都一致按照行抽样和按照 ID 抽样结果差异不大但如果你发现结果不稳定那就用另一个思路——自己写一个小循环按 ID 重抽样然后每次重新拟合 GEE、计算 RERI最后取分位数。这段代码不复杂但确实需要额外处理不是 interactionR 一键能搞定的。对论文严谨性要求高的朋友我建议在这个问题上多花点时间。5. 常见问题与排查技巧实录这部分是我踩坑三年的“血泪总结”挑最容易被卡住的几个问题说透。5.1 暴露变量的编码问题为什么 0/1 和 1/2 的 RERI 差好远太常见了。有人用as.factor创建变量对照组是 1、暴露组是 2。这种 1/2 编码放进模型后两个主效应的参照就变了乘积项的系数也不再代表“0/1”编码下的交互强度。我见过有人用 1/2 编码跑出来 RERI 为负以为存在拮抗作用后来改成 0/1 就变成正的显著交互了。这纯粹是编码尺度导致的假象。解决方法是在拟合模型之前显式把所有二分类暴露变量转成数值 0/1。data$X1 - ifelse(data$X1 exposed, 1, 0) data$X2 - ifelse(data$X2 exposed, 1, 0)转完之后跑一次table(data$X1, data$X2)确认四个格子的数值都符合预期。再没有效率的做法也比结果错了强。5.2 连续型暴露变量能算相加交互吗理论上可以但实际非常不建议。RERI/AP/SI 这几个指标是建立在“二分类暴露”的 2×2 四格表框架上的。如果你把连续变量丢进去interactionR会把它当成单个暴露水平来算输出的 RERI 没有一个“暴露 vs 不暴露”的实际解释。很多论文在做交互时用的是 BMI 分成两组比如 ≥24 vs 24、吸烟包年数分组比如 ≥20 vs 20正是为了迁就这个框架。所以你在设计分析时从开始就要想清楚暴露变量的切分方式最好在统计分析计划里提前写出来免得事后做多种切分被审稿人批评多重比较。如果实在想保留连续变量的连续属性建议用rms包的回归样条加上张量积光滑项画三维交互效应图——那套方法跟 RERI 就不是一个路线了。5.3 Bootstrap 卡住或者报错常见原因与排除顺序互动R包进行 bootstrap 时最常碰到的报错有这几类Error: object X1 not found说明你的数据框里有缺失值或者模型公式里用了$符号比如data$X1interactionR 在重抽样后找不到对应列。解决办法统一用数据框列名写公式拟合前na.omit。Error in model.frame.default模型对象本身有缺失或者因子水平不平衡interactionR 没法从模型对象里正确提取设计矩阵。解决办法用model.matrix看看设计矩阵是否完整确认每个暴露组合都有足够的观测数。convergence warning from glmer前面说过了换优化器或者简化随机效应结构先保证原模型收敛再谈 bootstrap。长时间不出结果如果样本量大boot_p TRUE时 bootstrap 可能要跑几分钟。这算是正常的。但如果你感觉异常慢先看是不是组数特别多或者模型太复杂。可以按 ESC 中断然后在interactionR调用里把boot参数调低默认是 1000可以写成boot 500或者用ci.fit TRUEboot_p FALSE先拿 delta 置信区间顶一版。5.4 审稿人要 RERI/AP/SI但你在多层逻辑回归里怎么解释很多审稿人看到 RERI 的“OR”是从混合模型GLMM里来的就会追问你这 OR 是条件 ORconditional OR还是边际 ORmarginal ORRERI 的理论框架最早是针对普通逻辑回归的 OR 推导的。在 GLMM 里模型给的是给定随机效应水平下的条件 OR不是人群平均的边际 OR。严格来讲把条件 OR 代入 RERI 公式数学上是缺乏严格支撑的。那该怎么办我自己通常是这么做如果主要分析选了 GLMM我会在论文的敏感性分析里同跑一个没有随机效应的普通逻辑回归或者选择 GEE 作为替代比较两种方法下 RERI 的方向和显著性是否一致。两边都指向同一个结论审稿人就不会再揪着条件 OR 和边际 OR 的区别不放。这个方法不是万能的但至少表现出了你的统计自觉。GEE 在这一点上反而有优势——它天然估计的是边际 OR所以 RERI 的公式能直接套用。如果你的数据结构和研究问题同时支持 GLMM 和 GEE为了跟 RERI 的理论框架保持一致我个人会优先推 GEE。这个优先偏好在学术上是站得住脚的建议你在方法部分解释一句GEE 提供边际估计与相加交互效应的尺度更匹配。5.5 怎么把结果优雅输出到论文表格最后一个实用技巧把结果转成表格。interactionR 的输出本身是列表我通常手动转成数据框再存成 CSV 或直接在 R Markdown 里用kable输出tidy_result - data.frame( Measure c(RERI, AP, SI), Estimate sapply(result_logit, function(x) x[1]), CI_lower sapply(result_logit, function(x) x[2]), CI_upper sapply(result_logit, function(x) x[3]) ) knitr::kable(tidy_result, digits 2)这样生成的三行表格对应许多医学期刊的 Table 2 格式。假如你还放了森林图那就是完整的“Measures of interaction on an additive scale”模块了。期刊投稿时我一般会在图注里写清楚The additive interaction was assessed by the relative excess risk due to interaction (RERI), the attributable proportion (AP), and the synergy index (SI); 95% confidence intervals were obtained by 1000 bootstrap resampling. 这句几乎不会引起审稿人的额外意见。做了这么多例子的经验告诉我相加交互分析不是难在算那几个数而是难在理解指标背后的尺度逻辑、编码问题、模型选择与置信区间方法。interactionR 包确实能把 5 分钟出结果这件事实现在绝大多数常见模型上但你要真正写好一篇论文、经得起审稿人反复追问还是要把上面这些细节弄明白。最后再分享一个我个人养成的习惯无论多简单的分析我都会用epiR的经典方法手动复算一遍核心三指标再用 interactionR 的结果对照。两套方法一致时我才会放心地把结果写进论文。多花这两分钟换来的是后续不需要返工修改的确定性。