##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 4 ESTIMATION OF AVERAGE TREATMENT EFFECT 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 data load(file="Chapter7_data_with_IPW.Rdata") #estimate the average treatment effect of logins per examinee library(survey) #define survey design designIPW <- svydesign(ids=~1,weights=~IPW3,data=data) #design with weights #estimate outcome model outcomeModel <- svyglm(formula= meanScale2014~logLoginsPerExaminee, design = designIPW, family=gaussian) summary(outcomeModel) #obtain a standardized coefficient coef(outcomeModel)[2]*sqrt(coef(svyvar(~logLoginsPerExaminee,designIPW)))/ sqrt(coef(svyvar(~meanScale2014,designIPW))) #estimate the average treatment effect of logins per examinee #with a model with polynomial terms outcomeModelNl <- svyglm(formula= meanScale2014~logLoginsPerExaminee+ I(logLoginsPerExaminee^2) + I(logLoginsPerExaminee^3), design = designIPW, family=gaussian) summary(outcomeModelNl) #compare linear and nonlinear models anova(outcomeModelNl, outcomeModel) #estimate the average treatment effect of logins per examinee #using the IPW and controlling for covariates outcomeModelDR <- svyglm(formula= meanScale2014~logLoginsPerExaminee + Charter+Magnet.+ Title.I.School.+locationRural+locationSize+Students.+ SeniorHigh+numOfStud2014+meanScale2012+lev1Perc2012+ lev5Perc2012+perc.free.lunch+perc.reduced.lunch, design = designIPW, family=gaussian) summary(outcomeModelDR) #obtain a standardized coefficient coef(outcomeModelDR)[2]*sqrt(coef(svyvar(~logLoginsPerExaminee,designIPW)))/ sqrt(coef(svyvar(~meanScale2014,designIPW))) #sensitivity analysis using Cinelli's method #article: Carlos Cinelli and Chad Hazlett (2020). #Making Sense of Sensitivity: Extending Omitted Variable Bias. #Journal of the Royal Statistical Society, Series B (Statistical Methodology). library(sensemakr) sensitivity = sensemakr( outcomeModelDR, treatment="logLoginsPerExaminee", q = 1, #proportion of change in effect deemed problematic (1 = 100%) alpha = 0.05) summary(sensitivity)