#code for Chapter 6 "Propensity Score Methods for Multiple Treatments" of book:
#Leite, W. L. (2017). Practical propensity score methods using R.
#Thousand Oaks, CA: Sage.
#
#PART 5 - ESTIMATE THE AVERAGE TREATMENT EFFECT FOR MULTIPLE TREATMENT VERSIONS
#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
#this code estimates the effects of assigning mentors
#of different areas to new teachers on the probability
#that they will continue in the teaching profession
#in the following year
#using data from the 1999--2000 School and Staffing Survey (SASS)
#and 2000--2001 Teacher Follow-up Survey (TFS)
#load data with generalized propensity scores
load("chapter6_dataset_with_mmws_weights.Rdata")
#=============================================================
#estimate ATE with propensity score weighting
#Fdefine the design
#school ids are provided to control for effects of clustering
require(survey)
#set up the survey design
design.IPTW <- svydesign(id = ~SCHCNTL, weights = ~mmwsFinal, data = imputedData)
#create replicate weights for bootstrapping
design.IPTW.boot <- as.svrepdesign(design.IPTW, type=c("bootstrap"),replicates=1000)
#obtain ATE as weighted proportions.
(weightedProportions <- svyby(formula=~leftTeaching,by=~Treat,design=design.IPTW.boot,
FUN=svymean,covmat=TRUE))
#obtain the weighted proportions with group names to be used in setting up constrats
print(weightedProportions)
#calculate pairwise weighted differences between treatment versions
pairwise.ATE <- svycontrast(weightedProportions,
contrasts=list(
sameArea.noMentor= c("sameArea:leftTeaching1"=1,"noMentor:leftTeaching1"=-1),
sameArea.otherArea=c("sameArea:leftTeaching1"=1,"otherArea:leftTeaching1"=-1),
otherArea.noMentor=c("otherArea:leftTeaching1"=1,"noMentor:leftTeaching1"=-1)))
#Estimate treatmetn effects with regression
model.IPTW <- svyglm("leftTeaching~Treat",design=design.IPTW,
family="quasibinomial")
summary(model.IPTW)
#get the coefficients
effects <- unclass(summary(model.IPTW))$coefficients
print(effects)