---
title: "Getting Started with the `multimediate` package"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Getting Started with the `multimediate` package}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---
## Multiple mediation analysis 
```{r, include = FALSE}
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)
```

```{r setup}
library(multimediate)
```
This function is based on the work presented in Jerolon et al 2020 and the R package mediation presented in Tingley et al 2014. Models supported by our R function are shown in the following table.

<table border="4" width="100%">
  <thead>
    <tr>
      <th colspan="5">Models supported by multimediate.</th>
    </tr>
  </thead>
  <tr>
    <td>&nbsp;</td>
    <td colspan="5" bgcolor="#333AFF">Outcome</td>
  </tr>
  <tr>
    <td bgcolor="#FF5733">Mediators </td><td bgcolor="#337AFF">Linear</td><td bgcolor="#337AFF">Binary</td><td bgcolor="#337AFF">Ordered categorical</td><td bgcolor="#337AFF">Non-ordered categorical</td>
  </tr>
  <tr>
    <td bgcolor="#FFA533">Linear </td><td>Ok</td><td>Ok</td><td>Ok</td><td>Not yet</td>
  </tr>
  <tr>
    <td bgcolor="#FFA533">Binary (Probit)</td><td>Ok</td><td>Ok</td><td>Ok</td><td>Not yet</td>
  </tr>
  <tr>
    <td bgcolor="#FFA533">Ordered categorical</td><td>Ok</td><td>Ok</td><td>Ok</td><td>Not yet</td>
  </tr>
  <tr>
    <td bgcolor="#FFA533">Non-ordered categorical</td><td>Not yet</td><td>Not yet</td><td>Not yet</td><td>Not yet</td>
  </tr>
</table>

 
## Example 1

<img src="data1.png">

We will do a mediation analysis with the following model.

- Binary treatment : $T$
- Linear Outcome : $Y$
- Linear mediators : $M_1,M_2,M_3$ 
- Binary covariables : $C_1,C_2,C_3$ 

In this first example the covariable $C_1$ is affecting all three mediators and the outcome $Y$; the covariable $C_2$ is affecting $M_1$ and the outcome $Y$; and the covariable $C_3$ is affecting $M_2$ and the outcome $Y$ as described in the diagramm.



```{r DATA1}
data(data1)
data1$Treatment=as.factor(data1$Treatment)
data1$C1=as.factor(data1$C1)
data1$C2=as.factor(data1$C2)
data1$C3=as.factor(data1$C3)
data1$M1=as.numeric(data1$M1)
data1$M2=as.numeric(data1$M2)
data1$M3=as.numeric(data1$M3)
data1$Outcome=as.numeric(data1$Outcome)
summary(data1)
```

The regression models are then as follows. 

```{r REG1}
M1reg=lm(M1~ Treatment + C1 + C2, data=data1)
M2reg=lm(M2~ Treatment + C1 + C3, data=data1)
M3reg=lm(M3~ Treatment + C1     , data=data1)

Yreg=lm(Outcome~ Treatment + M1 + M2 + M3 + C1 + C2 + C3, data=data1)
```

Once the regressions are done, we can then proceed to the multiple mediation analysis using the multimediate function. Then display the results with the summary function.

```{r MED1eval, eval=FALSE}
med.analysis=multimediate(lmodel.m=list(M1reg,M2reg,M3reg),correlated=TRUE,model.y=Yreg,
                        treat="Treatment",treat.value=1,control.value=0,J=1000,conf.level=0.95)
```

```{r MED1include, include=FALSE}
med.analysis=multimediate(lmodel.m=list(M1reg,M2reg,M3reg),correlated=TRUE,model.y=Yreg,
                        treat="Treatment",treat.value=1,control.value=0,J=1000,conf.level=0.95)
```



<table>
  <thead>
    <tr>
      <th colspan="6">True values of causal effects in data 1.</th>
    </tr>
  </thead>
    <tr>
      <td>Total Effect</td><td>Direct Effect</td><td>Joint Indirect Effect</td><td>Indirect Effect by M1</td><td>Indirect Effect by M2</td><td>Indirect Effect by M3</td>
    </tr>
    <tr>
      <td>92</td><td>10</td><td>82</td><td>10</td><td>18</td><td>54</td>
    </tr>
</table>

```{r}
summary(med.analysis,opt="avg")

#plot(med.analysis)
```

## Example 2 
<img src="data2.png">

- Binary treatment : $T$
- Binary Outcome : $Y$
- Linear mediators : $M_1,M_2,M_3$ 
- Binary covariables : $C_1,C_2,C_3$ 



```{r DATA2}
data(data2)
data2$Treatment=as.factor(data2$Treatment)
data2$C1=as.factor(data2$C1)
data2$C2=as.factor(data2$C2)
data2$C3=as.factor(data2$C3)
data2$M1=as.numeric(data2$M1)
data2$M2=as.numeric(data2$M2)
data2$M3=as.numeric(data2$M3)
data2$Outcome=as.factor(data2$Outcome)
summary(data2)
```


```{r REG2}
M1reg=lm(M1~ Treatment + C1, data=data2)
M2reg=lm(M2~ Treatment + C2, data=data2)
M3reg=lm(M3~ Treatment + C3, data=data2)

Yreg=glm(Outcome~ Treatment + M1 + M2 + M3 + C1 + C2 + C3, data=data2, family = binomial("logit"))
```



```{r MED2eval, eval=FALSE}
med.analysis=multimediate(lmodel.m=list(M1reg,M2reg,M3reg),correlated=TRUE,model.y=Yreg,
                        treat="Treatment",treat.value=1,control.value=0,J=1000,conf.level=0.95)

```

```{r MED2include, include=FALSE}
med.analysis=multimediate(lmodel.m=list(M1reg,M2reg,M3reg),correlated=TRUE,model.y=Yreg,
                        treat="Treatment",treat.value=1,control.value=0,J=1000,conf.level=0.95,fun=median)

```


<table>
  <thead>
    <tr>
      <th colspan="6">True values of causal effects in data 2 on log OR scale.</th>
    </tr>
  </thead>
    <tr>
      <td>Total Effect</td><td>Direct Effect</td><td>Joint Indirect Effect</td><td>Indirect Effect by M1</td><td>Indirect Effect by M2</td><td>Indirect Effect by M3</td>
    </tr>
    <tr>
      <td>1.311</td><td>0.71</td><td>0.601</td><td>0.195</td><td>0.19</td><td>0.216</td>
    </tr>
</table>

```{r MED2res}
 summary(med.analysis,opt="avg",logit="all")

# summary(med.analysis,opt="avg",logit="effects")
# plot(med.analysis)


# summary(med.analysis,opt="avg",logit="OR")
#plot(med.analysis,logit = "OR")

#summary(med.analysis,opt="avg",logit="logOR")
#plot(med.analysis,logit = "logOR")

```

## Example 3

<img src="data3.png">

- Binary treatment $T$
- Polytomic Outcome $Y$
- Linear mediators $M_1,M_2,M_3$ 
- Binary covariables $C_1,C_2,C_3$ 


```{r DATA3}
data(data3)
data3$Treatment=as.factor(data3$Treatment)
data3$C1=as.factor(data3$C1)
data3$C2=as.factor(data3$C2)
data3$C3=as.factor(data3$C3)
data3$M1=as.numeric(data3$M1)
data3$M2=as.numeric(data3$M2)
data3$M3=as.numeric(data3$M3)
data3$Outcome=as.factor(data3$Outcome)

summary(data3)
```

```{r REG3}
M1reg=lm(M1~ Treatment + C1 + C3, data=data3)
M2reg=lm(M2~ Treatment + C1 + C2, data=data3)
M3reg=lm(M3~ Treatment + C2 + C3, data=data3)

library(MASS)
Yreg=polr(Outcome ~ Treatment + M1 + M2 + M3 + C1 + C2 + C3 , data = data3, method = "probit")

```



```{r MED3,include=FALSE}
med.analysis=multimediate(lmodel.m=list(M1reg,M2reg,M3reg),correlated=TRUE,model.y=Yreg,
                        treat="Treatment",treat.value=1,control.value=0,J=100,conf.level=0.95)
```

```{r MED3eval,eval=FALSE}
med.analysis=multimediate(lmodel.m=list(M1reg,M2reg,M3reg),correlated=TRUE,model.y=Yreg,
                        treat="Treatment",treat.value=1,control.value=0,J=100,conf.level=0.95)

summary(med.analysis,opt="avg")
```

<table>
  <thead>
    <tr>
      <th colspan="6">Values of causal effects in data 3.</th>
    </tr>
  </thead>
    <tr>
      <td>Total Effect</td><td>Direct Effect</td><td>Joint Indirect Effect</td><td>Indirect Effect by M1</td><td>Indirect Effect by M2</td><td>Indirect Effect by M3</td>
    </tr>
    <tr>
      <td>2.21</td><td>0.5</td><td>1.71</td><td>0.32</td><td>0.44</td><td>0.60</td>
    </tr>
</table>

```{r MED3res}
summary(med.analysis,opt="avg")
```

## Example 4
<img src="data4.png">

- Binary treatment $T$
- Linear Outcome $Y$
- Linear mediator $M_1$
- Binary mediator $,M_2$
- Polytomic mediator $M_3$ 
- Binary covariables $C_1,C_2,C_3$ 




```{r DATA4}
data(data4)
data4$Treatment=as.factor(data4$Treatment)
data4$C1=as.factor(data4$C1)
data4$C2=as.factor(data4$C2)
data4$C3=as.factor(data4$C3)
data4$M1=as.numeric(data4$M1)
data4$M3=as.factor(data4$M3)
data4$M2=as.factor(data4$M2)
data4$Outcome=as.numeric(data4$Outcome)
summary(data4)
```

```{r REG4}
M1reg=lm(M1~  Treatment + C1 + C2 + C3, data = data4)
M2reg=glm(M2~ Treatment + C1 + C3, data = data4, family = binomial("probit"))
M3reg=polr(M3~Treatment + C2 + C3     , data = data4, method = "probit")

Yreg=lm(Outcome~ Treatment + M1 + M2 + M3 + C1 + C2 + C3, data=data4)
```



```{r MED4include, include=FALSE}
med.analysis=multimediate(lmodel.m=list(M1reg,M2reg,M3reg),correlated=TRUE,model.y=Yreg,
                        treat="Treatment",treat.value=1,control.value=0,J=1000,conf.level=0.95)


```

```{r MED4eval, eval=FALSE}
med.analysis=multimediate(lmodel.m=list(M1reg,M2reg,M3reg),correlated=TRUE,model.y=Yreg,
                        treat="Treatment",treat.value=1,control.value=0,J=1000,conf.level=0.95)
```

<table>
  <thead>
    <tr>
      <th colspan="6">Values of causal effects in data 4.</th>
    </tr>
  </thead>
    <tr>
      <td>Total Effect</td><td>Direct Effect</td><td>Joint Indirect Effect</td><td>Indirect Effect by M1</td><td>Indirect Effect by M2</td><td>Indirect Effect by M3</td>
    </tr>
    <tr>
      <td>57.68</td><td>10</td><td>47.48</td><td>34.2</td><td>0</td><td>13.48</td>
    </tr>
</table>

```{r MED4res}
summary(med.analysis,opt="avg")
```


## Example 5


- Binary treatment $T$
- Survival Outcome $Y$
- Linear mediators $M_1,M_2,M_3$


```{r DATA5}
data(data5)
# data5$Exposure=as.factor(data5$Exposure)
# data5$M1=as.numeric(data5$M1)
# data5$M3=as.numeric(data5$M3)
# data5$M2=as.numeric(data5$M2)
# data5$Outcome=as.numeric(data5$Outcome)
summary(data5)
```

```{r REG5}
modM1=lm(M1 ~ Exposure, data = data5)    
modM2=lm(M2 ~ Exposure, data = data5)
modM3=lm(M3 ~ Exposure, data = data5)
lmodel.m=list(modM1,modM2, modM3)

library(timereg)
model.y=aalen(Surv(surv_time, event) ~ const(Exposure) + const(M1) + const(M2) + const(M3), data = data5, robust=T)
```



```{r MED5include, include=FALSE}
multi.media=multimediate(lmodel.m,correlated=TRUE,model.y,treat='Exposure',treat.value=1,control.value=0,J=1000,conf.level=0.95,data=data5)
```

```{r MED5eval, eval=FALSE}
multi.media=multimediate(lmodel.m,correlated=TRUE,model.y,treat='Exposure',treat.value=1,control.value=0,J=1000,conf.level=0.95,data=data5)
```


```{r MED5res}
summary(multi.media, opt='avg')
```