% setwd("d:/dev/twang/{t/doc")
% Sweave("twang.rnw"); system("texify twang.tex"); system("c:\\MiKTeX\\miktex\\bin\\yap.exe twang.dvi",wait=FALSE)

\SweaveOpts{prefix.string=twang}

\documentclass{article}
\bibliographystyle{plain}
\usepackage[active]{srcltx}
\usepackage{url}
\addtolength{\topmargin}{-0.5in}
\addtolength{\textheight}{0.75in}
\addtolength{\oddsidemargin}{-0.5in}
\addtolength{\textwidth}{1in}
\newcommand{\EV}{\mathrm{E}}
\newcommand{\Var}{\mathrm{Var}}
\newcommand{\aRule}{\begin{center} \rule{5in}{1mm} \end{center}}

 %% Tracking changes and notes code added by Dan Mc 6-30-2014 %%
\usepackage[normalem]{ulem} % for "tracked changes"
\usepackage{color} % for "tracked changes"
\newcommand{\note}[1]{{\color{blue} #1}}
\newcommand{\dele}[1]{\textcolor{red}{\sout{#1}}} % delet in an equation
\newcommand{\ins}[1]{\textcolor{blue}{\bf #1}} % insert

\title{Propensity scores for repeated treatments:\\A
tutorial for the \texttt{iptw} function in the \texttt{twang} package}

\author{Lane Burgette, Beth Ann Griffin and Dan McCaffrey 
  \footnote{The development of this software and tutorial was funded
    by National Institute of Drug Abuse grants number 1R01DA015697 (PI:
    McCaffrey) and 1R01DA034065 (PIs: Griffin/McCaffrey).} 
\\RAND Corporation}

%\VignetteIndexEntry{Propensity scores for time-varying treatments: A tutorial for the iptw function}
%\VignetteDepends{gbm,survey,xtable}
%\VignetteKeywords{propensity score, nonresponse weighting}
%\VignettePackage{twang}

 \newcommand{\mathgbf}[1]{{\mbox{\boldmath$#1$\unboldmath}}}

\begin{document}
\SweaveOpts{concordance=TRUE}

\maketitle

\section{Introduction}

While standard propensity score methods attempt to answer the question
of how expected outcomes change if a group of individuals received one
treatment instead of another, researchers are often
interested in understanding how {\em sequences} of treatments
impact outcomes of interest. In this case, time-varying confounders
may be impacted by prior treatments. Consequently, simply controlling
for the time-varying confounders in standard regression models can
yield biased results. Instead, it is possible to perform weighted
regressions that account for time-varying confounders via {\em
  marginal structural models} (MSMs; Robins et al., 2000). In this
method, observations are weighted by the 
inverse of the estimated probability of 
receiving the observed sequence of treatments the individual actually
received, referred to as an {\em inverse probability of treatment
  weight} (IPTW). It has been proposed to use nonparametric
models to estimate IPTWs (Griffin et al.,
2014). Accordingly, we refer 
to the function in \texttt{twang} that performs this weighting as
\texttt{iptw}, for inverse probability of treatment weighting. 

For binary treatments, the \texttt{iptw} methods and syntax build
directly on the \texttt{ps} 
functionality; users are encouraged to study that tutorial before
using \texttt{iptw}. For treatment regimes with more than two
categories, the \texttt{iptw} methods build on the \texttt{mnps}
methods and software. For more background on marginal structural
models, see e.g., Robins et al.\ (2000) and Cole and Hern\'an (2008). 


\section{An IPTW example}

<<echo=FALSE>>=
options(width=80)
@


For the sake of illustration, we simulated data to demonstrate
the functionality of the \texttt{iptw} command. For time-varying treatment
data, one can either imagine 
a ``wide'' dataset, with one row per subject, or a ``long'' dataset
with one row for each subject/time point combination. Our artificial
data include time-invariant characteristics \texttt{gender}, and
\texttt{age} at time of 
study enrollment. Conceptually, we have a substance use index that is
measured four times: at baseline, after the first treatment period,
after the second treatment period, and after the third treatment
period, which concludes the study and is the outcome of interest.  In
the ``wide'' version of the 
dataset called \texttt{iptwExWide}, we have the \texttt{outcome},
baseline and intermediate 
measures, \texttt{use0}, \texttt{use1}, and \texttt{use2}. The
treatment indicators are, in chronological order, \texttt{tx1},
\texttt{tx2}, and \texttt{tx3}.  Our goal is to estimate the average
effect of each additional dose of treatment on substance use measured
at the end of the study (which is recorded in \texttt{outcome}). 

The ``long'' format data have a somewhat different form, and are
included in the data object \texttt{iptwExLong}. For the long format,
the outcomes are split from the covariates, and are available as
\texttt{iptwExLong\$outcome}. Similarly, the covariates and treatment
indicators are available in \texttt{covariates}, which 
includes data elements \texttt{gender}, \texttt{age},  
\texttt{use}, and \texttt{tx}; these include the same information as 
the wide version. Additionally, the long version contains elements 
\texttt{ID} (an individual-level identifier) and \texttt{time}, which 
corresponds to the study period.   

One of the benefits of GBM for applied researchers is the automatic 
handling of missing data. For MSMs, however, this does not extend to 
partially observed treatment histories. We assume throughout that 
missingness exists only in the covariates. 

\subsection{Fitting the model} 

To begin, we will work with the ``wide'' version of the data, which
are available after loading the twang package:

<<>>=
library(twang)
data(iptwExWide)
@

Next, we can fit the model using the \texttt{iptw} function. Unlike
for the standard \texttt{ps} function, we are only able to use a
single stop.method at a time. The treatment assignment models are
specified as a list of formulas, starting at the earliest time
period. For coding parsimony, terms that should appear in all of the
formulas can be specified once via a one-sided formula using the
\texttt{timeInvariant} argument.  Similarly, including treatment
indicators from previous models is achieved by setting
\texttt{priorTreatment = TRUE}.  Finally, if all terms included at
period $t$ should be included in the period $t+1$ model (as is
typically the case in MSM models), setting
\texttt{cumulative = TRUE} automatically includes all elements on the
right-hand side of previous models.

Thus, the model

<<eval=FALSE>>=
iptw.Ex <- iptw(list(tx1 ~ use0 + gender + age, 
                     tx2 ~ use1 + use0 + tx1 + gender + age, 
                     tx3 ~ use2 + use1 + use0 + tx2 + tx1 + gender + age),
                 timeInvariant ~ gender + age,
                 data = iptwExWide, 
                 cumulative = FALSE,
                 priorTreatment = FALSE,
                 verbose = FALSE, 
                 stop.method = "es.max", 
                 n.trees = 5000)
@ 

can be specified more simply as:

<<>>=
iptw.Ex <- iptw(list(tx1 ~ use0, tx2 ~ use1, tx3 ~ use2),
                 timeInvariant ~ gender + age,
                 data = iptwExWide, 
                 cumulative = TRUE,
                 priorTreatment = TRUE,
                 verbose = FALSE, 
                 stop.method = "es.max", 
                 n.trees = 5000)
@ 

After having fit the \texttt{iptw} object, the diagnostic checks are
similar to those specified for \texttt{ps} objects. 

First, we check to make sure that the GBM models were allowed to run
long enough (i.e., \texttt{n.trees} is sufficiently large). 


\begin{center}
<<fig=TRUE, echo=TRUE, include = TRUE, height = 7, width = 10.5>>= 
    plot(iptw.Ex, plots = 1)
@
\end{center}

Next, we can get a quick sense of the balance at each timepoint via
<<>>=
summary(iptw.Ex)
@
Further detail regarding the model at, e.g., the third time period is available using 
<<>>=
bal.table(iptw.Ex, timePeriods = 3)
@ 


Next, we can examine propensity score overlap at each time point: 

\begin{center}
<<fig=TRUE, echo=TRUE, include = TRUE, height = 7, width = 10.5>>= 
    plot(iptw.Ex, plots = 2)
@
\end{center}

These figures can focus on the results from particular time periods
using the \texttt{timePeriods} argument:

\begin{center}
<<fig=TRUE, echo=TRUE, include = TRUE, height = 7, width = 10.5>>= 
    plot(iptw.Ex, plots = 2, timePeriods = 2:3)
@
\end{center}

Next, we can check balance as measured by standardized mean
differences between the treated and control samples at each of the
time points by specifying 
\texttt{plots = 3}. As with other
TWANG figures, we can specify \texttt{color = FALSE} to produce black
and white 
figures. 

\begin{center}
<<fig=TRUE, echo=TRUE, include = TRUE, height = 7, width = 10.5>>= 
    plot(iptw.Ex, plots = 3, color = FALSE)
@
\end{center}

Finally, we can compare the differences between the treated and
control samples using $t$-test and KS $p$-values by specifying
\texttt{plots = 4} and \texttt{plots = 5}, respectively. 

\begin{center}
<<fig=TRUE, echo=TRUE, include = TRUE, height = 7, width = 10.5>>= 
    plot(iptw.Ex, plots = 4)
@
\end{center}

\begin{center}
<<fig=TRUE, echo=TRUE, include = TRUE, height = 7, width = 10.5>>= 
    plot(iptw.Ex, plots = 5)
@
\end{center}

Further, \texttt{iptw} can accommodate treatments with more than two
levels (McCaffrey et al., 2013). An example can be explored in the
following call, though we do 
not discuss it further in this vignette. See the \texttt{mnps}
vignette for more information on the diagnostic plots. 

<<eval = FALSE>>=
data(mnIptwExWide)
mniptw.Ex <- iptw(list(tx1 ~ use0, tx2 ~ use1, tx3 ~ use2),
                 timeInvariant ~ gender + age,
                 data = mnIptwExWide, 
                 cumulative = TRUE,
                 priorTreatment = TRUE,
                 verbose = FALSE, 
                 stop.method = "es.max", 
                 n.trees = 5000)
@ 

\section{Estimating treatment effects}

After having estimated the relevant propensity scores, the final step
is translating them into analytic weights and estimating treatment
effects. Twang provides several functions to facilitate this
process. For this analysis, we assume an additive treatment model,
where the mean change in outcomes depends on the number of periods of
treatment. Because the weights often have substantial variation, the
weights are commonly {\em stabilized} where the standard inverse
probability of treatment weights are multiplied by the estimated
probability of receiving the treatment that each individual received,
conditioning only on previous periods' treatment indicators. 

 To begin, we calculate {\em unstablilized} weights. These are
 computed as the inverse probability of treatment weight, and are
 available as 
<<>>=
unstabWt1 <- get.weights.unstab(iptw.Ex)
@ 


 We can estimate the treatment effect using these unstabilized weights
as follows. the number of periods of treatment for each individual
<<>>=
library(survey)
nTx <- with(iptwExWide, tx1 + tx2 + tx3)
outDatUnstab <- data.frame(outcome = iptwExWide$outcome,
                     nTx, 
                     wt = unstabWt1$es.max.ATE)
sv1unstab <- svydesign(~1, weights = ~wt, data = outDatUnstab)
@ 

We can then calculate the point estimate and 95\% confidence interval using the
unstabilized weights as
<<>>=
fitUnstab <- svyglm(outcome ~ nTx, sv1unstab)
coef(fitUnstab)
confint(fitUnstab)
@ 

To calculate the stabilized weights, we additionally calculate a
stabilizing factor that depends on on the marginal probabilities of
treatment. This can be done via 
<<>>=
fitList <- list(glm(tx1 ~ 1, family = binomial, data = iptwExWide),
                glm(tx2 ~ tx1, family = binomial, data = iptwExWide), 
                glm(tx3 ~ tx1 + tx2, family = binomial, data = iptwExWide))
numWt <- get.weights.num(iptw.Ex, fitList)
stabWt1 <- unstabWt1 * numWt
outDatStab <- data.frame(outcome = iptwExWide$outcome,
                     nTx, 
                     wt = stabWt1$es.max.ATE)
sv1stab <- svydesign(~1, weights = ~wt, data = outDatStab)
@ 

As before, we can then estimate the treatment effect and associated
confidence interval
<<>>=
fitStab <- svyglm(outcome ~ nTx, sv1stab)
coef(fitStab)
confint(fitStab)
@ 

Since these are simulated data, we know that the true treatment effect
is -0.1. We can see that both of the propensity score-weighted
estimates cover the true treatment effect. 

For comparison, we examine the unadjusted effect estimate, which we
see does not include the true value:
<<>>=
confint(lm(iptwExWide$outcome ~ nTx))
@ 

\section{Conclusion}

Frequently, researchers are interested in treatments that may vary
period-by-period. Twang's \texttt{iptw} function provides a
nonparametric method for calculating inverse probability of treatment
weights for marginal structural models. The function can accommodate
treatments with two or more levels. The diagnostic figures and tables
build on of the \texttt{mnps} and \texttt{ps} commands, with
additional features to help manage the numerous possible comparisons. 

\section{References}

Cole, S.R.\ and Hern\'an (2008). ``Constructing inverse probability
weights for marginal structural models.'' {\em American Journal of
  Epidemiology}, 168(6), 656-664.

Griffin, B.A., R.\ Ramchand, D.\ Almirall, M.E.\ Slaughter, L.F.\
Burgette, and D.F.\ McCaffrey (2014). ``Estimating the causal effects
of cumulative treatment episodes for adolescents using marginal
structural models and inverse probability of treatment weighting. {\em
  Drug and Alcohol Dependence}, 136(1), 69--78. 

McCaffrey, D.F., B.A.\ Griffin, D.\ Almirall, M.E.\ Slaughter, R.\
Ramchand, and L.F.\ Burgette (2013). ``A tutorial on propensity score
estimation for multiple treatments using generalized boosted models.''
{\em Statistics in Medicine}, 32(19): 3388-3414. 

Robins, J.M., M.\'A. Hern\'an, and B.\ Brumback (2000). ``Marginal structural
models and causal inference in epidemiology.'' {\em Epidemiology},
11(5), 550--560. 

\end{document}