clear
set seed 9055797
set obs 100000
gen x=runiform()<.5
gen c=rnormal()
gen m=1.5*x+c+rnormal()
gen y=x+.5*m+2*c+rnormal()*2

* Direct effect of x on y, not mediated by m is 1
* Indirect effect of x on y, via m is 0.75 (1.5*.5)
* Total effect of x on y is 1+.75=1.75 and there is no interaction between m and x

* Use BK
reg y x
scalar total_effect_est=_b[x]
reg y m
reg y m x
scalar biased_direct_effect_est=_b[x]

* Check that you could get the right answer if you controlled for c
reg y m x c
scalar direct_effect_est=_b[x]
dis direct_effect_est
dis in red "indirect effect estimated w/o control for c is: " round(total_effect_est-biased_direct_effect_est,.001)
dis in red "indirect effect estimated w/ control for c is: " round(total_effect_est-direct_effect_est,.001)
dis in red "total effect is: " round(total_effect_est,.001)

** Check it w/ paramed (you may need to install)
paramed y ,avar(x) mvar(m) a0(0) a1(1) m(0) yreg(linear) mreg(linear) cvars(c)

* Add some measurement error to see what happens
gen m_measured=m+rnormal()*.5
sum

* Redo above
* Use BK
reg y x
reg y m_measured
reg y m_measured x
reg y m_measured x c




/*
(1) make an identifier for your data (in stata, "gen id=_n")
(2) make 3 copies of every observation (in stata, use "expand 3"); now you have 2 fake copies of each observation and one real copy.
(3) for the first  "fake" copy of each observation, set x to 0, m to 0 and y to . 
(4) for the second "fake" copy of each observation, set x to 1, m to 0 and y to .
(5) estimate a regression model predicting the outcome as a function of exposure, mediator, the interaction of the exposure and mediator, 
     and the mediator-outcome confounder (C), using only the real observations.
(6) for the first fake copy of each observation, use the predict statement to predict the counterfactual value of y 
    setting x to 0 and m to 0
(7) for the second fake copy of each observation, use the predict statement to predict the counterfactual value of y
    setting x to 1 and m to 0
(8) estimate the controlled direct effect of x on y, setting m to 0

My naming convention is that cf_y_x0_m0 is the counterfactual value of y setting x to 0 and m to 0.
and cf_y_x0_cf_m_x1 is the counterfactual vlaue of y setting x to 0 and m to the value it would take if x were set to 1
*/
preserve
*sample  1000 , count
gen id=_n

expand 3
sort id
by id: gen copy=_n
replace m=0 if copy==2
replace x=0 if copy==2
replace y=. if copy==2
replace m=0 if copy==3
replace x=1 if copy==3
replace y=. if copy==3
*reg m x c if copy==1
*predict m if copy>1
gen mx=m*x
reg y x c m mx if copy==1
predict cf_y_x0_m0 if copy==2
predict cf_y_x1_m0 if copy==3
* What's the average value of y if x is set to 0 and m is set to 0?
sum cf_y_x0_m0 if copy==2
scalar mean_cf_y_x0_m0=r(mean)
* What the average potential outcome for y if x is set to 1 and m is set to 0?
sum cf_y_x1_m0 if copy==3
scalar mean_cf_y_x1_m0=r(mean)
* or you can just summarize
sum
reg y x 
dis "Estimated direct effect of x, setting m to 0, is: " round(mean_cf_y_x1_m0-mean_cf_y_x0_m0,.001)

predict cf_y_x0 if copy==2
predict cf_y_x1 if copy==3
sum cf_y_x0
scalar mean_cf_y_x0=r(mean)
sum cf_y_x1
scalar mean_cf_y_x1=r(mean)
dis "Estimated total effect of x on y is: " round(mean_cf_y_x1-mean_cf_y_x0,.001)

* compare to paramed

/* Bonus hw if you're having fun. Go back to your original data
(1) make an identifier for your data (in stata, "gen id=_n")
(2) make 3 copies of every observation (in stata, use "expand 3"); now you have 2 fake copies of each observation and one real copy.
(3) for the first "fake" copy of each observation, set x to 0 and m to . and y to . 
(4) for the second "fake" copy of each observation, set x to 1, m to . and y to .
(5) estimate a regression model predicting the mediator as a function of x and c, using only the real observations
(6) predict cf_m_x0 in the first fake copy
(7) predict cf_m_x1 in the second fake copy
(8) estimate a regression model predicting the outcome as a function of exposure, mediator, the interaction of the exposure and mediator, 
     and the mediator-outcome confounder (C), usin gonly the real observation
(9) in the first fake copy, set m to cf_m_x0 and set x to 1
(10) in the second fake copy, set m to cf_m_x1 
(11) in the first fake copy,  predict y based on the regression model in (8), to estimate cf_y_x1_cf_m_x0
(12) in the second fake copy, predict y based on the regression model in (8), to estimate cf_y_x1_cf_m_x1
(13) estimate the natural indirect effect of x on y, mediated by m
*/
restore
preserve
gen id=_n
expand 3
sort id
by id: gen copy=_n
replace m=. if copy==2
replace x=0 if copy==2
replace y=. if copy==2
replace m=. if copy==3
replace x=1 if copy==3
replace y=. if copy==3

gen mx=m*x if copy==1
reg m x c if copy==1
predict cf_m_x0 if copy==2
predict cf_m_x1 if copy==3

replace m=cf_m_x0 if copy==2
replace m=cf_m_x1 if copy==3
replace x=1 if copy==2
replace mx=m*x

reg y x m mx c if copy==1
predict cf_y_x1_cf_m_x0 if copy==2
predict cf_y_x1_cf_m_x1 if copy==3
sum cf_y_x1_cf_m_x0 if copy==2
scalar mean_cf_y_x1_cf_m_x0=r(mean)
sum cf_y_x1_cf_m_x1 if copy==3
scalar mean_cf_y_x1_cf_m_x1=r(mean)
dis "the natural indirect effect of x on y, mediated by m is: " round(mean_cf_y_x1_cf_m_x1-mean_cf_y_x1_cf_m_x0,.001)

sort copy
by copy: sum