/////generating random unmeasured confounder with the values a touch make plausible//////

drawnorm U_norm, n(5000) means(0) sds(4) clear

generate U=U_norm

////generating random error term for exposure X, outcome Y,////
generate error_Y = rnormal(0,1)
generate error_X = rnormal(0,1)

////generating X and Y (dependent on X and U)
generate X = 12 + U + error_X*2
generate Y = 80 + U - X*3.0 + error_Y*1.5
generate Y_ind = 80 + U + error_Y*1.5

////generating random bilevel indicator variables for S2 dependent on U/////
generate S2 =1
generate S2_90=1
generate S2_runi= runiform()
generate S2_temp = S2_runi + U*0.10

sort S2_temp
replace S2 = 0 in 2501/5000
replace S2_90 = 0 in 4501/5000

///generating random bilevel indicator variable for S1 dependent on X and Y////
generate S1 = 1
generate S1_90 = 1

generate S1_runi= runiform()
generate S1_temp = S1_runi + 3*X + 4*Y
generate S1_temp_indep = S1_runi + 3*X + 4*Y_ind

sort S1_temp

replace S1 = 0 in 2501/5000
replace S1_90 = 0 in 4501/5000

generate S1_indep = 1
generate S1_90_indep = 1

sort S1_temp_indep

replace S1_indep = 0 in 2501/5000
replace S1_90_indep = 0 in 4501/5000


///running regression models to confirm that the sample reasonably approximates my input parameters and IT DOES!!///
regress Y X

regress Y X U

generate S_total_5050 = S1+S2
regress Y X if S_total_5050==2

generate S_total_5090 = S1+S2_90
regress Y X if S_total_5090==2

generate S_total_9050 = S1_90+S2
regress Y X if S_total_9050==2

generate S_total_9090 = S1_90+S2_90
regress Y X if S_total_9090==2

generate S_total_5050_ind = S1_indep+S2
regress Y_ind X if S_total_5050_ind==2

generate S_total_5090_ind = S1_indep+S2_90
regress Y_ind X if S_total_5090_ind==2

generate S_total_9050_ind = S1_90_indep+S2
regress Y_ind X if S_total_9050_ind==2

generate S_total_9090_ind = S1_90_indep+S2_90
regress Y_ind X if S_total_9090_ind==2


/////running bootstrap option in stata for repeated sample size estimation/////

bootstrap, reps(1000) size(500): regress Y_ind X if S_total_5050_ind==2
bootstrap, reps(1000) size(500): regress Y_ind X if S_total_5090_ind==2
bootstrap, reps(1000) size(500): regress Y_ind X if S_total_9050_ind==2
bootstrap, reps(1000) size(500): regress Y_ind X if S_total_9090_ind==2