02-effmod.Rmd

```{r chunkSetup2, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
require(interactionR)
require(mosaic)
require(epiR)
require(jtools)
require(Publish)
```


# Effect modification

- Multiplicative odds ratio scale are very popular for examining effect modification. 
- However, often multiplicative scale is not helpful for public health related interpretations, and is considered as inadequate for such assessment. 
- Additive probability scale should be additionally used for such assessment.

We will discuss the assessment of effect modification using one example data, and show how we can calculate the measures using R or calculate them by hand. 

## Estimates from Tabulation

```{r}
require(interactionR)
data(OCdata)
outcome = "oc"
ex = "alc"
dataset = OCdata
```

Now we will be working with stratified data (stratification by levels of effect modifier)

### First strata: Smoking = 1

```{r}
str = "smk"
str.lab = 1
dataset1 <- dataset[dataset[[str]] == str.lab,]
M <- table(dataset1[[ex]], dataset1[[outcome]])
M3 <- matrix(c(M[2,2],M[2,1],M[1,2],M[1,1]), 
             nrow = 2, 
             byrow = TRUE)
res1 <- epi.2by2(dat = M3, 
                 method = "cross.sectional",
         conf.level = 0.95, 
         units = 1, 
         interpret = FALSE, 
         outcome = "as.columns")
res1
```


Check the results yourself
```{r}
p1.s1 <- as.numeric(strsplit(as.character(res1$tab$`              Prevalence *`), " ")[[1]][1])
p0.s1 <- as.numeric(strsplit(as.character(res1$tab$`              Prevalence *`), " ")[[2]][1])
RR.s1 <- p1.s1/p0.s1 
RR.s1
RD.s1 <- p1.s1 - p0.s1 
RD.s1
OR.s1a <- (p1.s1/(1-p1.s1))
OR.s1b <- (p0.s1/(1-p0.s1))
OR.s1 <- OR.s1a/OR.s1b
OR.s1
```

### second strata: : Smoking = 0

```{r}
str.lab = 0
dataset1 <- dataset[dataset[[str]] == str.lab,]
M <- table(dataset1[[ex]], dataset1[[outcome]])
M3 <- matrix(c(M[2,2],M[2,1],M[1,2],M[1,1]), 
             nrow = 2, byrow = TRUE)
res0 <- epi.2by2(dat = M3, 
                 method = "cross.sectional",
         conf.level = 0.95, 
         units = 1, 
         interpret = FALSE, 
         outcome = "as.columns")
res0
```

Check the results yourself
```{r}
p1.s0 <- as.numeric(strsplit(as.character(res0$tab$`              Prevalence *`), " ")[[1]][1])
p0.s0 <- as.numeric(strsplit(as.character(res0$tab$`              Prevalence *`), " ")[[2]][1])
RR.s0 <- p1.s0/p0.s0 
RR.s0
RD.s0 <- p1.s0 - p0.s0 
RD.s0
OR.s0a <- (p1.s0/(1-p1.s0))
OR.s0b <- (p0.s0/(1-p0.s0))
OR.s0 <- OR.s0a/OR.s0b
OR.s0
```


### Does effect modification exist?

OR: `r OR.s1` vs. `r OR.s0`

Ratio of the OR in two groups (null value = 1)
```{r}
OR.s1 / OR.s0
```

RD: `r RD.s1` vs. `r RD.s0`

Difference of the RD in two groups (null value = 0). This is also known as interaction contrast (IC)
```{r}
RD.s1 - RD.s0
```

RR: `r RR.s1` vs. `r RR.s0`

Ratio of the RR in two groups (null value = 1)
```{r}
RR.s1 / RR.s0
```

Are they equivalent in both stratified groups?


#### Are the lines parallel?

We can utilize the probability and log-odds calculated before to plot these lines:

```{r}
png('images/pem.png')
par(mar=c(5.1, 4.1, 4.1, 8.1), xpd=TRUE)
plot(c(1,0),c(p1.s1,p0.s1), 
     ylim = range(c(p1.s1,p1.s0, p0.s1,p0.s0)), 
     xlim=c(0,1), 
     xlab = "exposure", 
     ylab = "Probability of outcome", 
     lty = 1, 
     col = "green", 
     type = "l", 
     main = "Probability scale")
lines(c(1,0),c(p1.s0,p0.s0), 
      lty = 2, 
      col = "red")
legend("bottomright", 
       inset=c(-0.2,0), 
       legend=c("strata 1","strata 0"), 
       pch=c(1,3), 
       title="effect modifier", 
       col = c("green", "red"))
dev.off()
```

```{r dagemp, echo=FALSE, out.width = '50%'}
knitr::include_graphics("images/pem.png")
```


```{r}
png('images/loem.png')
par(mar=c(5.1, 4.1, 4.1, 8.1), xpd=TRUE)
plot(c(1,0),c(log(OR.s1a), log(OR.s1b)), 
     ylim = range(c(log(OR.s1a), log(OR.s0a),log(OR.s1b),log(OR.s0b))), 
     xlim=c(0,1), xlab = "exposure", 
     ylab = "log-odds of outcome", 
     lty = 1, 
     col = "green", 
     type = "l", 
     main = "Log-odds scale")
lines(c(1,0),c(log(OR.s0a),log(OR.s0b)), 
      lty = 2, 
      col = "red")
legend("bottomright", inset=c(-0.2,0), 
       legend=c("strata 1","strata 0"), 
       pch=c(1,3), 
       title="effect modifier", 
       col = c("green", "red"))
dev.off()
```

```{r dagemlo, echo=FALSE, out.width = '50%'}
knitr::include_graphics("images/loem.png")
```



## Estimates from Regression

- Add interaction term in the regression
- check the coef of the interaction terms, and compare that number with the numbers calculated at the previous section (rations and differences): they should be equivalent if not equal

Using regression could help accommodate the following

- covariates could be added
- exposure variable could be continuous covariate 


### OR: multiplicative scale

```{r}
OCdata$alc <- as.factor(OCdata$alc)
OCdata$smk <- as.factor(OCdata$smk)
reg1 = glm(oc~alc*smk, 
           family = binomial(link = 'logit'), 
           data = OCdata)
require(jtools)
reg.res1 <- summ(reg1, exp = TRUE)
reg.res1$coeftable
```

P-values are used to determine significance of the interaction term.

Let us estimate the ORs by effect modifier

```{r}
require(Publish)
saveres <- publish(reg1)
```

We can also plot this to check the overlap of CIs visually

```{r}
png('images/em.png')
plot(c(0,1), saveres$rawTable[3:4,"OddsRatio"], 
     xlab="effect modifier", ylab="OR", 
     main="OR by effect modifier", xlim = c(0,1), 
     ylim = range(saveres$rawTable[3:4,4:5]))
lines(c(0,0), saveres$rawTable[3,c("Lower","Upper")], 
      col="blue", lty=2)
lines(c(1,1), saveres$rawTable[4,c("Lower","Upper")], 
      col="blue", lty=2)
dev.off()
```

```{r dagemx, echo=FALSE, out.width = '50%'}
knitr::include_graphics("images/em.png")
```


<!---
### RD: additive scale

```{r}
reg2 = glm(oc~alc*smk, 
           family = binomial(link = 'identity'), 
           data = OCdata)
reg.res2 <- summ(reg2, 
                 exp = FALSE, 
                 confint = TRUE)
reg.res2$coeftable
```

Convergence is often an issue, and hence the following can be used instead. But for the following, robust standard error estimation is necessary. One disadvantage is that the estimated probabilities could be beyond 0 and 1.

```{r}
reg3 = glm(oc~alc*smk, 
           family = gaussian, 
           data = OCdata)
reg.res3 <- summ(reg3, 
                 exp = FALSE, 
                 robust = TRUE, 
                 confint = TRUE)
reg.res3$coeftable
```

### RR: multiplicative scale

```{r}
reg4 = glm(oc~alc*smk, 
           family = binomial(link = 'log'), 
           data = OCdata)
reg.res4 <- summ(reg4, 
                 exp = TRUE, 
                 robust = TRUE)
reg.res4$coeftable
```

or, alternatively, the following could be used.

```{r}
reg5 = glm(oc~alc*smk, 
           family = poisson(link = 'log'), 
           data = OCdata)
reg.res5 <- summ(reg5, 
                 exp = TRUE, 
                 robust = TRUE)
reg.res5$coeftable
```



#### Prevalence estimation

```{r, warning=FALSE, message=FALSE}
## library(devtools)
## install_github("Raydonal/prLogistic")
## require(prLogistic)
## OCdata$alc <- as.numeric(as.character(OCdata$alc))
## OCdata$smk <- as.numeric(as.character(OCdata$smk))
## OCdata$oc <- as.numeric(as.character(OCdata$oc))
## reg1x = glm(oc~alc*smk, family = binomial(link = 'logit'), data = OCdata)
##pr10 <- prLogisticBootCond(reg1, data = OCdata)
##pr10
## d1 <- subset(OCdata, smk == 1)
## d0 <- subset(OCdata, smk == 0)
## d1$alc <- as.numeric(as.character(d1$alc))
## d0$alc <- as.numeric(as.character(d0$alc))
## reg1s1 = glm(oc~alc, family = binomial(link = 'logit'), data = d1)
## reg1s0 = glm(oc~alc, family = binomial(link = 'logit'), data = d0)
##  prLogisticBootMarg(reg1s1, data = d1)
##  prLogisticBootMarg(reg1s0, data = d0)
```
--->

### Converted estimates from OR

```{r, warning=FALSE, message=FALSE}
require(interactionR)
reg1 = glm(oc~alc*smk, 
           family = binomial(link = 'logit'), 
           data = OCdata)
table_object = interactionR(reg1, 
                            exposure_names = c("smk1", "alc1"), 
                            ci.type = "mover", 
                            ci.level = 0.95, 
                            em=TRUE, 
                            recode = FALSE)
table_object$dframe
interactionR_table(table_object)
```

Compare with previously obtained results, and also calculate by hand

```{r}
beta1 <- as.numeric(reg1$coefficients[2])
OR10 <- exp(beta1)
OR10
beta2 <- as.numeric(reg1$coefficients[3])
OR01 <- exp(beta2)
OR01
beta3 <- as.numeric(reg1$coefficients[4])
OR11 <- exp(beta1 + beta2 + beta3)
OR11
```


```{r}
RERI <- OR11 - OR10 - OR01 + 1
RERI
```

<!---
```{r}
AP <- RERI / OR11
AP
SI <- (OR11 - 1)/ (OR10 - 1 + OR01 - 1)
SI
```
--->


### RR and RD

Estimation approach is same as interaction (adding an interaction term in the corresponding regression), we are just interested in estimates from each level of the effect modifier (hence less estimates are needed compared to those estimated in interaction assessment).