clear set seed 1234 set obs 10000 /* x=septic shock c=confounders m=mediator y=PICU LOS */ gen x=runiform()<.5 gen c=rnormal() gen m=2*x+c+rnormal() gen y=x+.25*m+5*c+rnormal()*2 * Direct effect of x on y, not mediated by m is 1 when c=0 * Indirect effect of x on y, via m is 0.5 (2*.25) * Total effect of x on y is 1+.5=1.5 and there is no interaction between m and x * Baron-Kenny decomposition to estimate the direct and indirect effects of the exposure on the outcome. /* Total effect of x on y, including mediation by m, not controlling for C is 1.6*/ reg y x scalar total_effect_est=_b[x] /*the effect of m on y is 2.1*/ reg y m /* Biased direct effect of x on y, not mediated by m, not controlling for C is -4 */ reg y m x scalar biased_direct_effect_est=_b[x] * Check that you could get the right answer if you controlled for c /* Direct effect of x on y, not mediated by m, controlling for C is 1*/ reg y m x c scalar direct_effect_est=_b[x] dis direct_effect_est /* Direct effect of x on y, not mediated by m, controlling for C is 1*/ dis in red "indirect effect estimated w/o control for c is: " round(total_effect_est-biased_direct_effect_est,.001) /* indirect effect estimated w/o control for c is: 5.637*/ dis in red "indirect effect estimated w/ control for c is: " round(total_effect_est-direct_effect_est,.001) /*indirect effect estimated w/ control for c is: .631*/ dis in red "total effect is: " round(total_effect_est,.001) /*total effect is: 1.64*/ ** 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) /* controlled direct effect= 1.0110056 nde:natural direct effect = 1.0110515 nie:natural indirect effect= .50238486 mte:marginal total effect = 1.5134364 */ * Add some measurement error to see what happens gen m_measured=m+rnormal()*.5 sum * Baron-Kenny decomposition to estimate the direct and indirect effects of the exposure on the outcome with the noisy mediator. /* Total effect of x on y, including mediation by m, not controlling for C is 1.6*/ reg y x /*the effect of m on y is 1.959073 */ reg y m_measured /* Biased direct effect of x on y, not mediated by m, not controlling for C is -3.340545 */ reg y m_measured x /* Direct effect of x on y, not mediated by m, controlling for C is 1.115947 */ reg y m_measured x c ***Estimating the CDE /* 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*/ *(1) make an identifier for your data (in stata, "gen id=_n") preserve *sample 1000 , count 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. expand 3 *(3) for the first "fake" copy of each observation, set x to 0, m to 0 and y to . sort id by id: gen copy=_n replace m=0 if copy==2 replace x=0 if copy==2 replace y=. if copy==2 *(4) for the second "fake" copy of each observation, set x to 1, m to 0 and y to . replace m=0 if copy==3 replace x=1 if copy==3 replace y=. if copy==3 *(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. *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 /*Direct Effect of x on Y = 1.011006*/ *(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 predict cf_y_x0_m0 if copy==2 * What's the average value of y if x is set to 0 and m is set to 0? = -.120016 sum cf_y_x0_m0 if copy==2 scalar mean_cf_y_x0_m0=r(mean) *(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 predict cf_y_x1_m0 if copy==3 * What the average potential outcome for y if x is set to 1 and m is set to 0? = .8909896 sum cf_y_x1_m0 if copy==3 scalar mean_cf_y_x1_m0=r(mean) * or you can just summarize sum /* cf_y_x0_m0 | 10,000 -.120016 5.034749 -17.338 16.76364 cf_y_x1_m0 | 10,000 .8909896 5.034749 -16.32699 17.77464*/ *(8) estimate the controlled direct effect of x on y, setting m to 0 *Total, biased effect = 1.64 reg y x *Estimated direct effect of x, setting m to 0, is: 1.011 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 * the counterfactual value of y setting x to 0 = -.1909555 sum cf_y_x0 scalar mean_cf_y_x0=r(mean) * the counterfactual value of y setting x to 1 = 1.449322 sum cf_y_x1 scalar mean_cf_y_x1=r(mean) *Estimated total effect of x on y is: 1.64 dis "Estimated total effect of x on y is: " round(mean_cf_y_x1-mean_cf_y_x0,.001) *Check it with paramed paramed y ,avar(x) mvar(m) a0(0) a1(1) m(0) yreg(linear) mreg(linear) cvars(c) /* cde:controlled direct effect= 1.0110056 nde:natural direct effect = 1.011025 nie:natural indirect effect= .16823403 mte:marginal total effect = 1.179259 */ /* Bonus hw. */ *(1) Go back to your original data; make an identifier for your data (in stata, "gen id=_n") restore preserve 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. expand 3 sort id by id: gen copy=_n *(3) for the first "fake" copy of each observation, set x to 0 and m to . and y to . replace m=. if copy==2 replace x=0 if copy==2 replace y=. if copy==2 *(4) for the second "fake" copy of each observation, set x to 1, m to . and y to . replace m=. if copy==3 replace x=1 if copy==3 replace y=. if copy==3 *(5) estimate a regression model predicting the mediator as a function of x and c, using only the real observations gen mx=m*x if copy==1 reg m x c if copy==1 *(6) predict cf_m_x0 in the first fake copy predict cf_m_x0 if copy==2 *(7) predict cf_m_x1 in the second fake copy predict cf_m_x1 if copy==3 *(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), using only the real observation replace mx=m*x reg y x m mx c if copy==1 *(9) in the first fake copy, set m to cf_m_x0 and set x to 1 replace m=cf_m_x0 if copy==2 replace x=1 if copy==2 *(10) in the second fake copy, set m to cf_m_x1 replace m=cf_m_x1 if copy==3 *(11) in the first fake copy, predict y based on the regression model in (8), to estimate cf_y_x1_cf_m_x0 (the counterfactual value of y setting x to 1 and m to the value it would take if x were set to 0)= .8831828 predict cf_y_x1_cf_m_x0 if copy==2 *(12) in the second fake copy, predict y based on the regression model in (8), to estimate cf_y_x1_cf_m_x1 (the counterfactual value of y setting x to 1 and m to the value it would take if x were set to 1)= 1.385568 predict cf_y_x1_cf_m_x1 if copy==3 *(13) the natural indirect effect of x on y, mediated by m is: .502*/ 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 /*