#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
#Estimate the effect of high school student participation in
#career academies on future income
#using data from the Education Longitudinal Study (ELS)
#=====================================================================
#Doubly-robust (DR) estimation with regression estimation of
#propensity score weighted means using the propensity score as a covariate
#load required package
#load required survey library
require(survey)
options(survey.lonely.psu = "adjust")
#load data
load(file="Chapter_3_treatment_effect_estimates.Rdata")
#define survey design for treated observations
surveyDesign2006T <- subset(surveyDesign2006Boot, treat==1)
#define survey design for untreated observations
surveyDesign2006U <- subset(surveyDesign2006Boot, treat==0)
#fit model for treated and untreated observations
modelT <- svyglm(F2ERN5P2~pScores+I(pScores^2)+I(pScores^3),surveyDesign2006T)
modelU <- svyglm(F2ERN5P2~pScores+I(pScores^2)+I(pScores^3),surveyDesign2006U)
#obtain predicted outcomes with treated observations
Yt1 <- predict(modelT, newdata=data.frame(pScores=with(ELS.data.imputed,pScores[treat==1])),
vcov=TRUE, type="response")
Yt0 <- predict(modelU, newdata=data.frame(pScores=with(ELS.data.imputed,pScores[treat==1])),
vcov=TRUE, type="response")
#obtain the difference between predicted outcomes
diff <- Yt1 - Yt0
ATT2006DR <- svycontrast(diff, contrasts=rep(1/length(diff),length(diff)))
#save all objects
save(list=ls(), file="Chapter_3_treatment_effect_estimates.Rdata",compress=T)