企业官网建设流程全解析

交互效应这东西，做流行病学、临床研究和社科数据分析的人几乎天天碰到。但大多数教程讲的是“乘积交互项”，也就是在模型里加 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 核心原理：从拟合模型到计算指标

整个分析的原理可以压缩成四步：

1. 拟合回归模型，包含因素 A、因素 B 以及 A×B 乘积项；
2. 从模型中提取两个主效应项的回归系数和乘积项的回归系数；
3. 根据公式计算三个相加交互指标的估计值；
4. 计算置信区间，目前主流方法有 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 的时候，需要自己动手写一个计算函数，提取固定效应系数和方差协方差矩阵，再做估计。这部分代码也不难，后面我会给出可直接复用的版本。

所以我的工具选型策略是：

2.3 四类模型的适用场景差异

为什么需要四种模型分别处理？核心是数据结构和研究设计不同：

• 逻辑回归：结局是二分类，数据独立，最经典的病例对照或横断面研究。
• Cox 回归：结局是生存时间（含删失），用于队列研究和临床试验的生存分析。
• GLMM：数据有层级结构（如多中心、重复测量），需要在模型里加入随机效应。这里的“相加交互”指标要用固定效应部分的系数计算，不能把随机效应的方差算进去。
• GEE：同样是处理相关数据，但更强调“群体平均效应”，常用于纵向数据或多层聚类数据。GEE 给的是稳健标准误，计算置信区间时直接用稳健方差矩阵即可。

在实际操作中，如果你只是验证一下逻辑回归的交互结果，epiR一把梭没问题；但如果你的数据是纵向的，忽略相关性直接跑普通逻辑回归，很可能导致标准误偏小、假阳性上升，换 GLMM 或 GEE 是更严谨的选择。

3. 实战：5分钟跑通完整的相加交互分析

3.1 数据准备与变量编码

先说编码问题，这一步错了后面全白搭。计算 RERI 时，两个暴露因素必须是二分变量，而且编码方向要和你的研究假设一致。我习惯把“暴露组”编码为 1，“非暴露组”编码为 0。

以吸烟和高血压为例：

• smoke：1 = 吸烟，0 = 不吸烟；
• htn：1 = 有高血压，0 = 无高血压；
• y：1 = 发生心血管事件，0 = 未发生；
• time：生存时间（Cox 用）；
• id：个体 ID（GLMM/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：是否输出暴露-混杂四格表形式的估计结果，建议设为TRUE；
• ci：是否计算置信区间，建议TRUE；
• conf.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.82（95% CI: 0.15, 1.49），AP = 0.31（95%CI: 0.10, 0.52），S = 1.75（95%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 有两种思路：

1. 对个体（聚类单元）做 Bootstrap，每次重抽样后重新拟合 GEE；
2. 用基于渐近正态的 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自己画——把四个 OR（OR11、OR10、OR01、参照组 1）和 RERI、AP、S 的点估计加置信区间放到一张图里，审稿人看了会非常直观。

三是和亚组分析结合。比如按性别分层，分别算男性和女性的 RERI，再比较两个 RERI 的差异。这个玩法在临床研究里很受欢迎，因为它能回答“协同效应是否在某个亚组里更强”。

最后再分享一个小技巧：不管你用哪种模型，跑完后顺手把sessionInfo()里epiR和lme4的版本记录下来，写论文方法部分时直接引用。审稿人问到具体算法时，你可以明确回答用的是哪个包的哪个版本，这种细节在返修时非常加分。我自己的习惯是所有交互分析代码都放在同一个 R 脚本里，注释写好数据和模型类型，三个月后回来看还能直接复现——这才是“几分钟搞定”的底气所在。