【交互作用】02. 加法交互 & 乘法交互 (R包 interactionR)

生物学交互作用的评价应该基于是否有相加交互作用， 而流行病学研究中常运用logistic和Cox等广义线性模型， 并纳入乘积项分析因素间交互作用，其是否有意义仅反映相乘交互作用， 并不能反映两因素间相加或生物学交互作用的有无。上篇【交互作用】01. 加法交互 & 乘法交互 (R包 epiR) 介绍了交互作用的基本概念、三个相加交互作用评价指标（RERI、AP和S）和 epiR R包的应用等。本篇内容主要介绍实现交互作用的另一个好用的R包 interactionR。

1. 原理和方法回顾

以两因素两水平为例。假设两暴露因素分别为A、B，1表示因素存在，0表示因素不存在，因变量为疾病发生与否。

OR_A0B0表示A、Ｂ都不存在时发病的OR值，分析时作为参照组；OR_A1B0表示仅A存在、Ｂ不存在时发病的OR值；OR_A0B1表示A不存在、仅Ｂ存在时发病的OR值；OR_A1B1表示A、Ｂ共同存在时发病的OR值。

评价相加交互作用的三个指标：
① 相对超危险度比：RERI = OR_A1B1 - OR_A0B1 - OR_A1B0 + 1；
② 归因比：AP = RERI / OR_A1B1；
③ 交互作用指数：SI = (OR_A1B1 - 1) / ［(OR_A0B1 - 1) + (OR_A1B0 - 1)］。

判定准则：如果两因素无相加交互作用，则RERI和AＰ的可信区间应包含0，SI的可信区间应包含1。

交互作用指标的点估计：
可通过以下两种方法，建立logistic回归模型计算OR_A1B1、OR_A0B1和OR_A1B0，代入交互作用指标的计算公式。
(1) 用两因素A、B及乘积项A×B构建logistic回归模型1：
ln[p/(1-p)] = β0 + β1A + β2B + β3(A×B)；
OR_A1B0 = exp(β1)，
OR_A0B1 = exp(β2)，
OR_A1B1 = exp(β1 + β2 + β3)。
(2) 根据两因素A、B，建立新的交互作用哑变量A_B，构建logistic回归模型2：
A0B0表示A=0且B=0，分析时作为参照组；A0B1表示A=0且B=1，A0B1表示A=0且B=1，A1B1表示A=1且B=1。
ln[p/(1-p)] = β0 + β1A1B0 + β2A0B1 + β3A1B1；
OR_A1B0 = exp(β1)，
OR_A0B1 = exp(β2)，
OR_A1B1 = exp(β3)；
模型1和2中的β1、β2相同，而模型2中的β3等于模型1中的β1 + β2 + β3。

注意：一般以高风险的一类作为暴露组，尤其是在保护因素时，应当将无暴露设置为1，有暴露设置为0，以避免解释上混乱。当暴露变量为多分类或连续变量时，置信区间（CI）估计的delta方法以及Andersson编制的Excel法均不适用。而Bootstrap方法可以估计所有解释变量类型的交互情况的CI。

2. interactionR R包介绍

interactionR 可直接导出出版级别的Word，结果包括联合效应、暴露效应和交互效应，乘法交互，加法交互，以及分层分析。

R包内置3种置信区间的估计方法：

• delta method (Hosmer and Lemeshow (1992), [doi:10.1097/00001648-199209000-00012]),
• variance recovery method (Zou (2008), [doi:10.1093/aje/kwn104]),
• percentile bootstrapping (Assmann et al. (1996), [doi:10.1097/00001648-199605000-00012]).

R包的安装：

devtools::install_github("epi-zen/interactionR")
library(interactionR) 

参数介绍：

• model: 包含交互项的回归模型。可以是logistic glm(formula, family = binomial(link = "logit"), data)、条件logistic clogit() 或cox coxph()回归模型。模型可以包含适当的协变量。
• exposure_names: 模型中可能存在交互作用的两个二分类变量（乘积项）。
• ci.type: 加法交互作用的置信区间估计方法 (“delta” 或 “mover”) ，默认为delta方法。
• ci.level: 置信区间水平。
• em: TRUE, for effect modification assessment. FALSE, for interaction.
• recode: If TRUE, recodes the exposures - if at least one of the exposures is protective - such that the stratum with the lowest risk becomes the new reference category when the two exposures are considered jointly (See Knol et al (2011) [doi: 10.1007/s10654-011-9554-9]).

示例数据代码实现：

已报道的饮酒和吸烟对口腔癌的联合作用数据。包括两个二分类暴露因素，即饮酒 (alc) 和吸烟 (smk)，结局为二分类变量，即口腔癌 (oc)。

data (OCdata)

## fit the interaction model
model.glm <- glm(oc ~ alc*smk, family = binomial(link = "logit"), data = OCdata)
summary(model.glm)

## 1) analysis----
table_object = interactionR(model.glm, exposure_names = c("alc", "smk"), ci.type = "mover", ci.level = 0.95, em = F, recode = F)
table_object$dframe 
# Measures Estimates       CI.ll      CI.ul
# 1                        OR00 1.0000000          NA         NA
# 2                        OR01 2.9629630   0.6800459 12.9096430
# 3                        OR10 3.3333333   0.7006066 15.8592734
# 4                        OR11 9.0361446   2.6413389 30.9130753
# 5  OR(smk on outcome [alc==0] 2.9629630   0.6800459 12.9096430
# 6  OR(smk on outcome [alc==1] 2.7108434   0.9969750  7.3709689
# 7  OR(alc on outcome [smk==0] 3.3333333   0.7006066 15.8592734
# 8  OR(alc on outcome [smk==1] 3.0496988   1.2948765  7.1826638
# 9        Multiplicative scale 0.9149096   0.1543611  5.4227369
# 10                       RERI 3.7398483 -11.4297248 21.8721579
# 11                         AP 0.4138765  -0.3775073  0.8113231
# 12                         SI 1.8704819   0.6460433  5.4155854

# 乘法：
# 饮酒作用：OR10/OR00=3.33
# 吸烟作用：OR01/OR00=2.96
# 烟酒联合作用：OR11/OR00=9.04
# 相乘作用：OR11/(OR01*OR10)=0.9149096(Multiplicative scale)=exp(-0.08893)(alc:smk)=exp(2.2012)/exp(1.0862)/exp(1.2040)
# P interaction= 0.92197(即alc:smk的P值)

# 加法：
# 饮酒作用：OR10-OR00=3.33-1=2.33
# 吸烟作用：OR01-OR00=2.96-1=1.96
# 烟酒联合作用：OR11-OR00=8.0361446
# 相加作用：OR11-OR01-OR10+1=9.04-3.33-2.96+1=3.7398483(RERI)

## 2) output word----
interactionR_table(table_object) #将分层分析与乘法和加法交互作用的结果输出到Word。

## CI estimating methods: 
CI1 <- interactionR_mover(model.glm, exposure_names = c("alc", "smk"), ci.level = 0.95, em = F, recode = F);CI1$dframe
CI2 <- interactionR_delta(model.glm, exposure_names = c("alc", "smk"), ci.level = 0.95, em = F, recode = F);CI2$dframe 
CI3 <- interactionR_boot(model.glm, ci.level = 0.95, em = F, recode = F, seed = 12345, s = 1000);CI3$dframe 
# interactionR_boot() #A fitted model object of class glm. Requires that the two binary exposure variables are listed first in the call formula.

输出到Word中的表格如下：

参考阅读：
[1] 许敏锐,强德仁,周义红,石素逸,秦晶,陶源.应用R软件进行logistic回归模型的交互作用分析[J].中国卫生统计,2017,34(04):670-672+675.
[2] Babatunde Alli (2021). interactionR: Full Reporting of Interaction Analyses. R package version 0.1.2. https://CRAN.R-project.org/package=interactionR.