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