高分SCI中几乎必备有交互效应分析表，通常为文章的表三或者表五，上一篇文章我们已经讲过怎么使用SPSS做出交互效应分析，今天我们R手把手实战做出下图这篇20分的文章的交互效应表。

我们先要看他是怎么分析的，他把血小板进行4分组和2分组后，在进行血小板和叶酸的交互作用分析，从各个层面分析了血小板对叶酸的影响，所以别人能发高分不是没有原因的。在模型中，他还使用了调整模型和未调整模型进行比较，支持自己的证据力度。我们还是沿用上次的数据来进行分析
把数据导入并查看数据
library(foreign)
library(survival)
bc <- read.spss(“E:/r/Breast cancer survival agec.sav”,
use.value.labels=F, to.data.frame=T)
bc <- na.omit(bc)
head(bc)

假如我们研究的是乳腺癌各个pr(孕激素受体状态)和死亡的关系，我们想知道pr(孕激素受体状态)和年龄是否产生交互作用。
先使用CUT函数进行分割

age2<-cut(bc$age,breaks=c(0,40,60,100),include.lowest=T,
labels = c(1,2,3))#把age划分为0-20，20-60，60到100这样3个区间
dd<-cbind(bc,age2)#把变量加入表格

刚才是3分组的，再来一个2分组的

Age3<-cut(bc$age,breaks=c(0,40, 100),include.lowest=T,
          labels = c(1,2))
df<-cbind(bc,age2)

然后通过函数subset和if函数设定子集，每个分类变量都要设定一个子集
我们先来做3分类变量的age2,它的3个标签是，1,2,3，简单设置一下
目前age2是因子变量，我们把它转换成数字的

n.age2<-as.numeric(as.character(df$age2))
df<-cbind(df,n.age2)

然后再对它取其中是1的子集

agec1<-subset(df,df$n.age2<=1) #就是age2中取1的子集

得出这个子集后就对他进行分析
先来个广义线性方程得个回归系数，

f1<-glm(status ~pr,family = binomial，)
summary(f1)

计算OR和95%CI

exp(confint(f1))
exp(coef(f1))

由此已经算出年龄0-40岁的亚组的OR，95%CI和P值，同理可以算出其他的组
现在我们来做2分组的交互效应
交互效应就不能拿子集了，要用原来的数据

f2<-glm(status ~pr+age3,family = binomial,data =df)
summary(f2)

求交互效应其实很简单，加入一个a*b就可以了,箭头所示就是交互的P值了

f3<-glm(status ~pr+age3+pr*age3,family = binomial,data =df)

调整模型就是在模型中加入其他变量就可以了,我这里就调整了er，其实做出来用到R的代码并不难，关键是你要思路清晰，明白它是怎么做出来的,

动动小手关注一下吧，更多精彩文章尽在零基础说科研

交互作用效应(p for Interaction)在SCI文章中可以算是一个必杀技，几乎在高分的SCI中必出现，因为把人群分为亚组后再进行统计可以增强文章结果的可靠性，进行可视化后可以清晰的表明变量之间的关系。不仅如此，交互作用还可以使用来进行数据挖掘。在既往文章中，我们已经介绍了怎么使用R语言和SPSS对logistic回归亚组交互效应（交互作用）进行可视化分析（见下图）


后台有粉丝问能不能进行COX回归交互效应可视化分析和对连续变量进行交互效应可视化分析，都是可以的，下面我们来一一进行演示。
COX回归交互效应可视化分析和logistic回归交互项的可视化分析的步骤基本一样的，继续使用我们的乳腺癌数据（公众号回复：乳腺癌，可以获得数据）我们先导入数据和R包

library(foreign)
library(visreg)
library(survival)
bc <- read.spss("E:/r/test/Breast cancer survival agec.sav",
                use.value.labels=F, to.data.frame=T)

我们先来看看数据：
age表示年龄，pathsize表示病理肿瘤大小（厘米），lnpos表示腋窝淋巴结阳性，histgrad表示病理组织学等级，er表示雌激素受体状态，pr表示孕激素受体状态，status结局事件是否死亡，pathscat表示病理肿瘤大小类别（分组变量），ln_yesno表示是否有淋巴结肿大，time是生存时间，后面的agec是我们自己设定的，不用管它。

接下来删除缺失变量和把分类变量转成因子

bc <- na.omit(bc)
bc$er<-as.factor(bc$er)
bc$pr<-as.factor(bc$pr)
bc$ln_yesno<-as.factor(bc$ln_yesno)
bc$histgrad<-as.factor(bc$histgrad)

建立cox回归方程，我们假设age和histgrad有交互

f1<- coxph(Surv(time, status) ~ age + histgrad+
           pathsize+ er + pr + ln_yesno+age*histgrad, data=bc)
summary(f1)

这里发现age：histgrad的P值大于0.05，两者之间没有交互效应，我们不管他，这里主要是演示怎么做图

visreg(f1, "age", "histgrad", ylab="log(Hazard ratio)")

可以把3张图合并起来

visreg(f1, "age", "histgrad", ylab="log(Hazard ratio)", overlay = T)

还可以进一步修改

plot(visreg(f1,xvar = "age",by="histgrad",plot=F),xlab="age",ylab="log(Hazard ratio)",
     overlay = T,partial = F,rug=F,,line=list(lty=1:6))
legend("bottomright",
       c("histgrad1","histgrad2","histgrad3"),
       lty=c(1,1,1),
       col=c("red","green","blue"),
       lwd=c(2,1,1),
       bty="n")

也可以按箱线图来表示

plot(visreg(f1, "histgrad", by="age",plot=F),xlab="age",ylab="log(Hazard ratio)",
     overlay = T,partial = F,rug=F,,line=list(lty=1:6))
legend("bottomright",
       c("age38","age55","age74"),
       lty=c(1,1,1),
       col=c("red","green","blue"),
       lwd=c(2,1,1),
       bty="n")

从图上中数据的可信区间严重重合也可以看出来两者并没有交互。
下面来演示连续变量的交互可视化分析，使用我们的臭氧数据（公众号回复：臭氧可以获得数据）
我们先导入R包和臭氧数据

library(splines)
be <- read.spss("E:/r/test/ozone.sav",
                use.value.labels=F, to.data.frame=T)

数据中有七个变量，ozon每日臭氧水平为结局变量，Inversion base height（ibh）反转基准高度，Pressure gradient (mm Hg) 压力梯度（mm Hg），Visibility (miles) 能见度（英里），Temperature (degrees F) 温度（华氏度），Day of the year日期，vh我也不知道是什么，反正就是一参数，这里所有的变量都是连续的。
假设我们想知道反转基准高度和温度对臭氧浓度有无存在交互影响

fit <- lm(ozon ~ vh +ns(ibh, df=2)*ns(temp, df=2), data=be)
summary(fit)

visreg2d(fit, "ibh", "temp")

也可以自己更改想要的颜色

visreg2d(fit, "ibh", "temp", color.palette=colorRampPalette(c("black", "white", "purple")))

也可以绘制3D立体图

visreg2d(fit,"ibh", "temp", plot.type="persp")

也可以绘制成这种

plot(visreg(fit,xvar = "ibh",by="temp",plot=F),xlab="age",ylab="log(Probability)",
     overlay = T,partial = F,rug=F,,line=list(lty=1:6))
legend("topright",
       c("temp41","temp62","temp81"),
       lty=c(1,1,1),
       col=c("red","green","blue"),
       lwd=c(2,1,1),
       bty="n")

1. 交互作用概述

(1) 概念：
某一因素的真实效应（单独效应）随着另一因素水平的改变而改变。当两种或两种以上暴露因素同时存在时所致的效应不等于它们单个作用相联合的效应时，则称因素之间存在交互作用。
① 因素A的效应在因素B的不同水平上存在差异，则认为因素A、B之间存在交互作用。
② 因素A、B的联合效应不等于两因素独立效应之和或之积。

(2) 识别：
①统计表：所有可能交叉组合情况下的结局指标；
②统计图：观察另一因素各水平下，某因素与结局指标关系图是否平行（Addictive interation）。

(3) 理解交互作用：
在不同B因素水平，A因素的效应有统计学差异；同理，在不同A因素水平，B因素的效应也有差异。交互作用就是看二者合在一起的效应，与二者单独效应之和（相加）/之乘（相乘），相等还是不相等。不相等，那就是有交互。Mediation是解释暴露因素how导致的结局，而interaction是说明who是高风险目标人群。
首先需明确，“交互作用”是指统计学上的交互作用，不能直接说它们具有生物学意义上的交互作用。因为统计模型不止一种，你会看到有些研究中，相乘模型有统计学意义，相加模型没有统计学意义，或者相乘模型为拮抗作用，相加模型为协同作用；所以，这才有了谁更合理的争论。根据“A因素和B因素在相加尺度上具有交互作用”，可以推测有生物学交互，而不能直接写“具有生物交互作用”。同样，没检测到统计学上的交互作用，也不能说明没有生物学交互作用。⚠️注意：有无统计学交互作用很大程度上取决于所选择的模型，所以在报告有无交互作用时，需要说明分析所用的模型。

(4) 相加交互与相乘交互的差异：
统计建模中一般线性模型交互项反映的是因素间相加交互作用（additive interactions，INT_A），而logistic和Cox等广义线性模型则反映因素间统计学上相乘交互作用（multiplicative interaction，INT_M），即logistic乘积项仅反映统计学上的交互，而只有在生物学机制上病因因素间存在相加交互作用才可解释为相加交互作用。

(5) 评价交互作用的目的和意义：
①在资源有限的情况下，识别对某干预受益最大的人群（亚组）。
②在某一主要暴露因素不容易施加干预的情况下，识别最有可能施加干预的互作协变量以降低主要暴露因素的效应。
③揭示暴露影响疾病发生的机制。
④提高评价某暴露因素对结局影响总效应的把握度。
⑤仅从统计学角度考虑，通常包含交互作用项的模型拟合数据更好。

2. 二分类解释变量交互作用

2.1 相乘交互作用

(1) 相乘交互作用的定义：假设研究多风险因素中交互作用的两暴露因素为AB，则OR₀₀表示AB均无暴露，即OR₀₀＝1；OR₁₀表示A暴露、B无暴露，OR₀₁表示A无暴露、B暴露，OR₁₁表示A、B均暴露。则相乘交互作用INT_M＝OR_A×B＝OR₁₁/（OR₁₀×OR₀₁）。

(2) 相乘交互作用的判定：logistic 等回归乘积项95%CI不包含1，表明有相乘交互作用。交互项得到的OR值<1，拮抗作用；交互项 OR>1，协同作用。

2.2 相加交互作用

(1) 定量评价流行病学研究中暴露因素间及暴露因素与基因间相加交互作用需要3项指标：交互对比度（interaction contrast ratio，ICR）又称交互作用超额相对危险度（relative excess risk due to interaction，RERI），交互作用归因比（attributable proportion due to interaction，AP）和协同指数（the synergy index，S）。

(2) 相加交互作用3项指标定义为：
RERI = ICR = OR₁₁ - OR₁₀ - OR₀₁ + 1
AP = ICR/OR₁₁
S = (OR₁₁ - 1)/(OR₁₀ + OR₀₁)=(OR₁₁ - 1)/(OR₁₁ + 1 - ICR)

此处公式有误，之前一直懒得去核实，现在做个修改。感谢x小x小姐姐提供的文献参考，如有需要，请自行下载原文。—— Revised: 2022-04-01
Knol MJ, VanderWeele TJ, Groenwold RH, Klungel OH, Rovers MM, Grobbee DE. Estimating measures of interaction on an additive scale for preventive exposures. Eur J Epidemiol. 2011 Jun;26(6):433-8.

在二分类logistic回归模型中，ln[P/(1-P)] = β0 + β1A + β2B + β3A×B，OR₁₀=exp(β1)，OR₀₁=exp(β2)，OR₁₁=INT_A=exp(β1 + β2 + β3)

(3) 相加交互作用的判定：如果两因素有相加交互作用，则RERI 95%CI、AP 95%CI应不包含0，S 95%CI应不包含1。超额相对危险度RERI以及归因比AP均>0，且可信区间不包括0， S>1且可信区间不包括1，表示存在交互作用，且为协同作用。RER1以及AP均<0，S<1，表示存在交互作用，且为拮抗作用。

2.3 计算置信区间CI的方法

Delta法、Wald法、轮廓似然置信区间（profile likelihood confidence intervals，PL）法、variance recovery method和percentile bootstrapping等。

3. R 代码示例

epiR R包可实现将两个二元解释变量作为交互项纳入logistic回归，计算其乘法和加法交互作用（RERI、AP、S）。置信区间的计算基于Hosmer和Lemeshow(1992)描述的delta方法。

(1) R 包安装及数据介绍：

install.packages("epiR")
library(epiR) 
## Data from Rothman and Keller (1972) evaluating the effect of joint exposure to alcohol and tabacco on risk of cancer of the mouth and pharynx (cited in Hosmer and Lemeshow, 1992):
can <- c(rep(1, times = 231), rep(0, times = 178), rep(1, times = 11), 
         rep(0, times = 38))
smk <- c(rep(1, times = 225), rep(0, times = 6), rep(1, times = 166), 
         rep(0, times = 12), rep(1, times = 8), rep(0, times = 3), rep(1, times = 18), 
         rep(0, times = 20))
alc <- c(rep(1, times = 409), rep(0, times = 49))
dat <- data.frame(alc, smk, can)

#因子化前后，回归分析结果一致
dat$smk <- factor(dat$smk)
dat$alc <- factor(dat$alc)
dat$can <- factor(dat$can)
summary(dat)

(2) 相乘交互作用及二元Logist回归：

# 1) 相乘交互作用及二元Logist回归：
fit <- glm(can ~ alc + smk + alc:smk, family = binomial(link = "logit"), data = dat)
summary(fit)$coefficients
coef <- summary(fit)$coefficients[,1]
se <- summary(fit)$coefficients[,2]
Results <- cbind(exp(coef),exp(coef-1.96*se),exp(coef+1.96*se))
exp(confint(fit)) #另外一种形式显示可信区间
P <- summary(fit)$coefficients[4,4]
dimnames(Results)[[2]] <- c("OR", "lower","upper")
Results

#                   OR      lower      upper
#(Intercept) 0.1500000 0.04457283  0.5047918
#alc1        3.3333333 0.70058649 15.8597280
#smk1        2.9629630 0.68002749 12.9099922
#alc1:smk1   0.9149096 0.15435605  5.4229143

从结果可知，饮酒与吸烟的乘法交互效应为0.91<1，交互P值为0.92，不显著。

(3) 相加交互作用及二元Logist回归：

# 2) 相加交互作用及二元Logist回归：
## Table 2 of Hosmer and Lemeshow (1992):
dat.glm01 <- glm(can ~ alc + smk + alc:smk, family = binomial, data = dat)
summary(dat.glm01) #P interaction= 0.92197

## What is the measure of effect modification on the additive scale?
## RERI
epi.interaction(model = dat.glm01, param = "product", coef = c(2,3,4), type = "RERI", conf.level = 0.95)
## Measure of interaction on the additive scale: RERI 3.73 
## (95% CI -1.84 -- 9.32), page 453 of Hosmer and Lemeshow (1992).

#AP 和 S
epi.interaction(model = dat.glm01, param = "product", coef = c(2,3,4), type = "APAB", conf.level = 0.95)
#        est       lower     upper
#1 0.4138765 -0.07306308 0.9008162
epi.interaction(model = dat.glm01, param = "product", coef = c(2,3,4), type = "S", conf.level = 0.95)
#     est     lower    upper
#1 1.870482 0.6460433 5.415585

结果表明，加法交互的3个指标RERI、AP和S均不显著，推测可能不存在生物学交互作用。

(4) 哑变量实现交互作用：

# 3) 哑变量实现交互作用
## Rothman defines an alternative coding scheme to be employed for parameterising an interaction term. Using this approach, instead of using two risk factors and one product term to represent the interaction (as above) the risk factors are combined into one variable with (in this case)
## four levels:
## a.neg b.neg: 0 0 0
## a.pos b.neg: 1 0 0
## a.neg b.pos: 0 1 0
## a.pos b.pos: 0 0 1

dat$d <- rep(NA, times = nrow(dat))
dat$d[dat$alc == 0 & dat$smk == 0] <- 0
dat$d[dat$alc == 1 & dat$smk == 0] <- 1
dat$d[dat$alc == 0 & dat$smk == 1] <- 2
dat$d[dat$alc == 1 & dat$smk == 1] <- 3
dat$d <- factor(dat$d)

## Table 3 of Hosmer and Lemeshow (1992):
dat.glm02 <- glm(can ~ d, family = binomial, data = dat)
summary(dat.glm02)

# 1> 乘法交互尺度上的效应修饰作用
## What is the measure of effect modification on the multiplicative scale?
## See VanderWeele and Knol (2014) page 36 and Knol and Vanderweele (2012) for details.
beta1 <- as.numeric(dat.glm02$coefficients[2])
beta2 <- as.numeric(dat.glm02$coefficients[3])
beta3 <- as.numeric(dat.glm02$coefficients[4])
exp(beta3) / (exp(beta1) * exp(beta2))
## Measure of interaction on the multiplicative scale: 0.92.

# 2> 加法交互
## What is the measure of effect modification on the additive scale?
# coef: a vector listing the positions of the coefficients of the interaction terms in the model. 
epi.interaction(model = dat.glm02, param = "dummy", coef = c(2,3,4), type = "RERI", conf.level = 0.95) #超额相对危险度
## Measure of interaction on the additive scale: RERI 3.73 
## (95% CI -1.84 -- 9.32), page 455 of Hosmer and Lemeshow (1992).

## 计算AP、S指标
epi.interaction(model = dat.glm02, param = "dummy", coef = c(2,3,4), type = "APAB", conf.level = 0.95) #归因比
epi.interaction(model = dat.glm02, param = "dummy", coef = c(2,3,4), type = "S", conf.level = 0.95) #协同指数
# Skrondal (2003) advocates for use of the synergy index as a summary measure of additive interaction, showing that when regression models adjust for the effect of confounding variables (as in the majority of cases) RERI and AP may be biased, while S remains unbiased.

(5) 结果可视化

# 4) 可视化
# 1> 加法交互效应
ORalc <- 3.3333333
ORsmk <- 2.9629630
RERI <- 3.739848
bar_d <- matrix(c(1, 1, 1, 1,
                  ORalc-1, 0, ORalc-1, 0,
                  ORsmk-1 ,ORsmk-1, 0, 0,
                  RERI, 0, 0, 0),
                c(4,4), byrow = T,
                dimnames = list(c('U','alc','smk1','alc1 & smk1'),c("OR_A1B1","OR_A1B0","OR_A0B1","OR_A0B0")))

plot <- barplot(bar_d, legend=rownames(bar_d))

该图直观地显示了相加交互作用的大小：

# 2> 加法交互是否显著
library(visreg)
plot(visreg(fit,xvar = "alc",by="smk",plot=F),xlab="alc",ylab="predict",overlay = T,partial = F,rug=F)
legend("topleft", c("Smoke: No","Smoke: Yes"), lty=c(1,1), col=c("red","blue"), lwd=c(1,1),bty="n")

上述分析显示，RERI、AP的可信区间包含0，S的可信区间包含1，说明加性交互作用不显著。图示如下：

4. 实例分析

(1) 问题描述及数据：
以2004-2014年广东省先天性心脏病（先心）监测网数据库病例对照资料为例[10, 11]，分析孕母是否被动吸烟（自报围孕期家庭、工作任一环境中接触吸烟平均时长＞15 min/d且持续1周以上）与家庭人均月收入有无相乘和相加交互作用（表 1），并比对采用Andersson分析方法的结果。

(2) 结果解释：
INT_M 95%CI＞1有相乘交互作用，A暴露（如，孕期被动吸烟）与B暴露（如，家庭人均月收入低）有正向相乘交互作用；Delta法、Wald法、PL法估计 ICR 95%CI＞0、AP 95%CI＞0，S 95%CI＞1，A暴露与B暴露有协同相加交互作用，即同时暴露于A、B者，其结局（如，胎儿患先心）的发生风险显著增加。ICR和S的意义相同，AP表示全部病例中可归因于两因素交互作用的病例所占比例，如，AP＝0.324（95%CI：0.153～0.519），说明全部先心病例中归因于被动吸烟与家庭人均月收入低的相加交互作用所引起的病例占32.4%。

(3) 总结：
统计学交互作用的线性模型一般为加法模型，乘积项表示有无相加交互作用；而广义线性属于乘法模型，乘积项表示有无相乘交互作用。生物学交互作用是多风险因素在发病的生物机制上的定性概念，两因素皆为病因前提下生物学机制的相互联系，包括协同和拮抗，生物学交互不同于统计模型中乘积项的分析。故相乘交互项OR值95%CI＜1为负相乘，OR值95%CI＞1为正相乘，相加交互项ICR值95%CI＞0为协同作用，＜0为拮抗作用。有时相乘与相加交互方向相悖，此时更宜运用生物学交互。在置信区间估计方法的选择上，Andersson等使用Delta法估计，结果不够稳健易造成假阴性结果。在风险比不对称时，应采用PL法，因Wald法的结果不稳定。

参考阅读：
[1] 聂志强, 欧艳秋, 庄建, 曲艳吉, 麦劲壮, 陈寄梅, 刘小清. 实现logistic与Cox回归相乘相加交互作用的临床实践宏程序. 中华流行病学杂志, 2016, 37(5): 737-740.
[2] 交互作用: 相加交互，相乘交互
[3] epiR: Tools for the Analysis of Epidemiological Data. R package version 2.0.19.
[4] https://max.book118.com/html/2016/1031/60871118.shtm
[5] https://max.book118.com/html/2017/0216/92080186.shtm