bootstrapLavaan works for single-level mediation but not multilevel mediation

[daqing liu]

Hi everyone,

I used lavaan to build a 1-1-1 multilevel mediation model, whcih worked ok. But when I try to use bootstrapLavaan() to get the CIs with a function to extract the unstandardised coefficients (a*b of particular interest), the new CIs produced by bootstraplavaan() did not match the estimated coefficients, that is, the coefficients of the model were not in the middle of the CIs, and sometimes even out of the CIs. 

I tried to build a single level mediation model with the same data, using the cluster-robust SEs, and then used bootstrapLavaan() to get the CIs for the unstandardised coefficients. The results seemed right; all the CIs matched the coefficients. This makes me belive my bootstrap code was fine (I followed some suggestions found here), but I couldn't understand why the bootstrapLavaan() works for single-level mediation, but not multilevel mediation. 

Is this because the bootsrapLavaan() can't simulate multilevel data well, or is there any particular parameters I need to specify when I fit the multilevel model or bootsrap it?

The next question is not about lavaan, but I'd appreciate it if anyone could share some thoughts. I've see suggestion of using the Monte Carlo method to get the CIs for the indirect path for 1-1-1 mediation model ( http://www.quantpsy.org/medmc/medmc111.htm). But because my model only has random intercepts, this method, which requires level-2 covariance between the a_j and b_j slopes, is not applicable. I wonder if I could use the single level Monte Carlo method (http://www.quantpsy.org/medmc/medmc.htm) to get the CIs for my multilevel mediation model, since the slopes of the coefficients are not varying anyway?

Thanks！

Daqing

[Terrence Jorgensen]

Is this because the bootsrapLavaan() can't simulate multilevel data well

It doesn't simulate, it resamples.  I'm not sure it is designed to work with multilevel data.  What would you resample?

 

I've see suggestion of using the Monte Carlo method to get the CIs for the indirect path for 1-1-1 mediation model

That's what I would suggest, even for a single-level model.  If you read the simulation studies closely (even among those who recommend bootstrapping), the benefits of bootstrapping have always been negligible, only manifesting minor power advantages compared to less-intensive methods when samples are so small that the Type II error rate is huge regardless (so who really cares if you have 32% power instead of 27% power?).

 

because my model only has random intercepts, this method, which requires level-2 covariance between the a_j and b_j slopes, is not applicable

Well, if your slopes don't vary, you should be able to simply set their (co)variances to 0.

 

I wonder if I could use the single level Monte Carlo method (http://www.quantpsy.org/medmc/medmc.htm) to get the CIs for my multilevel mediation model, since the slopes of the coefficients are not varying anyway?

You should be able to pass your lavaan model to semTools::monteCarloCI(), and it will give your MC CIs for any user-defined parameters in your model.

Terrence D. Jorgensen

Assistant Professor, Methods and Statistics

Research Institute for Child Development and Education, the University of Amsterdam

http://www.uva.nl/profile/t.d.jorgensen

[daqing liu]

Hi Terrence,

Thank you so much for your reply.  I'll defintely look up semTools::monteCarloCI(). 

I did find bootstrapping and the MC method give very similar results to the model's estimated CIs.

 Just a side note, the interactive page (MC CIs for multilevel mediation, http://www.quantpsy.org/medmc/medmc111.htm) doesn't allow the covariance between a and b slopes to be 0. But it's an informative tool to use, if anyone is interested.

Best wishes,

Daqing

[Terrence Jorgensen]

the interactive page (MC CIs for multilevel mediation, http://www.quantpsy.org/medmc/medmc111.htm) doesn't allow the covariance between a and b slopes to be 0. 

Yes it does.  I just plugged in some arbitrary value for other parameters, highlighting the covab=0 in the generated syntax below.

Terrence D. Jorgensen    (he, him, his)

Assistant Professor, Methods and Statistics
Research Institute for Child Development and Education, the University of Amsterdam
http://www.uva.nl/profile/t.d.jorgensen

################################################
# This code can be edited in this window and   #
# submitted to Rweb, or for faster performance #
# and a nicer looking histogram, submit        #
# directly to R.                               #
################################################
a=5
b=10
covajbj=1
vara=1
varb=1
covab=0
varcovajbj=5
rep=20000
conf=95
dvec=rnorm(rep)
avec=dvec*sqrt(vara)+a
bvec=dvec*covab/sqrt(vara)+sqrt(varb)*rnorm(rep,sd=sqrt(1-(covab^2)/(vara*varb)))+b
cvec=rnorm(rep)*sqrt(varcovajbj)+covajbj
ab=avec*bvec+cvec
low=(1-conf/100)/2
upp=((1-conf/100)/2)+(conf/100)
LL=quantile(ab,low)
UL=quantile(ab,upp)
LL4=format(LL,digits=4)
UL4=format(UL,digits=4)
################################################
# The number of columns in the histogram can   #
# be changed by replacing 'FD' below with      #
# an integer value.                            #
################################################
hist(ab,breaks='FD',col='skyblue',xlab=paste(conf,'% Confidence Interval ','LL',LL4,'  UL',UL4), main='Distribution of Indirect Effect')