#code for Chapter 3 "Propensity Score Weighting" of book:
#Leite, W. L. (2017). Practical propensity score methods using R. 
#Thousand Oaks, CA: Sage. 
# 
#this is the code that was used to generate the example results in the book
#As the R software and the R packages used in this example are updated frequently
#some incompatibilities between the current code and new R versions or package versions
#may appear
#Any updates to the code will be posted at:
# http://www.practicalpropensityscore.com

#======================================================
#STRATEGIES TO DEAL WITH EXTREME PROPENSITYS CORE WEIGHTS WEIGHTS

#LOAD DATA
load(file="Chapter3_ELS_data_imputed_with_weights.Rdata")

#load required package
require(survey)
#change options in the package
options(survey.lonely.psu = "adjust")

#check the 99th percentile weight:
with(ELS.data.imputed,quantile(weightATE, 0.99))
#check how many observatios will have truncated weights
with(ELS.data.imputed,sum(weightATE>quantile(weightATE, 0.99)))
#truncate weights at the 99th percentile:
ELS.data.imputed$weightATETruncated = with(ELS.data.imputed,
                                           ifelse(weightATE > quantile(weightATE, 0.99), 
                                                  quantile(weightATE, 0.99),weightATE))    
#check truncated weights
with(ELS.data.imputed, by(weightATETruncated,treat,summary))


#obtain stabilized weights for the estimation of ATE
ELS.data.imputed$C <- with(ELS.data.imputed,ifelse(treat==1,pScores,1-pScores))
surveyDesign <- svydesign(ids=~psu, strata=~STRAT_ID, weights=~bystuwt,
                         data = ELS.data.imputed, nest=T) #recreate survey design
(constants <- svyby(~C, by=~treat, design=surveyDesign, FUN=svymean))
ELS.data.imputed$stabilizedWeightATE <- ifelse(ELS.data.imputed$treat==1, 
                                              constants[1,2]/ELS.data.imputed$C, 
                                              constants[2,2]/ELS.data.imputed$C)

#check stabilized weigts for extremeness
with(ELS.data.imputed, by(stabilizedWeightATE,treat,summary))

#compare the distributions of weights and stabilized weights for the ATE
tiff("Chapter3_figure3-1.tif", res=600, compression = "lzw", height=6, width=15, units="in")
with(ELS.data.imputed,
     hist(weightATE[treat==1], density = 10, angle = 45, main="Propensity Score Weights for the ATE",
             xlab="Shaded = IPTW | Gray = Stabilized IPTW") )
with(ELS.data.imputed,
     hist(stabilizedWeightATE[treat==1], col=gray(0.4,0.25), add=T) )
dev.off()

#save data with truncated and stabilized weights
save(ELS.data.imputed, file="Chapter3_ELS_data_imputed_with_weights.Rdata", compress=T)