##code for Chapter 7 "Propensity Score Methods for Continuous Treatment Doses" of book: #Leite, W. L. (2016). 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 3 EVALUATION OF COVARIATE BALANCE WITH INVERSE PROBABILITY WEIGHTS #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_with_IPW.Rdata") #define covariates covariateNames <- c("Charter", #dummy indicator of charter school "Magnet.", #dummy indicator of magnet school "Title.I.School.", #dummy indicator of title 1 school "locationRural", #rural location "locationSize", #size of location (large, midsize, small) "Students.", #total number of students in 2012 "SeniorHigh" , #dummy indicator of whether it is a senior high school (the reference group is high school "numOfStud2014", #number of test takers in 2014 "meanScale2012", #mean scaled scores in 2012 "lev1Perc2012",#percent achieving level 1 in 2012 "lev5Perc2012", #percent achieving level 5 in 2012 "perc.free.lunch", #percent free lunch in 2012 "perc.reduced.lunch") #percent reduced lunch in 2012 #set survey design library(survey) designIPW <- svydesign(ids=~1,weights=~IPW,data=data) #design with weights designBase <- svydesign(ids=~1,weights=~1,data=data) #create empty balance table balanceTableIPW <- data.frame() #storage #loop through covariates for (var in 1:length(covariateNames)) { #formula with treatment dose as a function of a single covariate balanceFormula <- paste("logLoginsPerExaminee~",covariateNames[var],sep="") #regress dose on covariate without weights maxEffBaseline <- max(abs(coef(svyglm(balanceFormula,designBase))[-1])) #regress dose on covariate with IPW maxEffIPW <- max(abs(coef(svyglm(balanceFormula,designIPW))[-1])) balanceTableIPW <- rbind(balanceTableIPW,c(var,maxEffBaseline,maxEffIPW)) } #close loop #put variable names on table names(balanceTableIPW) <- c("variable","coefBaseline","coefIPW") balanceTableIPW$variable <- covariateNames #standardize coefficients with respect to sd of outcome balanceTableIPW$coefBaseline <- balanceTableIPW$coefBaseline/ sqrt(coef(svyvar(~logLoginsPerExaminee,designBase))) balanceTableIPW$coefIPW <- balanceTableIPW$coefIPW/ sqrt(coef(svyvar(~logLoginsPerExaminee,designIPW))) #save balance table write.csv(balanceTableIPW,file="Chapter7_table72.csv",row.names=F) #====================================== #Evaluate covariate balance with the CBPS package library(CBPS) #the object cpbs was created in Part 2 of the code balance.table.cbps <- balance(cbps) #re-structure the list into a single data frame balance.table.cbps <- data.frame(original=balance.table.cbps$unweighted, weighted=balance.table.cbps$balanced) names(balance.table.cbps)<- c("Pearson r baseline", "Person r IPW") summary(balance.table.cbps) #================================================ #use the package cobalt to evaluate covariate balance library(cobalt) #obtain covariate balance for weights obtained with the linear model for treatment dose #and Gaussian density balance.table = bal.tab(logLoginsPerExaminee~Charter+Magnet.+ Title.I.School.+locationRural+locationSize+Students.+ SeniorHigh+numOfStud2014+meanScale2012+lev1Perc2012+ lev5Perc2012+perc.free.lunch+perc.reduced.lunch, treat = "logLoginsPerExaminee", data = data, weights = data$IPW) balance.table <- balance.table$Balance #covariate balance for weights obtained with CBPS balance.table.cbps2 = bal.tab(cbps, r.threshold=0.1)$Balance #covariate balance for weights obtained with linear model for treatment dose #and kernel density estimation balance.table.2 = bal.tab(weight.model, r.threshold=0.1)$Balance #covariate balance for weights obtained Bayesian Additive Regression Trees (BART) #and kernel density estimation balance.table.3 = bal.tab(weight.model2, r.threshold=0.1)$Balance