##code for Chapter 9 "Weighting Methods for Time-Varying Treatments" 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

#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) #check distribution of weights

#add stabilized weights to daset
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) #check distribution of weights

data.long$bsiptw <- bsiptw$ipw.weights #add basic stabilized 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()

#================================================
#COVARIATE BALANCE EVALUATION

#obtain final weight with stabilized 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)

#========================================================
#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)