multimediate

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.5494 79.8501 83.0861     0
#> 2      PM.joint     0.8859  0.8781  0.8942     0
#> 3       ACME.M1    10.3988  9.6326 11.1910     0
#> 4         PM.M1     0.1130  0.1049  0.1218     0
#> 5       ACME.M2    17.5109 16.6270 18.3698     0
#> 6         PM.M2     0.1902  0.1804  0.2003     0
#> 7       ACME.M3    53.6397 52.3675 54.9432     0
#> 8         PM.M3     0.5828  0.5658  0.5995     0
#> 9           ADE    10.5010  9.7397 11.2122     0
#> 10 Total Effect    92.0504 90.5099 93.5028     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.0770  0.0550 0.0990 0.000   OR ACME.joint     1.9387
#> 2      PM.joint     0.4474  0.3292 0.6112 0.000     OR PM.joint     0.6286
#> 3       ACME.M1     0.0245 -0.0005 0.0580 0.054      OR ACME.M1     1.2400
#> 4         PM.M1     0.1429 -0.0031 0.3694 0.054        OR PM.M1     0.2482
#> 5       ACME.M2     0.0210  0.0030 0.0490 0.024      OR ACME.M2     1.2043
#> 6         PM.M2     0.1232  0.0190 0.3012 0.024        OR PM.M2     0.2198
#> 7       ACME.M3     0.0280  0.0125 0.0475 0.000      OR ACME.M3     1.2733
#> 8         PM.M3     0.1628  0.0738 0.3019 0.000        OR PM.M3     0.2777
#> 9           ADE     0.0950  0.0520 0.1395 0.000          OR ADE     2.2671
#> 10 Total Effect     0.1720  0.1240 0.2240 0.000 OR Total Effect     4.4320
#>     IC.inf IC.sup P.val                  . Estimation  IC.inf IC.sup P.val
#> 1   1.6156 2.3564 0.000   logOR ACME.joint     0.6620  0.4797 0.8572 0.000
#> 2   0.5163 0.7429 0.000     logOR PM.joint     0.4474  0.3312 0.6040 0.000
#> 3   0.9925 1.6333 0.060      logOR ACME.M1     0.2151 -0.0076 0.4906 0.060
#> 4  -0.0102 0.5193 0.060        logOR PM.M1     0.1445 -0.0055 0.3662 0.060
#> 5   1.0357 1.4992 0.032      logOR ACME.M2     0.1859  0.0351 0.4049 0.032
#> 6   0.0477 0.4329 0.032        logOR PM.M2     0.1248  0.0247 0.2842 0.032
#> 7   1.1359 1.4911 0.000      logOR ACME.M3     0.2416  0.1274 0.3995 0.000
#> 8   0.1528 0.4410 0.000        logOR PM.M3     0.1623  0.0808 0.2966 0.000
#> 9   1.5484 3.2267 0.000          logOR ADE     0.8185  0.4372 1.1715 0.000
#> 10  2.9130 6.6957 0.000 logOR Total Effect     1.4889  1.0692 1.9015 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.5281 0.8425 1.9638  0.00
#> 2      PM.joint     0.6993 0.5295 0.8447  0.00
#> 3       ACME.M1     0.4036 0.1178 0.6961  0.00
#> 4         PM.M1     0.1905 0.0504 0.3798  0.00
#> 5       ACME.M2     0.3570 0.0509 0.6855  0.04
#> 6         PM.M2     0.1689 0.0222 0.3522  0.04
#> 7       ACME.M3     0.4568 0.2180 0.7421  0.00
#> 8         PM.M3     0.2153 0.0940 0.3852  0.00
#> 9           ADE     0.6520 0.3229 1.0757  0.00
#> 10 Total Effect     2.1800 1.4735 2.5175  0.00

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.9586 46.2903 49.5585 0.000
#> 2      PM.joint     0.8272  0.8061  0.8469 0.000
#> 3       ACME.M1    33.8316 33.0873 34.6129 0.000
#> 4         PM.M1     0.5837  0.5671  0.6016 0.000
#> 5       ACME.M2     0.0008 -0.0026  0.0056 0.644
#> 6         PM.M2     0.0000  0.0000  0.0001 0.644
#> 7       ACME.M3    14.1262 12.7203 15.5276 0.000
#> 8         PM.M3     0.2437  0.2183  0.2694 0.000
#> 9           ADE    10.0152  8.9144 11.1702 0.000
#> 10 Total Effect    57.9738 56.9722 58.9797 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 113526.7278  66341.5781 163930.8039 0.000
#> 2      PM.joint      0.5494      0.2540      1.1237 0.000
#> 3       ACME.M1   2103.6958  -4308.7569   9567.5636 0.514
#> 4         PM.M1      0.0104     -0.0210      0.0510 0.514
#> 5       ACME.M2  48168.9636  -1098.7624  94776.2611 0.056
#> 6         PM.M2      0.2322     -0.0043      0.5559 0.056
#> 7       ACME.M3  63254.0684  23719.4878 103248.6928 0.000
#> 8         PM.M3      0.3047      0.0965      0.6380 0.000
#> 9           ADE 112991.2513 -13482.2225 233612.2017 0.086
#> 10 Total Effect 226517.9791 112864.3347 342020.5498 0.000