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)
# 1) 相乘交互作用及二元Logist回归:
fit <- glm(can ~ alc + smk + alc:smk, family = binomial, 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
# 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) 哑变量实现交互作用
## 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.
# 4) 可视化
bar_d <- matrix(c(1,1,1,1,
3.3333333,0,3.3333333,0,
2.9629630 ,2.9629630 ,0,0,
3.739848, 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))
联系客服