#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 4 - OBTAIN MARGINAL MEAN WEIGHTS THROUGH STRATIFICATION AND CHECK COVARIATE BALANCE #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_generalized_propensity_score_estimation_results.Rdata") #======================================================= ##FUNCTION TO OBTAIN MARGINAL MEAN WEIGHTS through sratification #Arguments: #id = variable with id numbers for dataset #treat = treatment indicator (must be a factor) #ps = data.frame of generalized propensity scores with variable names equal to treatment levels. #numb.strata = number of strata desired mmws.weights <- function(id,treat, ps,numb.strata) { all.weights <- data.frame(id) #loop through levels of treatment for (t in levels(treat)) { #create strata strata <- cut(x=ps[,t],breaks=quantile(ps[,t], prob = seq(0, 1, 1/numb.strata)), labels=1:numb.strata,include.lowest=T) #create treatment by strata table treatment <- as.numeric(treat==t) treat.by.strata <- data.frame(xtabs(~strata+treatment)) #create strata table strata.table <- data.frame(table(strata)) names(strata.table)[2] <- "Freq.strata" #merge two tables treat.by.strata <- merge(treat.by.strata,strata.table) print(treat.by.strata) #create marginal mean weights treat.by.strata$mm.weight <- ifelse(treat.by.strata$treatment==1, mean(treatment)*treat.by.strata$Freq.strata/treat.by.strata$Freq, (1-mean(treatment))*treat.by.strata$Freq.strata/treat.by.strata$Freq) #create full dataset data.weight <- data.frame(id,treatment,strata) data.weight <- merge(treat.by.strata[,c(1,2,5)],data.weight) #convert to zero the weights of individuals that did not receive this level of treatment data.weight$mm.weight[data.weight$treatment==0] <- 0 data.weight <- data.weight[,3:4] names(data.weight)[1] <- paste("W.",t,sep="") all.weights <- merge(all.weights,data.weight) } #close loop through levels of treatment #create a single weight all.weights$mmws <- apply(all.weights[,-1],1,sum) return(all.weights[,c("id","mmws")]) } #close function #========================================================= #run function to obtain MMWS mmws <- mmws.weights(imputedData$CNTLNUM, imputedData$Treat, ps2, numb.strata=5) names(mmws)[1] <- "CNTLNUM" imputedData <- merge(imputedData,mmws) by(imputedData$mmws, imputedData$Treat, summary) #multiply by sampling weight and normalize imputedData$mmwsFinal <- with(imputedData,(mmws*TFNLWGT)/mean(mmws*TFNLWGT)) by(imputedData$mmwsFinal, imputedData$Treat, summary) #======================================================== #evaluate covariate balance with MMWS #function to evaluate covariate balance pairwise pairwise.balance <- function(condition1, condition2, data, iptw){ require(twang) data <- subset(data, Treat == condition1 | Treat == condition2) data$Treat <- as.numeric(data$Treat==condition1) balance.iptw <- bal.stat(data=data,#dataset is defined vars=covariateNames, #list of the covariates that are used in the treatment assignment model treat.var ="Treat" , #treatment assignment indicator w.all = data[,iptw], #This time, PS weight times Sampling weight is defined get.ks=F, #Avoid estimating KS statistic sampw = data$TFNLWGT, #indicate the sampling weight estimand="ATE", #Define proper estimand multinom=F) #Indicate that we make pairwise comparisons, return( balance.iptw$results) } #close function #=========================================== #obtain balance for all groups pairwise balance.mmws <- list(pair12 <- pairwise.balance("noMentor", "sameArea", imputedData, "mmwsFinal"), pair13 <- pairwise.balance("noMentor", "otherArea", imputedData, "mmwsFinal"), pair23 <- pairwise.balance("sameArea", "otherArea", imputedData, "mmwsFinal")) #obtain a dataset with just standardized effect sizes std.eff.mmws <- data.frame(abs(balance.mmws[[1]][5]), abs(balance.mmws[[2]][5]), abs(balance.mmws[[3]][5])) summary(std.eff.mmws) #identify covariates/categories with standardized mean difference above 0.25, which will be #used in covariance adjustment unbalanced.covariates.mmws <- row.names(std.eff.mmws)[as.logical(apply(std.eff.mmws>0.25,1,max))] unbalanced.covariates.mmws <- strsplit(unbalanced.covariates.mmws, ":") unbalanced.covariates.mmws <- unique(sapply(unbalanced.covariates.mmws, "[", 1)) print(unbalanced.covariates.mmws) #save dataset with weights save(imputedData, file="chapter6_dataset_with_mmws_weights.Rdata", compress=T)