方法一
基于R包epiR: Tools for the Analysis of Epidemiological Data
https://cran.r-project.org/web/packages/epiR/index.html
一、R包加载，以及示例数据构建
rm(list = ls())
#install.packages("epiR")
library(epiR)
## EXAMPLE 1:
## 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)
二、相乘交互作用
fit <- glm(can ~ alc + smk + alc:smk, family = binomial, data = dat)
summary(fit)$coefficients
coef <- summary(fit)$coefficients[,1]
se <- summary(fit)$coefficients[,2]
#CI <-exp(confint(fit))#另外一种形式显示可信区间
Results <- cbind(exp(coef),exp(coef-1.96*se),exp(coef+1.96*se))
P <- summary(fit)$coefficients[4,4]#0.9219744
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，不显著
三、相加交互作用
## Table 2 of Hosmer and Lemeshow (1992):
dat.glm01 <- glm(can ~ alc + smk + alc:smk, family = binomial, data = dat.df01)
summary(dat.glm01)
## What is the measure of effect modification on the additive scale?
epi.interaction(model = dat.glm01, param = "product", coef = c(2,3,4),
conf.level = 0.95)
## Measure of interaction on the additive scale: RERI 3.73
## (95% CI -1.84 to 9.32), page 453 of Hosmer and Lemeshow (1992).
## 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.glm01$coefficients[2])
beta2 <- as.numeric(dat.glm01$coefficients[3])
beta3 <- as.numeric(dat.glm01$coefficients[4])
exp(beta3) / (exp(beta1) * exp(beta2))
## Measure of interaction on the multiplicative scale: 0.093.
四、哑变量实现交互作用
## EXAMPLE 2:
## 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 comprised of
## (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.df01$d <- rep(NA, times = nrow(dat.df01))
dat.df01$d[dat.df01$alc == 0 & dat.df01$smk == 0] <- 0
dat.df01$d[dat.df01$alc == 1 & dat.df01$smk == 0] <- 1
dat.df01$d[dat.df01$alc == 0 & dat.df01$smk == 1] <- 2
dat.df01$d[dat.df01$alc == 1 & dat.df01$smk == 1] <- 3
dat.df01$d <- factor(dat.df01$d)
## Table 3 of Hosmer and Lemeshow (1992):
dat.glm02 <- glm(can ~ d, family = binomial, data = dat.df01)
summary(dat.glm02)
## What is the measure of effect modification on the additive scale?
epi.interaction(model = dat.glm02, param = "dummy", coef = c(2,3,4),
conf.level = 0.95)
## Measure of interaction on the additive scale: RERI 3.74
## (95% CI -1.84 to 9.32), page 455 of Hosmer and Lemeshow (1992).
image.png
image.png
五、结果可视化
加法交互效应，堆叠图
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),
args.legend=c(x=4,y=10)#图示的位置
)
加法交互效应，堆叠图
方法二
基于interactionR
https://github.com/epi-zen/interactionR
rm(list = ls())
#devtools::install_github("epi-zen/interactionR")
library(interactionR)
#Example: The joint effect of alcohol and smoking on oral cancer.
data (OCdata)
## fit the interaction model
model.glm <- glm(oc ~ alc*smk, family = binomial(link = "logit"),
data = OCdata)
table_object = interactionR(model.glm,
exposure_names = c("alc", "smk"),
ci.type = "mover", ci.level = 0.95,
em = F, recode = F)
interactionR_table(table_object)#可以输出成Word格式
输出结果
#Example 2: Effect of sports participation and smoking on herniated lumbar disc.
data(HDiscdata)
m2 = glm(h ~ ns*smk, family = binomial(link = 'logit'),
data = HDiscdata)
table_object2 = interactionR_boot(m2, ci.level = 0.95,
em = F, recode = F,
seed = 1234, s = 1000)
interactionR_table(table_object2)
hist(table_object2$bootstrap)
输出结果
image.png
背景知识学习
https://max.book118.com/html/2018/0607/171193334.shtm
https://www.docin.com/p-467347939.html
https://zhuanlan.zhihu.com/p/138018634
https://max.book118.com/html/2016/1031/60871118.shtm
https://max.book118.com/html/2017/0216/92080186.shtm
https://www.sohu.com/a/155275085_170798
