【交互作用】02. 加法交互 & 乘法交互 [R包 interactionR]

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

1. 原理和方法回顾

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

ORA0B0表示A、B都不存在时发病的OR值,分析时作为参照组;ORA1B0表示仅A存在、B不存在时发病的OR值;ORA0B1表示A不存在、仅B存在时发病的OR值;ORA1B1表示A、B共同存在时发病的OR值。

评价相加交互作用的三个指标:
① 相对超危险度比:RERI = ORA1B1 - ORA0B1 - ORA1B0 + 1;
② 归因比:AP = RERI / ORA1B1;
③ 交互作用指数:SI = (ORA1B1 - 1) / [(ORA0B1 - 1) + (ORA1B0 - 1)]。

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

交互作用指标的点估计:
可通过以下两种方法,建立logistic回归模型计算ORA1B1、ORA0B1和ORA1B0,代入交互作用指标的计算公式。
(1) 用两因素A、B及乘积项A×B构建logistic回归模型1:
ln[p/(1-p)] = β0 + β1A + β2B + β3(A×B);
ORA1B0 = exp(β1),
ORA0B1 = exp(β2),
ORA1B1 = 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;
ORA1B0 = exp(β1),
ORA0B1 = exp(β2),
ORA1B1 = 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.

(0)

相关推荐