##code for Chapter 9 "Weighting Methods for Time-Varying Treatments" of book:
#Leite, W. L. (2016). Practical propensity score analysis 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
#this example estimates the effect of self-employment on job satisfaction
#using data from the National Longitudinal Survey of Youth 1979 (NLSY79)
#load the dataset
load("chapter_9_long_data_imputed.Rdata")
#calculate inverse probability of treatment weights
require(ipw)
#estimate unstabilized IPTW
iptw <- ipwtm(exposure = selfEmploy,timevar=timeRecoded,
family = "binomial", link="logit", weights=data.long$normWeight,
denominator = ~ cumSE + age + female + mngr.m + rotter.index + white + degree +
divorce + experience + married + urban +
age:female + mngr.m:female + rotter.index:female +
white:female + degree:female + divorce:female +
experience:female + married:female + urban:female,
id = id, type = "all", data = data.long)
quantile(iptw$ipw.weights)
#add iptw to dataset
data.long$iptw <- iptw$ipw.weights #add weights to dataset
#obtain stabilized weights using an indicator of self empolyment in the previous year
#and the number of measurements with self employment in previous years
#in the numerator
siptw <- ipwtm(exposure = selfEmploy,timevar=timeRecoded, weights=data.long$normWeight,
family = "binomial", link="logit", numerator= ~1 + cumSE,
denominator = ~ cumSE + age + female + mngr.m + rotter.index + white + degree +
divorce + experience + I(experience^2) + experience:cumSE + I(experience^2):cumSE + married + urban +
age:female + mngr.m:female + rotter.index:female +
white:female + degree:female + divorce:female +
experience:female + married:female + urban:female,
id = id, type = "all", data = data.long)
quantile(siptw$ipw.weights)
data.long$siptw <- siptw$ipw.weights #add stablized weights to dataset
#obtain basic stabilized weights using only the proportion treated in the numerator
bsiptw <- ipwtm(exposure = selfEmploy,timevar=timeRecoded, weights=data.long$normWeight,
family = "binomial", link="logit", numerator= ~1,
denominator = ~ cumSE+age + female + mngr.m + rotter.index + white + degree +
divorce + experience + married + urban +
age:female + mngr.m:female + rotter.index:female +
white:female + degree:female + divorce:female +
experience:female + married:female + urban:female,
id = id, type = "all", data = data.long)
quantile(bsiptw$ipw.weights)
data.long$bsiptw <- bsiptw$ipw.weights #add basic stablized weights to dataset
# #plot inverse probability weights
# require(lattice)
# tiff("Chapter9_figure9-1.tif", res=600, compression = "lzw", height=6, width=15, units="in")
# bwplot(iptw~selfEmploy | time, data = data.long,
# ylab = "IPTW", auto.key = TRUE)
#dev.off()
#
#plot stabilized inverse probability weights
require(lattice)
tiff("Chapter9_figure9-1.tif", res=600, compression = "lzw", height=6, width=15, units="in")
bwplot(siptw~time | selfEmploy, data = data.long,
ylab = "Stabilized IPTW", auto.key = TRUE)
dev.off()
#
# #plot basic stabilized inverse probability weights
# require(lattice)
# tiff("Chapter9_figure9-3.tif", res=600, compression = "lzw", height=6, width=15, units="in")
# bwplot(bsiptw~selfEmploy | time, data = data.long,
# ylab = "Basic Stabilized IPTW", auto.key = TRUE)
#dev.off()
#obtain final weight
data.long$finalWeight <- with(data.long,siptw*normWeight)
#check balance at every wave with stabilized IPTW.
require(twang)
balance <- bal.stat(data=data.long,#dataset
var=colnames(data.long)[c(3:5,11:13,15:16)],#list of the variables to be included in the balance check
treat.var="selfEmploy",#treatment assignment variable
w.all=data.long$finalWeight,#the weight
sampw=data.long$normWeight, ,#define the sampling weights
get.means=T, #compute means and variances
get.ks=F,#overwrite the default for obtaining KS statistic
estimand="ATE",#specify the estimand of interest
multinom=F)#underline that multinomial ps are not used
summary(balance$results$std.eff.sz)
#obtain final weight with basic stabilized weight
data.long$finalWeight2 <- with(data.long,normWeight*bsiptw)
#evaluate balance with the basic stabilized weights
balance2=bal.stat(data=data.long,#dataset
var=colnames(data.long)[c(3:5,11:13,15:16)],#list of the variables to be included in the balance check
treat.var="selfEmploy",#treatment assignment variable
w.all=data.long$finalWeight2,#the weight
sampw=1,#define the sampling weights
get.means=T, #compute means and variances
get.ks=F,#overwrite the default for obtaining KS statistic
estimand="ATE",#specify the estimand of interest
multinom=F)#underline that multinomial ps are not used
summary(balance2$results$std.eff.sz)
#=====================================================
# #estimate inverse probability of treatment weights with covariate balancing propensity scores
# require(CBPS)
# CBPS.formula = "selfEmploy ~ age + female + mngr.m + rotter.index + white + degree + divorce + experience + married + urban"
# F96.12 = formula("selfEmploy ~ age + female + mngr.m + rotter.index + white + degree +
# divorce + experience + married + urban")
#
# cbiptw = CBMSM(formula = list(CBPS.formula,CBPS.formula,CBPS.formula,CBPS.formula,CBPS.formula,
# CBPS.formula,CBPS.formula,CBPS.formula,CBPS.formula,CBPS.formula),
# id=data.long$id, time=data.long$time, data=data.long, type="MSM", twostep = TRUE,
# msm.variance = "approx", time.vary = T)
#========================================================
#grand mean center the continuous covariates that will be included in the outcome model
#to preserve interpretation of the effect as the ATE
data.long$cumSE <- with(data.long, cumSE - mean(cumSE))
data.long$experience <- with(data.long, experience - mean(experience))
#create the design object for IPTW weights
require(survey)
design.iptw <- svydesign(ids=~id,nest=F,weights =~finalWeight,data=data.long)
#obtain boostrapping of standard errors
# surveyDesignBoot = as.svrepdesign(design.iptw, type=c("bootstrap"),replicates=1000)
# #obtain the group means and standard errors
# (means.iptw=svyby(formula=~genSatis,by=~selfEmploy,design=design.iptw, FUN=svymean))
# #obtain group mean difference by contrast statement
# (contrast.iptw<- svycontrast(means.iptw, contrasts=c(-1,1)))#not self employed group is contrasted by -1
# (contrast.iptw.t=coef(contrast.iptw)/SE(contrast.iptw))#T statistic is estimated as coefficient over standard error
# (contrast.iptw.p=2*(1-pt(abs(contrast.iptw.t),nrow(data.long)-2)))#appropriate p value is obtained from the student's t distribution
#
# #==================================================
#estimate the treatment effect with weighted least squares regression with
#cluster effect adjustment of standard errors through Taylor series linearization
glm.iptw <- svyglm(genSatis~selfEmploy+cumSE+experience+timeRecoded+
selfEmploy:cumSE + selfEmploy:experience+selfEmploy:timeRecoded,design.iptw)
summary(glm.iptw)
#=========================================================
#generalized estimating equations
library(geepack)
gee.iptw <- geeglm(genSatis~selfEmploy+cumSE+experience+timeRecoded+
selfEmploy:cumSE + selfEmploy:experience+selfEmploy:timeRecoded,
data = data.long, weights = data.long$finalWeight, family = gaussian,
id = data.long$id, waves=data.long$time, corstr = "ar1")
summary(gee.iptw)
#save all results
#save(list=ls(), file="Chapter9_all_results.Rdata",compress=T)