mediation analysis in R with incomplete data I Multiple Imputation

[Xibit]

I want to conduct a mediation analysis in R with incomplete data for my master's thesis. For the missings I use MICE multiple imputation according to van Buuren & Groothuis-Oudshoorn (2011).

For the mediation analysis I have considered to refer to Tingley et al. (2014)

Here, you can find his code: https://cran.r-project.org/web/packages/mediation/mediation.pdf

# Hypothetical example
## To use mediations, must make list of multiple datasets. Then, ## must also repeat the treatment assignment list as many times ## as you have data sets.
datasets <- list(D1=D1, D2=D2) # list of multiply imputed data sets
mediators <- c("M1")
outcome <- c("Ycont1")
treatment <- c("T1","T1") # note how the treatment indicator is repeated
covariates <- c("X1+X2")
olsols <- mediations(datasets, treatment, mediators, outcome, covariates,
    families=c("gaussian","gaussian"), interaction=FALSE, conf.level=.90, sims=1000)
output <- amelidiate(olsols) # summary(output)
plot(output)

For my master thesis I want to find out if the therapeutic alliance between an eCoach and a patient mediates the relationship of the number of completed modules of an internet-based depression intervention and the severity of depression.

Number of completed modules (ModuleEr) -> therapeutic alliance (WAI_C) -> severity of depression (QIDS_t1)

My variables are numeric. I am especially interested in the standardized coefficients.I tried Tingley’s Code aswell and adapted to my variables:

df <-data
library(mice)
library(mediation)

simple imputation for testing reasons

imp <- mice(df, m=5, seed = 1234)

I convert every imputed data set

imp1 <- complete(imp,1)
imp2 <- complete(imp,2)
imp3 <- complete(imp,3)
imp4 <- complete(imp,4)
imp5 <- complete(imp,5)

I standardize my data sets, only numeric variables

imp1z <- scale(imp1[,1:13])
imp2z <- scale(imp2[,1:13])
imp3z <- scale(imp3[,1:13])
imp4z <- scale(imp4[,1:13])
imp5z <- scale(imp5[,1:13])

I make two lists of multiply imputed data sets: one list with standardized variables, one with unstandardized variables

datasets1 <- list(imp1=imp1, imp2=imp2, imp3=imp3, imp4=imp4, imp5=imp5)
datasets2<- list(imp1z=imp1z, imp2z=imp2z, imp3z=imp3z, imp4z=imp4z, imp5z=imp5z)

then I define my mediator, outcome, treatment an covarites

mediators <- c("WAI_C")
outcome <- c("QIDS_t1")
treatment <- c("ModuleEr", "ModuleEr", "ModuleEr", "ModuleEr", "ModuleEr" ) # is it correct to repeat the treatment as many times as you have data sets? I hav read it here: https://rdrr.io/cran/mediation/src/R/amelidiate.R
covariates <- c("QIDS_t0")

next I use the mediations function with unstandardized variables. I was wondering how to define the families=c("gaussian","gaussian“). My variables are not normally distributed…

olsols1 <- mediations(datasets1, treatment, mediators, outcome, covariates, interaction=FALSE, conf.level=.95, sims=10)
output1 <- amelidiate(olsols1)
summary(output1)
plot(output1) # the plot function doesn’t work: Error in plot.new() : figure margins too large

My output 1 looks like this

If I repeat the mediations function with the standardized data sets (datasets2), I I get an error messages:

olsols2 <- mediations(datasets2, treatment, mediators, outcome, covariates, interaction=FALSE, conf.level=.95, sims=10)
output2 <- amelidiate(olsols2)
summary(output2)

Error in eval(predvars, data, env) : Object 'WAI_C' not found
Additional: Warning:
In dataarg$weight <- weight <- rep(1, nrow(dataarg)) :
Convert left side to a list

What is wrong? Why can WAI_C not been found? I checked all 5 standardized data sets, they all include WAI_C. Why is especially the mediator variable affected?

And how can I convert left side to a list?

Since the output1 only shows ACME,ADE, total Effect.. How do I obtain the associated regression coefficients? Is it enough to simply set up the regression equations?

Med1 <- with(data = imp, exp = lm(scale(WAI_C) ~ scale(ModuleEr) + scale(QIDS_t0)))
Med1.pooled <- pool(Med1)
summary(Med1.pooled)

Med2 <- with(data = imp, exp = lm(scale(QIDS_t1) ~ scale(ModuleEr) + scale(WAI_P) + scale(QIDS_t0)))
Med2.pooled <- pool(Med2)
summary(Med2.pooled)