bootstrapping standard errors in simulation

[Jessica Lum]

Hi all,

I'd like to run a simulation of bootstrapped standard errors estimated from inverse probability of treatment weighting. The code below is just an example of how I've organized the simulation:

            
clear all
set more off
set seed 123456

capture program drop test01222020
program define test01222020, rclass

version 14.2 
syntax [, nobs(integer 1000)]

preserve
set obs `nobs'
g x1 = rnormal()
g x2 = rnormal() 
g x3 = rnormal()

g xb2 = .3*(x1 + x2 + x3)
g xb3 = -.6*(x1 + x2 + x3) 

g p1 = 1 / (1 + exp(xb2) + exp(xb3))
g p2 = exp(xb2) / (1 + exp(xb2) + exp(xb3))
g p3 = exp(xb3) / (1 + exp(xb2) + exp(xb3))

g prob = runiform()
g treat = 3 
replace treat = 2 if prob < p1 + p2 
replace treat = 1 if prob < p1 

g trueY1 = -4 + .2*(x1 + x2 + x3) 
g trueY2 = -3 + .2*(x1 + x2 + x3) 
g trueY3 = -1 + .2*(x1 + x2 + x3) 

g trueY = cond(treat == 1, trueY1, cond(treat == 2, trueY2, cond(treat == 3, trueY3, .)))

mlogit treat x1 x2 x3
predict ps1 ps2 ps3

g weight = .
replace weight = 1/ps1 if treat == 1 
replace weight = 1/ps2 if treat == 2
replace weight = 1/ps3 if treat == 3

regress trueY ib2.treat [w=weight] if inlist(treat, 1, 2)
return scalar iptw_ATE1vs2_b = _b[1.treat]
return scalar iptw_ATE1vs2_se = _se[1.treat] 

simulate iptw_ATE12b = r(iptw_ATE12b), reps(100): myboot
sum iptw_ATE12b
return scalar iptw_ATE12bstrap_se = r(sd)

restore 
end 


capture program drop myboot 
program define myboot, rclass
bsample 
regress trueY ib2.treat [pw=weight] if inlist(treat, 1, 2)
return scalar iptw_ATE12b = _b[1.treat]
end 

simulate ///
iptw_ATE1vs2_b = r(iptw_ATE1vs2_b) ///
iptw_ATE1vs2_se = r(iptw_ATE1vs2_se) ///
iptw_ATE12bstrap_se = r(iptw_ATE12bstrap_se) ///
, reps(5) seed(123456): test01222020, nobs(1000)

Could someone confirm whether the bootstrapped SEs are correctly estimated or not? If not, I'd appreciate any assistance on this.

Thanks!
Jessica

[FernandoRios]

Hi Jessica,
i dont think your Bootstrapped SE are correct. First, you need a preserve and restore statement in the bootstrap program
Second, you should also estimate the IPW within the bootstrap program.

capture program drop myboot
program define myboot, rclass
  preserve
  bsample
  mlogit treat x1 x2 x3
  predict ps1 ps2 ps3
  g weight2 = .
  replace weight2 = 1/ps1 if treat == 1
  replace weight2 = 1/ps2 if treat == 2
  replace weight2 = 1/ps3 if treat == 3

  regress trueY ib2.treat [pw=weight2] if inlist(treat, 1, 2)  
  return scalar iptw_ATE12b = _b[1.treat]
restore
end

HTH
Fernando

[Jessica Lum]

Thanks for the input, Fernando!