Chapter 2 Interaction

We will be using three measures of effect to explain the concept and necessary implementation ideas

2.1 Modelling to obtain OR

odds ratio (OR), Multiplicative scale
- Logit[Pr(Y=1)] = \(\alpha_0\) + \(\alpha_A\) A + \(\alpha_B\) B + \(\alpha_{AB}\) (A \(\times\) B)
- Condition for presence of interaction (“departure from multiplicativity”)
  - OR10 \(\times\) OR01 \(\ne\) OR11
  - \(\exp(\alpha_A)\) \(\times\) \(\exp(\alpha_B)\) \(\ne\) \(\exp(\alpha_A+ \alpha_B + \alpha_{AB})\)
- synergism / super-multiplicative
  - OR11 > OR10 \(\times\) OR01
- antagonism / sub-multiplicative
  - OR11 < OR10 \(\times\) OR01

Note: The above assumes that both exposures are risk factors for the outcome (OR > 1), not protective factors. If protective, estimates of RERI and AP (defined later) will be invalid, although the estimate of SI is not affected by this condition.

2.1.1 Base

2.1.1.1 no alcohol, no smoking for OR(smk1 on outcome [alc1==0] and OR(alc1 on outcome [smk1==0]

Model-based SE estimates from logistic regression model (Naimi and Whitcomb 2020):

require(interactionR)
data(OCdata)
Obs.Data <- OCdata
Obs.Data$smk <- as.factor(Obs.Data$smk)
Obs.Data$smk <- relevel(Obs.Data$smk, ref = "0")
Obs.Data$alc <- as.factor(Obs.Data$alc)
Obs.Data$alc <- relevel(Obs.Data$alc, ref = "0")
fit.or.int11 <- glm(oc ~ alc + smk + alc:smk, 
                    family = binomial(link = 'logit'), 
                    data = Obs.Data)
require(jtools)
results.or.model11 <- summ(fit.or.int11, 
                           model.info = FALSE, 
                           model.fit = FALSE, 
                           exp = TRUE, 
                           confint = TRUE)
results.or.model11

	exp(Est.)	2.5%	97.5%	z val.	p
(Intercept)	0.15	0.04	0.50	-3.06	0.00
alc1	3.33	0.70	15.86	1.51	0.13
smk1	2.96	0.68	12.91	1.45	0.15
alc1:smk1	0.91	0.15	5.42	-0.10	0.92
Standard errors: MLE

2.1.2 Estimates of different ORs

2.1.2.0.1 Alcohol

## OR00 = 1
OR10 <- exp(sum(summary(fit.or.int11)$coef[c('alc1'),'Estimate']))
OR10 ## OR_A=1

## [1] 3.333333

\(\alpha_A\) = 1.2
OR_{A=1} = OR10 = \(\exp(\alpha_A)\) = 3.33

2.1.2.0.2 Smoking

OR01 <- exp(sum(summary(fit.or.int11)$coef[c('smk1'),'Estimate']))
OR01 ## OR_B=1

## [1] 2.962963

\(\alpha_B\) = 1.09
OR_{B=1} = OR01 = \(\exp(\alpha_B)\) = 2.96

2.1.2.0.3 Both alcohol and smoking

OR11 <- exp(sum(summary(fit.or.int11)$coef[c('alc1','smk1','alc1:smk1'),'Estimate'])) 
OR11 ##OR_A=1,B=1

## [1] 9.036145

\(\alpha_{AB}\) = -0.09
OR_{A=1,B=1} = OR11 = 9.04

2.1.2.0.4 Joint effect?

Is OR11 \(\ne\) OR10 * OR01

OR10 * OR01

## [1] 9.876543

2.1.2.0.5 Additive scales

Relative excess risk due to interaction (RERI)
- ranges from 0 to \(\infty\)
- 0 means no interaction
- +ve means positive interaction
Proportion of outcome among those with both exposures that is attributable to the interaction (AP)
- ranges from -1 to 1
- 0 means no interaction
Synergy index (SI)
- ratio of the combined effects and the individual effects
- ranges from 0 to \(\infty\)
- greater than 1 means positive interaction

RERI <- OR11 - OR10 - OR01 + 1
RERI

## [1] 3.739848

AP <- RERI / OR11
AP

## [1] 0.4138765

SI <- (OR11 - 1)/ (OR10 - 1 + OR01 - 1)
SI

## [1] 1.870482

2.1.2.0.6 Recoding to obtain different combinations

2.1.2.0.6.1 Base: alcohol drinker, no smoking for OR(smk1 on outcome [alc1==1]

Obs.Data <- OCdata
Obs.Data$smk <- as.factor(Obs.Data$smk)
Obs.Data$smk <- relevel(Obs.Data$smk, ref = "0")
Obs.Data$alc <- as.factor(Obs.Data$alc)
Obs.Data$alc <- relevel(Obs.Data$alc, ref = "1")
fit.or.int01 <- glm(oc ~ alc + smk + alc:smk, 
                    family = binomial(link = 'logit'), 
                    data = Obs.Data)
results.or.model01 <- summ(fit.or.int01, 
                           model.info = FALSE, 
                           model.fit = FALSE, 
                           exp = TRUE, 
                           confint = TRUE)
results.or.model01

	exp(Est.)	2.5%	97.5%	z val.	p
(Intercept)	0.50	0.19	1.33	-1.39	0.17
alc0	0.30	0.06	1.43	-1.51	0.13
smk1	2.71	1.00	7.37	1.95	0.05
alc0:smk1	1.09	0.18	6.48	0.10	0.92
Standard errors: MLE

2.1.2.0.6.2 Base: Smoker, no alcohol for OR(alc1 on outcome [smk1==1]

Obs.Data <- OCdata
Obs.Data$smk <- as.factor(Obs.Data$smk)
Obs.Data$smk <- relevel(Obs.Data$smk, ref = "1")
Obs.Data$alc <- as.factor(Obs.Data$alc)
Obs.Data$alc <- relevel(Obs.Data$alc, ref = "0")
fit.or.int10 <- glm(oc ~ alc + smk + alc:smk, 
                    family = binomial(link = 'logit'), 
                    data = Obs.Data)
results.or.model10 <- summ(fit.or.int10, 
                           model.info = FALSE, 
                           model.fit = FALSE, 
                           exp = TRUE,  
                           confint = TRUE)
results.or.model10

	exp(Est.)	2.5%	97.5%	z val.	p
(Intercept)	0.44	0.19	1.02	-1.91	0.06
alc1	3.05	1.29	7.18	2.55	0.01
smk0	0.34	0.08	1.47	-1.45	0.15
alc1:smk0	1.09	0.18	6.48	0.10	0.92
Standard errors: MLE

OR10b <- exp(sum(summary(fit.or.int10)$coef[c('alc1'),'Estimate']))
OR10b

## [1] 3.049699

2.1.2.1 Multiplicative scale

Multiplicative.scale <- exp(sum(summary(fit.or.int11)$coef[c('alc1:smk1'),'Estimate'])) 
Multiplicative.scale

## [1] 0.9149096

Multiplicative.scale = \(\exp(\alpha_{AB})\) = 0.91

# estimates of alcohol after changing smoking ref levels
OR10b / OR10

## [1] 0.9149096

2.1.3 Reporting guideline

Ref: Knol and VanderWeele (2012), Alli (2021)

int.object <- interactionR(fit.or.int11, 
                            exposure_names = c("alc1", "smk1"), 
                            ci.type = "mover", ci.level = 0.95, 
                            em=FALSE, recode = FALSE)
kable(int.object$dframe[,1:4], digits = 2)

Measures	Estimates	CI.ll	CI.ul
OR00	1.00	NA	NA
OR01	2.96	0.68	12.91
OR10	3.33	0.70	15.86
OR11	9.04	2.64	30.91
OR(smk1 on outcome [alc1==0]	2.96	0.68	12.91
OR(smk1 on outcome [alc1==1]	2.71	1.00	7.37
OR(alc1 on outcome [smk1==0]	3.33	0.70	15.86
OR(alc1 on outcome [smk1==1]	3.05	1.29	7.18
Multiplicative scale	0.91	0.15	5.42
RERI	3.74	-11.43	21.87
AP	0.41	-0.38	0.81
SI	1.87	0.65	5.42

2.1.4 Other available functions

Ref: Mathur and VanderWeele (2018)

## Mathur MB & VanderWeele TJ (2018)
source("Additive interactions function.R") 
additive_interactions(fit.or.int11)

## Loading required package: msm

##        Stat       Est       CI.lo     CI.hi    p.val.0 p.val.epi
## 1      RERI 3.7398483 -1.83667500 9.3163716 0.09435021 0.2704345
## 2        AP 0.4138765 -0.07306308 0.9008162 0.04786863        NA
## 3      alc1 0.2903548 -0.13508981 0.7157995 0.09050950        NA
## 4      smk1 0.2442668 -0.08791751 0.5764510 0.07475981        NA
## 5 alc1:smk1 0.4653784 -0.10296969 1.0337265 0.05426119        NA
##   p.val.suff.cause
## 1        0.1677822
## 2               NA
## 3               NA
## 4               NA
## 5               NA

2.1.5 Interaction definitions for OR

Summary of interaction definition for OR.
Groups and Conditions	Risk	Odds Ratio
Baseline (exposed to none)	\(R_{A = 0, B = 0} = P[Y^{a=0, b=0} \| L]\)	Reference
Exposed to \(A\) only	\(R_{A = 1, B = 0} = P[Y^{a=1, b=0} \| L]\)	\(OR_{A = 1} = [R_{A = 1, B = 0} / (1-R_{A = 1, B = 0})] / [R_{A = 0, B = 0} / (1 - R_{A = 0, B = 0})]\)
Exposed to \(B\) only	\(R_{A = 0, B = 1} = P[Y^{a=0, b=1} \| L]\)	\(OR_{B = 1} = [R_{A = 0, B = 1} / (1-R_{A = 0, B = 1})] / [R_{A = 0, B = 0} / (1 - R_{A = 0, B = 0})]\)
Exposed to both \(A\) and \(B\)	\(R_{A = 1, B = 1} = P[Y^{a=1, b=1} \| L]\)	\(OR_{A = 1, B = 1} = [R_{A = 1, B = 1} / (1-R_{A = 1, B = 1})] / [R_{A = 0, B = 0} / (1 - R_{A = 0, B = 0})]\)
Condition for interaction		\(OR_{A = 1, B = 1} \ne OR_{A = 1} \times OR_{B = 1}\)
Synergism		\(OR_{A = 1, B = 1} > OR_{A = 1} \times OR_{B = 1}\)
Antagonism		\(OR_{A = 1, B = 1} < OR_{A = 1} \times OR_{B = 1}\)

2.2 Modelling to obtain RR

risk ratio (RR), Multiplicative scale
- Very similar ideas as in OR

Obs.Data <- OCdata
Obs.Data$smk <- as.factor(Obs.Data$smk)
Obs.Data$smk <- relevel(Obs.Data$smk, ref = "0")
Obs.Data$alc <- as.factor(Obs.Data$alc)
Obs.Data$alc <- relevel(Obs.Data$alc, ref = "0")
fit.rr.int11 <- glm(oc ~ alc + smk + alc:smk, 
                    family = poisson(link = 'log'), 
                    data = Obs.Data)
results.rr.model11 <- summ(fit.rr.int11, 
                           model.info = FALSE, 
                           model.fit = FALSE, 
                           exp = TRUE, 
                           robust = "HC1", 
                           confint = TRUE)
results.rr.model11

	exp(Est.)	2.5%	97.5%	z val.	p
(Intercept)	0.13	0.04	0.40	-3.53	0.00
alc1	2.56	0.64	10.22	1.33	0.18
smk1	2.36	0.63	8.89	1.27	0.20
alc1:smk1	0.73	0.15	3.46	-0.39	0.69
Standard errors: Robust, type = HC1

Ref: Naimi and Whitcomb (2020)

Alternate formulation using log-binomial is possible (using model based SE, but confidence intervals somewhat differ from those reported from poisson(log)):

fit.rr.int11lb <- glm(oc ~ alc + smk + alc:smk, 
                      family = binomial(link = 'log'), 
                      data = Obs.Data)
results.rr.model11lb <- summ(fit.rr.int11lb, 
                             model.info = FALSE, 
                             model.fit = FALSE, 
                             exp = TRUE, 
                             confint = TRUE)
results.rr.model11lb

	exp(Est.)	2.5%	97.5%	z val.	p
(Intercept)	0.13	0.05	0.37	-3.78	0.00
alc1	2.56	0.74	8.84	1.48	0.14
smk1	2.36	0.71	7.85	1.40	0.16
alc1:smk1	0.73	0.19	2.88	-0.45	0.66
Standard errors: MLE

2.2.1 Estimates of different RRs

2.2.1.1 Alcohol

## RR00 = 1
RR10 <- exp(sum(summary(fit.rr.int11)$coef[c('alc1'),'Estimate']))
RR10 ## RR_A=1

## [1] 2.555555

\(\alpha_A\) = 0.94
RR_{A=1} = RR10 = \(\exp(\alpha_A)\) = 2.56

2.2.1.2 Smoking

RR01 <- exp(sum(summary(fit.rr.int11)$coef[c('smk1'),'Estimate']))
RR01 ## RR_B=1

## [1] 2.358974

\(\alpha_B\) = 0.86
RR_{B=1} = RR01 = \(\exp(\alpha_B)\) = 2.36

2.2.1.3 Both alcohol and smoking

RR11 <- exp(sum(summary(fit.rr.int11)$coef[c('alc1','smk1','alc1:smk1'),'Estimate'])) 
RR11 ##RR_A=1,B=1

## [1] 4.411764

\(\alpha_{AB}\) = -0.31
RR_{A=1,B=1} = RR11 = 4.41

2.2.1.4 Joint effect?

Is RR11 \(\ne\) RR10 * RR01

RR10 * RR01

## [1] 6.028489

2.2.1.5 Multiplicative scale

Multiplicative.scale <- exp(sum(summary(fit.rr.int11)$coef[c('alc1:smk1'),'Estimate'])) 
Multiplicative.scale

## [1] 0.7318192

Multiplicative.scale = \(\exp(\alpha_{AB})\) = 0.73

2.2.1.6 Additive scales

RERI <- RR11 - RR10 - RR01 + 1
RERI

## [1] 0.4972348

AP <- RERI / RR11
AP

## [1] 0.1127066

SI <- (RR11 - 1)/ (RR10 - 1 + RR01 - 1)
SI

## [1] 1.170606

2.2.1.7 Reporting guideline

int.object <- interactionR(fit.rr.int11, 
                            exposure_names = c("alc1", "smk1"), 
                            ci.type = "mover", ci.level = 0.95, 
                            em=FALSE, recode = FALSE)
kable(int.object$dframe[,1:4], digits = 2)

Measures	Estimates	CI.ll	CI.ul
OR00	1.00	NA	NA
OR01	2.36	0.63	8.89
OR10	2.56	0.64	10.22
OR11	4.41	1.41	13.78
OR(smk1 on outcome [alc1==0]	2.36	0.63	8.89
OR(smk1 on outcome [alc1==1]	1.73	0.77	3.88
OR(alc1 on outcome [smk1==0]	2.56	0.64	10.22
OR(alc1 on outcome [smk1==1]	1.87	0.92	3.79
Multiplicative scale	0.73	0.15	3.46
RERI	0.50	-9.97	7.01
AP	0.11	-0.82	0.75
SI	1.17	0.42	3.23

Ref: Alli (2021)

2.2.1.8 Interaction definitions for RR

Summary of interaction definition for RR.
Groups and Conditions	Risk	Risk Ratio
Baseline (exposed to none)	\(R_{A = 0, B = 0} = P[Y^{a=0, b=0} \| L]\)	Reference
Exposed to \(A\) only	\(R_{A = 1, B = 0} = P[Y^{a=1, b=0} \| L]\)	\(RR_{A = 1} = R_{A = 1, B = 0} / R_{A = 0, B = 0}\)
Exposed to \(B\) only	\(R_{A = 0, B = 1} = P[Y^{a=0, b=1} \| L]\)	\(RR_{B = 1} = R_{A = 0, B = 1} / R_{A = 0, B = 0}\)
Exposed to both \(A\) and \(B\)	\(R_{A = 1, B = 1} = P[Y^{a=1, b=1} \| L]\)	\(RR_{A = 1, B = 1} = R_{A = 1, B = 1} / R_{A = 0, B = 0}\)
Condition for interaction		\(RR_{A = 1, B = 1} \ne RR_{A = 1} \times RR_{B = 1}\)
Synergism		\(RR_{A = 1, B = 1} > RR_{A = 1} \times RR_{B = 1}\)
Antagonism		\(RR_{A = 1, B = 1} < RR_{A = 1} \times RR_{B = 1}\)

2.3 Modelling to obtain RD

risk difference (RD), Additive scale
- risk(Y) = \(\alpha_0\) + \(\alpha_A\) A + \(\alpha_B\) B + \(\alpha_{AB}\) (A \(\times\) B)
  - RD10 \(\times\) RD01 \(\ne\) RD11
  - \(\alpha_A\) + \(\alpha_B\) \(\ne\) (\(\alpha_A\) + \(\alpha_B\) + \(\alpha_{AB}\))
- synergism / super-additive
  - RD11 > RD10 + RD01
- antagonism / sub-additive
  - RD11 < RD10 + RD01

Obs.Data <- OCdata
Obs.Data$smk <- as.factor(Obs.Data$smk)
Obs.Data$smk <- relevel(Obs.Data$smk, ref = "0")
Obs.Data$alc <- as.factor(Obs.Data$alc)
Obs.Data$alc <- relevel(Obs.Data$alc, ref = "0")
fit.rd.int11 <- glm(oc ~ alc + smk + alc:smk, 
                    family = gaussian(link = 'identity'), 
                    data = Obs.Data)
results.rd.model11 <- summ(fit.rd.int11, 
                           model.info = FALSE, 
                           model.fit = FALSE, 
                           exp = FALSE, 
                           robust = "HC1", 
                           confint = TRUE)
results.rd.model11

	Est.	2.5%	97.5%	t val.	p
(Intercept)	0.13	-0.07	0.33	1.28	0.20
alc1	0.20	-0.10	0.50	1.32	0.19
smk1	0.18	-0.10	0.45	1.27	0.20
alc1:smk1	0.06	-0.29	0.42	0.36	0.72
Standard errors: Robust, type = HC1

Ref: Naimi and Whitcomb (2020)

Alternate formulation using binomial(identity) is possible (using model based SE, but confidence intervals slightly differ from those reported from gaussian(identity)):

fit.rd.int11bi <- glm(oc ~ alc + smk + alc:smk, 
                      family = binomial(link = 'identity'), 
                      data = Obs.Data)
results.rd.model11bi <- summ(fit.rd.int11bi, 
                             model.info = FALSE, 
                             model.fit = FALSE, 
                             exp = FALSE, 
                             confint = TRUE)
results.rd.model11bi

	Est.	2.5%	97.5%	z val.	p
(Intercept)	0.13	-0.01	0.27	1.86	0.06
alc1	0.20	-0.05	0.46	1.54	0.12
smk1	0.18	-0.05	0.40	1.55	0.12
alc1:smk1	0.06	-0.25	0.38	0.40	0.69
Standard errors: MLE

2.3.1 Base: Smoker, no alcohol for RD(alc1 on outcome [smk1==1]

Alternate combination of smoking:

Obs.Data <- OCdata
Obs.Data$smk <- as.factor(Obs.Data$smk)
Obs.Data$smk <- relevel(Obs.Data$smk, ref = "1")
Obs.Data$alc <- as.factor(Obs.Data$alc)
Obs.Data$alc <- relevel(Obs.Data$alc, ref = "0")
fit.rd.int01 <- glm(oc ~ alc + smk + alc:smk, 
                    family = gaussian(link = 'identity'), 
                    data = Obs.Data)
results.rd.model01 <- summ(fit.rd.int01, 
                           model.info = FALSE, 
                           model.fit = FALSE, 
                           exp = FALSE, 
                           robust = "HC1", 
                           confint = TRUE)
results.rd.model01

	Est.	2.5%	97.5%	t val.	p
(Intercept)	0.31	0.12	0.49	3.22	0.00
alc1	0.27	0.07	0.46	2.71	0.01
smk0	-0.18	-0.45	0.10	-1.27	0.20
alc1:smk0	-0.06	-0.42	0.29	-0.36	0.72
Standard errors: Robust, type = HC1

2.3.2 Estimates of different RDs

2.3.2.1 Alcohol

## RD00 = 0
RD10 <- sum(summary(fit.rd.int11)$coef[c('alc1'),'Estimate'])
RD10 ## RD_A=1

## [1] 0.2028986

RD_{A=1} = RD10 = \(\alpha_A\) = 0.2

From alternate combination of smoking

RD10b <- sum(summary(fit.rd.int01)$coef[c('alc1'),'Estimate'])
RD10b ## RD_A=1

## [1] 0.2677553

2.3.2.2 Smoking

RD01 <- sum(summary(fit.rd.int11)$coef[c('smk1'),'Estimate'])
RD01 ## RD_B=1

## [1] 0.1772575

RD_{B=1} = RD01 = \(\alpha_B\) = 0.18

2.3.2.3 Both alcohol and smoking

RD11 <- sum(summary(fit.rd.int11)$coef[c('alc1','smk1','alc1:smk1'),'Estimate'])
RD11 ##RD_A=1,B=1

## [1] 0.4450128

RD_{A=1,B=1} = RD11 = \(\alpha_{AB}\) = 0.06

2.3.2.4 Joint effect?

Is RD11 \(\ne\) RD10 + RD01

RD10 + RD01

## [1] 0.3801561

2.3.2.5 Additive scale (interaction contrast)

IC <- RD10b - RD10
IC

## [1] 0.06485671

Current version of interactionR do not support RD, although the documentation says so: Alli (2021)

2.3.3 Interaction definitions for RD

Summary of interaction definition for RD.
Groups and Conditions	Risk	Risk Difference
Baseline (exposed to none)	\(R_{A = 0, B = 0} = P[Y^{a=0, b=0} \| L]\)	Reference
Exposed to \(A\) only	\(R_{A = 1, B = 0} = P[Y^{a=1, b=0} \| L]\)	\(RD_{A = 1} = R_{A = 1, B = 0} - R_{A = 0, B = 0}\)
Exposed to \(B\) only	\(R_{A = 0, B = 1} = P[Y^{a=0, b=1} \| L]\)	\(RD_{B = 1} = R_{A = 0, B = 1} - R_{A = 0, B = 0}\)
Exposed to both \(A\) and \(B\)	\(R_{A = 1, B = 1} = P[Y^{a=1, b=1} \| L]\)	\(RD_{A = 1, B = 1} = R_{A = 1, B = 1} - R_{A = 0, B = 0}\)
Condition for interaction		\(RD_{A = 1, B = 1} \ne RD_{A = 1} + RD_{B = 1}\)
Synergism		\(RD_{A = 1, B = 1} > RD_{A = 1} + RD_{B = 1}\)
Antagonism		\(RD_{A = 1, B = 1} < RD_{A = 1} + RD_{B = 1}\)

2.4 Reporting guidelines

How to determine whether interaction should be reported in additive or multiplicative scale?
- Not via empirical assessment (analysis of data at hand)
- Study objectives should be clear
- data generating mechanism should be considered carefully
What if results differ?
- Very possible to have additive interaction present, but not at a multiplicative scale and vice versa
- When conclusions differ (interaction is present in additive, but absent in multiplicative or vice versa), both should be reported.
- When both are concluding the same (interaction is present in both scales or absent), questions related to magnitude remains

Reference

Alli, Babatunde Y. 2021. “InteractionR: An r Package for Full Reporting of Effect Modification and Interaction.” Software Impacts, 100147.

Knol, Mirjam J, and Tyler J VanderWeele. 2012. “Recommendations for Presenting Analyses of Effect Modification and Interaction.” International Journal of Epidemiology 41 (2): 514–20.

Mathur, Maya B, and Tyler J VanderWeele. 2018. “R Function for Additive Interaction Measures.” Epidemiology (Cambridge, Mass.) 29 (1): e5.

Naimi, Ashley I, and Brian W Whitcomb. 2020. “Estimating Risk Ratios and Risk Differences Using Regression.” American Journal of Epidemiology 189 (6): 508–10.