Mediation Analysis in a Multilevel Model: How to bootstrap indirect effects after gsem

[Sule Yaylaci]

Dear Statalist users,
I am using Stata 14 SE. My data is composed of two-wave observations at the district-level and districts are nested in provinces. Data example is below. I would like to see if the independent variables (X1 X2 X3) have indirect effects on Y via Z.
Given the multilevel nature of data I am using gsem and am trying to do mediation analysis.

I use the command:

gsem (Y<- X1 i.X2 X3 Z P[province] C[province>district_no]) ///
(Z<-X1 i.X2 X3 P[province] C[province>district_no]), nocapslatent latent (P C)

Because "estat teffects" is not supported after gsem, I manually calculated the indirect effects:

nlcom _b[Z:X1]*_b[Y:Z]
nlcom _b[Z:1.X2]*_b[Y:Z]
nlcom _b[Z:2.X2]*_b[Y:Z]
nlcom _b[Z:X3]*_b[Y:Z]

In the mediation literature, there is an argument for using bootstrapped confidence intervals when calculating indirect effects (e.g. Preacher and Hayes 2004)
I saw guidelines about how to do that after sem :

mat bi = r(indirect)
mat bd = r(direct)
mat bt = r(total)
return scalar indir = el(bi,1,3)
return scalar direct = el(bd,1,3)
return scalar total = el(bt,1,3)

But I cannot figure how to apply this to my nlcom results i.e. how to make a matrix after manually calculated indirect effects to run bootstrapped confidence intervals on the beta coefficients of the indirect effects.
I tried running the gsem with bootstrap prefix but after 3 hours it was still running. I read a bit about 'boottest' but I could not figure how to do use it for my purposes.

Any help with getting bootstrapped CIs for indirect effects after a two-level gsem model would be much appreciated.

* Example generated by -dataex-. To install: ssc install dataex
clear
input int district_no byte(wave province) double (Y Z) float X1 byte X2 float X3
 7 1  1 .23223635003739715  .869795996674554  0 0 0
 4 1  1  .4701195219123506 .8039459375485475  0 0 0
 8 1  1 .42910860429108605 .8222902932232359  0 0 0
 2 1  1 .25843545684528696 .8517262864682598  0 0 0
 3 1  1  .2314806361495061 .8671302188663923  0 0 0
 1 1  1  .5008916001877053 .8836196467457053  0 0 0
 9 1  1  .4536984981126014 .8757469606429013  0 0 0
11 1  1    .36644963615473 .8615101724805369  0 0 0
15 1  1  .3377397403399747 .8314157170778479  0 0 0
 6 1  1  .5529871977240398 .9015807301467821  0 0 0
10 1  1  .4836112708453134 .7833479404031551  0 0 0
 5 1  1 .45914704629798625 .8386024120303737  0 0 0
13 1  1 .35363741339491916 .8177992041628406  0 0 0
14 1  1 .32956786802940646 .8485232696897375  0 0 0
12 1  1  .2582014753593243 .8491196205853588  0 0 0
 8 2  1  .5013144023806029 .8456993069130976 13 1 6
 2 2  1 .33096601673721704 .8558748788826122 13 2 6
 5 2  1  .5308816595945309 .8622581288649511 13 0 6
 3 2  1 .28064805600966664 .8749831214675643 13 2 6
12 2  1  .3228475641790513 .8528333602230218 13 2 6
 7 2  1  .2716732739920712 .8649493081680801 13 0 6
15 2  1 .42843334243252795 .8406463605192804 13 1 6
 9 2  1  .5270582609388699 .8802646998000965 13 0 6
 1 2  1  .5701738334858188  .891889725917615 13 0 6
11 2  1  .4484892121448794 .8800895139308706 13 1 6
 4 2  1  .5468765275197967  .825685903500473 13 0 6
 6 2  1  .5972691441441441 .9105439383410197 13 0 6
10 2  1  .5462431001464458 .8043728423475259 13 1 6
13 2  1  .4731265652090156 .8134437035333794 13 0 6
14 2  1 .37584912406149446 .8606675244077803 13 0 6
16 1  2  .5449415852219232 .8457498530453939  0 2 0
17 1  2   .524244480400856 .7929118002416432  0 0 0
24 1  2  .5623674911660778 .8395382395382396  0 0 0
22 1  2  .7121102248005802 .7838338895068595  0 0 0
18 1  2  .5755950385517935 .8811443932411674  0 0 0
21 1  2   .685989894350023 .7928177975148384  0 0 0
20 1  2  .4849688681767829 .8204195205479452  0 2 0
23 1  2  .9071300179748353 .8856816985436041  0 0 0
19 1  2  .7091292483254775 .6773241515002459  0 0 0
23 2  2  .9593705293276109 .9326021581461171  2 0 2
24 2  2  .6620594333102972 .8542088516054382  2 2 2
20 2  2  .5811804708578187 .8416793893129771  2 0 2
21 2  2  .8277890608586036  .828457731311777  2 2 2
22 2  2  .8619329388560157  .872473077649726  2 0 2
17 2  2  .6770883478172743 .8163368642780467  2 0 2
19 2  2    .85121412803532 .7589862514493954  2 0 2
18 2  2  .6560495938435229 .9165220744168112  2 0 2
16 2  2   .641944955764613 .8683048852266039  2 0 2
42 1  3 .24974731232197003 .8973921874433282  0 0 0
29 1  3  .4544952285283777 .8783187717363644  0 0 0
27 1  3  .5145569620253164  .868237347294939  0 0 0
30 1  3  .5304798962386511 .9205705009276438  0 0 0
36 1  3  .5650262617035853 .9233479726279236  0 0 0
25 1  3  .5348067182412929 .9103541429696387  0 0 0
28 1  3  .6256125821524903 .8756493401735875  0 0 0
26 1  3  .4125722543352601 .8892910634048926  0 0 0
40 1  3   .569474921630094 .8865721434528774  0 0 0
38 1  3 .43885714285714283 .8626991565135895  0 0 0
31 1  3 .43370756482224004 .9102380952380953  0 0 0
35 1  3   .559327566508895 .8404325464855598  0 0 0
39 1  3  .5349692529496572 .8505993873465352  0 0 0
34 1  3  .4162415833503367 .8761123713139068  0 0 0
37 1  3  .6467490520994242 .9294468787705594  0 0 0
33 1  3  .4869785664899747  .778960223307746  0 0 0
32 1  3  .4417435328386157 .8575885377549252  0 0 0
41 1  3  .5593326906149139 .8820998278829604  0 0 0
33 2  3  .5746084480303749 .7624526498389209  2 0 0
27 2  3    .62026913372582 .8601643069393463  2 0 0
34 2  3 .47278770253427505 .8665964542741794  2 0 0
35 2  3   .660734327400994 .8594414893617022  2 0 0
38 2  3  .4796839729119639 .8629751290473956  2 0 0
36 2  3  .7066111111111111 .9333808336302102  2 0 0
26 2  3  .4892944388561575  .900830606594513  2 0 0
42 2  3 .30645011600928074 .9088443737344518  2 1 0
37 2  3  .7727925586485193 .9328652917946467  2 0 0
41 2  3  .6688803780964798 .8722279220266751  2 0 0
25 2  3  .6653963139734789  .915298976671581  2 1 0
40 2  3  .6965041965041965 .8881137465949106  2 0 0
39 2  3  .6348095224320963 .8626212058616248  2 0 0
31 2  3  .5171763437963087 .9005827090022595  2 0 0
30 2  3  .6460984702403908 .9102711397058824  2 0 0
32 2  3  .5016402405686168 .8635209235209235  2 1 0
28 2  3  .7158580413297394 .8784253184098804  2 0 0
29 2  3  .5412363492612542 .8783167145512929  2 0 0
52 1 68  .5795023847696735 .8339957901642929  0 0 0
55 1 68  .5809735921094495 .7684285375681441  0 0 0
57 1  5  .5786761791518034  .702775532201563  0 0 0
51 1 68   .640272373540856 .8042553191489362  0 0 0
54 1 68              .5541  .774018944519621  0 0 0
53 1 68  .6289149686802505  .869839519784393  0 0 0
56 1 68  .6004657351962741 .7849462365591398  0 0 0
56 2 68    .71238570241064 .7945223149023992  2 0 0
57 2  5  .6419597989949749 .6763807937829587  1 0 3
52 2 68  .7150396119644301 .8562599887634934  2 1 0
51 2 68   .761794723666474 .8159802560542958  2 0 0
55 2 68  .7036920659858602 .7909467023606717  2 0 0
54 2 68  .6571784550507955 .8054483541430193  2 1 0
53 2 68   .777793237790455 .8943317859760567  2 0 0
end

[Sule Yaylaci]

I actually found a solution, thanks to the help of a Stata mathematician, and want to share it here:

capture program drop
boot1 program boot1, rclass
gsem (Y<- X1 i.X2 X3 Z P[province] C[province>district_no]) ///
(Z<-X1 i.X2 X3 P[province] C[province>district_no]), nocapslatent latent (P C)
return scalar X1_indirect=(_b[Z:X1]*_b[Y:Z])
return scalar X21_indirect=(_b[Z:1.X2]*_b[Y:Z])
return scalar X22_indirect=(_b[Z:2.X2]*_b[Y:Z])
return scalar X3_indirect=(_b[Z:X3]*_b[Y:Z])
end
bootstrap r(X1_indirect) r(X21_indirect) r(X22_indirect) r(X3_indirect), seed(12345) reps(200): boot1 estat boot, bc percentile

[Chris Martin]

I found a typo in the code above (line break at wrong place in first line). In case it's useful to anyone else, here is the fixed version:

capture program drop boot1
program boot1, rclass
gsem (Y<- X1 i.X2 X3 Z P[province] C[province>district_no]) ///
(Z<-X1 i.X2 X3 P[province] C[province>district_no]), nocapslatent latent (P C)
return scalar X1_indirect=(_b[Z:X1]*_b[Y:Z])
return scalar X21_indirect=(_b[Z:1.X2]*_b[Y:Z])
return scalar X22_indirect=(_b[Z:2.X2]*_b[Y:Z])
return scalar X3_indirect=(_b[Z:X3]*_b[Y:Z])
end
bootstrap r(X1_indirect) r(X21_indirect) r(X22_indirect) r(X3_indirect), seed(12345) reps(200): boot1 estat boot, bc percentile

[gi peters]

Hello Sule Yaylaci and Chris Martin, thanks for posting this solution. I want to make sure that I understand the code properly because Im also trying to use Hayes' approach for multiple mediation but using STATA. So if I am not using a nested data and would like to derive the bootstrapped CI of a multiple mediation model, should my code be

   
 capture program drop boot1 program boot1, rclass gsem (x m1 m2->y)(x->m1)(x->m2) return scalar indirectm1=(_b[m1:x]*_b[y:m1]) return scalar indirectm2=(_b[m2:x]*_b[y:m2]) return scalar totindirect=(_b[m1:x]*_b[y:m1])+(_b[m2:x]*_b[y:m2]) return scalar toteffect=_b[y:x]+((_b[m1:x]*_b[y:m1])+(_b[m2:x]*_b[y:m2])) end bootstrap r(indirectm1) r(indirectm2) r(totindirect) r(toteffect), seed(12345) reps(500): boot1 estat boot, bc percentile

where m1 = mediator 1 and m2=mediator 2
indirectm1 = indirect effect through m1
indirectm2 = indirect effect through m2
totindirect = total indirect effect
toteffect = total effect of x on y

[Mike Zyphur]

Hi Sule
Hmm, you should be very careful making inference with the parametric bootstrapping approach using clustered/nested data. The problem is that the data involve two sampling schemes (at two levels), but the bootstrapping won't be sensitive to this if it is only resampling the L1 observations. How should resampling be done in this case, using L1, L2, or a combination (L2 then L1)? To be correct, the resampling should match the actual sampling design. This is a topic in statistics literature that has not been fully resolved to my knowledge. Whatever the case, a Monte Carlo (parametric) bootstrap or one that relies on Bayes credible intervals would be preferred in this case.

In case it's relevant, Instats is offering a series of seminars on Stata's sem command. All seminars are combined in a Structured Course in case it's of interest.

Best wishes
Mike