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.

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.

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.5751 79.9442 83.2013     0
#> 2      PM.joint     0.8860  0.8785  0.8942     0
#> 3       ACME.M1    10.3920  9.6338 11.1993     0
#> 4         PM.M1     0.1129  0.1045  0.1216     0
#> 5       ACME.M2    17.5171 16.6219 18.4503     0
#> 6         PM.M2     0.1903  0.1800  0.2013     0
#> 7       ACME.M3    53.6660 52.3990 54.9648     0
#> 8         PM.M3     0.5829  0.5665  0.6007     0
#> 9           ADE    10.4924  9.7525 11.1828     0
#> 10 Total Effect    92.0675 90.5881 93.5012     0

#plot(med.analysis)

Example 2

Binary treatment : \(T\)
Binary Outcome : \(Y\)
Linear mediators : \(M_1,M_2,M_3\)
Binary covariables : \(C_1,C_2,C_3\)

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.0550 0.0995 0.000   OR ACME.joint     1.9354
#> 2      PM.joint     0.4518 0.3151 0.6121 0.000     OR PM.joint     0.6312
#> 3       ACME.M1     0.0250 0.0000 0.0590 0.048      OR ACME.M1     1.2455
#> 4         PM.M1     0.1448 0.0000 0.3762 0.048        OR PM.M1     0.2574
#> 5       ACME.M2     0.0210 0.0030 0.0485 0.014      OR ACME.M2     1.2075
#> 6         PM.M2     0.1232 0.0172 0.3074 0.014        OR PM.M2     0.2228
#> 7       ACME.M3     0.0275 0.0120 0.0480 0.000      OR ACME.M3     1.2690
#> 8         PM.M3     0.1613 0.0659 0.3194 0.000        OR PM.M3     0.2780
#> 9           ADE     0.0940 0.0510 0.1405 0.000          OR ADE     2.2714
#> 10 Total Effect     0.1710 0.1190 0.2220 0.000 OR Total Effect     4.4184
#>    IC.inf IC.sup P.val                  . Estimation IC.inf IC.sup P.val
#> 1  1.6094 2.3513 0.000   logOR ACME.joint     0.6603 0.4759 0.8550 0.000
#> 2  0.4989 0.7458 0.000     logOR PM.joint     0.4516 0.3251 0.6153 0.000
#> 3  1.0023 1.6607 0.050      logOR ACME.M1     0.2195 0.0023 0.5072 0.050
#> 4  0.0029 0.5222 0.050        logOR PM.M1     0.1478 0.0014 0.3646 0.050
#> 5  1.0325 1.4767 0.016      logOR ACME.M2     0.1886 0.0320 0.3898 0.016
#> 6  0.0414 0.4420 0.016        logOR PM.M2     0.1269 0.0222 0.2977 0.016
#> 7  1.1264 1.4915 0.000      logOR ACME.M3     0.2382 0.1191 0.3998 0.000
#> 8  0.1472 0.4529 0.000        logOR PM.M3     0.1603 0.0779 0.3157 0.000
#> 9  1.5477 3.3132 0.000          logOR ADE     0.8204 0.4368 1.1979 0.000
#> 10 2.8186 6.7064 0.000 logOR Total Effect     1.4858 1.0362 1.9031 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

Binary treatment \(T\)
Polytomic Outcome \(Y\)
Linear mediators \(M_1,M_2,M_3\)
Binary covariables \(C_1,C_2,C_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.6047 1.1323 1.9730  0.00
#> 2      PM.joint     0.7198 0.5873 0.8651  0.00
#> 3       ACME.M1     0.3571 0.0292 0.7153  0.02
#> 4         PM.M1     0.1648 0.0140 0.3641  0.02
#> 5       ACME.M2     0.4230 0.1086 0.8363  0.00
#> 6         PM.M2     0.1907 0.0508 0.3847  0.00
#> 7       ACME.M3     0.4936 0.2290 0.7829  0.00
#> 8         PM.M3     0.2281 0.0956 0.4302  0.00
#> 9           ADE     0.6242 0.2862 1.0028  0.00
#> 10 Total Effect     2.2289 1.8156 2.5385  0.00

Example 4

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\)

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.9719 46.4846 49.5663 0.000
#> 2      PM.joint     0.8270  0.8056  0.8469 0.000
#> 3       ACME.M1    33.8248 33.0611 34.6076 0.000
#> 4         PM.M1     0.5832  0.5650  0.5998 0.000
#> 5       ACME.M2     0.0008 -0.0026  0.0059 0.624
#> 6         PM.M2     0.0000  0.0000  0.0001 0.624
#> 7       ACME.M3    14.1464 12.7963 15.5253 0.000
#> 8         PM.M3     0.2439  0.2199  0.2686 0.000
#> 9           ADE    10.0358  8.8892 11.2862 0.000
#> 10 Total Effect    58.0077 57.0576 59.1006 0.000

Example 5

Binary treatment \(T\)
Survival Outcome \(Y\)
Linear mediators \(M_1,M_2,M_3\)

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 110747.9799  62646.1987 160435.6382 0.000
#> 2      PM.joint      0.5318      0.2478      1.0495 0.000
#> 3       ACME.M1   2645.2348  -3850.4118  10732.0961 0.404
#> 4         PM.M1      0.0126     -0.0185      0.0559 0.404
#> 5       ACME.M2  47252.5758   2643.0232  94616.5612 0.042
#> 6         PM.M2      0.2285      0.0128      0.5896 0.042
#> 7       ACME.M3  60850.1693  21385.8671  97639.5324 0.002
#> 8         PM.M3      0.2914      0.0924      0.6254 0.002
#> 9           ADE 114408.2405  -6823.4424 231804.7809 0.072
#> 10 Total Effect 225156.2204 115043.4609 337720.8040 0.000