The hardware and bandwidth for this mirror is donated by dogado GmbH, the Webhosting and Full Service-Cloud Provider. Check out our Wordpress Tutorial.
If you wish to report a bug, or if you are interested in having us mirror your free-software or open-source project, please feel free to contact us at mirror[@]dogado.de.
R package mma is used for general third-variable effect analysis (Yu and Li 2022). The third-variable effect (TVE) refers to the effect conveyed by intervening variables to an observed relationship between an exposure and a response variable (outcome). In this package, the exposure is also called the predictor, the intervening variables are called mediators or confounders (MC). The TVE include total effect, direct effect, and indirect effects from different third variables. The effects are defined in Yu et al. (2014). The algorithms on how to make inferences on TVEs and how to explain the effects are described in Yu et al. (2014).
To use the R package mma, we first install the package in R
(install.packages("mma")
) and load it.
library(mma)
#> Loading required package: gbm
#> Loaded gbm 2.1.8.1
#> Loading required package: splines
#> Loading required package: survival
#> Loading required package: car
#> Loading required package: carData
#> Loading required package: gplots
#>
#> Attaching package: 'gplots'
#> The following object is masked from 'package:stats':
#>
#> lowess
#source('C:/Users/qyu/Desktop/Research/mma/version current/R/mma.r')
We use weight_behavior data as an example to show how to do a preliminary data analysis to identify potential MCs and covariates. The function data.org is used for this purpose. The list of arguments for data.org are defined in Yu and Li (2017). In the following example, the predictor, “pred”, is gender, and the outcome, “y”, is a binary indicator of overweight (1) or not (0). The reference group of the exposure variable is the male as defined by “predref=”M”“. There are 12 variables in”x”, which includes all potential MCs and covariates. As an argument in the function data.org, “mediator=5:12” indicates that the variables from column 5 to column 12 in x should be tested for potential MCs. Other columns in “x” are to be treated as covariates. The “mediator” argument can also be a character vector with names of the potential MCs in x. Two tests are needed to identify a potential MC: first, the variable is significantly related with the predictor adjusting for other covariates. The other covariates here is defined by “cova”. If cova is a data frame by itself, all covariates in cova are for all potential MCs. If cova is a list, the first item will be the covariate data frame and the second item list the names or column numbers of MCs in “x” that should adjust for the covariates. The significance level is set by “alpha2”, whose default value is 0.1. Second, the variable has to be significantly related with the outcome adjusting (testtype=1, by default) or not adjusting (testtype=2) for the predictor and other variables. The significance level is set by “alpha”. In the following example, p-value 1 shows the results for the second test. and p-value 2 are the results for the first test. A variable that pass the second test but not the first test is considered as covariates. Variables do not pass the second test are not adopted in further analysis. Variables in “x” but not in “mediator” are forced in further analysis as covariates.
“jointm” is a list with the first item identify the number of groups of MCs whose joint effects are of interesting, where each group is defined in an item of “jointm”. In the above example, the joint effect of Variables 5, 7, 9 is of interesting. These variables are forced into analysis as potential MCs. A variable can be shown in “jointm” multiple times. The joint effect of each group of variables will be reported. A trick to force in variables as potential MCs is to assign them in a separate group of “jointm”.
data("weight_behavior")
#binary predictor x
#binary y
x=weight_behavior[,c(2,4:14)]
pred=weight_behavior[,"sex"]
y=weight_behavior[,"overweigh"]
data.b.b.1<-data.org(x,y,mediator=5:12,jointm=list(n=1,j1=7:9),
pred=pred,predref="M", alpha=0.4,alpha2=0.4)
summary(data.b.b.1)
#> Identified as mediators:
#> [1] "tvhours" "cmpthours" "cellhours" "exercises" "sweat" "sports"
#> Selected as covariates:
#> [1] "age" "race" "numpeople" "car"
#> Tests:
#> P-Value 1.y P-Value 2.pred
#> age 0.836 NA
#> race 0.850 NA
#> numpeople 0.561 NA
#> car 0.052 NA
#> gotosch 0.710 NA
#> snack 0.724 NA
#> tvhours - 0.830 0.484
#> cmpthours - 0.826 0.760
#> cellhours - 0.067 0.688
#> sports * 0.000 0.003
#> exercises * 0.176 0.203
#> sweat * 0.181 0.046
#> pred 0.178 NA
#> ----
#> *:mediator,-:joint mediator
#> P-Value 1:Type-3 tests in the full model (data.org) or estimated coefficients (data.org.big)
#> when testtype=1, univariate relationship test with the outcome when testtype=2
#> P-Value 2:Tests of relationship with the Predictor
The results from “data.org” function is summarized as above.
A potential MC defined in “mediator” argument is identified as categorical or binary if it has only two unique values or if the variable is a vector factors or characters. By default, the reference group is the first level of the variable. If the reference group needs to be changed, the categorical variable should be defined in binmed or catmed, and then the corresponding reference group in binref and catref. Continuous MCs can be defined by contmed. The following example did the same analysis as above, only that it defines that potential continuous mediators are columns 7,8,9,11, and 12 of x, binary mediators are columns 6 and 10, and categorical mediator is column 5 of x with 1 to be the reference group for all categorical and binary mediators.
The functions in mma can deal with binary, categorical, continuous, or time-to-event outcomes. If the outcome is time-to-event, it should be defined by the “Surv” function in the survival package. An example of application on breast cancer survival can be found in Yu et al. (2019).
In addition, the package can handle multivariate and multicategorical predictors. If the predictor is multicategorical of k levels, “data.org” first transform the predictor to k-1 binary predictors. If a variable significantly relates with any of the k-1 predictors, the variable passes the first test described above. P-value 2 is shown for each predictor. In the following example, the predictor is the race with six levels.
#multivariate predictor
x=weight_behavior[,c(2:3,5:14)]
pred=weight_behavior[,4]
y=weight_behavior[,15]
data.mb.b <-data.org(x,y,mediator=5:12,jointm=list(n=1,j1=c(5,7,9)),
pred=pred,predref="OTHER", alpha=0.4,alpha2=0.4)
summary(data.mb.b)
#> Identified as mediators:
#> [1] "tvhours" "cellhours" "sweat" "sports" "gotosch"
#> Selected as covariates:
#> [1] "age" "sex" "numpeople" "car" "exercises"
#> Tests:
#> P-Value 1.y P-Value 2.pred P-Value 2.pred.AFRICAN
#> age 0.836 NA NA
#> sex 0.178 NA NA
#> numpeople 0.561 NA NA
#> car 0.052 NA NA
#> gotosch - 0.710 0.112 0.796
#> snack 0.724 NA NA
#> tvhours - 0.830 0.748 0.535
#> cmpthours 0.826 NA NA
#> cellhours - 0.067 0.880 0.994
#> sports * 0.000 0.587 0.163
#> exercises 0.176 0.629 0.731
#> sweat * 0.181 0.647 0.174
#> pred 0.543 NA NA
#> pred.AFRICAN 0.663 NA NA
#> pred.CAUCASIAN 0.242 NA NA
#> pred.INDIAN 0.890 NA NA
#> pred.MIXED 0.782 NA NA
#> P-Value 2.pred.CAUCASIAN P-Value 2.pred.INDIAN
#> age NA NA
#> sex NA NA
#> numpeople NA NA
#> car NA NA
#> gotosch - 0.996 0.097
#> snack NA NA
#> tvhours - 0.334 0.916
#> cmpthours NA NA
#> cellhours - 0.707 0.555
#> sports * 0.816 0.453
#> exercises 0.723 0.912
#> sweat * 0.214 0.374
#> pred NA NA
#> pred.AFRICAN NA NA
#> pred.CAUCASIAN NA NA
#> pred.INDIAN NA NA
#> pred.MIXED NA NA
#> P-Value 2.pred.MIXED
#> age NA
#> sex NA
#> numpeople NA
#> car NA
#> gotosch - 0.033
#> snack NA
#> tvhours - 0.092
#> cmpthours NA
#> cellhours - 0.383
#> sports * 0.091
#> exercises 0.881
#> sweat * 0.203
#> pred NA
#> pred.AFRICAN NA
#> pred.CAUCASIAN NA
#> pred.INDIAN NA
#> pred.MIXED NA
#> ----
#> *:mediator,-:joint mediator
#> P-Value 1:Type-3 tests in the full model (data.org) or estimated coefficients (data.org.big)
#> when testtype=1, univariate relationship test with the outcome when testtype=2
#> P-Value 2:Tests of relationship with the Predictor
Similarly, the package can deal with multivariate outcomes. The following code deals with multivariate predictors and multivariate responses. If the variable is significantly related with any one of the outcomes, it passes the second test described above. The results from “data.org” are summarized for each combination of the exposure-outcome relationship.
#multivariate responses
x=weight_behavior[,c(2:3,5:14)]
pred=weight_behavior[,4]
y=weight_behavior[,c(1,15)]
data.mb.mb<-data.org(x,y,mediator=5:12,jointm=list(n=1,j1=c(5,7,9)),
pred=pred,predref="OTHER", alpha=0.4,alpha2=0.4)
summary(data.mb.mb)
#> Identified as mediators:
#> [1] "tvhours" "cellhours" "sweat" "sports" "gotosch"
#> Selected as covariates:
#> [1] "age" "sex" "numpeople" "car" "exercises"
#> Tests:
#> P-Value 1.bmi P-Value 1.overweigh P-Value 2.pred
#> age 0.458 0.836 NA
#> sex 0.003 0.178 NA
#> numpeople 0.282 0.561 NA
#> car 0.059 0.052 NA
#> gotosch - 0.527 0.710 0.112
#> snack 0.830 0.724 NA
#> tvhours - 0.505 0.830 0.748
#> cmpthours 0.676 0.826 NA
#> cellhours - 0.084 0.067 0.880
#> sports * 0.001 0.000 0.587
#> exercises 0.089 0.176 0.629
#> sweat * 0.078 0.181 0.647
#> pred 0.344 0.543 NA
#> pred.AFRICAN 0.900 0.663 NA
#> pred.CAUCASIAN 0.032 0.242 NA
#> pred.INDIAN 0.446 0.890 NA
#> pred.MIXED 0.833 0.782 NA
#> P-Value 2.pred.AFRICAN P-Value 2.pred.CAUCASIAN
#> age NA NA
#> sex NA NA
#> numpeople NA NA
#> car NA NA
#> gotosch - 0.796 0.996
#> snack NA NA
#> tvhours - 0.535 0.334
#> cmpthours NA NA
#> cellhours - 0.994 0.707
#> sports * 0.163 0.816
#> exercises 0.731 0.723
#> sweat * 0.174 0.214
#> pred NA NA
#> pred.AFRICAN NA NA
#> pred.CAUCASIAN NA NA
#> pred.INDIAN NA NA
#> pred.MIXED NA NA
#> P-Value 2.pred.INDIAN P-Value 2.pred.MIXED
#> age NA NA
#> sex NA NA
#> numpeople NA NA
#> car NA NA
#> gotosch - 0.097 0.033
#> snack NA NA
#> tvhours - 0.916 0.092
#> cmpthours NA NA
#> cellhours - 0.555 0.383
#> sports * 0.453 0.091
#> exercises 0.912 0.881
#> sweat * 0.374 0.203
#> pred NA NA
#> pred.AFRICAN NA NA
#> pred.CAUCASIAN NA NA
#> pred.INDIAN NA NA
#> pred.MIXED NA NA
#> ----
#> *:mediator,-:joint mediator
#> P-Value 1:Type-3 tests in the full model (data.org) or estimated coefficients (data.org.big)
#> when testtype=1, univariate relationship test with the outcome when testtype=2
#> P-Value 2:Tests of relationship with the Predictor
##Third-Variable Effect Analysis Next, we use med function to estimate the TVE using the results from data.org function. The default approach is using a generalized linear model to fit the final full model. If “nonlinear=TRUE”“, multiple Additive Regression Trees (MART) will be used to fit the relationship between outcome(s) and other variables, and smoothing splines used to fit the relationship between the potential MC and exposure variables(s).
Using the linear method and to show the results:
med.b.b.1<-med(data=data.b.b.2,n=2)
med.b.b.1
#>
#>
#> For the predictor pred :
#> The estimated total effect:[1] 0.6156
#>
#> The estimated indirect effect:
#> y1.all y1.tvhours y1.cmpthours y1.cellhours y1.exercises y1.sweat
#> 0.2158 -0.0007 0.0207 0.0177 -0.0177 0.0080
#> y1.sports y1.j1
#> 0.1633 0.0192
We can plot the marginal effect of the MC on the outcome, and the marginal effect of the predictor on the MC. If the predictor is binary, plot.med draws a histogram or boxplot of the marginal density of the MC at each different value of the predictor.
For survival outcome, the default option is to fit the final full model using Cox proportional hazard model. Note that when use type=“lp”, the linear predictor is the type of predicted value as in predict.coxph, the results is similar to those from MART. Readers are referred to Yu et al. (2019) for more details.
boot.med function is used to make inferences on the TVE. Results return from data.org function is input as the first argument.
bootmed.b.b.1<-boot.med(data=data.b.b.2,n=2,n2=50)
#> [1] 1
#> [1] 2
#> [1] 3
#> [1] 4
#> [1] 5
#> [1] 6
#> [1] 7
#> [1] 8
#> [1] 9
#> [1] 10
#> [1] 11
#> [1] 12
#> [1] 13
#> [1] 14
#> [1] 15
#> [1] 16
#> [1] 17
#> [1] 18
#> [1] 19
#> [1] 20
#> [1] 21
#> [1] 22
#> [1] 23
#> [1] 24
#> [1] 25
#> [1] 26
#> [1] 27
#> [1] 28
#> [1] 29
#> [1] 30
#> [1] 31
#> [1] 32
#> [1] 33
#> [1] 34
#> [1] 35
#> [1] 36
#> [1] 37
#> [1] 38
#> [1] 39
#> [1] 40
#> [1] 41
#> [1] 42
#> [1] 43
#> [1] 44
#> [1] 45
#> [1] 46
#> [1] 47
#> [1] 48
#> [1] 49
#> [1] 50
summary(bootmed.b.b.1)
#> MMA Analysis: Estimated Mediation Effects Using GLM
#> For Predictor/Moderator at pred
#> $total.effect
#> est mean sd upbd lwbd upbd_q lwbd_q upbd_bcbi
#> 0.517 0.660 0.421 1.485 -0.165 1.480 -0.094 1.165
#> lwbd_bcbi upbd_b lwbd_b upbd_win lwbd_win p_norm p_quan
#> -0.212 1.675 -0.235 1.439 -0.111 0.117 0.120
#>
#> $direct.effect
#> est mean sd upbd lwbd upbd_q lwbd_q upbd_bcbi
#> 0.369 0.403 0.419 1.224 -0.418 1.120 -0.341 1.120
#> lwbd_bcbi upbd_b lwbd_b upbd_win lwbd_win p_norm p_quan
#> -0.341 1.538 -0.372 1.154 -0.283 0.336 0.360
#>
#> $indirect.effect
#> y1.all y1.tvhours y1.cmpthours y1.cellhours y1.exercises y1.sweat
#> est 0.148 -0.007 0.051 -0.007 -0.013 0.000
#> mean 0.257 -0.013 0.036 0.027 -0.016 0.011
#> sd 0.149 0.036 0.059 0.053 0.028 0.025
#> upbd 0.549 0.057 0.152 0.130 0.038 0.061
#> lwbd -0.035 -0.083 -0.080 -0.076 -0.071 -0.039
#> upbd_q 0.528 0.052 0.182 0.139 0.020 0.065
#> lwbd_q -0.074 -0.085 -0.041 -0.051 -0.086 -0.033
#> upbd_bcbi 0.371 0.078 0.192 0.045 0.010 0.042
#> lwbd_bcbi -0.126 -0.062 -0.036 -0.069 -0.089 -0.039
#> upbd_b 0.529 0.086 0.194 0.145 0.021 0.068
#> lwbd_b -0.127 -0.087 -0.050 -0.070 -0.089 -0.040
#> upbd_win 0.541 0.052 0.187 0.161 0.013 0.067
#> lwbd_win -0.074 -0.087 -0.040 -0.036 -0.089 -0.026
#> p_norm 0.085 0.719 0.543 0.610 0.554 0.661
#> p_quan 0.080 0.680 0.520 0.600 0.600 0.800
#> y1.sports y1.j1
#> est 0.124 0.011
#> mean 0.169 0.072
#> sd 0.084 0.095
#> upbd 0.334 0.257
#> lwbd 0.003 -0.114
#> upbd_q 0.367 0.246
#> lwbd_q 0.002 -0.094
#> upbd_bcbi 0.172 0.138
#> lwbd_bcbi -0.031 -0.111
#> upbd_b 0.387 0.308
#> lwbd_b -0.032 -0.111
#> upbd_win 0.369 0.250
#> lwbd_win 0.022 -0.093
#> p_norm 0.045 0.449
#> p_quan 0.080 0.440
Similarly, the default approach is using a generalized linear model to fit the final full model. If the argument nonlinear=TRUE, multiple Additive Regression Trees (MART) will be used to fit the model in estimating the outcome. Applications for nonlinear Mediation Analysis are described in Yu et al. (2017, 2018).
summary(bootmed.b.b.2)
#> MMA Analysis: Estimated Mediation Effects Using MART
#> For Predictor/Moderator at pred
#> $total.effect
#> est mean sd upbd lwbd upbd_q lwbd_q upbd_bcbi
#> 0.164 0.194 0.084 0.359 0.030 0.287 0.123 0.287
#> lwbd_bcbi upbd_b lwbd_b upbd_win lwbd_win p_norm p_quan
#> 0.123 0.291 0.126 0.287 0.123 0.020 0.000
#>
#> $direct.effect
#> est mean sd upbd lwbd upbd_q lwbd_q upbd_bcbi
#> 0.070 0.021 0.074 0.167 -0.125 0.122 -0.020 0.132
#> lwbd_bcbi upbd_b lwbd_b upbd_win lwbd_win p_norm p_quan
#> -0.018 0.132 -0.017 0.122 -0.018 0.776 0.500
#>
#> $indirect.effect
#> y1.all y1.tvhours y1.cmpthours y1.cellhours y1.exercises y1.sweat
#> est 0.094 0.002 0.003 0.010 -0.003 0.006
#> mean 0.173 0.015 0.028 0.008 -0.016 0.000
#> sd 0.050 0.025 0.014 0.007 0.021 0.006
#> upbd 0.272 0.063 0.056 0.022 0.024 0.011
#> lwbd 0.075 -0.034 0.000 -0.006 -0.057 -0.011
#> upbd_q 0.241 0.035 0.041 0.017 0.005 0.006
#> lwbd_q 0.143 -0.014 0.011 0.002 -0.042 -0.005
#> upbd_bcbi 0.143 0.035 0.010 0.018 0.006 0.007
#> lwbd_bcbi 0.143 -0.015 0.010 0.003 -0.022 -0.005
#> upbd_b 0.248 0.036 0.038 0.018 0.006 0.007
#> lwbd_b 0.143 0.003 0.010 0.004 -0.044 -0.005
#> upbd_win 0.241 0.035 0.041 0.017 0.005 0.006
#> lwbd_win 0.143 -0.014 0.011 0.003 -0.042 -0.005
#> p_norm 0.001 0.556 0.050 0.276 0.430 0.992
#> p_quan 0.000 0.500 0.000 0.000 0.500 1.000
#> y1.sports y1.j1
#> est 0.078 0.015
#> mean 0.147 0.029
#> sd 0.055 0.029
#> upbd 0.255 0.086
#> lwbd 0.038 -0.028
#> upbd_q 0.194 0.061
#> lwbd_q 0.075 -0.003
#> upbd_bcbi 0.175 0.045
#> lwbd_bcbi 0.068 -0.004
#> upbd_b 0.196 0.063
#> lwbd_b 0.068 0.017
#> upbd_win 0.194 0.061
#> lwbd_win 0.068 -0.003
#> p_norm 0.008 0.321
#> p_quan 0.000 0.500
plot.mma() plots the marginal effect of the selected variable on the outcome, and the marginal effect of the predictor on the selected variable. Interpretations for similar plots can be found in Yu et al. (2017).
The mma function is a combined function that automatically identify potential MCs, based on which to make statistical inference on the mediation effects.
x=weight_behavior[,c(2,4:14)]
pred=weight_behavior[,3]
y=weight_behavior[,15]
mma.b.b.glm<-mma(x,y,pred=pred,contmed=c(7:9,11:12),binmed=c(6,10),binref=c(1,1), catmed=5,catref=1,predref="M",alpha=0.4,alpha2=0.4,n=2,n2=2)
#> [1] 1
#> [1] 2
summary(mma.b.b.glm)
#> MMA Analysis: Estimated Mediation Effects Using GLM
#> For Predictor/Moderator at pred
#> $total.effect
#> est mean sd upbd lwbd upbd_q lwbd_q upbd_bcbi
#> 0.572 0.956 0.304 1.551 0.361 1.160 0.752 0.742
#> lwbd_bcbi upbd_b lwbd_b upbd_win lwbd_win p_norm p_quan
#> 0.742 1.171 1.171 1.160 0.752 0.002 0.000
#>
#> $direct.effect
#> est mean sd upbd lwbd upbd_q lwbd_q upbd_bcbi
#> 0.417 0.785 0.304 1.381 0.189 0.989 0.581 0.570
#> lwbd_bcbi upbd_b lwbd_b upbd_win lwbd_win p_norm p_quan
#> 0.570 1.000 1.000 0.989 0.581 0.010 0.000
#>
#> $indirect.effect
#> y1.all y1.exercises y1.sweat y1.sports
#> est 0.155 0.014 0.000 0.091
#> mean 0.171 -0.024 0.008 0.174
#> sd 0.001 0.011 0.019 0.003
#> upbd 0.173 -0.002 0.045 0.181
#> lwbd 0.170 -0.046 -0.029 0.168
#> upbd_q 0.172 -0.016 0.020 0.177
#> lwbd_q 0.171 -0.031 -0.005 0.172
#> upbd_bcbi 0.171 -0.016 0.020 0.172
#> lwbd_bcbi 0.171 -0.016 -0.005 0.172
#> upbd_b 0.172 0.172 0.172 0.172
#> lwbd_b -0.016 -0.016 -0.016 -0.016
#> upbd_win 0.172 -0.016 0.020 0.177
#> lwbd_win 0.171 -0.031 -0.005 0.172
#> p_norm 0.000 0.034 0.677 0.000
#> p_quan 0.000 0.000 1.000 0.000
When the argument nonlinear=TRUE, multiple Additive Regression Trees (MART) will be used to fit the final full model in estimating the outcome.
mma.b.b.mart<-mma(x,y,pred=pred,contmed=c(7:9,11:12),binmed=c(6,10),binref=c(1,1), catmed=5,catref=1,predref="M",alpha=0.4,alpha2=0.4,nonlinear=TRUE,n=2,n2=5)
summary(mma.b.b.mart)
summary(mma.b.b.mart)
#> MMA Analysis: Estimated Mediation Effects Using MART
#> For Predictor/Moderator at pred
#> $total.effect
#> est mean sd upbd lwbd upbd_q lwbd_q upbd_bcbi
#> 0.123 0.362 0.075 0.509 0.214 0.448 0.257 0.246
#> lwbd_bcbi upbd_b lwbd_b upbd_win lwbd_win p_norm p_quan
#> 0.246 0.387 0.246 0.453 0.249 0.000 0.000
#>
#> $direct.effect
#> est mean sd upbd lwbd upbd_q lwbd_q upbd_bcbi
#> 0.067 0.262 0.104 0.467 0.058 0.322 0.098 0.078
#> lwbd_bcbi upbd_b lwbd_b upbd_win lwbd_win p_norm p_quan
#> 0.078 0.322 0.078 0.322 0.078 0.012 0.000
#>
#> $indirect.effect
#> y1.all y1.exercises y1.sweat y1.sports
#> est 0.056 -0.001 0.003 0.065
#> mean 0.099 0.002 0.009 0.110
#> sd 0.067 0.015 0.018 0.046
#> upbd 0.231 0.032 0.044 0.201
#> lwbd -0.032 -0.028 -0.026 0.020
#> upbd_q 0.175 0.021 0.033 0.156
#> lwbd_q 0.036 -0.016 -0.010 0.045
#> upbd_bcbi 0.173 0.021 0.030 0.125
#> lwbd_bcbi 0.035 -0.013 -0.011 0.038
#> upbd_b 0.168 0.014 0.035 0.157
#> lwbd_b 0.034 -0.017 -0.011 0.038
#> upbd_win 0.176 0.021 0.035 0.157
#> lwbd_win 0.040 -0.013 -0.008 0.040
#> p_norm 0.139 0.873 0.602 0.017
#> p_quan 0.000 0.800 0.400 0.000
These binaries (installable software) and packages are in development.
They may not be fully stable and should be used with caution. We make no claims about them.
Health stats visible at Monitor.