##code for Chapter 7 "Propensity Score Methods for Continuous Treatment Doses" 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

#PART 2 INVERSE PROBABILITY WEIGHTING

#this example estimates the effect of school participation
#in the Algebra Nation virtual learning environment
#on the mean student scores on Florida's Algebra I End-of-Course (EOC) Assessment.
#load data
load(file = "Chapter7_data_Algebra_Nation_with_GPS.Rdata")
#====================================================
# Estimation of continuous treatment effects
#with inverse probability weights

# numerator of weighting equation
data$numerator <- with(data, dnorm(logLoginsPerExaminee,
                 mean=mean(logLoginsPerExaminee),
                 sd=sd(logLoginsPerExaminee)))

# calculate the IPW
data$IPW <- with(data, numerator/GPS)

summary(data$IPW)
#normalize weights
data$IPW <- data$IPW/mean(data$IPW)
summary(data$IPW)
#=============================================
#calculation of IPW using the ipw package
library(ipw)

#obtain the IPW
IPW2 <- ipwpoint(exposure=logLoginsPerExaminee,
                family="gaussian", numerator= ~ 1,
              denominator= ~ Charter+Magnet.+Title.I.School.+
                  locationRural+locationSize+Students.+
                  SeniorHigh+numOfStud2014+meanScale2012+
                  lev1Perc2012+lev5Perc2012+
                  perc.free.lunch+perc.reduced.lunch,
              data= data)

data$IPW2 <- IPW2$ipw.weights
summary(data$IPW2)
#normalize weights
data$IPW2 <- data$IPW2/mean(data$IPW2)
summary(data$IPW2)

with(data,cor(IPW,IPW2)) #compare the two IPW

#====================================================
#calculation of IPW using the Covariate Balance Propensity Score
#obtain the covariate-balancing propensity score
library(CBPS)

#define model for treatment dose
formulaDose <- formula("logLoginsPerExaminee~Charter+Magnet.+
                           Title.I.School.+locationRural+locationSize+Students.+
                        SeniorHigh+numOfStud2014+meanScale2012+lev1Perc2012+
                        lev5Perc2012+perc.free.lunch+perc.reduced.lunch")


cbps <- CBPS(formulaDose, data=data, standardize =F, method="exact")
summary(cbps)

data$IPW3 <- cbps$weights
summary(data$IPW3)
#normalize weights
data$IPW3 = data$IPW3/mean(data$IPW3)
summary(data$IPW3)

#========================================================
#estimate weights with a linear model for the treatment dose
#and with kerney density estimation for obtaining the generalized propensity score
#see https://en.wikipedia.org/wiki/Kernel_density_estimation

#using the weightit package
library(WeightIt)
weight.model = weightit(formula("logLoginsPerExaminee~Charter+Magnet.+
                           Title.I.School.+locationRural+locationSize+Students.+
                        SeniorHigh+numOfStud2014+meanScale2012+lev1Perc2012+
                        lev5Perc2012+perc.free.lunch+perc.reduced.lunch"), data = data,
                        method = "ps", use.kernel=T)

#check weights
summary(weight.model$weights)

#add weights to data
data$IPW4 <- weight.model$weights

#=========================================================
#estimate weights by modeling the treatment dose
#with Bayesian additive regression trees using
#kerney density estimation for obtaining the generalized propensity score

#using the weightit package
library(WeightIt)
weight.model2 = weightit(formula("logLoginsPerExaminee~Charter+Magnet.+
                           Title.I.School.+locationRural+locationSize+Students.+
                        SeniorHigh+numOfStud2014+meanScale2012+lev1Perc2012+
                        lev5Perc2012+perc.free.lunch+perc.reduced.lunch"), data = data,
                        method = "bart", use.kernel=T)

#check weights
summary(weight.model2$weights)

#add weights to data
data$IPW5 <- weight.model2$weights


#check correlation betwen IPW versions
with(data,cor(cbind(IPW,IPW2,IPW3,IPW4, IPW5)))

#save results
save(list=ls(),file="Chapter7_data_with_IPW.Rdata", compress=T)