*Inverse Probability Weighting
clear all 
set seed 22019
set obs 10000

*********************************************************************
*Let's review the case of the conditionally randomized experiment ***
*Let's compare standardization and IPWs for getting the correct *****
*Causal estimate (risk ratio or risk difference) ********************
*********************************************************************
*gen unique_id = _n
gen female = runiform()<0.60
gen treat = runiform()<0.60 if female==1
replace treat = runiform()<0.40 if female==0
tab female treat, row
gen y1 = int(runiform()<0.45 - 0.15*treat) if female==1 /*treatment protective for women*/
replace y1 = int(runiform()<0.45 + 0.15*treat) if female==0 /*harms men*/

*counterfactual outcome values
gen y1_fem1 = int(runiform()<0.45 - 0.15*1) if female==1
gen y1_fem0 = int(runiform()<0.45 - 0.15*0) if female==1

*causal risk ratio for women
di 0.30/0.45

*causal risk difference for women
di 0.30 - 0.45

gen y1_male1 = int(runiform()<0.45 + 0.15*1) if female==0
gen y1_male0 = int(runiform()<0.45 + 0.15*0) if female==0

*causal risk ratio for men 
di 0.60/0.45

*causal risk difference for men 
di 0.60-0.45

*What happens if we tried to estimate the causal risk ratios via regression
*for women
glm y1 treat if female==1, family(binomial) link(log) eform

*for men
glm y1 treat if female==0, family(binomial) link(log) eform

*overall
glm y1 treat, family(binomial) link(log) eform

*adding female in to the regression model
glm y1 treat female, family(binomial) link(log) eform

*Let's use standardization to get overall causal risk ratio
*Without standardization
di(0.30 + 0.60)/(0.45 + 0.45)

*With standardization
*Denominator = Risk of Y=1|A=0
di(0.45*0.6038) + (0.45*.3962)

*Numerator = Risk of Y=1|A=1
di(0.30*0.6038) + (0.60*.3962)

*overall causal risk ratio should be 
di 0.41 / 0.45

*Standardization via regression
logistic y1 i.treat##i.female
margins treat, post
nlcom(risk_ratio: _b[1.treat]/_b[0.treat])

*Inverse probability of treatment weighting
logistic treat female
predict treat_p

twoway (histogram treat_p if treat==1, start(0) color(green)) ///
       (histogram treat_p if treat==0, start(0) ///
	   fcolor(none) lcolor(black)), legend(order(1 "Treated" 2 "Not Treated" ))

gen wt1 = 1/treat_p if treat==1
replace wt1 = 1/1-treat_p if treat==0
sum wt1, detail 
*pay attention to the ratio of the smallest to the largest!

glm y1 treat [pw = wt1], family(binomial) link(log) eform 
/*robust standard errors were estimated above, but below are two ways to add to double-check*/
glm y1 treat [pw = wt1], family(binomial) link(log) eform cluster(unique_id) 
glm y1 treat [pw = wt1], family(binomial) link(log) eform robust