Performing Bootstrap Using Macros and Foreach Loops

[Roman Le]

Dear Statalists,

I would like to consult with you about using Macros and Foreach Loops to perform Bootstrap for mediation analysis, in STATA 17. My goal is to estimate the indirect effects of X2 on Y through the channel X1. My code is:

encode CompanyName, gen(id)
xtset id Year
capture program drop channel
program define channel, rclass
    xtreg X1 X2 CONTROLS  i.Year, vce(cluster id2)
    scalar a = _b[X2]
    xtreg Y X1 X2 CONTROLS  i.Year, vce(cluster id2)
    scalar b = _b[X1]
    return scalar indirect = a*b
end
bootstrap r(indirect), reps(1000) idcluster(id) cluster(CompanyName) seed(12345): channel

Now, my problem arises when I try to generalize the codes because my goal now is to estimate a set of indirect effects of X2 on each of Y1, Y2, Y3. Therefore, I consider using Macros and Foreach loop. Here is what I code:

capture program drop channel
program define channel, rclass
    xtreg X1 X2 CONTROLS  i.Year, vce(cluster id2)
    scalar a = _b[X2]
    xtreg `Y' X1 X2 CONTROLS  i.Year, vce(cluster id2)
    scalar b = _b[X1]
    return scalar indirect = a*b
end
foreach Y of varlist Y1 Y2 Y3{
bootstrap r(indirect), reps(1000) idcluster(id) cluster(CompanyName) seed(12345): channel}

However, STATA returns as "invalid syntax"

I wish to listen to your expertise. Meanwhile, How does the code look like if X2 becomes a set of variables X21, X22, X23,...?
Thanks so much!

[Andrew Musau]

Probably, some renaming trick would work here, assuming that the referenced variables in the program do not exist in the dataset. You should find a way to save the intermediate results, otherwise they will be overwritten.

encode CompanyName, gen(id)
xtset id Year
capture program drop channel
program define channel, rclass
    xtreg X1 X2 CONTROLS  i.Year, vce(cluster id2)
    scalar a = _b[X2]
    xtreg Y X1 X2 CONTROLS  i.Year, vce(cluster id2)
    scalar b = _b[X1]
    return scalar indirect = a*b
end

forval i=1/3{
    rename (Y`i' X2`i') (Y X2)      
    bootstrap r(indirect), reps(1000) idcluster(id) cluster(CompanyName) seed(12345): channel
    rename (Y X2) (Y`i' X2`i')
}

[Manh Hoang Ba]

The problem has been solved, and I just want to further optimize the code.

• First, I want to provide a varlist for the channel command, so you can specify different dependent variables.
• Then, I removed the vce option in the xtreg command, because we only need bootstrap for the parameter, so adjusting the standard error in the xtreg command would only waste time without any benefit

My recommended code:

encode CompanyName, gen(id)
xtset id Year
capture program drop channel
program define channel, rclass
    syntax varlist(ts fv)
    gettoken Y rest : varlist
    gettoken X1 rest : rest
    gettoken X2 ctrl : rest
    
    tempname a b
    xtreg `X1' `X2' `ctrl'
    scalar `a' = _b[`X2']
    xtreg `Y' `X1' `X2' `ctrl'
    scalar `b' = _b[`X1']
    return scalar indirect = `a'*`b'
end

foreach Y of varlist Y1 Y2 Y3 {
    bootstrap r(indirect), reps(1000) seed(12345) ///
        idcluster(id) cluster(CompanyName) group(id) : ///
        channel `Y' X1 X2 CONTROLS  i.Year
}

[Nick Cox]

To spell out some more details about #1

The program will not work as intended unless you define local macro Y. But that is never done. How did you suppose that Stata would know what it meant?

However, not defining a local macro is not itself an error. Stata ignores references to local macros that don't exist, or more precisely it replaces references to them with empty strings. But then

 
    xtreg `Y' X1 X2 CONTROLS  i.Year, vce(cluster id2)

would be interpreted as

 
    xtreg X1 X2 CONTROLS  i.Year, vce(cluster id2)

which is just a command used earlier in the code.

In short, not defining a local macro is not only not an error, it doesn't in this case imply an error of different kind either.

It follows that your code was in error for a different reason. What was that? Perhaps it was that you referred to CONTROLS. which is not a variable in your dataset.

You have so far as I can tell good positive answers at #2 and #3, but in #1 you didn't specify data that would allow anyone to test the code directly.

For your next question, please read the FAQ Advice first, all the way down to #18.

[Joseph Coveney]

My goal is to estimate . . . a set of indirect effects of X2 on each of Y1, Y2, Y3.

Wouldn't those three outcome variables be correlated? If so, then your estimating each separately with xtreg doesn't capture that key feature of the data-generating process. My understanding is that the use of structural equation modeling (SEM) greatly eases estimating indirect effects; moreover, it readily accommodates multiple (multivariate) outcomes such as you have.

So sem would be an obvious choice, and it allows vce(bootstrap), but it doesn't conveniently allow random effects, i.e., your CompanyName variable. You can reshape wide and fit something akin to a latent growth-curve model, but if the number of years is large, then that approach might prove cumbersome.

On the other hand, gsem does allow specifying random effects directly, and although it doesn't allow vce(bootstrap), you can nevertheless bootstrap all three of your estimates of indirect effects readily and in one line of code—see below. (Begin at the "Begin here" comment; the stuff above is to create a dataset for illustration inasmuch as you neglected to furnish yours, as Nick has pointed out. I use lower-case variable names and shorten the name of the so-called controls variable.)

version 19

clear *

// seedem
set seed 97147652

// Companies
tempname Corr
matrix define `Corr' = J(4, 4, 0.5) + I(4) * 0.5
quietly drawnorm e0 e1 e2 e3, double corr(`Corr') n(250)
generate `c(obs_t)' cid = _n

// Predictors, one endogenous
generate double ctr = runiform(-0.5, 0.5) 
generate double x2 = runiform(-0.5, 0.5)

quietly expand 5
bysort cid: generate byte tim = _n
generate double x1 = rnormal(x2 + e0, 1)

// Three correlated outcomes
forvalues i = 1/3 {
    generate double y`i' = rnormal(x1 + x2 + ctr + e`i', 1)
}

*
* Begin here
*
#delimit ;
bootstrap 
    (y1: indirect = (_b[x1:x2] * _b[y1:x1]))
    (y2: indirect = (_b[x1:x2] * _b[y2:x1]))
    (y3: indirect = (_b[x1:x2] * _b[y3:x1])), 
        cluster(cid) idcluster(aid) reps(40): 
    gsem 
        (x1 <- c.(x2 ctr) i.tim M[aid])
        (y1 y2 y3 <- c.(x1 x2 ctr) i.tim M[aid]);
#delimit cr

exit

I've attached the complete do-file and log file below, but I believe all of these commands will work also with your Stata Release 17.