#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)