How to estimate direct/indirect effects by bootstrap with -mi- in a program?

[Senor Massao]

Hi,
I am trying to run the following program to estimate direct and indirect effects using the approach proposed by Eric Tchetgen.
Ref: Tchetgen Tchetgen EJ. Inverse odds ratio-weighted estimation for causal mediation analysis. Stat Med. 2013;32(26):4567–4580.

I first imputed data using –mi-. Than here is the stata commands I am using to estimate total, direct, and indirect effects;

program myprog , rclass
capture drop lo pp io wt

mi estimate, saving(miest): logit x m1 m2 c1 c2

mi predict lo using miest
mi xeq: gen pp = exp(lo)/(1+exp(lo))
mi xeq: gen io = ((1-pp)/pp)

mi xeq: gen wt = 1 if x==0
mi xeq: replace wt = io if x==1

mi estimate, saving(miest) eform post: glm y x c1 c2, fam(poisson) link(log) vce(robust)

matrix bb_total= e(b)
scalar b_total=(bb_total[1,1])
return scalar b_total=bb_total[1,1]

mi estimate, saving(miest) eform post: glm y x c1 c2 [pweight=wt] , fam(poisson) link(log) vce(robust)
matrix bb_direct = e(b)
scalar b_direct=(bb_direct[1,1])
return scalar b_direct=bb_direct[1,1]


return scalar b_indirect = b_total-b_direct

end


bootstrap exp(r(b_indirect)) exp(r(b_direct)) exp(r(b_total)), seed(32222) reps(50): myprog



The first error I get is:
-myprog already defined-
I change the program name to something else, but it would be nice not having to do that.
The second problem (after giving a new name to program) is that the total observations in the results indicate that the analysis does not use imputed datasets. So the bootstrapping is only giving results (estimates for b_indirect, b_direct, and b_total) for the complete cases.
The third issue is that the confidence intervals do not correspond with the p value. The CI include 1, but the p value is <0.01.

I'm not a programming expert so any feedback is highly appreciated.

Thankfully,
Massao

[Clyde Schechter]

First a disclaimer: I know nothing about the Tchetgen reference. But there are some programming errors here, independent of that.

1. Once you define a program, it stays defined and in memory until it is explicitly removed. I'm betting you didn't get the -myprog already defined- error the very first time you tried this. But you get it every time thereafter. To drop a program you either have to close Stata and start over again, issue a -program drop myprog- command, or issue a -clear*- command. It's generally good practice to put clear* before the "business" part of any Stata do file: that way you can be sure that your do file doesn't depend on something that happened to be lying around once, but won't be when you try again. It also assures that if your do-file contains a -program define- command,that has been revised since last use, that he old version of that program will not be lingering in memory and blocking you. When you are developing a program within a do-file, a simple way of avoiding this problem is to put -capture drop myprog- before the -program define myprog- statement.

2. When you run your -mi estimate- commands, the default is for e(b) to contain the results of the estimation in the last data run, not the results for the combined estimation with Rubin's rules. You will find the Rubin's rules combined estimates in e(b_mi) instead.

3. The p-values you get from -bootstrap- test the null hypothesis that the parameter being estimated is zero, not one. -bootstrap- does not understand that 1 is a more reasonable null value for an exponentiated expression. So you basically need to ignore the p-values here. If you really need pvalues, then run the -bootstrap- command with the -saving- option specified. Then in the saved bootstrap dataset, you can -ttest var = 1- to get pvalues for each of the variables therein.

[Clyde Schechter]

By the way, I am somewhat surprised you even got this to run as far as it did. Each of your -mi estimate- commands specifies -saving(miest)-. But when it reaches the second such command, I would expect Stata to break and complain that miest.ster already exists since you didn't specify the -replace- suboption.

Also, I'm really not sure if this can work properly even with all of these things corrected. I don't know how things play out when you try to bootstrap with a multiply imputed data set. It may depend on how you -mi set- the style. If your data is -mi set long-, I don't think -bootstrap- knows to pick corresponding observations in each imputation, and if it doesn't, I don't know what that implies for the distributions of the bootstrapped results. If it is -mi set wide- I imagine things might work out. (Of course, you can always -mi convert wide- to be sure.)

[Senor Massao]

Hi,
Thank you for the detailed reply :-)