Getting Started with the multimediate package

Multiple mediation analysis

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.

Models supported by multimediate.
  Outcome
Mediators Linear Binary Ordered categorical Non-ordered categorical
Linear Ok Ok Ok Not yet
Binary (Probit) Ok Ok Ok Not yet
Ordered categorical Ok Ok Ok Not yet
Non-ordered categorical Not yet Not yet Not yet Not yet

Example 1

We will do a mediation analysis with the following model.

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.

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)
#>  Treatment C1      C2      C3            M1               M2         
#>  0:480     0:501   0:417   0:596   Min.   :-1.410   Min.   :-0.8455  
#>  1:520     1:499   1:583   1:404   1st Qu.: 2.850   1st Qu.: 4.7567  
#>                                    Median : 5.809   Median : 7.2685  
#>                                    Mean   : 5.431   Mean   : 7.2333  
#>                                    3rd Qu.: 7.921   3rd Qu.: 9.7710  
#>                                    Max.   :12.054   Max.   :15.6520  
#>        M3             Outcome      
#>  Min.   :-0.8298   Min.   : 12.88  
#>  1st Qu.: 3.5296   1st Qu.:100.61  
#>  Median : 6.8410   Median :145.90  
#>  Mean   : 6.6268   Mean   :145.54  
#>  3rd Qu.: 9.7148   3rd Qu.:190.45  
#>  Max.   :13.5039   Max.   :295.79

The regression models are then as follows.

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.

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)
True values of causal effects in data 1.
Total Effect Direct Effect Joint Indirect Effect Indirect Effect by M1 Indirect Effect by M2 Indirect Effect by M3
92 10 82 10 18 54 
summary(med.analysis,opt="avg")
#>               . Estimation  IC.inf  IC.sup P.val
#> 1    ACME.joint    81.5825 79.8114 83.2367     0
#> 2      PM.joint     0.8861  0.8780  0.8946     0
#> 3       ACME.M1    10.4082  9.6572 11.1916     0
#> 4         PM.M1     0.1131  0.1047  0.1218     0
#> 5       ACME.M2    17.5295 16.6181 18.4545     0
#> 6         PM.M2     0.1904  0.1802  0.2012     0
#> 7       ACME.M3    53.6448 52.2872 54.9504     0
#> 8         PM.M3     0.5827  0.5662  0.6011     0
#> 9           ADE    10.4811  9.7115 11.2080     0
#> 10 Total Effect    92.0637 90.5538 93.5323     0

#plot(med.analysis)

Example 2

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)
#>  Treatment C1      C2      C3            M1                M2          
#>  0:514     0:500   0:430   0:610   Min.   :-2.7798   Min.   :-2.19561  
#>  1:486     1:500   1:570   1:390   1st Qu.:-0.1127   1st Qu.:-0.04618  
#>                                    Median : 0.5923   Median : 0.61047  
#>                                    Mean   : 0.5699   Mean   : 0.60425  
#>                                    3rd Qu.: 1.2644   3rd Qu.: 1.23342  
#>                                    Max.   : 3.8948   Max.   : 3.75137  
#>        M3           Outcome   
#>  Min.   :-3.3517   FALSE:757  
#>  1st Qu.:-0.2396   TRUE :243  
#>  Median : 0.4367              
#>  Mean   : 0.4572              
#>  3rd Qu.: 1.1357              
#>  Max.   : 4.0804
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"))
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)
True values of causal effects in data 2 on log OR scale.
Total Effect Direct Effect Joint Indirect Effect Indirect Effect by M1 Indirect Effect by M2 Indirect Effect by M3
1.311 0.71 0.601 0.195 0.19 0.216 
 summary(med.analysis,opt="avg",logit="all")
#> Warning in summary.mm(med.analysis, opt = "avg", logit = "all"): The
#> proportions mediated on the OR scale can be considered if the outcome is rare,
#> otherwise the proportions mediated on effects scale and/or logOR scale have to
#> be considered.
#>               . Estimation IC.inf IC.sup P.val               . Estimation
#> 1    ACME.joint     0.0765 0.0545 0.0990 0.000   OR ACME.joint     1.9407
#> 2      PM.joint     0.4485 0.3288 0.5824 0.000     OR PM.joint     0.6313
#> 3       ACME.M1     0.0245 0.0000 0.0540 0.046      OR ACME.M1     1.2448
#> 4         PM.M1     0.1436 0.0000 0.3611 0.046        OR PM.M1     0.2543
#> 5       ACME.M2     0.0205 0.0035 0.0480 0.014      OR ACME.M2     1.1990
#> 6         PM.M2     0.1216 0.0219 0.3009 0.014        OR PM.M2     0.2164
#> 7       ACME.M3     0.0270 0.0120 0.0470 0.000      OR ACME.M3     1.2749
#> 8         PM.M3     0.1610 0.0706 0.3135 0.000        OR PM.M3     0.2805
#> 9           ADE     0.0955 0.0535 0.1380 0.000          OR ADE     2.2984
#> 10 Total Effect     0.1710 0.1180 0.2190 0.000 OR Total Effect     4.4488
#>     IC.inf IC.sup P.val                  . Estimation  IC.inf IC.sup P.val
#> 1   1.6064 2.3588 0.000   logOR ACME.joint     0.6630  0.4740 0.8581 0.000
#> 2   0.5038 0.7323 0.000     logOR PM.joint     0.4481  0.3328 0.5852 0.000
#> 3   0.9962 1.6349 0.058      logOR ACME.M1     0.2190 -0.0038 0.4916 0.058
#> 4  -0.0048 0.5065 0.058        logOR PM.M1     0.1467 -0.0023 0.3584 0.058
#> 5   1.0376 1.4833 0.016      logOR ACME.M2     0.1815  0.0369 0.3943 0.016
#> 6   0.0491 0.4411 0.016        logOR PM.M2     0.1216  0.0261 0.2930 0.016
#> 7   1.1242 1.4975 0.000      logOR ACME.M3     0.2429  0.1171 0.4038 0.000
#> 8   0.1385 0.4366 0.000        logOR PM.M3     0.1627  0.0755 0.2978 0.000
#> 9   1.5848 3.2651 0.000          logOR ADE     0.8322  0.4605 1.1833 0.000
#> 10  2.9467 6.5280 0.000 logOR Total Effect     1.4926  1.0807 1.8761 0.000

# 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

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)
#>  Treatment C1      C2      C3            M1               M2         
#>  0:487     0:487   0:399   0:590   Min.   :-1.005   Min.   :-0.3562  
#>  1:513     1:513   1:601   1:410   1st Qu.: 5.971   1st Qu.: 7.4968  
#>                                    Median : 8.068   Median :12.3671  
#>                                    Mean   : 8.640   Mean   :11.2222  
#>                                    3rd Qu.:11.451   3rd Qu.:14.2538  
#>                                    Max.   :18.107   Max.   :21.2199  
#>        M3          Outcome
#>  Min.   :-0.7423   0:516  
#>  1st Qu.: 5.5540   1: 64  
#>  Median : 9.6019   2:102  
#>  Mean   : 9.5298   3:318  
#>  3rd Qu.:13.3874          
#>  Max.   :19.1836
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")
#> Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
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")
Values of causal effects in data 3.
Total Effect Direct Effect Joint Indirect Effect Indirect Effect by M1 Indirect Effect by M2 Indirect Effect by M3
2.21 0.5 1.71 0.32 0.44 0.60 
summary(med.analysis,opt="avg")
#>               . Estimation IC.inf IC.sup P.val
#> 1    ACME.joint     1.5421 1.0425 1.9567     0
#> 2      PM.joint     0.6977 0.5774 0.8146     0
#> 3       ACME.M1     0.4165 0.1204 0.7676     0
#> 4         PM.M1     0.1929 0.0519 0.3616     0
#> 5       ACME.M2     0.3601 0.0345 0.7152     0
#> 6         PM.M2     0.1654 0.0153 0.3423     0
#> 7       ACME.M3     0.4564 0.2353 0.6824     0
#> 8         PM.M3     0.2104 0.1049 0.3440     0
#> 9           ADE     0.6629 0.4111 0.9683     0
#> 10 Total Effect     2.2050 1.6413 2.5277     0

Example 4

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)
#>  Treatment C1      C2      C3            M1         M2      M3     
#>  0:530     0:501   0:399   0:605   Min.   : 7.944   0:841   0:533  
#>  1:470     1:499   1:601   1:395   1st Qu.:16.721   1:159   1:224  
#>                                    Median :21.915           2:243  
#>                                    Mean   :21.988                  
#>                                    3rd Qu.:26.679                  
#>                                    Max.   :37.638                  
#>     Outcome      
#>  Min.   : 38.22  
#>  1st Qu.: 79.09  
#>  Median :109.25  
#>  Mean   :111.24  
#>  3rd Qu.:140.45  
#>  Max.   :194.37
M1reg=lm(M1~  Treatment + C1 + C2 + C3, data = data4)
M2reg=glm(M2~ Treatment + C1 + C3, data = data4, family = binomial("probit"))
#> Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
M3reg=polr(M3~Treatment + C2 + C3     , data = data4, method = "probit")
#> Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

Yreg=lm(Outcome~ Treatment + M1 + M2 + M3 + C1 + C2 + C3, data=data4)
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)
Values of causal effects in data 4.
Total Effect Direct Effect Joint Indirect Effect Indirect Effect by M1 Indirect Effect by M2 Indirect Effect by M3
57.68 10 47.48 34.2 0 13.48 
summary(med.analysis,opt="avg")
#>               . Estimation  IC.inf  IC.sup P.val
#> 1    ACME.joint    47.9419 46.1842 49.6012 0.000
#> 2      PM.joint     0.8270  0.8028  0.8496 0.000
#> 3       ACME.M1    33.8436 33.0623 34.5909 0.000
#> 4         PM.M1     0.5840  0.5662  0.5997 0.000
#> 5       ACME.M2     0.0009 -0.0028  0.0071 0.608
#> 6         PM.M2     0.0000  0.0000  0.0001 0.608
#> 7       ACME.M3    14.0974 12.5787 15.5737 0.000
#> 8         PM.M3     0.2432  0.2159  0.2695 0.000
#> 9           ADE    10.0248  8.7977 11.3318 0.000
#> 10 Total Effect    57.9667 56.9544 58.9418 0.000

Example 5

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)
#>     Exposure            M1               M2               M3        
#>  Min.   :0.0000   Min.   : 4.395   Min.   : 3.663   Min.   : 6.050  
#>  1st Qu.:0.0000   1st Qu.: 7.332   1st Qu.: 6.623   1st Qu.: 9.563  
#>  Median :1.0000   Median : 8.006   Median : 7.301   Median :10.269  
#>  Mean   :0.5097   Mean   : 8.004   Mean   : 7.324   Mean   :10.292  
#>  3rd Qu.:1.0000   3rd Qu.: 8.668   3rd Qu.: 8.031   3rd Qu.:11.011  
#>  Max.   :1.0000   Max.   :11.361   Max.   :11.176   Max.   :13.883  
#>      event     surv_time      
#>  Min.   :1   Min.   :0.00000  
#>  1st Qu.:1   1st Qu.:0.02000  
#>  Median :1   Median :0.04000  
#>  Mean   :1   Mean   :0.05808  
#>  3rd Qu.:1   3rd Qu.:0.08000  
#>  Max.   :1   Max.   :0.57000
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)
#> Loading required package: survival
model.y=aalen(Surv(surv_time, event) ~ const(Exposure) + const(M1) + const(M2) + const(M3), data = data5, robust=T)
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)
summary(multi.media, opt='avg')
#>               .  Estimation      IC.inf      IC.sup P.val
#> 1    ACME.joint 115205.6896  67042.0465 162811.6616 0.000
#> 2      PM.joint      0.5476      0.2625      1.1254 0.000
#> 3       ACME.M1   2348.2280  -4745.8431  10745.8561 0.458
#> 4         PM.M1      0.0110     -0.0213      0.0518 0.458
#> 5       ACME.M2  47750.3556   2288.7470  96271.4029 0.040
#> 6         PM.M2      0.2258      0.0104      0.5552 0.040
#> 7       ACME.M3  65107.1059  26265.4908 105670.6870 0.000
#> 8         PM.M3      0.3083      0.1100      0.6783 0.000
#> 9           ADE 114531.7341 -20309.6088 234776.5196 0.082
#> 10 Total Effect 229737.4237 113347.7146 344473.1502 0.000