Monte carlo analysis for moderated mediation, what's wrong for the code?

[Fred Lee]

Refering this post, https://stats.oarc.ucla.edu/stata/fa...irect-effects/
I wrote the code below, what's wrong for the code since it cannnot work well.

sysuse auto,clear

set seed 10000
capture program drop ModMed
program define ModMed, rclass
quietly summarize gear_ratio
   return list
   global m=r(mean)
   global sd=r(sd)
   reg price mpg c.length##c.gear_ratio
   generate IV = rnormal(_b[c.length],_se[c.length])
   generate Mod = rnormal(_b[c.length#c.gear_ratio],_se[c.length#c.gear_ratio])


   reg weight price mpg c.length
   generate Med = rnormal(_b[price],_se[price])

end

simulate indL=(((IV+Mod*($m-$sd)))*(Med)) indM=(((IV+Mod*($m)))*(Med)) indH=(((IV+Mod*($m+$sd)))*(Med)), reps(50) seed(1234): ModMed
ci means indL indM indH

[Fred Lee]

Can anyone help me out, thanks so much!

[Maarten Buis]

The statement that some code "does not work" tells us exactly nothing:you are asking a question here, so it would be surprising if everything worked fine...

You told us what you did, that is a good start. You also need to tell us what Stata showed you, what you expected, and why you think that what Stata showed does not correspond to what you expected.

[Fred Lee]

Thanks Maarten Buis! Here are the results of Stata:

      Command: ModMed
         indL: ((IV+Mod*(3.014864866797988-.4562870967076303)))*(Med)
         indM: ((IV+Mod*(3.014864866797988)))*(Med)
         indH: ((IV+Mod*(3.014864866797988+.4562870967076303)))*(Med)

Simulations (50)
----+--- 1 ---+--- 2 ---+--- 3 ---+--- 4 ---+--- 5
.xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx    50

. ci means indL indM indH

    Variable |        Obs        Mean    Std. err.       [95% conf. interval]
-------------+---------------------------------------------------------------
        indL |          1   -4.128647           .               .           .
        indM |          1   -7.442639           .               .           .
        indH |          1   -10.75663           .               .           .

I expect the Stata can output the confidence interval

Thanks a ton!

[William Lisowski]

The problem is that you expect that for each of the 50 simulations the values of $m and $sd set by ModMed will be substituted into the three expressions on the simulate command line.

That is not what happens.

      Command: ModMed
         indL: ((IV+Mod*(3.014864866797988-.4562870967076303)))*(Med)
         indM: ((IV+Mod*(3.014864866797988)))*(Med)
         indH: ((IV+Mod*(3.014864866797988+.4562870967076303)))*(Med)

tells you that the values of $m (3.014864866797988) and $sd (.4562870967076303) at the time the simulate command was processed by Stata were substituted into the expressions. These values were probably created by an earlier attempt to run the simulation.

Further, on the second and successive simulations, ModMed failed (all the x's) because the three variables IV, Mod, and Med already existed.

Simulations (50)
----+--- 1 ---+--- 2 ---+--- 3 ---+--- 4 ---+--- 5
.xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx    50

Hence there was only a single observation from the simulation, no way to create a confidence interval.

. ci means indL indM indH

    Variable |        Obs        Mean    Std. err.       [95% conf. interval]
-------------+---------------------------------------------------------------
        indL |          1   -4.128647           .               .           .
        indM |          1   -7.442639           .               .           .
        indH |          1   -10.75663           .               .           .

You have correctly defined ModMed as an rclass program, so you need to return r(m) and r(sd) and use those in your expressions.

sysuse auto,clear
set seed 10000

capture program drop ModMed
program define ModMed, rclass
   summarize gear_ratio
   return scalar m=r(mean)
   return scalar sd=r(sd)
   reg price mpg c.length##c.gear_ratio
   capture drop IV Mod Med
   generate IV = rnormal(_b[c.length],_se[c.length])
   generate Mod = rnormal(_b[c.length#c.gear_ratio],_se[c.length#c.gear_ratio])
   reg weight price mpg c.length
   generate Med = rnormal(_b[price],_se[price])
end

simulate indL=(((IV+Mod*(r(m)-r(sd))))*(Med)) indM=(((IV+Mod*(r(m))))*(Med)) indH=(((IV+Mod*(r(m)+r(sd))))*(Med)), reps(50) seed(1234): ModMed
ci means indL indM indH

. simulate indL=(((IV+Mod*(r(m)-r(sd))))*(Med)) indM=(((IV+Mod*(r(m))))*(Med)) indH=(((IV+Mod*(r(m)+r(sd))))*(Med)), reps(50) seed(1234): ModMed

      Command: ModMed
         indL: ((IV+Mod*(r(m)-r(sd))))*(Med)
         indM: ((IV+Mod*(r(m))))*(Med)
         indH: ((IV+Mod*(r(m)+r(sd))))*(Med)

Simulations (50)
----+--- 1 ---+--- 2 ---+--- 3 ---+--- 4 ---+--- 5
..................................................    50

. ci means indL indM indH

    Variable |        Obs        Mean    Std. err.       [95% conf. interval]
-------------+---------------------------------------------------------------
        indL |         50    2.365787    .8399577         .677829    4.053745
        indM |         50    .4486443    .9061497       -1.372332     2.26962
        indH |         50   -1.468499    .9823586       -3.442622    .5056248

.

[Fred Lee]

Thanks William Lisowski, I found just adding the code below will output confidence interval.

capture drop IV Mod Med

There is no difference whether to add

return scalar m=r(mean)
   return scalar sd=r(sd)

since r(m), r(sd) will be replaced to exact values.

[Fred Lee]

By the way, why the resutls are different from yours? William Lisowski

sysuse auto,clear
set seed 10000

capture program drop ModMed
program define ModMed, rclass
   summarize gear_ratio
   return scalar m=r(mean)
   return scalar sd=r(sd)
   reg price mpg c.length##c.gear_ratio
   return scalar IV = rnormal(_b[c.length],_se[c.length])
   return scalar Mod = rnormal(_b[c.length#c.gear_ratio],_se[c.length#c.gear_ratio])
   reg weight price mpg c.length
   return scalar Med = rnormal(_b[price],_se[price])
end

simulate indL=(((r(IV)+r(Mod)*(r(m)-r(sd))))*(r(Med))) indM=(((r(IV)+r(Mod)*(r(m))))*(r(Med))) indH=(((r(IV)+r(Mod)*(r(m)+r(sd))))*(r(Med))), reps(50) seed(1234): ModMed
ci means indL indM indH

      Command: ModMed
         indL: ((r(IV)+r(Mod)*(r(m)-r(sd))))*(r(Med))
         indM: ((r(IV)+r(Mod)*(r(m))))*(r(Med))
         indH: ((r(IV)+r(Mod)*(r(m)+r(sd))))*(r(Med))

Simulations (50)
----+--- 1 ---+--- 2 ---+--- 3 ---+--- 4 ---+--- 5
..................................................    50

. ci means indL indM indH

    Variable |        Obs        Mean    Std. err.       [95% conf. interval]
-------------+---------------------------------------------------------------
        indL |         50     1.93562    .7815475        .3650412    3.506198
        indM |         50    .0573611    .8679277       -1.686805    1.801527
        indH |         50   -1.820897    .9627525       -3.755621    .1138262

[William Lisowski]

The results in post #7 differ from the results in post #5 because your code in post #7 is different than your code in post #1 which I copied in post #5. In particular the code in post #1 apparently generates 74 observations of 3 variables, using 222 random values, while the code in post #7 generates 3 random values.

[Fred Lee]

Thanks William Lisowski, Can you explain how do you know post #5 generates 74 observations, using 222 radom values, I couldn't find such details.
And if I want to examine the moderated mediation, regarding legenth as IV, gear_ratio as moderator, price as mediator and weight as DV, which setting between post #5 and post #7 is correct?
Here is official guide to examine mediation effects using Monte Carlo:https://stats.oarc.ucla.edu/stata/fa...irect-effects/
Thanks a ton!

[William Lisowski]

The autos dataset contains 74 observations.

The code for ModMed in post #5 generates 3 new variables, each of which will have 74 observations.

Each of the new variables generates 1 random normal value.

74 x 3 x 1 = 224.

Regarding moderated mediation, I have no knowledge of that process. Your question suggests that either the code in post #5 is correct or the code in post #7 is correct. Perhaps neither is; I lack the background to make that judgement. And I note that the simulate command accomplishes nothing that couldn't be accomplished with a few, non-looping, lines of Stata code, following the guidance in the document you linked to.

sysuse auto,clear
set seed 10000
set obs 5000

   summarize gear_ratio
   scalar m=r(mean)
   scalar sd=r(sd)
   reg price mpg c.length##c.gear_ratio
   generate IV = rnormal(_b[c.length],_se[c.length])
   generate Mod = rnormal(_b[c.length#c.gear_ratio],_se[c.length#c.gear_ratio])
   regress weight price mpg c.length
   generate Med = rnormal(_b[price],_se[price])

generate indL=(IV+Mod*(m-sd)) * Med 
generate indM=(IV+Mod*(m))    * Med
generate indH=(IV+Mod*(m+sd)) * Med

ci means indL indM indH

. ci means indL indM indH

    Variable |        Obs        Mean    Std. err.       [95% conf. interval]
-------------+---------------------------------------------------------------
        indL |      5,000    2.276629    .0728486        2.133814    2.419445
        indM |      5,000    .4235689    .0780027        .2706494    .5764883
        indH |      5,000   -1.429492    .0842831       -1.594723    -1.26426

.

[Fred Lee]

Thanks William Lisowski, I thought #10 and #5 is not correct for examining the moderated mediation, since setting the observation number I thought might be wrong. Here is the way of bootstraping method to examine moderated mediation.
https://stats.oarc.ucla.edu/stata/fa...tion-in-stata/

According to the guide above, do you thnk #7 is correct?
For #10, the number of Monte carlo iterations is 1,
for #5, the number of Monte carlo iteration is 50 and the number of variables is 74, the random values is 222
for #7 the number of Monte carlo iteration is 50, the random value is 3,
are these judgments correct?

[William Lisowski]

I have no more to add to this discussion. Post #5 corrected the Stata coding mistakes in your post #1; that was what I set out to accomplish. I have neither the ability nor the interest to evaluate methodological correctness.

[Fred Lee]

OK, thank you.