/* 
    Code for matching effects ra, ipw, ipwra, and aipw
	
	Code could be more efficient but efficient code is harder to read 

	Marcelo Coca Perraillon - mcoca@uchicago.edu 
*/ 

*cd ""

log using "matching teffects code.txt", replace text

webuse set "https://perraillon.com/s/"
webuse "help_1_stata12.dta", clear

/// -- Recode vars, etc
rename i1 ndrinks
rename homeless intervention
label var intervention "1 if received intervention"
recode intervention (0 =1)(1=0)
* keep non-missing for all variables 
qui reg pcs intervention age female ndrinks drugrisk
keep if e(sample)

/// --- teffect ra -- regression adjustment following definition of causal effects
* See Wooldridge (2010) Chapter 21

qui reg pcs age female ndrinks drugrisk if intervention == 1
predict double yhat_t 
qui reg pcs age female ndrinks drugrisk if intervention == 0
predict double yhat_c 

sum yhat_t
local pom_t = r(mean)
sum yhat_c
local pom_c = r(mean)
di `pom_t' - `pom_c'

* same as teffects RA
teffects ra (pcs age female ndrinks drugrisk) (intervention)

* same as fully interacted model, averaging over other vars 
qui reg pcs i.intervention##(c.age i.female c.ndrinks c.drugrisk)
margins, dydx(intervention)

* ATET 
qui reg pcs age female ndrinks drugrisk if intervention == 0
predict double yhat_atet if intervention ==1

* compare POMs
sum pcs if intervention ==1
local pom_t = r(mean)
sum yhat_atet
local pom_c = r(mean)
di `pom_t' - `pom_c'

* same as teffects ra with ATET option
teffects ra (pcs age female ndrinks drugrisk) (intervention), atet 


** Logistic
preserve
 
use https://www.stata-press.com/data/r16/cattaneo2, clear

logit lbweight  mage medu mrace if mbsmoke ==1 
predict double phat_t 
logit lbweight  mage medu mrace if mbsmoke ==0 
predict double phat_c 

sum phat_t
local pom_t = r(mean)
sum phat_c
local pom_c = r(mean)
di `pom_t' - `pom_c'

teffects ra (lbweight  mage medu mrace, logit) ///
             (mbsmoke)

restore 

/// --- teffects IPW
capture drop yhat_t yhat_c

* Propensity score 
qui logit intervention ndrinks age female drugrisk
predict pscore if e(sample)

* ipw weights
gen     ipw = 1/pscore if intervention == 1
replace ipw = 1/(1-pscore) if intervention ==0

* Regression weighted by ipw, no covariates
reg pcs i.intervention [pweight= ipw], robust

* same teffects ipw
teffects ipw (pcs) (intervention age female ndrinks drugrisk) 
 
/// --- teffects IPWRA 
capture drop pom_t pom_c
reg pcs intervention age female ndrinks drugrisk [pweight= ipw] if intervention ==1
predict double pom_t
reg pcs intervention age female ndrinks drugrisk [pweight= ipw] if intervention ==0
predict double pom_c
mean pom_c pom_t

* Matches teffects ipwra
teffects ipwra (pcs age female ndrinks drugrisk) ///
             (intervention ndrinks  age female drugrisk), pom

/// --- teffects AIPW
* Need two separate weights
gen double ipw0 = 0.intervention/(1-pscore)
gen double ipw1 = 1.intervention/pscore 

* No weights in outcome model 
qui reg pcs age female ndrinks drugrisk if intervention == 1
predict double pom1
* Use weight to correct prediction (note observed outcome there)
replace pom1 = pom1 + ipw1*(pcs - pom1)

qui reg pcs age female ndrinks drugrisk if intervention == 0
predict double pom0
replace pom0 = pom0 + ipw0*(pcs - pom0)

sum pom0 pom1

* matches teffects aipw
teffects aipw (pcs age female ndrinks drugrisk) ///
    (intervention ndrinks age female drugrisk), pom 

	
log close