\documentclass[nojss]{jss}
%\VignetteEngine{knitr::knitr}
%\VignetteIndexEntry{Time-to-event outcomes with BART}
%\usepackage{thumbpdf,lmodern}
\usepackage{verbatim}
\usepackage{statex2}
\usepackage[authoryear,round]{natbib}
\usepackage{rotating}
%\usepackage{pdfsync}
%\synctex=1
%% -- Article metainformation (author, title, ...) -----------------------------

%% - \author{} with primary affiliation
%% - \Plainauthor{} without affiliations
%% - Separate authors by \And or \AND (in \author) or by comma (in \Plainauthor).
%% - \AND starts a new line, \And does not.
%\author{Achim Zeileis\\Universit\"at Innsbruck
%   \And Second Author\\Plus Affiliation}
%\Plainauthor{Achim Zeileis, Second Author}
\author{Rodney Sparapani\\Medical College of Wisconsin
\And Robert McCulloch\\Arizona State University}
\Plainauthor{Rodney Sparapani, Robert McCulloch}

%% - \title{} in title case
%% - \Plaintitle{} without LaTeX markup (if any)
%% - \Shorttitle{} with LaTeX markup (if any), used as running title
%\title{A Short Demo Article: Regression Models for Count Data in \proglang{R}}
%\Plaintitle{A Short Demo Article: Regression Models for Count Data in R}
%\Shorttitle{A Short Demo Article in \proglang{R}}

\title{Time-to-event outcomes with BART}
\Plaintitle{Time-to-event outcomes with BART}
\Shorttitle{Time-to-event outcomes with BART}

%% - \Abstract{} almost as usual
\Abstract{
  This short article illustrates how to analyze time-to-event
outcomes with the {\bf BART} \proglang{R} package.
}

%% - \Keywords{} with LaTeX markup, at least one required
%% - \Plainkeywords{} without LaTeX markup (if necessary)
%% - Should be comma-separated and in sentence case.
\Keywords{Bayesian Additive Regression Trees} 
\Plainkeywords{Bayesian Additive Regression Trees}

%% - \Address{} of at least one author
%% - May contain multiple affiliations for each author
%%   (in extra lines, separated by \emph{and}\\).
%% - May contain multiple authors for the same affiliation
%%   (in the same first line, separated by comma).

\Address{
  Rodney Sparapani {rsparapa@mcw.edu}\\
  Division of Biostatistics, Institute for Health and Equity\\
  Medical College of Wisconsin, Milwaukee campus
%  8701 Watertown Plank Road\\
%  Milwaukee, WI\ \ 53226, USA\\
%  E-mail: {rsparapa@mcw.edu}
}

\begin{document}

\maketitle


\begin{comment}
Funding for this research was provided, in part, by the
Advancing Healthier Wisconsin Research and Education Program
under awards 9520277 and 9520364.
\end{comment}

\section{Time-to-event outcomes with BART}

The {\bf BART} package supports time-to-event outcomes including
survival analysis, competing risks and recurrent events.

\subsection{Background: survival analysis with the Cox proportional hazard model}

The inspiration for survival analysis with BART is the classic
standard: the Cox proportional hazard model \citep{Cox72}.  The data
is $(s_i, \delta_i, \bm{x}_{i})$ where $s_i$ is the time of an
absorbing event, $\delta_i=1$, or right censoring, $\delta_i=0$, and
$ \bm{x}_{i} $ is a vector of covariates (which can be time-dependent,
but, for simplicity, we assume that they are known at time zero).  We
construct a grid of the ordered distinct event times,
$0=t_{(0)}<\dots<t_{(J)}<\infty$, and we consider the following time
intervals:
$(0, t_{(1)}], (t_{(1)}, t_{(2)}], \. (t_{(J-1)}, t_{(J)}]$.  The
general form of the Cox proportional hazard model is the following:
$\lambda(t_{(j)},
\bm{x}_i)=\lambda_0(t_{(j)})\exp(\bm{\beta}'\bm{x}_i)$
where $\lambda(t_{(j)}, \bm{x}_i)$ is the hazard, $\lambda_0(t_{(j)})$
is a nonparametric baseline hazard defined at the grid of time points
and $\exp(\bm{\beta}'\bm{x}_i)$ is a parametric multiplier which we
call linear proportionality.  To perform estimation and inference of
$\bm{\beta}$, we utilize what is known as the partial likelihood:
$ \wrap{\bm{\beta}|\lambda_0(t)} = \prod_i
\frac{\e{\bm{\beta}'\bm{x}_{i}}}{\sum_{{j \in R(t_i)}}
  \e{\bm{\beta}'\bm{x}_{j}}} $
where $R(t_i)$ is the set of subjects at risk for an event at time
point $t_i$ (which for events is a grid point by definition).  The
cumulative baseline hazard can be estimated as
$\hat\Lambda_0(t_{(j)}) = \sum_{t_i\le t_{(j)} }
\frac{\delta_i}{\sum_{j \in R(t_i)} \e{\bm{\hat\beta}'\bm{x}_{i}}}$.
The baseline survival is
$\hat{S}_0(t_{(j)}) = \e{-\hat\Lambda_0(t_{(j)})} $ and the general
survival is
$\hat{S}(t_{(j)}, \bm{x}_{i}) =
\hat{S}_0(t_{(j)})^{\exp(\bm{\hat\beta}'\bm{x}_{i}) }$.
Notice that we don't directly estimate the survival; rather, we
estimate $\bm{\beta}$ and survival is a consequence of this estimate
by construction.  This feature, and the time grid, foreshadow elements
of survival analysis with BART.

\subsection{Survival analysis with BART}

Survival analysis with BART is provided by the \code{surv.bart}
function for serial computation and \code{mc.surv.bart} for parallel
computation.
%There are many potential approaches that could be taken to utilize
%BART in survival analysis. 
The complete details of our approach can be found in
\citet{SparLoga16} and a brief introduction follows.  We take an
approach that is tantamount to discrete-time survival analysis
\citep{Thom77,ArjaHaar87,FahrTutz94,Fahr98}.  Following the
capabilities of BART, we do not stipulate linearity nor proportional
hazards.  We use the same notation developed for the Cox proportional
hazards model above.

Now, consider event indicators $y_{ij}$ for each subject $i$ at each
distinct time $t_{(j)}$ up to and including the subject's last
observation time $t_i=t_{(n_i)}$ with
%$n_i=\#\{j:t_{(j)}\leq t_i\}$ or
$n_i=\arg \max_j\wrap{t_{(j)}\leq t_i}$.  This means $y_{ij}=0$ if
$j<n_i$ and $y_{in_i}=\delta_i$. Denote the probability of an event at
time $t_{(j)}$, conditional on no previous event, by $p_{ij}$.  Now, our
model for $y_{ij}$ is a nonparametric probit regression of
$y_{ij}$ on the time $t_{(j)}$ and the covariates $\bm{x}_i$.  We
utilize the \citet{AlbeChib93} truncated Normal latent variables
$z_{ij}$ to recast it as a continuous BART
model where the latents are the outcome.
We choose \citet{AlbeChib93} Normal latents as the
default for computational efficiency, but we also provide the
optional \citet{HolmHeld06} Logistic latents by specifying
\code{type='lbart'}.

So the model is
\begin{align*}%\label{eq:survbartdata}
y_{ij} & = \delta_i \I{s_i=t_{(j)}},\ j=1, \dots, n_i \\
y_{ij}|p_{ij} & ~  \B{p_{ij}} \\
p_{ij}|f & =  \Phi(\mu_{ij}),\ \mu_{ij} =\mu_0+f(t_{(j)}, \bm{x}_i)  \\
f & \prior  \mathrm{BART} \\
z_{ij}|y_{ij},f & ~  \begin{cases} 
\N{ \mu_{ij}}{1}\I{-\infty, 0} & \If y_{ij}=0 \\
\N{ \mu_{ij}}{1}\I{0, \infty} & \If y_{ij}=1 \\
\end{cases} 
\end{align*}
where $\Phi$ is the standard Normal cumulative distribution function.
% With the data vector $\bm{y}$ made up of $n$ independent
% sequences of $0$'s and $1$'s given $\bm{p}$ (the entire collection of
% $p_{ij}$'s), we have the following likelihood:
% $\wrap{\bm{y}|\bm{p}}\ = \prod_{i=1}^n\prod_{j=1}^{n_i} p_{ij}^{y_{ij}} 
% (1-p_{ij})^{1-y_{ij}}$\ . The product over $j$ is a result of the
% definition of $p_{ij}$'s as conditional probabilities, and not a consequence of
% an assumption of independence.

If the event indicators, $y_{ij}$, have already been computed, then
you can specify them with the \code{y.train} argument; this results in
the default $\mu_0=0$ (which you can over-ride with the
\code{binaryOffset} argument).  However, it is likely that the
indicators would need to be constructed, so for convenience, you can
specify $(s_i, \delta_i)$ by the arguments \code{times} and
\code{delta} respectively.  In this case, the default value of $\mu_0$
assumes that the times follow the Exponential distribution and the
covariates are not involved, i.e.,
$\mu_0=\Phi^{-1}\wrap{1-\exp\wrap[()]{\frac{\sum_i \delta_i}{\sum_i
      s_i}}}$
(which you can over-ride with the \code{binaryOffset} argument).  For
BART with continuous outcomes, typically the outcome is centered and
$\mu_0$ is taken to be $\b{y}$.  While centering can be helpful for
small samples with \citet{AlbeChib93}, it is unnecessary for moderate
to large samples because of the flexibility of $f$ (for
\citet{HolmHeld06} with Logistic latents which have heavier tails,
centering is unnecessary even for small samples so $\mu_0$ is fixed at
zero for \code{type='lbart'}).

So just like in the Cox model case, we have to construct quantities of
interest with BART for survival analysis.  In discrete-time survival
analysis, the probability of an event in an interval is the foundation
which essentially replaces
the instantaneous hazard from continuous-time survival analysis: 
$p(t,\bm{x})=\Phi(\mu_0+f(t,\bm{x}))$.  Now, the survival function
is constructed as follows: 
$ S(t_{(j)}|\bm{x})=Pr(T>t_{(j)}|\bm{x})=\prod_{l=1}^j (1-p(t_{(l)},\bm{x}))$.

Survival data pairs $(s, \delta)$ are converted to indicators by the
helper function \code{surv.pre.bart} which is called automatically by
\code{surv.bart} if \code{y.train} is not provided.
\code{surv.pre.bart} returns a list which contains \code{y.train} for
the indicators; \code{tx.train} for the covariates corresponding to
\code{y.train} for training $f(t, \bm{x})$ (which includes time in the
first column, and the rest of the covariates afterward, if any, i.e.,
rows of $\wrap{t, \bm{x}}$, hence the name \code{tx.train} to
distinguish it from the original \code{x.train}); \code{tx.test} for
the covariates to predict $f(t, \bm{x})$ rather than to train;
\code{times} which is the grid of ordered distinct time points;
\code{K} which is the length of \code{times}; and \code{binaryOffset}
which is $\mu_0$.  Here is a very simple example of a data set with
three observations and no covariates re-formatted for display 
(no covariates is an interesting special case but we
will discuss the more common case with covariates further below).
\begin{verbatim}
times <- c(2.5, 1.5, 3.0)
delta <- c(  1,   1,   0)
surv.pre.bart(times=times, delta=delta)

$y.train  $tx.train  $tx.test  $times    $K     $binaryOffset 
[1]              t          t  [1]       [1] 3  [1] -0.6791459
    0     [1,] 1.5   [1,] 1.5    1.5  
    1     [2,] 2.5   [2,] 2.5    2.5
    1     [3,] 1.5   [3,] 3.0    3.0
    0     [4,] 1.5
    0     [5,] 2.5
    0     [6,] 3.0
\end{verbatim}

Here is a schematic of the input and output for the \code{surv.pre.bart} function.\\
\code{pre <- surv.pre.bart(times, delta, x.train, x.test=x.train)}   \\
\begin{align*}
%\mbox{\code{pre} is a list as such} & \\
\mbox{\code{pre} is a list with the matrix \code{pre\$tx.train}\ } & \mbox{\& \code{pre\$y.train} which is a vector} \\
\wrap{\begin{array}{cc}
t_{(1)} & \bm{x}_1 \\
\vdots & \vdots \\
t_{({n_1})} & \bm{x}_1 \\
\vdots & \vdots \\
t_{(1)} & \bm{x}_N \\
\vdots & \vdots \\
t_{({n_N})} & \bm{x}_N \\
\end{array}} & \wrap{\begin{array}{c}
y_{11}=0 \\
\vdots \\
y_{1n_1}=\delta_1 \\
\vdots \\
y_{N1}=0 \\
\vdots \\
y_{Nn_N}= \delta_N \\
\end{array}}\\
\end{align*}
For \code{pre\$tx.test}, ${n_i}$ is replaced by $K$ which is very
helpful so that each subject contributes an equal number of settings
for programmatic convenience and noninformative estimation, i.e., if
high-risk subjects with earlier events did not appear beyond their
event, then estimates of survival for latter times would be biased
upward.  For other outcomes besides time-to-event, we provide two
matrices of covariates, \code{x.train} and \code{x.test}, where
\code{x.train} is for training and \code{x.test} is for validation.
However, due to the variable $n_i$ for time-to-event outcomes, we
generally provide two arguments as follows: \code{x.train,
  x.test=x.train} where the former matrix will be expanded by
\code{surv.pre.bart} to $\sum_{i=1}^N n_i$ rows for training
$f(t, \bm{x})$ while the latter matrix will be expanded to
$N \times K$ rows for $f(t, \bm{x})$ estimation only.  If you still
need to perform validation, then you can make a separate call to the
\code{predict} function.

N.B.\ the argument \code{ndpost=M} is the length of the chain to be
returned and the argument \code{keepevery} is used for thinning, i.e.,
return \code{M} observations where \code{keepevery} are culled in
between each returned value.  For BART with time-to-event outcomes,
the default is \code{keepevery=10} (rather than \code{keepevery=1} for
other outcomes) since the grid of time points creates data set
observations of order $N \times K$ which have a tendency towards higher
auto-correlation, therefore, making thinning more necessary.  To avoid
unnecessarily enlarged data sets, it is often prudent to coarsen the
time axis appropriately, i.e., re-scale from days to weeks or months.
You can coarsen automatically by supplying the optional \code{K}
argument to coarsen the times to a grid of time quantiles: \code{1/K},
\code{2/K}, ..., \code{K/K} (not to be confused with the \code{k}
argument which is a prior parameter for the distribution of the leaf
terminal values).

Here is a schematic of the input and output for the \code{surv.bart}
function for serial computation and \code{mc.surv.bart} for parallel
computation.\\  
\code{set.seed(99)}\\
\code{post=surv.bart(x.train, times=times, delta=delta, x.test=x.train, ndpost=M)} or \\
\code{post=mc.surv.bart(x.train, times=times, delta=delta, x.test=x.train, ndpost=M, mc.cores=C, seed=99)}  \\
\begin{align*}
\mbox{Input vector \code{times} with \code{K} distinct values and \code{x.train}:\ } &
\wrap{\begin{array}{c}
\bm{x}_{1} \\
\bm{x}_{2} \\
\vdots \\
\bm{x}_{N} \\
\end{array}} \mbox{\ or\ } \bm{x}_i \\
\mbox{Output \code{post} of type \code{survbart} which is essentially a list of } & \\
\mbox{objects including the matrix: \code{post\$surv.test}:\ } \hat{S}_m(t_{(j)}, \bm{x}_i) & \\
\wrap{\begin{array}{ccccccc}
\scriptstyle \hat{S}_1(t_{(1)}, \bm{x}_1)& \scriptstyle \. & \scriptstyle \hat{S}_1(t_{(K)}, \bm{x}_1) & \scriptstyle \.
& \scriptstyle \hat{S}_1(t_{(1)}, \bm{x}_N)& \scriptstyle \. & \scriptstyle \hat{S}_1(t_{(K)}, \bm{x}_N) \\
\vdots & \vdots & \vdots & \vdots & \vdots & \vdots & \vdots \\
\scriptstyle \hat{S}_M(t_{(1)}, \bm{x}_1)& \scriptstyle \. & \scriptstyle \hat{S}_M(t_{(K)}, \bm{x}_1) & \scriptstyle \.
& \scriptstyle \hat{S}_M(t_{(1)}, \bm{x}_N)& \scriptstyle \. & \scriptstyle \hat{S}_M(t_{(K)}, \bm{x}_N) \\
\end{array}} & \\
\end{align*}

Here is a schematic of the input and output for the \code{predict.survbart} function.\\
\code{pred <- predict(post, {pre\$tx.test}, mc.cores=C)} \\
\begin{align*}
\mbox{Input: \code{x.test}\ } &
\wrap{\begin{array}{c}
\bm{x}_{1} \\
\bm{x}_{2} \\
\vdots \\
\bm{x}_{Q} \\
\end{array}} \mbox{\ or\ } \bm{x}_i \\
\mbox{Output:\ } & \mbox{ \code{pred} of type \code{survbart} with \code{pred\$surv.test}:\ } 
{\hat{S}_m(t_{(j)}, \bm{x}_i)} \\
& \wrap{\begin{array}{ccccccc}
\scriptstyle \hat{S}_1(t_{(1)}, \bm{x}_1)& \scriptstyle \. & \scriptstyle \hat{S}_1(t_{(K)}, \bm{x}_1) & \scriptstyle \.
& \scriptstyle \hat{S}_1(t_{(1)}, \bm{x}_Q)& \scriptstyle \. & \scriptstyle \hat{S}_1(t_{(K)}, \bm{x}_Q) \\
\vdots & \vdots & \vdots & \vdots & \vdots & \vdots & \vdots \\
\scriptstyle \hat{S}_M(t_{(1)}, \bm{x}_1)& \scriptstyle \. & \scriptstyle \hat{S}_M(t_{(K)}, \bm{x}_1) & \scriptstyle \.
& \scriptstyle \hat{S}_M(t_{(1)}, \bm{x}_Q)& \scriptstyle \. & \scriptstyle \hat{S}_M(t_{(K)}, \bm{x}_Q) \\
\end{array}} \\
\end{align*}

As previously noted, BART does not directly provide a summary of the
effect of a single covariate, or a subset of covariates, on the
outcome.  For survival analysis, we use Friedman's partial dependence
function \citep{Frie01} with BART to summarize the marginal effect due
to a subset of the covariates, $(t, \bm{x}_S)$, by aggregating over
the complement of covariates, $\bm{x}_C$, i.e.,
$\bm{x} =\wrap{\bm{x}_S,\bm{x}_C}$.  This marginal dependence function
is defined by fixing $(t, \bm{x}_S)$ while aggregating over the
observed settings of the complement covariates in the cohort:
$f(t, \bm{x}_S)={N^{-1}}\sum_{i=1}^N f(t, \bm{x}_S,\bm{x}_{iC})$.  For
survival analysis, the $f$ function is often not directly of interest;
rather, the survival function is more readily interpretable:
$S(t, \bm{x}_S) = {N^{-1}} \sum_{i=1}^N S(t, \bm{x}_S,\bm{x}_{iC})$.
Other marginal functions can be obtained in a similar fashion.
Estimates can be derived via functions of
the posterior samples such as means, quantiles, e.g.,
$\hat{S}(t, \bm{x}_S) = {M^{-1}} {N^{-1}} \sum_{m=1}^M \sum_{i=1}^N
S_m(t, \bm{x}_S,\bm{x}_{iC})$
where $m$ indexes posterior samples.  Friedman's partial dependence
function is a concept which is very flexible.  So flexible that as of
yet, we are unable to provide abstract functional support in the {\bf
  BART} package; rather, we provide examples of the many practical
uses in the \code{demo} directory.

Here we present an example that is available in the {\bf BART} package.\\
\code{system.file('demo/lung.surv.bart.R', package='BART')}.  The
North Central Cancer\\ Treatment Group surveyed 228 advanced lung cancer
patients \citep{LoprLaur94}.  This data can be found in the
\code{lung} data set.  The study focused on prognostic variables.
Patient responses were paired with a few clinical variables.  We control
for age, gender and Karnofsky performance score as rated by their
physician.  We compare the survival for males and females with
Friedman's partial dependence function; see Figure~\ref{lung}.  We
also analyze this data set with Logistic
latents and the results are quite similar (not shown):
\code{system.file('demo/lung.surv.lbart.R', package='BART')}.
Furthermore, we perform convergence diagnostics on the chain:\\
\code{system.file('demo/geweke.lung.surv.bart.R', package='BART')}.
\begin{comment}
\begin{figure}[t!]
\begin{center}
\includegraphics[scale=0.45]{figures/geweke-lung-surv-bart.pdf}
\caption{Convergence diagnostics of the advanced lung cancer example.}
\end{center}
\end{figure}
\end{comment}

\subsection{Survival analysis and the concordance probability}

The concordance probability \citep{GoneHell05} is a measure of the
discriminatory ability of survival analysis analogous to the
area under the receiver operating characteristic curve for binary
outcomes.  A simple analytic example with the Exponential
distribution is as follows.
\begin{align*}
t_i & ~ \Exp{\lambda_i} \where i \in \{1, 2\} \\
\P{t_1<t_2} & = \int_0^{\infty} \int_0^{t_2} \lambda_2 \e{-\lambda_2 t_2} 
              \lambda_1 \e{-\lambda_1 t_1} \d{t_1}  \d{t_2} 
              = \frac{\lambda_1}{\lambda_1+\lambda_2} \\
1-\P{t_1>t_2} & = 1- \frac{\lambda_2}{\lambda_1+\lambda_2}
  = \frac{\lambda_1}{\lambda_1+\lambda_2}\\
\end{align*}
Notice that the concordance is symmetric with respect to $t_1$ and $t_2$.

We can make a similar calculation based on our BART survival analysis model.
First, we calculate $ \P{s_1 < s_2} $.
\begin{align*}
\P{s_1 < s_2}  = &\P{s_1=t_{(1)}, s_2>t_{(1)}}+\\
&\P{s_1=t_{(2)}, s_2>t_{(2)}|s_1>t_{(1)}, s_2>t_{(1)}}\P{s_1>t_{(1)}, s_2>t_{(1)}}+\.\\
= & \sum_{j=1}^K \P{s_1=t_{(j)}, s_2>t_{(j)}|s_1>t_{(j-1)}, s_2>t_{(j-1)}}
    \P{s_1>t_{(j-1)}, s_2>t_{(j-1)}} \\
= & \sum_{j=1}^K p_{1j} q_{2j} S_{1}(t_{(j-1)}) S_{2}(t_{(j-1)}) \\
\end{align*}
Now, we calculate the mirror image relationship.
\begin{align*}
1- \P{s_1 > s_2} = &  1-\sum_{j=1}^K q_{1j} p_{2j} S_{1}(t_{(j-1)}) S_{2}(t_{(j-1)}) \\
 = &1-\sum_{j=1}^K (1-p_{1j}) (1-q_{2j})  S_{1}(t_{(j-1)}) S_{2}(t_{(j-1)}) \\
 = & 1-\sum_{j=1}^K (1-p_{1j}-q_{2j}+p_{1j}q_{2j})  S_{1}(t_{(j-1)}) S_{2}(t_{(j-1)}) \\
= & 1-\sum_{j=1}^K p_{1j}q_{2j}  S_{1}(t_{(j-1)}) S_{2}(t_{(j-1)})
-\sum_{j=1}^K (q_{1j}-q_{2j})  S_{1}(t_{(j-1)}) S_{2}(t_{(j-1)}) \\
\end{align*}
However, note that these probabilities are not symmetric.  Therefore,
we enforce symmetry as follows.
\begin{align*}
  \kappa_{s_1, s_2} = & 0.5 \wrap[()]{\P{s_1 < s_2}+1- \P{s_1 > s_2}} 
   =  0.5\wrap{1  -\sum_{j=1}^K (q_{1j}-q_{2j}) S_{1}(t_{(j-1)}) S_{2}(t_{(j-1)})}  \\
\end{align*}
See the concordance probability example at 
\code{system.file('demo/concord.surv.bart.R', package='BART')}.

\subsection{Competing risks with BART}

Competing risks are supported by the function \code{crisk.bart} for
serial computation and\\ \code{mc.crisk.bart} for parallel computation.
% As has been described, dealing with dependent censoring in the
% recurrent events framework is challenging
% \citep{CookLawl97,GhosLin00,WangQin01,GhosLin03}.  Our approach is to
% combine the recurrent events and competing risks paradigms.
Typically, competing risks \citep{FineGray99,KalbPren02} deal with events
which are mutually exclusive, say, death from cardiovascular disease
vs.\ death from other causes, i.e., a patient experiencing one of the
events is incapable of experiencing another.  % Our application is
% slightly different in that we have two events, death and hospital
% admission, which are not mutually exclusive.  Suffering death prevents a patient from experiencing a
% future hospital admission, but the converse is not true, i.e., a
% hospital admission does not prevent a future death.  So, technically,
% we have what are termed semi-competing events, yet the competing events framework is sufficient
% for our needs.  For the simplicity of this exposition, we assume
% that the hospital admissions are instantaneous, i.e., hospital stays
% are of length zero. 
% (from here on we use the
% term competing risks loosely to represent either competing events or
% semi-competing events)
%  Nevertheless, competing risks offers us an
% approach to handle dependent censoring.
We adopt the subdistribution concept of \citet{FineGray99} for
competing risks.  Let's suppose we have two kinds of events: events of
kind 1, death from cause 1 which is the cause of interest, and events
of kind 2, death from cause 2 which is any other cause.  The
distribution function of an event time is
$F(t, \bm{x})=G_1(t, \bm{x})+G_2(t, \bm{x})$ where
$G_1(t, \bm{x})=p F_1(t, \bm{x})$ and
$G_2(t, \bm{x})=(1-p) F_2(t, \bm{x})$.  $F_1$ and $F_2$ are
distribution functions, i.e., $F_h(\infty, \bm{x})=1$.  While $G_1$
and $G_2$ are subdistribution functions, i.e.,
$G_h(\infty, \bm{x})<1$.  % integrate to $p$ and $1-p$ respectively
\citet{FineGray99} model the subdistribution functions rather than the
distribution functions; and we do the same.  But, here we part ways
with \citet{FineGray99} since they assume linear proportionality and
the Exponential distribution while we impose neither precarious
restrictive assumption.

To accomodate competing risks, we adapt our notation slightly:
$(s_i, \delta_i)$ where $\delta_i=1$ for kind 1 events, $\delta_i=2$
for kind 2 events, or $\delta_i=0$ for censoring times.  We create a
single grid of time points for the ordered distinct times based on
either kind of event or censoring:
$0=t_{(0)}< t_{(1)}<\cdots < t_{(J)} < \infty$.  % To accomodate
% competing risks, we adapt our notation slightly: $(s_i, \delta_i)$ are
% death, $\delta_i=1$, or censoring, $\delta_i=0$, time; $t_{ik}$ are
% hospital admissions; $i=1, \dots, n$ indexes patients;
% $j=1, \dots, n_i=\arg \max_j \wrap{ t_{(j)}\le s_i } $ indexes time
% points on the grid; and $k=1, \dots, N_i=N_i(s_i)$ indexes hospital
% admissions.  Now, an analysis of recurrent events proceeds as follows.
We model the probability for an event of kind 1,
$p_1(t_{(j)}, \bm{x}_i)$, and an event of kind 2 conditioned on subject $i$ being alive at time
$t_{(j)}$, $p_2(t_{(j)}, \bm{x}_i)$.  % So, we can estimate the
% survival function and the cumulative incidence functions as follows.
% \begin{align*}
% S(t, \bm{x}_{i}) & = 1-F(t, \bm{x}_{i}) = \prod_{j=1}^k (1-p_1(t_{(j)}, \bm{x}_i))(1-p_{2}(t_{(j)}, \bm{x}_i)) \where k=\arg \max_j \wrap{t_{(j)} \le t}  \\
% F_1(t, \bm{x}_{i}) & = \int_0^t S(u-, \bm{x}_{i}) \lambda_1(u, \bm{x}_{i}) \mathrm{d}u 
% = \sum_{j=1}^k S(t_{(j-1)}, \bm{x}_{i}) p_{1}(t_{(j)}, \bm{x}_i) \\
% F_2(t, \bm{x}_{i}) &  = \int_0^t S(u-, \bm{x}_{i}) \lambda_2(u, \bm{x}_{i}) \mathrm{d}u 
% = \sum_{j=1}^k S(t_{(j-1)}, \bm{x}_{i}) (1-p_{1}(t_{(j)}, \bm{x}_i)) p_{2}(t_{(j)}, \bm{x}_i) \\
% \end{align*}
Now, as before, we create event indicators.  
% First, we develop event 1, then event 2 
% which is necessarily conditioned on patient $i$ being alive at time $t_{(j)}$.
\begin{align*}
y_{1ij} & = \I{\delta_i=1} \I{j=n_i} \where j=1, \dots, n_i\\ %  \where y_{1i0}=0 \\
%y_{1ij}|(y_{1i\,j-1}=0, p_{1ij}) & ~ \B{p_{1ij}} \where p_{1i0}=0 \\
y_{1ij}|p_{1ij} & ~ \B{p_{1ij}} \\ % \where p_{1i0}=0 \\
 p_{1ij} & = \Phi(\mu_1+f_1(t_{(j)}, {\bm{x}}_{i}) \where f_1 \prior \mathrm{BART} \\
y_{2ij} & = \I{\delta_i=2} \I{j=n_i}  \where j=1, \dots, n_i-y_{1in_i} \\ % \where y_{2i0}=0 \\
y_{2ij}|p_{2ij} & ~ \B{p_{2ij}} \\ %\where p_{2i0}=0\\
%y_{2ij}|(y_{1i\,j-1}=y_{1ij}=0, p_{2ij}) & ~ \B{p_{2ij}} \where p_{2i0}=0\\
 p_{2ij} & = \Phi(\mu_2+f_2(t_{(j)}, {\bm{x}}_{i}) \where f_2 \prior \mathrm{BART} \\
\end{align*}

Based on this BART framework, we can estimate the survival function and
the cumulative incidence functions as follows.
\begin{align*}
S(t, \bm{x}_{i}) & = 1-F(t, \bm{x}_{i}) = \prod_{j=1}^k (1-p_{1ij})(1-p_{2ij}) \where k=\arg \max_j \wrap{t_{(j)} \le t}  \\
F_1(t, \bm{x}_{i}) & = \int_0^t S(u-, \bm{x}_{i}) \lambda_1(u, \bm{x}_{i}) \mathrm{d}u 
= \sum_{j=1}^k S(t_{(j-1)}, \bm{x}_{i}) p_{1ij} \\
F_2(t, \bm{x}_{i}) &  = \int_0^t S(u-, \bm{x}_{i}) \lambda_2(u, \bm{x}_{i}) \mathrm{d}u 
= \sum_{j=1}^k S(t_{(j-1)}, \bm{x}_{i}) (1-p_{1ij}) p_{2ij} \\
\end{align*}
% There are pros and cons to this approach.  We are modeling the
% subdistribution function rather than the distribution function;
% however, for applications like ours, inference with respect to
% the subdistribution function is probably adequate.
The returned object of type \code{criskbart} from \code{crisk.bart}
or \code{mc.crisk.bart} provides the cumulative incidence functions
and survival corresponding to \code{x.test} as follows: $F_1$ is
\code{cif.test}, $F_2$ is \code{cif.test2} and $S$ is \code{surv.test}.

Here, we present the Mayo Clinic liver transplant waiting list data
from 1990-1999 with $N=815$ patients.  During the study period, the
liver transplant organ allocation policy was flawed.  Blood type is an
important matching factor to avoid organ rejection.  Donor livers from
subjects with blood type O can be used by patients with A, B, AB or O
blood types; whereas a donor liver from the other types will only be
transplanted to a matching A, B or AB recipient.  Therefore, type O
subjects on the waiting list were at a disadvantage since the pool of
competitors was larger for type O donor livers.  This data is of
historical interest and provides a useful example of competing risks,
but it has little relevance today.  Current liver transplant policies
have evolved and now depend on each individual patient's risk/need
which are assessed and updated regularly while a patient is on the
waiting list.  However, there still remains an acute shortage of donor
livers today.  The \code{transplant} data set is provided by the {\bf
  BART} R package as is this example:
\code{system.file('demo/liver.crisk.bart.R', package='BART')}.  We
compare the nonparametric Aalen-Johansen competing risks estimator
with BART for the transplant event of type O patients which are in
general agreement; see Figure~\ref{liver-transplant}.

\begin{comment}
\begin{figure}[t!]
\begin{center}
\includegraphics[scale=0.45]{figures/liver-AJ.pdf}
\end{center}
\caption{Liver transplant competing risks for type O patients and the
Aalen-Johansen estimator available in {\bf survival} R package.}
\end{figure}
\end{comment}

\subsection{Recurrent events with BART}%\label{recur}

The {\bf BART} package supports recurrent events with
\code{recur.bart} for serial computation and \code{mc.recur.bart} for
parallel computation.  BART provides much desired flexibility in
modeling the dependence of recurrent events on covariates.  Consider
data in the form: $\delta_i, s_i, \bm{t}_{i}, \bm{u}_{i}, \bm{x}_i(t)$
where $i=1, \dots, n$ indexes subjects; $s_i$ is the end of the
observation period (death, $\delta_i=1$, or censoring, $\delta_i=0$);
$N_i$ is the number of events during the observation period;
$\bm{t}_{i}=\wrap{t_{i1}, \dots, t_{iN_i }}$ and $t_{ik}$ is the event
start time of the $k$th event (let $t_{i0}=0$);
$\bm{u}_{i}=\wrap{u_{i1}, \dots, u_{iN_i}}$ and $u_{ik}$ is the event
end time of the $k$th event (let $u_{i0}=0$); and $\bm{x}_i(t)$ is a
vector of time-dependent covariates.  Both start and end times of
events are necessary to define risk set eligibility for events of
stochastic duration like readmissions since patients currently
hospitalized cannot be readmitted.  For instantaneous events (or
roughly instantaneous events such as emergency department visits with
time measured in days), the end times can be simply ignored.

We denote the $J$ collectively distinct event start and end times for
all subjects by $0<t_{(1)}< \dots< t_{(J)}<\infty$ thus taking
$t_{(j)}$ to be the $j^{th}$ order statistic among distinct
observation times and, for convenience, $t_{(j')}=0 \where j' \le 0$
(note that $t_{(j)}$ are constructed from all event start/end times
for all subjects,
but they may be a censoring time for any given subject).  Now consider
binary event indicators $y_{ij}$ for each subject $i$ at each distinct
time $t_{(j)}$ up to %and including
the subject's last observation time $t_{(n_i)} \le s_{i}$ with
$n_i=\arg \max_j \wrap{ t_{(j)}\leq s_{i}}$, i.e.,
$y_{i1}, \dots, y_{in_i} \in \{0, 1\}$.  We then denote by $p_{ij}$
the probability of an event at time $t_{(j)}$ conditional on
%$\wrap[()]{t_{(j)}, N_i(t_{(j)}-), v_i(t_{(j)}), \bm{x}_i(t_{(j)})}$.
$\wrap[()]{t_{(j)}, \tilde{\bm{x}}_i(t_{(j)})}$ where
$\tilde{\bm{x}}_i(t_{(j)})=\wrap[()]{ N_i(t_{(j-1)}), v_i(t_{(j)}),
  \bm{x}_i(t_{(j)})}$.
Let $N_i(t-) \equiv \lim\limits_{s\uparrow t} N_i(s)$ be the counting
process of events for subject $i$ just prior to time $t$ and we also
note that $N_i=N_i(s_i)$.  Let $v_i(t)=t-u_{N_i(t-)}$ be the sojourn time
for subject $i$, i.e., time since last event, if any.  Notice that we
can replace $N_i(t_{(j)}-)$ with $N_i(t_{(j-1)})$ since, by construction,
the state of information available at time $t_{(j)}-$ is the same as
that available at $t_{(j-1)}$.
%For simplicity of presentation, we let
%$\tilde{\bm{x}}_i(t_{(j)})=\wrap[()]{ N_i(t_{(j-1)}), v_i(t_{(j)}),
%  \bm{x}_i(t_{(j)})}$.
Assuming a constant intensity and constant covariates,
$\tilde{\bm{x}}_i(t_{(j)})$, in the interval $(t_{(j-1)} , t_{(j)}]$, we
define the cumulative intensity process as:
\begin{align}\label{cum-int}
\Lambda(t_{(j)}, \tilde{\bm{x}}_i(t_{(j)})) &
= \int_0^{ t_{(j)}} \mathrm{d} \Lambda(t, \tilde{\bm{x}}_i(t)) 
% = \sum_{j'=1}^j \int_{t_{(j'-1)}}^{t_{(j')}} \Pr{\mathrm{d} N_i(t)=1\left|\right. t, \tilde{\bm{x}}_i(t)} \\
 = \sum_{j'=1}^j \Pr{N_i(t_{(j')})-N_i(t_{(j'-1)})=1\left|\right. t_{(j')}, \tilde{\bm{x}}_i(t_{(j')})} 
  =\sum_{j'=1}^j p_{ij'} 
\end{align}
where these $p_{ij}$ are currently unspecified and we 
provide their definition later \eqref{recur}.
Note: the terms ``intensity'' and ``hazard'' are generally interchangeable.
% (some use the terms ``intensity'' and ``hazard'' interchangeably, but
% we prefer ``intensity'' for recurrent events to avoid confusion with
% standard time to a single event survival analysis which is more often
% associated with ``hazard'').  
% So, in our framework, we estimate the cumulative intensity process as follows.
% \begin{align*}
% \widehat{\Lambda}(t_{(j)}, \tilde{\bm{x}}_i(t_{(j)})) =\sum_{i=1}^j \hat{p}_{ij} 
% \end{align*}

With absorbing events such as mortality there is no concern about the
conditional independence of future events because there will never be
any.  Conversely, with recurrent events, there is a valid concern.  Of
course, conditional independence can be satisfied by conditioning on
the entire event history, denoted by $N_i(s) \where 0\le s < t$.
However, conditioning on the entire event history is impractical.
Rather, we condition on both $N_i(t-)$ and $v_i(t)$ to satisfy any
concern of conditional independence.

We now write the model for $y_{ij}$ as a nonparametric probit
regression of $y_{ij}$ on
$\wrap[()]{t_{(j)}, \tilde{\bm{x}}_i(t_{(j)})}$ tantamount to
parametric models of discrete-time intensity
\citep{Thom77,ArjaHaar87,FahrTutz94,Fahr98}.  Specifically, with
temporal data converted from
$\delta_i, s_i, \bm{t}_{i}, \bm{u}_{i}, \bm{x}_i(t)$ to a sequence of
longitudinal binary events as follows:
$y_{ij}\ =\ \max_k \I{t_{ik}=t_{(j)}}$. %,\ j=1, \dots, n_i
% But, what is the range of $j$ in the equation above?  As we shall see,
% there is a caveat in translating this into the likelihood
% which we address by a vignette. 
%(the likelihood is presented below in equation~\eqref{L}).  
However, note that the definition of $j$ is currently unspecified.
To understand the impetus of the range of $j$, let's look at an example.

Suppose that we have two subjects with the following values:
\begin{align*}
N_1=2, s_1=9, {t}_{11}=3, {u}_{11}=7, {t}_{12}=8, {u}_{12}=8 & 
\Rightarrow y_{11}=1, y_{12}=y_{13}=0, y_{14}=1, y_{15}=0\ (2.3)\\
N_2=1, s_2=12, {t}_{21}=4, {u}_{21}=7 & 
\Rightarrow y_{21}=0, y_{22}=1, y_{23}=y_{24}=y_{25}=y_{26}=0 \notag
\end{align*}
which creates the grid of times $(3, 4, 7, 8, 9, 12)$.  For subject 1
(2), notice that $y_{12}=y_{13}=0$ ($y_{23}=0$) as it should be since
no event occurred at times 4 or 7 (7).  However, there were no events 
since their first event had not ended yet, i.e., these subjects are not
chronologically at risk for an event and, therefore, no corresponding random 
behavior contributed to the likelihood.  The {\bf BART} package provides the
\code{recur.pre.bart} function which you can use to construct these
data sets.  Here is a short demonstration of its capabilities adapted
from \code{demo/data.recur.pre.bart.R} (re-formatted for display purposes).
%\newpage
\begin{verbatim}
> library(BART)
> times <- matrix(c(3, 8, 9, 4, 12, 12), nrow=2, ncol=3, byrow=TRUE)
> tstop <- matrix(c(7, 8, 0, 7,  0,  0), nrow=2, ncol=3, byrow=TRUE)
> delta <- matrix(c(1, 1, 0, 1,  0,  0), nrow=2, ncol=3, byrow=TRUE)
> recur.pre.bart(times=times, delta=delta, tstop=tstop)
$K   $times   $y.train  $tx.train     $tx.test   
[1]  [1]      [1]             t v N         t v N
 6    3        1        [1,]  3 3 0   [1,]  3 3 0
      4        1        [2,]  8 5 1   [2,]  4 1 1
      7        0        [3,]  9 1 2   [3,]  7 4 1
      8        0        [4,]  3 3 0   [4,]  8 5 1
      9        1        [5,]  4 4 0   [5,]  9 1 2
     12        0        [6,]  8 4 1   [6,] 12 4 2
               0        [7,]  9 5 1   [7,]  3 3 0
               0        [8,] 12 8 1   [8,]  4 4 0
                                      [9,]  7 3 1
                                     [10,]  8 4 1
                                     [11,]  9 5 1
                                     [12,] 12 8 1
\end{verbatim}
Notice that \code{\$tx.test} is not limited to the same time points as
\code{\$tx.train}, i.e., we often want/need to estimate $f$
at counter-factual values not observed in the data so each subject
contributes an equal number of evaluations for estimation purposes.
 
It is now clear that the $y_{ij}$
which contribute to the likelihood are those that correspond to $j
\in R_i(t_{(j)})$ which is the risk set for subject
$i$ at time $t_{(j)}$ that either contains $j$
or is empty.  We formally define the risk set as $R_i(t_{(j)})
= \wrap[\{\}]{j : ( j \in \{1, \dots, n_i\} ) \cap
  \wrap[()]{\cap_{k=1}^{N_i} \{t_{(j)} \not\in (t_{ik}, u_{ik}) \} }
}$.
% and we denote the binary outcomes for subject $i$ by the
% longitudinally ordered vector
% $\bm{y}_i=\wrap{y_{ij}} \where j \in R_i(t_{(j)})$.

Putting it all together, we arrive at the following recurrent events
discrete-time model.
\begin{align}\label{recur}
y_{ij}|p_{ij} & ~  \B{p_{ij}} \where i=1, \dots, n; j \in R_i(t_{(j)}) \notag \\
p_{ij}|f & = \Phi(\mu_{ij}),\ \mu_{ij} =\mu_0+f(t_{(j)}, \tilde{\bm{x}}_i(t_{(j)})) \\ % & \If  j \in R_i \\ 
f & \prior \mathrm{BART} \notag
\end{align}
For computational efficiency, we carry out the probit regression via
truncated Normal latent variables $z_{ij}$ to reduce it to a
continuous outcome BART model like so \citep{AlbeChib93} (this 
default can be over-ridden to utilize \citet{HolmHeld06} Logistic
latents by specifying \code{type='lbart'}).
\begin{align*}
z_{ij}|y_{ij},f & ~ \begin{cases} 
\N{ \mu_{ij}}{1}\I{-\infty, 0} & \If y_{ij}=0 \\
\N{ \mu_{ij}}{1}\I{0, \infty} & \If y_{ij}=1 \\
\end{cases}  \\
\end{align*}
% Consequently, we have the following Bernoulli likelihood
% $\prod_{i=1}^n\prod_{j \in R_i} p_{ij}^{y_{ij}} (1-p_{ij})^{1-y_{ij}}$.
For binary data, $\mu_0=\Phi^{-1}({p_0})$ can be used for centering
the latents around the probability of an event $p_0$.
%where $\Phi^{-1}(.) $ is the inverse standard Normal
%cumulative distribution function.  
For recurrent event data, we can
similarly center the latents by assuming the times of recurrent
events follow the Exponential distribution and the covariates,
$\tilde{\bm{x}}$, have no impact, i.e.,
$\mu_0=\Phi^{-1}\!\wrap[()]{1-\exp\wrap[()]{\frac{\sum_i N_i}{\sum_i s_{i}}}}$.

With the data prepared as described in the above example, the BART
model for binary data treats the probability of an event within an
interval as a nonparametric function of time, $t$, and 
covariates, $\tilde{\bm{x}}(t)$. Conditioned on the data, BART
provides samples
from the posterior distribution of $f$. For any $t$ and $\tilde{\bm{x}}(t)$,
we obtain the posterior distribution of
$p(t, \tilde{\bm{x}}(t))  =\Phi(\mu_0+f(t, \tilde{\bm{x}}(t)))$.

For the purposes of recurrent events survival analysis, we are
typically interested in estimating the cumulative intensity function
as presented in formula \eqref{cum-int}.  With these estimates, one
can accomplish inference from the posterior via means, quantiles or
other functions of $ p(t, \tilde{\bm{x}}_i(t))$ or
$\Lambda(t,\tilde{\bm{x}}(t))$ as needed such as the relative
intensity, i.e.,
$RI(t,\tilde{\bm{x}}_n(t),\tilde{\bm{x}}_d(t))=\frac{p(t,\tilde{\bm{x}}_n(t))}{p(t,\tilde{\bm{x}}_d(t))}$\label{risk}
where $\tilde{\bm{x}}_n(t) $ and $\tilde{\bm{x}}_d(t) $ are two
settings we wish to compare like two
treatments. %treatment~1 vs.\ treatment~2.
% Also, the cumulative intensity could be calculated for arbitrary intervals
% on the grid, i.e., $\int_{t_{(a)}}^{ t_{(b)}} \mathrm{d} \Lambda(t, \tilde{\bm{x}}_i(t)) = \sum_{j=a}^b p_{ij}$.

An interesting example of recurrent events involves a clinical trial
conducted by the Veterans Administration Cooperative Urological
Research Group.  In this study, all patients had superficial bladder
tumors when they entered the trial.  These tumors were removed
transurethrally and patients were randomly assigned to one of three
treatments: placebo, thiotepa or pyridoxine (vitamin B6).  Many
patients had multiple recurrences of tumors during the study and new
tumors were removed at each visit.  For each patient, their recurrence
time, if any, was measured from the beginning of treatment.  There
were 118 patients enrolled but only 116 were followed beyond time zero
and contribute information.  This data set is loaded by
\code{data(bladder)} and the data frame of interest is
\code{bladder1}.  This data set is analyzed by
\code{system.file('demo/bladder.recur.bart.R', package='BART')}.  In
Figure~\ref{RI-Th-Pl}, notice that the relative intensity calculated
by Friedman's partial dependence function favors thiotepa over placebo
from roughly 3 to 18 months and afterward they are about equal, but
the 95\% credible intervals are wide throughout.  Similarly, the
relative intensity calculated by Friedman's partial dependence
function favors thiotepa over vitamin B6 from roughly 3 to 18 months
and afterward they are about equal, but the 95\% credible intervals
are wide throughout; see Figure~\ref{RI-Th-B6}.  And, finally, vitamin
B6 is no better than placebo, and possibly worse, but the 95\%
credible intervals are wide; see Figure~\ref{RI-B6-Pl}.

\bibliography{references}

\begin{figure}[t!]
\begin{center}
\includegraphics[scale=0.45]{figures/lung.pdf}
\caption{Advanced lung cancer example: Friedman's partial dependence
  function with 95\% credible intervals: males (blue) vs.\ females
  (red). \label{lung}}
\end{center}
\end{figure}

\begin{figure}[t!]
\begin{center}
\includegraphics[scale=0.45]{figures/liver-BART.pdf}
\end{center}
\caption{Liver transplant competing risks for type O patients
  estimated by BART and Aalen-Johansen. \label{liver-transplant}}
\end{figure}

\begin{figure}[t!]
\begin{center}
\includegraphics[scale=0.45]{figures/RI-Th-Pl.pdf}
\end{center}
\caption{Relative Intensity via Friedman's partial
  dependence function: Thiotepa vs.\ Placebo.  \label{RI-Th-Pl}}
\end{figure}

\begin{figure}[t!]
\begin{center}
\includegraphics[scale=0.45]{figures/RI-Th-B6.pdf}
\end{center}
\caption{Relative Intensity via Friedman's partial
  dependence function: Thiotepa vs.\ Vitamin B6.  \label{RI-Th-B6}}
\end{figure}

\begin{figure}[t!]
\begin{center}
\includegraphics[scale=0.45]{figures/RI-B6-Pl.pdf}
\end{center}
\caption{Relative Intensity via Friedman's partial
  dependence function: Vitamin B6 vs.\ Placebo.  \label{RI-B6-Pl}}
\end{figure}


\end{document}

\begin{comment}
\subsection{Case study: diabetes and recurrent hospital admissions}

We study a cohort of newly diagnosed diabetes patients from a single
health care system \citep{SparRein17}.  
% A roughly 20\% random sample
% of patients is available in the {\bf BART} package; see \code{xdm20.train}
% for the covariates (including time, number of previous events and the
% sojourn time) and \code{ydm20.train} for binary indicators.
In this study, we collected the electronic health
records (EHR) for these patients from 2007-2012: prior records may, or
may not, be available.  The EHR are an omnibus of digital health care
information. We focus on 82 covariates: patient demographics, health
insurance, health care charges, diagnoses, procedures, pharmaceutical
therapy, laboratory values and vital signs.  By its very nature, EHR
data is fundamentally time-varying.  EHR covariates are occasionally
missing even when carrying the last value forward since no earlier
value exists.  We imputed 15 continuous variables with Sequential BART
\citep{XuDani16} which is provided by the {\bf sbart} package
\citep{DaniSing17}.

488 patients were followed for 5 years from 2008-2012.  The survival
rate was high, 0.939, suggesting that censoring on death is
noninformative.  Yet this cohort experienced a high rate of hospital
admissions: 525 total: 63\% with no admissions and 37\% with 1 or
more.  In order to properly manage patients suffering diabetes in the
outpatient clinical setting and prevent future admissions, we have
several questions.  Which covariates increase the risk of admission?
What about the number of previous admissions and/or an acutely recent
admission?  What are the functional forms of the covariates, i.e.,
linear, quadratic, logarithm, etc.?  Are the covariate effects
additive or multiplicative?  Are there interactions? Are these effects
constant with respect to time, i.e., do we really want to assume
proportionality?  To let the data speak for itself, we want to avoid
precarious restrictive assumptions.  Hence, we choose to use Bayesian
Additive Regression Trees (BART).

To determine which variables are risk agonists for admission, we used
the variable selection technique known as Decoupling, Shrinkage and
Selection (DSS) \citep{HahnCarv15,McCuCarv15}.  We divided the cohort
at random into training and validation sets of roughly equal size and
performed DSS on each.  Based on the intersection of the findings, we
identified {peripheral vascular disease (PVD)}, {insulin treatment},
and the number of previous admissions, $N_i(t-)$, as risk agonists.


\begin{table}[t!]
\begin{center}
\caption{Diabetes and recurrent hospital admissions.}
\begin{tabular}{lrcrc}\hline
&\multicolumn{2}{c}{Patients}&\multicolumn{2}{c}{Admissions}\\ \hline
Number of Admissions&488&       & 525&\\
  0                 &308& (63.0)&   0& \\
  1                 & 79& (16.2)&  79& (15.0)\\
  2-3               & 50& (10.3)& 115& (21.9)\\
{4-16}              & 51& {(10.5)}& 331& {(63.1)} \\ \hline
Gender              &488&                       &525\\
  M                 &216& (44.3)                &228& (43.4)\\
  F                 &272& (55.7)                &297& (56.6)\\ \hline
Race                &488&                       &525\\
  Black             &174& (35.7)                &265& (50.5)\\
  White             &314& (64.3)                &260& (49.5)\\ \hline
Age                 &488&                       &525\\
  Mean, SD         & 60.9&15.0                & 60.3&15.7\\ \hline
% 21-44             & 66& (13.5)                & 87& (16.6)\\
% 45-54             & 94& (19.3)                & 99& (18.9)\\
% 55-64             &137& (28.1)                &115& (21.9)\\
% 65-74             & 92& (18.9)                &100& (19.0)\\
% 75+               & 99& (20.3)                &124& (23.6)\\ \hline
ZIP area           &488&                       &525\\
  urban         &378& (77.5)                &454& (86.5)\\
  suburb        &110& (22.5)                & 71& (13.5)\\ \hline
Insurance and Age          &488&                       &525\\
  Government 65+    &191& (39.1)                &224& (42.7)\\
  Government $<$65  &138& (28.3)                &208& (39.6)\\
  Commercial $<$65  &143& (29.3)                & 71& (13.5)\\
  Other $<$65       & 16& ( 3.3)                & 22& ( 4.2)\\ \hline
\end{tabular}
\end{center}
\end{table}

\begin{table}[t!]
\begin{center}
\caption{Diabetes and recurrent hospital admissions.}
\begin{tabular}{lrcrccc}\hline
&\multicolumn{2}{c}{Patients}&\multicolumn{2}{c}{Admissions}&Relative Intensity&95\% Credible Interval\\ \hline
% A$_{1c}$(\%)         &416&                       &488\\
%   Mean, SD         &  7.1&\multicolumn{1}{c}{1.6}                 &  6.9&\multicolumn{1}{c}{1.7}\\
%    Missing           & 72&                       & 37\\ \hline
% % Median, IQR      &  6.6&1.8                 &  6.4&1.3\\
% % Q1, Q3            &  6.0& 7.8                 &  5.9& 7.3\\
%   Min, Max          &  4.5& 15.6                &  4.1& 16.4\\
%    $<$ 9            &366& (88.0)                &437  & (89.5)\\
%  $\ge$ 9            & 50& (12.0)                & 51  & (10.5)\\
%   Missing           & 72&                       & 37\\ \hline
{Insulin}             &488&                       &525&            &{2.39} & 1.56, 3.25\\
  Yes               &206& (42.2)                &391& (74.5)\\
  No                &282& (57.8)                &134& (25.5)\\ \hline
{PVD}                 &488&                       &525&            &{2.90} &2.00, 3.89\\
  Yes               &272& (55.7)                &488& (93.0)\\
  No                &216& (44.3)                & 37& ( 7.0)\\ \hline
&                           & & &                              &\multicolumn{2}{c}{Friedman's partial dependence function} \\
%&                           & & &                              &\multicolumn{2}{c}{of BART relative intensities}\\
\end{tabular}
\end{center}
\end{table}

\begin{table}[t!]
\caption{Hospital admission risk profiles.}
\begin{center}
\begin{tabular}{l|cc|cccccc} \hline
& & & \multicolumn{6}{|l}{$N_i(t)$ with time in months} \\ 
Risk & Insulin & PVD & 0 &12 &24 &36 &48 & \multicolumn{1}{l}{60} \\ \hline
Low     & 0 & 0  & 0 & 0 & 0 & 0 & 0 & 0  \\
{Medium}  & 1 & 0  & 0 & 0 & 1 & 1 & 1 & 1  \\
{High}    & 1 & 1  & 0 & 1 & 2 & 3 & 4 & 4  \\ \hline
\end{tabular}
\end{center}
\end{table}

% \begin{center}
% \begin{tabular}{l|cc|cccccc}
% & & & \multicolumn{6}{|l}{$N_i(t)$ with time in months} \\
% Risk & Insulin & PVD & 0 &12 &24 &36 &48 & \multicolumn{1}{l}{60} \\ \hline
% Low     & 0 (60.9\%)& 0 (64.8\%) & 0 & 0 & 0 & 0 & 0 & 0 (63.0\%) \\
% Medium  & 1 (39.1\%)& 0 (64.8\%) & 0 & 0 & 1 & 1 & 1 & 1 (16.2\%) \\
% High    & 1 (39.1\%)& 1 (35.2\%) & 0 & 1 & 2 & 3 & 4 & 4 (10.5\%) \\
% \end{tabular}
% \end{center}

\begin{figure}[t!]
\begin{center}
\includegraphics[scale=0.45]{figures/count-pred-I-color.pdf}
\end{center}
\caption{Risk profiles: Cumulative Intensity via Friedman's partial dependence function.}
\end{figure}

\begin{figure}[t!]
\begin{center}
\includegraphics[scale=0.45]{figures/count-pred-RI-color.pdf}
\end{center}
\caption{Risk profiles: Relative Intensity and 95\% Credible Intervals
  via Friedman's partial dependence function.}
\end{figure}

\begin{figure}[t!]
\begin{center}
\includegraphics[scale=0.45]{figures/count-pred-PI-color.pdf}
\end{center}
\caption{Risk profiles: Relative Intensity and 95\% Prediction Intervals
  via Friedman's partial dependence function.}
\end{figure}

This study has important health policy implications.  The health care
of patients suffering diabetes should be carefully orchestrated to
ensure the delivery of quality clinical care which maximizes healthy
outcomes while preventing adverse events and costly unnecessary
hospital admissions.  Specifically, some diabetes patients are at a
high risk for hospital admission: those diagnosed with PVD, prescribed
insulin therapy, with a recent hospital admission and/or several
previous hospital admissions.  The {\bf BART} package contains a
roughly 20\% random sample of the data from this study
and a worked example in \code{demo/dm.recur.bart.R}.  50 patients
from training: \code{ydm20.train} for the binary indicators and
\code{xdm20.train} for the covariates (including time, number of
previous events and sojourn time).  50 patients from validation:
\code{xdm20.test}.  The complete data set is available online at
\url{http://www.mcw.edu/FileLibrary/Groups/Biostatistics/TechReports/TechReports5175/tr064.tar}.

Acknowledgment and Disclaimer.
\begin{quote}
  This research was supported, in part, by the Advancing Healthier
  Wisconsin Research and Education Program under award 9520277 and the
  National Center for Advancing Translational Sciences, National
  Institutes of Health, through Grant Number UL1TR001436. Its contents
  are solely the responsibility of the authors and do not necessarily
  represent the official views of the NIH.  The data have been
  supplied by the Clinical and Translational Science Institute of
  Southeast Wisconsin's Clinical Research Data Warehouse at the
  Medical College of Wisconsin.  The interpretation and reporting of
  these data are the responsibility of the author(s) and in no way
  should be seen as an official policy of, or interpretation by, the
  Medical College of Wisconsin.
\end{quote}
\end{comment}

%Local Variables:
%TeX-command-default: "weavePDF"
%End:
