\documentclass[nojss]{jss}
%\VignetteEngine{knitr::knitr}
%\VignetteIndexEntry{The BART R package}
%\usepackage{thumbpdf,lmodern}
\usepackage[]{graphicx}
\usepackage[]{color}
\usepackage{tikz}
\usepackage{upquote}
\usepackage{Sweave}
%% maxwidth is the original width if it is less than linewidth
%% otherwise use linewidth (to make sure the graphics do not exceed the margin)

% \makeatletter
% \def\maxwidth{ %
%   \ifdim\Gin@nat@width>\linewidth
%     \linewidth
%   \else
%     \Gin@nat@width
%   \fi
% }
% \makeatother

% \definecolor{fgcolor}{rgb}{0.345, 0.345, 0.345}
% \newcommand{\hlnum}[1]{\textcolor[rgb]{0.686,0.059,0.569}{#1}}%
% \newcommand{\hlstr}[1]{\textcolor[rgb]{0.192,0.494,0.8}{#1}}%
% \newcommand{\hlcom}[1]{\textcolor[rgb]{0.678,0.584,0.686}{\textit{#1}}}%
% \newcommand{\hlopt}[1]{\textcolor[rgb]{0,0,0}{#1}}%
% \newcommand{\hlstd}[1]{\textcolor[rgb]{0.345,0.345,0.345}{#1}}%
% \newcommand{\hlkwa}[1]{\textcolor[rgb]{0.161,0.373,0.58}{\textbf{#1}}}%
% \newcommand{\hlkwb}[1]{\textcolor[rgb]{0.69,0.353,0.396}{#1}}%
% \newcommand{\hlkwc}[1]{\textcolor[rgb]{0.333,0.667,0.333}{#1}}%
% \newcommand{\hlkwd}[1]{\textcolor[rgb]{0.737,0.353,0.396}{\textbf{#1}}}%
% \let\hlipl\hlkwb

% \usepackage{framed}
% \makeatletter
% \newenvironment{kframe}{%
%  \def\at@end@of@kframe{}%
%  \ifinner\ifhmode%
%   \def\at@end@of@kframe{\end{minipage}}%
%   \begin{minipage}{\columnwidth}%
%  \fi\fi%
%  \def\FrameCommand##1{\hskip\@totalleftmargin \hskip-\fboxsep
%  \colorbox{shadecolor}{##1}\hskip-\fboxsep
%      % There is no \\@totalrightmargin, so:
%      \hskip-\linewidth \hskip-\@totalleftmargin \hskip\columnwidth}%
%  \MakeFramed {\advance\hsize-\width
%    \@totalleftmargin\z@ \linewidth\hsize
%    \@setminipage}}%
%  {\par\unskip\endMakeFramed%
%  \at@end@of@kframe}
% \makeatother

% \definecolor{shadecolor}{rgb}{.97, .97, .97}
% \definecolor{messagecolor}{rgb}{0, 0, 0}
% \definecolor{warningcolor}{rgb}{1, 0, 1}
% \definecolor{errorcolor}{rgb}{1, 0, 0}
%\newenvironment{knitrout}{}{} % an empty environment to be redefined in TeX

\usepackage{alltt}
%\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 Charles Spanbauer\\Medical College of Wisconsin
\And Robert McCulloch\\Arizona State University}
\Plainauthor{Rodney Sparapani, Robert McCulloch, Charles Spanbauer}

%% - \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{The \pkg{BART} \proglang{R} package} 

\Plaintitle{The BART R package}

\Shorttitle{The \pkg{BART} package}

%% - \Abstract{} almost as usual
\Abstract{ In this article, we introduce the \pkg{BART} \proglang{R}
  package which is an acronym for Bayesian Additive Regression Trees.
  BART is a Bayesian nonparametric, machine learning, ensemble
  predictive modeling method for continuous, binary, categorical and
  time-to-event outcomes.  Furthermore, BART is a tree-based,
  black-box method which fits the outcome to an arbitrary random
  function, $f$, of the covariates.  The BART technique is relatively
  computationally efficient as compared to its competitors, but large
  sample sizes can be demanding.  Therefore, the \pkg{BART} package
  includes efficient state-of-the-art implementations for continuous,
  binary, categorical and time-to-event outcomes that can take
  advantage of modern off-the-shelf hardware and software
  multi-threading technology.  The \pkg{BART} package is written in
  \proglang{C++} for both programmer and execution efficiency.  The
  \pkg{BART} package takes advantage of multi-threading via forking as
  provided by the \pkg{parallel} package and OpenMP when available and
  supported by the platform.  The ensemble of binary trees produced by
  a BART fit can be stored and re-used later via the \proglang{R}
  \code{predict} function.  In addition to being an \proglang{R}
  package, the installed BART routines can be called directly from
  \proglang{C++}.  The \pkg{BART} package provides the tools for your
  BART toolbox.  }

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

\Keywords{binary trees, black-box, categorical, competing risks,
  continuous, ensemble predictive model, forking, multinomial,
  multi-threading, OpenMP, recurrent events, survival analysis}

\Plainkeywords{binary trees, black-box, categorical, competing risks,
  continuous, ensemble predictive model, forking, multinomial,
  multi-threading, OpenMP, recurrent events, survival analysis}


%% - \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}
}
%\IfFileExists{upquote.sty}{\usepackage{upquote}}{}
\begin{document}

\maketitle

N.B.\ This vignette is largely based on our previous
work \citep{SparSpan19}.  However, erratum and new
developments will by necessity appear here only.

\clearpage

\section{Introduction}

Bayesian Additive Regression Trees (BART) arose out of earlier
research on Bayesian model fitting of an outcome to a single tree
\citep{ChipGeor98}.  In this era from 1996 to 2001, the excellent
predictive performance of ensemble models became apparent
\citep{Brei96,KrogSoli97,FreuScha97,Brei01,Frie01,BaldBrun01}.
Instead of making a single prediction from a complex model, ensemble
models make a single prediction which is the summary of the
predictions from many simple models.  Generally, ensemble models have
desirable properties, e.g., they do not suffer from over-fitting
\citep{KuhnJohn13}.  Like bagging \citep{Brei96}, boosting
\citep{FreuScha97,Frie01} and random forests \citep{Brei01}, BART
relies on an ensemble of trees to predict the outcome; and, although,
there are similarities, there are also differences between these
approaches.

BART is a Bayesian nonparametric, sum of trees method for continuous,
dichotomous, categorical and time-to-event outcomes.  Furthermore,
BART is a black-box, machine learning method which fits the outcome
via an arbitrary random function, $f$, of the covariates.  So-called
black-box models generate functions of the covariates which are so
complex that interpreting the internal details of the fitted model is
generally abandoned in favor of assessment via evaluations of the
fitted function, $f$, at chosen values of the covariates.  As shown by
\citet{ChipGeor10}, BART's out-of-sample predictive performance is
generally equivalent to, or exceeds that, of alternatives like lasso
with L1 regularization \citep{EfroHast04} or black-box models such as
gradient boosting \citep{FreuScha97,Frie01}, neural nets with one
hidden layer \citep{Ripl07,VenaRipl13} and random forests
\citep{Brei01}.  Over-fitting is the tendency to overly fit a model to
an in-sample training data set at the expense of poor predictive
performance for unseen out-of-sample data.  Typically, BART does not
over-fit to the training data due to the regularization tree-branching
penalty of the BART prior, i.e., generally, each tree has few branches
and plays a small part in the overall fit.  So, the resulting fit from
the ensemble of trees as a whole is generally a good fit that does not
over-fit.  Essentially, BART is a Bayesian nonlinear model with all
the advantages of the Bayesian paradigm such as posterior inference
including point and interval estimation.  Conveniently, BART naturally
scales to large numbers of covariates and facilitates variable
selection; it does not require the covariates to be rescaled; neither
does it require the covariate functional relationship, nor the
interactions considered, to be pre-specified.

In this article, we give an overview of data analysis with BART and
the \pkg{BART} \proglang{R} package.  In Section~\ref{cont}, we
describe the \proglang{R} functions provided by the \pkg{BART} package
for analyzing continuous outcomes with BART.  In Section~\ref{boston},
we demonstrate the typical usage of BART via the classic example of
Boston housing values.  In Section~\ref{bincat}, we describe how BART
can be used to analyze binary and categorical outcomes.  In
Section~\ref{surv}, we describe how BART can be used to analyze
time-to-event outcomes with censoring including competing risks and
recurrent events.  In Appendix Section~\ref{BART}, we describe how to
get and install the \pkg{BART} package.  In Appendix
Section~\ref{trees}, we describe the basis of BART on binary trees
along with the details of the BART prior.  In Appendix
Section~\ref{post}, we briefly describe the posterior computations
required to use BART.  In Appendix Section~\ref{efficient}, we
describe how to perform the BART computations efficiently by resorting
to parallel processing with multi-threading 
% added to address Windows multi-threading concerns...
(N.B.\ by default, the Microsoft Windows operating system does not
provide the multi-threading interfaces employed by the \proglang{R}
environment; in lieu of the provision, the \pkg{BART} package is
single-threaded on Windows, yet, otherwise, completely functional; see
Appendix~\ref{efficient} for more details).

\section{Continuous outcomes with BART}\label{cont}

In this section, we document the analysis of continuous outcomes with
the \pkg{BART} \proglang{R} package.  We provide two functions for
continuous outcomes: 1)~\code{wbart} named for weighted BART; and
2)~\code{gbart} named for generic, or generalized, BART.  Both functions
have roughly the same functionality.  \code{wbart} has a verbose
interface while \code{gbart} is streamlined.  Also, \code{wbart}
is for continuous outcomes only whereas \code{gbart} also supports 
binary outcomes.

Typically, when calling the \code{wbart} and \code{gbart} functions,
many of the arguments can be omitted since the default values are
adequate for most purposes.  However, there are certain common
arguments which are either always needed or frequently provided.  The
\code{wbart} (\code{mc.wbart}) and \code{gbart} (\code{mc.gbart})
functions are for serial (parallel) computation; for more details on
parallel computation see the Appendix Section~\ref{efficient}.  The
outcome \code{y.train} is a vector of numeric values.  The covariates
for training (validation, if any) are \code{x.train} (\code{x.test})
which can be matrices or data frames containing factors; in the
display below, we assume matrices for simplicity.  N.B.\ throughout we
denote integer constants by upper case letters, e.g., in the following
display: $M$~for the number of posterior samples, $B$~for the number
of threads (generally, $B=1$ for Windows), $N$~for the number of
observations in the training set, and $Q$~for the number of
observations in the test set.
\begin{Sinput}
R> set.seed(99)
R> post <- wbart(x.train, y.train, x.test, ndpost=M)
R> post <- mc.wbart(x.train, y.train, x.test, ndpost=M, mc.cores=B, seed=99)
R> post <- gbart(x.train, y.train, x.test, ndpost=M)
R> post <- mc.gbart(x.train, y.train, x.test, ndpost=M, mc.cores=B, seed=99)
\end{Sinput}
\begin{align*}
\mbox{Input matrices, \code{x.train} and, optionally, \code{x.test}:\ } & 
\wrap{\begin{array}{c}
\bm{x}_{1} \\
\bm{x}_{2} \\
\vdots \\
\bm{x}_{N} \\
\end{array}} \mbox{made up of $\bm{x}_{i}$ as row vectors}\\
\mbox{\code{post}, of type \code{wbart}, which is essentially a list} & \\
\mbox{\code{post\$yhat.train} and \code{post\$yhat.test}:\ } &
\wrap{\begin{array}{ccc}
\hat{y}_{11}& \dots & \hat{y}_{N1} \\
\vdots & \ddots & \vdots \\
\hat{y}_{1M}& \dots & \hat{y}_{NM} \\
\end{array}} \begin{array}{l} \hat{y}_{im}=\mu_0+f_m(\bm{x}_i) \\
\mbox{\ $m$th posterior draw} \end{array} 
\end{align*}
The columns of \code{post\$yhat.train} and \code{post\$yhat.test}
represent different covariate settings and the rows, the \code{M}
draws from the posterior.

Often it is impractical to provide \code{x.test} in the call to
\code{wbart}/\code{gbart} due to the large number of predictions
considered, or all of the settings to be evaluated are not known at
that time.  To allow for this common problem, the \pkg{BART} package
returns the trees encoded in an ASCII string, \code{treedraws\$trees},
and provides a \code{predict} function to generate any predictions
needed (more details on trees and, the string representation of trees,
can be found in Appendix Section~\ref{trees}).  Note that if you need
to perform the prediction in some later \proglang{R} instance, then
you can save the \code{wbart} object returned and reload it when
needed, e.g., save with \code{saveRDS(post, "post.rds")} and reload,
\code{post <- readRDS("post.rds")}\ .  The \code{x.test} input can be
a matrix or a data frame; for
simplicity, we assume a matrix below.\\
For serial computation
\begin{Sinput}
R> pred <- predict(post, x.test)
\end{Sinput}
For parallel computation
\begin{Sinput}
R> pred <- predict(post, x.test, mc.cores=B)
\end{Sinput}
\begin{align*}
\mbox{Input: \code{x.test}:\ }  &
\wrap{\begin{array}{c}
\bm{x}_{1} \\
\bm{x}_{2} \\
\vdots \\
\bm{x}_{Q} \\
\end{array}} \mbox{made up of $\bm{x}_{h}$ as row vectors} \\
\mbox{\code{pred} is a matrix:} & 
%\mbox{\code{pred}:\ } &
\wrap{\begin{array}{ccc}
\hat{y}_{11}& \dots & \hat{y}_{Q1} \\
\vdots & \ddots & \vdots \\
\hat{y}_{1M}& \dots & \hat{y}_{QM} \\
\end{array}} \where \hat{y}_{hm}=\mu_0+f_m(\bm{x}_h)  
\end{align*}

\subsection{Posterior samples returned}

The number of MCMC samples discarded for burn-in is specified by the
\code{nskip} argument and the default is 100.  The number of MCMC
samples returned is specified by the \code{ndpost} argument and the
default is 1000.  Returning every $l^{th}$ value, or thinning, can be
specified by the \code{keepevery} argument which defaults to 1, i.e.,
no thinning.  Some, but not all, returned values can be thinned.  The
following arguments are available with \code{wbart} and default to
\code{ndpost}, but can be over-ridden as needed (with \code{gbart},
\code{ndpost} draws are always returned and can't be over-ridden).

\begin{itemize}
\item \code{nkeeptrain} : number of $f$ draws to return corresponding to \code{x.train}
\item \code{nkeeptest} :  number of $f$ draws to return corresponding to \code{x.test}
\item \code{nkeeptestmeam} : number of $f$ draws to use in computing \code{yhat.test.mean}
\item \code{nkeeptreedraws} : number of tree ensemble draws to return for use
with \code{predict}
\end{itemize}

Members of the object returned (which is essentially a list) include
\code{varprob} and \code{varcount} which correspond to the variable
selection probability and the observed counts in the ensemble of
trees.  When \code{sparse=TRUE}, \code{varprob} is the random variable
selection probability, $s_j$; otherwise, it is the fixed constant
$s_j=P^{-1}$.  Besides the posterior samples, also the mean over the
posterior is provided as \code{varprob.mean} and \code{varcount.mean}.
% Objects returned when the sparsity prior is being used include a
% matrix of the splitting probabilities at each iteration in
% \code{varprobs} and a matrix of splitting counts at each iteration in
% \code{varcounts}. In particular, the latter can be useful to assess
% how the ensemble is changing throughout the posterior sampling.

\section{The Boston housing values example}\label{boston} 

Now, let's examine the classic Boston housing values example
\citep{HarrRubi78}.  This data is from the 1970 US Census where each
observation represents a Census tract in the Boston Standard
Metropolitan Statistical Area.  For each tract, there was a localized
air pollution estimate, the concentration of nitrogen oxides, based on
a meteorological model that was calibrated to monitoring data.
Restricted to tracts with owner-occupied homes, there are 506
observations.  We'll predict the median value of owner-occupied homes
(in thousands of dollars truncated at 50), \code{y=mdev}, from two
covariates: \code{rm} and \code{lstat}.  \code{rm} is the number of
rooms defined as the average number of rooms for owner-occupied homes.
\code{lstat} is the percent of population that is lower status defined
as the average of the proportion of adults without any high school
education and the proportion of male workers classified as laborers.
Below, we present several observations of the data and scatter plots
in Figure~\ref{boston1}.
\begin{Sinput}
R> library("MASS")
R> x = Boston[, c(6, 13)] 
R> y = Boston$medv     
R> head(cbind(x, y))
\end{Sinput} 
%$
\begin{Soutput}
     rm lstat    y
1 6.575  4.98 24.0
2 6.421  9.14 21.6
3 7.185  4.03 34.7
4 6.998  2.94 33.4
5 7.147  5.33 36.2
6 6.430  5.21 28.7
\end{Soutput}
\begin{Sinput}
R> par(mfrow=c(2, 2))
R> plot(x[, 1], y, xlab="x1=rm", ylab="y=mdev")
R> plot(x[, 2], y, xlab="x2=lstat", ylab="y=mdev")
R> plot(x[, 1], x[, 2], xlab="x1=rm", ylab="x2=lstat")
R> par(mfrow=c(1, 1))
\end{Sinput}
\begin{figure}
\begin{center}
\includegraphics{figures/boston1.pdf}
\end{center}
\caption{\label{boston1}The Boston housing data was compiled from the
  1970 US Census where each observation represents a Census tract in
  Boston with owner-occupied homes. For each tract, we have the median
  value of owner-occupied homes (in thousands of dollars truncated at
  50), \code{y=mdev}, the average number of rooms, \code{x1=rm}, and
  the percent of the population that is lower status, \code{x2=lstat}.
  Here, we show scatter plots of the data.  }
\end{figure}

\subsection{wbart for continuous outcomes}

In this example, we fit the following BART model for continuous outcomes:
\begin{align*}
y_i &= \mu_0 + f(x_i) + \epsilon_i \where \epsilon_i \sim N(0,\sigma^2) \\
(f, \sd^2) & \prior \mathrm{BART}
\end{align*}
with $i$ indexing subjects; $i=1, \dots, N$.  We use Markov chain
Monte Carlo (MCMC) to get draws from the posterior distribution of the
parameter $(f,\sigma^2)$. 

\begin{Sinput}
R> library("BART") 
R> set.seed(99)  
R> nd = 200       
R> burn = 50     
R> post = wbart(x, y, nskip=burn, ndpost=nd)
\end{Sinput}
\begin{Soutput}
*****Into main of wbart
*****Data:
data:n,p,np: 506, 2, 0
y1,yn: 1.467194, -10.632806
x1,x[n*p]: 6.575000, 7.880000
*****Number of Trees: 200
*****Number of Cut Points: 100 ... 100
*****burn and ndpost: 50, 200
*****Prior:beta,alpha,tau,nu,lambda: 2.000000,0.950000,0.795495,3.000000,5.979017
*****sigma: 5.540257
*****w (weights): 1.000000 ... 1.000000
*****Dirichlet:sparse,a,b,rho,augment: 0,0.5,1,2,0
*****nkeeptrain,nkeeptest,nkeeptestme,nkeeptreedraws: 200,200,200,200
*****printevery: 100
*****skiptr,skipte,skipteme,skiptreedraws: 1,1,1,1

MCMC
done 0 (out of 250)
done 100 (out of 250)
done 200 (out of 250)
time: 1s
check counts
trcnt,tecnt,temecnt,treedrawscnt: 200,0,0,200
\end{Soutput}

\subsection{Results returned from wbart}
 
We returned the results of running \code{wbart} in the object
\code{post} of type \code{wbart} which is essentially a list.

\begin{Sinput}
R> names(post)
\end{Sinput}
\begin{Soutput}
 [1] "sigma"           "yhat.train.mean" "yhat.train"      "yhat.test.mean" 
 [5] "yhat.test"       "varcount"        "varprob"         "treedraws"      
 [9] "mu"              "varcount.mean"   "varprob.mean"    "rm.const"       
\end{Soutput}
\begin{Sinput}
R> length(post$sigma)
\end{Sinput}
%$
\begin{Soutput}
[1] 250
\end{Soutput}
\begin{Sinput}
R> length(post$yhat.train.mean)
\end{Sinput}
%$
\begin{Soutput}
[1] 506
\end{Soutput}
\begin{Sinput}
R> dim(post$yhat.train)
\end{Sinput}
%$
\begin{Soutput}
[1] 200 506
\end{Soutput}
%$
%\textcolor{red}
{Remember}, the training data has $n=$
506 observations, we had \code{burn}=50 
burn-in discarded draws and \code{nd}=$M$=200 draws kept.
Let's look at a couple of the key list components.\\
\code{\$sigma}: both the 50 burn-in and 250 draws are kept for $\sigma$;
burn-in are kept only for this parameter.\\
\code{\$yhat.train}: the $m$th row and $i$th column is $f_m(x_i)$ (the $m^{th}$ kept MCMC draw
evaluated at the $i^{th}$ training observation).\\
\code{\$yhat.train.mean}: the posterior estimate of $f(x_i)$, i.e.,
$M^{-1} \sum_m f_m(x_i)$.

\subsection{Assessing convergence with wbart} 
%{ with $\sigma$ Draws}

As with any high-dimensional MCMC, assessing convergence may be
non-trivial.  Posterior convergence diagnostics are recommended for
BART especially with large data sets and/or a large number of
covariates.  Besides diagnostics, routine counter-measures such as
longer chains, thinning and multiple chains may be warranted.  For
continuous outcomes, the simplest thing to look at are the draws of
$\sigma$.  See Section~\ref{geweke} for a primer on other convergence
diagnostic options for binary and categorical outcomes that are also
applicable for continuous outcomes.

Assessing convergence in this example, the parameter $\sigma$ is the
only identified parameter in the model and, of course, it is
indicative of the size of the errors.

\begin{Sinput}
R> plot(post$sigma, type="l")
R> abline(v=burn, lwd=2, col="red")
\end{Sinput}
%$
\begin{figure}
\begin{center}
\includegraphics{figures/boston2.pdf}
\end{center}
\caption{\label{boston2}The Boston housing data was compiled from the
  1970 US Census where each observation represents a Census tract in
  Boston with owner-occupied homes. For each tract, we have the median
  value of owner-occupied homes (in thousands of dollars truncated at
  50), \code{y=mdev}, the average number of rooms, \code{x1=rm}, and
  the percent of the population that is lower status, \code{x2=lstat}.
  With BART, we predict \code{y=mdev} from \code{rm} and \code{lstat}.
  Here, we show a trace plot of the error variance, $\sd$, that
  demonstrates convergence for BART rather quickly, i.e., by 50
  iterations or earlier.  }
\end{figure}

In Figure~\ref{boston2}, you can see that BART burned in very quickly.  
Just one initial draw looking a bit bigger than the rest.
Apparently, subsequent variation is legitimate posterior variation.  
In a more difficult problem you may see the $\sigma$ draws initially 
declining as the MCMC searches for a good fit.

\subsection{wbart and linear regression compared}

Let's look at the in-sample BART fit (\code{yhat.train.mean})  
and compare it to \code{y=medv} fits from a multiple linear regression.  

\begin{Sinput}
R> lmf = lm(y~., data.frame(x, y))
R> fitmat = cbind(y, post$yhat.train.mean, lmf$fitted.values)
R> colnames(fitmat) = c("y", "BART", "Linear")
R> cor(fitmat)
\end{Sinput}
\begin{Soutput}
               y      BART    Linear
y      1.0000000 0.9051200 0.7991005
BART   0.9051200 1.0000000 0.8978003
Linear 0.7991005 0.8978003 1.0000000
\end{Soutput}
\begin{Sinput}
R> pairs(fitmat)
\end{Sinput}
\begin{figure}
\begin{center}
\includegraphics{figures/boston3.pdf}
\end{center}
\caption{\label{boston3}The Boston housing data was compiled from the
  1970 US Census where each observation represents a Census tract in
  Boston with owner-occupied homes. For each tract, we have the median
  value of owner-occupied homes (in thousands of dollars truncated at
  50), \code{y=mdev}, the average number of rooms, \code{x1=rm}, and
  the percent of the population that is lower status, \code{x2=lstat}.
  With BART, we predict \code{y=mdev} from \code{rm} and \code{lstat}.
  Here, we show scatter plots comparing \code{y=mdev}, the BART fit
  (``BART'') and multiple linear regression (``Linear'').  }
\end{figure}
In Figure~\ref{boston3}, we present scatter plots between
\code{mdev}, the BART fit and the multiple linear regression.
The BART fit is noticeably different from the linear fit.

\subsection{Prediction and uncertainty with wbart}

In Figure~\ref{boston4}, we order the observations by the fitted house
value (\code{yhat.train.mean}) and then use boxplots to display the
draws of $f(x)$ in each column of \code{yhat.train}.
\begin{Sinput}
R> i = order(post$yhat.train.mean) 
R> boxplot(post$yhat.train[, i])   
\end{Sinput}
\begin{figure}
\begin{center}
\includegraphics{figures/boston4.pdf}
\end{center}
\caption{\label{boston4}The Boston housing data was compiled from the
  1970 US Census where each observation represents a Census tract in
  Boston with owner-occupied homes. For each tract, we have the median
  value of owner-occupied homes (in thousands of dollars truncated at
  50), \code{mdev}, the average number of rooms, \code{rm}, and the
  percent of the population that is lower status, \code{lstat}.  With
  BART, we predict \code{y=mdev} from \code{rm} and \code{lstat}.
  Here, we show boxplots of the posterior samples of predictions (on
  the y-axis) ordered by the average predicted home value per tract
  (on the x-axis).  }
\end{figure}

Substantial predictive uncertainty, but you can still be fairly certain that
some houses should cost more than other.

\subsection{Using the predict function with wbart}

We can get out of sample predictions in two ways. 
First, we can can just ask for them when we call \code{wbart} by supplying 
a matrix or data frame of test $x$ values. 
Second, we can call a \code{predict} method.
Now, let's split our data into train and test subsets.  

\begin{Sinput}
R> n = length(y)   
R> set.seed(14) 
R> i = sample(1:n, floor(0.75*n)) 
R> x.train = x[i, ]; y.train=y[i]  
R> x.test = x[-i, ]; y.test=y[-i] 
R> cat("training sample size = ", length(y.train), "\n")
R> cat("testing sample size = ", length(y.test), "\n")
\end{Sinput}
\begin{Soutput}
training sample size = 379  
testing sample size = 127 
\end{Soutput}

And now we can run \code{wbart} using the training data to 
learn and predict at \code{x.test}.
First, we'll just pass \code{x.test} to the \code{wbart} call.

\begin{Sinput}
R> set.seed(99)
R> post1 = wbart(x.train, y.train, x.test) 
R> dim(post1$yhat.test)
\end{Sinput}
%$
\begin{Soutput}
[1] 1000  127
\end{Soutput}
\begin{Sinput}
R> length(post1$yhat.test.mean)
\end{Sinput}
%$
\begin{Soutput}
[1] 127
\end{Soutput}

The testing data is handled similarly to the training data.\\
\code{\$yhat.test}: the $m$th row and $h$th column is $f_m(x_h)$ (the $m^{th}$ kept MCMC draw
evaluated at the $h^{th}$ testing observation).\\
\code{\$yhat.test.mean}: the posterior estimate of $f(x_h)$, i.e.,
$Q^{-1} \sum_m f_m(x_h)$.
% Now, \code{yhat.test}: the $i,j$ value is the $i^{th}$ kept MCMC draw of 
% $f(x_j)$ where $x_j$ is the $j^{th}$ row of \code{xtest}.\\
% \code{yhat.test.mean}: the $j^{th}$ value is the posterior mean of 
% $f(x_j)$, i.e., $f$ evaluated at the $j^{th}$ row of \code{xtest}.

Alternatively, we could run \code{wbart} saving all the MCMC results
and then call \code{predict}.
%\code{predict.wbart}.
\begin{Sinput}
R> set.seed(99)
R> post2 = wbart(x.train, y.train)
R> yhat = predict(post2, x.test)
\end{Sinput}
\begin{Soutput}
*****In main of C++ for bart prediction
tc (threadcount): 1
number of bart draws: 1000
number of trees in bart sum: 200
number of x columns: 2
from x,np,p: 2, 127
***using serial code
\end{Soutput}
\begin{Sinput}
R> dim(yhat)
\end{Sinput}
[1] 1000  127
\begin{Sinput}
R> summary(as.double(yhat - post1$yhat.test))
\end{Sinput}
%$
\begin{Soutput}
      Min.    1st Qu.     Median       Mean    3rd Qu.       Max. 
-9.091e-09 -1.186e-09  2.484e-11  2.288e-12  1.188e-09  6.790e-09 
\end{Soutput}
So \code{yhat} and \code{post1$yhat.test} are practically identical.%$

\subsection{wbart and thinning}

In our simple example of the Boston housing data, \code{wbart} runs
pretty fast.  But with more data and/or longer runs, you may want to
speed things up by saving fewer samples and then using \code{predict}.
Let's just keep a thinned subset of 200 tree ensemble draws.

\begin{Sinput}
R> set.seed(4) 
R> post3 = wbart(x.train, y.train, nskip=1000, ndpost=10000,
+    nkeeptrain=0, nkeeptest=0, nkeeptestmean=0, nkeeptreedraws=200)
R> yhatthin = predict(post3, x.test)
\end{Sinput}
\begin{Soutput}
*****In main of C++ for bart prediction
tc (threadcount): 1
number of bart draws: 200
number of trees in bart sum: 200
number of x columns: 2
from x,np,p: 2, 127
***using serial code
\end{Soutput}
\begin{Sinput}
R> dim(post3$yhat.train)
\end{Sinput}
%$
\begin{Soutput}
[1]   0 379
\end{Soutput}
\begin{Sinput}
R> dim(yhatthin)
\end{Sinput}
\begin{Soutput}
[1] 200 127
\end{Soutput}
Now, there are no kept draws of $f(x)$ for training $x$, and we have
200 tree ensemble draws to use with \code{predict}.  Of course, if we keep
200 out of 10000, then every 50th draw is kept.
%\code{predict.wbart}.  

% The thinning arguments for \code{wbart} are the following.\\
% \code{nkeeptrain} : number of $f(x)$ draws to save for training $x$\\
% \code{nkeeptest} : number of $f(x)$ draws to save for test $x$\\
% \code{nkeeptestmeam} : number of draws to use in computing \code{yhat.test.mean}\\
% \code{nkeeptreedraws} : number of tree ensemble draws to keep for
% \code{predict}

% The thinning arguments for \code{wbart}.
% \begin{itemize}
% \item \code{nkeeptrain} : number of $f(x)$ draws to save for training $x$
% \item \code{nkeeptest} : number of $f(x)$ draws to save for test $x$
% \item \code{nkeeptestmeam} : number of draws to use in computing 
% \code{yhat.test.mean}
% \item \code{nkeeptreedraws} : number of tree ensemble draws to keep
% for \code{predict}
% \end{itemize}

The default values are to keep all the draws (e.g.,
\code{nkeeptrain=ndpost}).
Now, let's have a look at the predictions.
\begin{Sinput}
R> fmat = cbind(y.test, post1$yhat.test.mean, apply(yhatthin, 2, mean)) 
R> colnames(fmat) = c("y", "yhat", "yhatThin")
R> pairs(fmat)
\end{Sinput}
%$
\begin{figure}
\begin{center}
\includegraphics{figures/boston5.pdf}
\end{center}
\caption{\label{boston5}The Boston housing data was compiled from the
  1970 US Census where each observation represents a Census tract in
  Boston with owner-occupied homes. For each tract, we have the median
  value of owner-occupied homes (in thousands of dollars truncated at
  50), \code{mdev}, the average number of rooms, \code{rm}, and the
  percent of the population that is lower status, \code{lstat}.  With
  BART, we predict \code{y=mdev} by \code{rm} and \code{lstat}.  The
  predictions labeled ``yhat'' are from a BART run with \code{seed=99}
  and all default values.  The predictions labeled ``yhatThin'' are
  thinned by 50 (after 1000 burnin discarded, 200 kept out of 10000
  draws) with \code{seed=4}.  It is very interesting how similar they
  are!  }
\end{figure}

In Figure~\ref{boston5}, we present scatter plots between \code{mdev},
``yhat'' and ``yhatThin''.  Recall, the predictions labeled ``yhat'' are
from a BART run with \code{seed=99} and all default values.  The
predictions labeled ``yhatThin'' are thinned by 50 (after 1000 burnin
discarded, 200 kept out of 10000 draws) with \code{seed=4}.  It is
very interesting how similar they are!

\subsection{wbart and Friedman's partial dependence function}
\label{Friedman-partial}

BART does not directly provide a summary of the effect of a single
covariate, or a subset of covariates, on the outcome.  This is also
the case for black-box, or nonparametric regression, models in general
that need to deal with this same issue.  Developed for such complex
models, Friedman's
partial dependence function \citep{Frie01} can be employed with BART
to summarize the marginal effect due to a subset of the covariates.
Friedman's partial dependence function is a concept that is
very flexible.  So flexible that we are unable to provide abstract
functional support in the \pkg{BART} package; rather, we provide
examples of the many practical uses in the \code{demo} directory.

We use $S$ to denote the indices of the covariates in the subset and
the collection itself, i.e., define the row vector for test setting
$h$ as $\bm{x}_{hS}=\wrap{{x}_{hj}} \where j \in S$.  Similarly, we
denote the complement of the subset as $C$ with $S \cup C$ spanning
all covariates.  The complement row vector for training
observation~$i$ is $\bm{x}_{iC}=\wrap{{x}_{ij}} \where j \in C$.  The
marginal dependence function is defined by fixing the subset at a test
setting while aggregating over the training observations of the
complement covariates:
$f(\bm{x}_{hS})={N^{-1}}\sum_{i=1}^N f(\bm{x}_{hS},\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, etc., e.g.,
$\hat{f}(\bm{x}_{hS}) = {M^{-1}} {N^{-1}} \sum_{m=1}^M \sum_{i=1}^N
f_m(\bm{x}_{hS}, \bm{x}_{iC})$ where $m$ indexes posterior samples.  
However, care must be taken in the interpretation of the marginal 
effect as estimated by Friedman's partial dependence function.  
If there are strong relationships among the covariates, it may be
unrealistic to assume that individual covariates can be manipulated
independently. 

For example, suppose that we want to summarize the median home value,
\code{medv} (variable 14 of the \code{Boston} data frame), by the
percent of the population with lower status, \code{lstat} (variable
13), while aggregating over the other twelve covariates in the Boston
housing data.  In Figure~\ref{boston6}, we demonstrate the marginal
estimate and its 95\% credible interval.
\begin{Sinput}
R> x.train = as.matrix(Boston[i, -14])
R> set.seed(12) 
R> post4 = wbart(x.train, y.train)
R> H = floor(0.75*length(y.train))
R> L = 41
R> x = seq(min(x.train[, 13]), max(x.train[, 13]), length.out=L)
R> x.test = cbind(x.train[, -13], x[1])
R> for(j in 2:L)
+    x.test = rbind(x.test, cbind(x.train[, -13], x[j]))
R> pred = predict(post4, x.test)
R> partial = matrix(nrow=1000, ncol=L)
R> for(j in 1:L) {
R>   h = (j - 1) * H + 1:H
R>   partial[, j] = apply(pred[, h], 1, mean)
R> }
R> plot(x, apply(partial, 2, mean), type='l',
+    xlab='percent lower status', ylab='median home value',
+    ylim=c(10, 50))
R> lines(x, apply(partial, 2, quantile, probs=0.025), lty=2)
R> lines(x, apply(partial, 2, quantile, probs=0.975), lty=2)
\end{Sinput}
\begin{figure}
\begin{center}
\includegraphics{figures/boston6.pdf}
\end{center}
\caption{\label{boston6}The Boston housing data was compiled from the
  1970 US Census where each observation represents a Census tract in
  Boston with owner-occupied homes. For each tract, we have the median
  value of owner-occupied homes (in thousands of dollars truncated at
  50), \code{mdev}, and the percent of the population that is lower
  status, \code{lstat}, along with eleven other covariates.  We
  summarize the marginal effect of \code{lstat} on \code{mdev} while
  aggregating over the other covariates with Friedman's partial
  dependence function.  The marginal estimate and its 95\% credible
  interval are shown. }
\end{figure}
Besides the marginal effect, we can 
define the conditional effect of $x_1|x_2$ as
  $\frac{f(x_1+\delta, {x}_2)-f(x_1, {x}_2)}{\delta}$.  However,
  BART is not fitting simple linear functions.  For example, suppose
  the data follows a sufficiently complex function like so:
  $f(x_1, x_2)=b_1 x_1 + b_2 x_1^2 + b_3 x_1 x_2$.  Then the
  conditional effect that BART is likely to fit is approximately
  $b_1 + 2 b_2 x_1 + b_2\delta+ b_3 {x}_2$.  This function is not
so easy to characterize (as the marginal effect) since it involves
$x_1$, $x_2$ and $\delta$.  Nevertheless, these functions can be
estimated by BART if these inputs are provided.  But, these 
functions have the same limitations as Friedman's partial dependence
function and, perhaps, even moreso.  See the conditional effect
example at the end of \code{demo("boston.R", package="BART")}.

\section{Binary and categorical outcomes with BART}\label{bincat}

The \pkg{BART} package supports binary outcomes via probit BART with
Normal latents and logit BART with Logistic latents.  Categorical
outcomes are supported with Multinomial BART which defaults to probit
for computational efficiency, but logit is available as an option.
Convergence diagnostics are provided and variable selection as well.

\subsection{Probit BART for binary outcomes}

Probit BART for binary outcomes is provided by the \pkg{BART} package
as the \code{pbart} and \code{gbart} functions.  In this case, the
outcome, \code{y.train}, is an integer with values of 0 or 1.  The
model is as follows with $i$ indexing subjects: $i=1, \dots, N$.
% where $\Phi$ is the standard Normal cumulative distribution function.
\begin{align*}
y_{i}|p_{i} & \ind \B{p_{i}} 
\where \B{.} \mbox{is the Bernoulli distribution} \\
p_{i} & =  \Phi({\mu}_0+f(\bm{x}_i)) \where f \prior \mathrm{BART} 
\mbox{\ and\ } \Phi(.) \mbox{\ is the standard Normal cdf}
\end{align*}
This setup leads to the following likelihood:
$\wrap{\bm{y}|f} = \prod_{i=1}^N p_{i}^{y_i}(1-p_{i})^{1-y_i}$. 

% \begin{comment}
% The BART function, $f$, is centered around a known constant, $\mu_0$,
% which is analogous to centering the probabilities, $p_i$, around
% $p_0=\Phi(\mu_0)$.  The default value of $\mu_0$ is
% $\Phi^{-1}(\bar{y})$ (which you can over-ride with the
% \code{binaryOffset} argument).
% \end{comment}

To extend BART to binary outcomes, we employ the technique of
\cite{AlbeChib93} that assumes there is an unobserved latent, $z_i$,
where $y_i=\I{z_i>0}$ and $i=1, \dots, n$ indexes subjects. Given
$y_i$, we generate the truncated Normal latents, $z_i$; these
auxiliary latents are efficiently sampled \citep{Robe95} and recast as
the outcome for a continuous BART with unit variance as follows.
\begin{align*}
z_{i}|y_{i},f & ~ \N{{\mu}_0+ f(\bm{x}_i)}{1} \begin{cases}
\I{-\infty, 0} & \If y_{i}=0 \\
\I{0, \infty} & \If y_{i}=1 \\
\end{cases}  
\end{align*}

Centering the latent $z_i$ around the constant $\mu_0$ is
  analogous to quasi-centering the probabilities, $p_i$, at
  $p_0=\Phi(\mu_0)$, i.e., $\E{p_i}$ is approximately equal
  to $p_0$ which is all that is necessary for inference to be
  performed.  The default value of $\mu_0$ is
$\Phi^{-1}(\bar{y})$ (which you can over-ride with the
\code{binaryOffset} argument).

% If $\mu_0=0$, which is the default, then the $p_i$
% are centered around 0.5; to center at a different value, say 0.05, pass
% the argument \code{binaryOffset=-1.645} in the \code{pbart} call.  
% The key
% insight into the probit BART technique is that the Gibbs conditional
% $ {f|z_i,y_i} {\;\stackrel{d}{=} f|z_i} $, i.e., given $z_i$, $y_i$ is
% unnecessary.  

%In the following, we assume that \code{binaryOffset=0} for convenience.
%(which is the default).  
The \code{pbart} (\code{mc.pbart}) and \code{gbart} (\code{mc.gbart})
functions are for serial (parallel) computation.  The outcome
\code{y.train} is a vector containing zeros and ones.  The covariates
for training (validation, if any) are \code{x.train} (\code{x.test})
which can be matrices or data frames containing factors; in the
display below, we assume matrices for simplicity. 
Notation: $M$~for the number of posterior samples, $B$~for the number
of threads (generally, $B=1$ for Windows), $N$~for the number of
observations in the training set, and $Q$~for the number of
observations in the test set.
\begin{Sinput}
R> set.seed(99)
R> post <- pbart(x.train, y.train, x.test, ndpost=M)
R> post <- mc.pbart(x.train, y.train, x.test, ndpost=M, mc.cores=B, seed=99)
R> post <- gbart(x.train, y.train, x.test, type='pbart', ndpost=M)
R> post <- mc.gbart(x.train, y.train, x.test, type='pbart', ndpost=M, ...
\end{Sinput}
N.B.\ for \code{pbart}, the thinning argument, \code{keepevery}
defaults to 1 while for \code{gbart} with\\ \code{type='pbart'},
\code{keepevery} defaults to 10.  

\begin{align*}
\mbox{Input matrices: \code{x.train} and, optionally, \code{x.test}:\ } & 
\wrap{\begin{array}{c}
\bm{x}_{1} \\
\bm{x}_{2} \\
\vdots \\
\bm{x}_{N} \\
\end{array}} \mbox{\ or\ } \bm{x}_{i} \mbox{\ as row vectors} \\
\mbox{\code{post}, of type \code{pbart}, which is essentially a list} & \\
\mbox{\code{post\$yhat.train} and \code{post\$yhat.test}:\ } &
\wrap{\begin{array}{ccc}
\hat{y}_{11}& \dots & \hat{y}_{N1} \\
\vdots & \ddots & \vdots \\
\hat{y}_{1M}& \dots & \hat{y}_{NM} \\
\end{array}} \hat{y}_{im}=\mu_0+f_m(\bm{x}_i) 
\end{align*}

The columns of \code{post\$yhat.train} and \code{post\$yhat.test}
represent different covariate settings and the rows, the \code{M}
draws from the posterior.  \code{post\$yhat.train} and
\code{post\$yhat.test}, when requested, are returned, although,
\code{post\$prob.train} and \code{post\$prob.test} are generally of
more interest (and \code{post\$prob.train.mean} and
\code{post\$prob.test.mean} which are the means of the posterior
sample columns, not shown).

\begin{align*}
\mbox{\code{post\$prob.train} and \code{post\$prob.test}:\ } &
\wrap{\begin{array}{ccc}
\hat{p}_{11}& \dots & \hat{p}_{N1} \\
\vdots & \ddots & \vdots \\
\hat{p}_{1M}& \dots & \hat{p}_{NM} \\
\end{array}} \where  \hat{p}_{im}=\Phi(\hat{y}_{im})  %\Phi(f_m(\bm{x}_i)) \\
\end{align*}

Often it is impractical to provide \code{x.test} in the call to
\code{pbart} due to the number of predictions considered or all the
settings to evaluate are simply not known at that time.  To allow for
this common problem, the \pkg{BART} package returns the trees encoded
in an ASCII string, \code{treedraws\$trees}, and provides a
\code{predict} function to generate any predictions needed.  Note that
if you need to perform the prediction in some later \proglang{R}
instance, then you can save the \code{pbart} object returned and
reload it when needed, e.g., save with \code{saveRDS(post,
  "post.rds")} and reload, \code{post <- readRDS("post.rds")}\ .
The \code{x.test} input can be a matrix or a data frame; for
simplicity, we assume a matrix below.\\
\begin{Sinput}
R> pred <- predict(post, x.test, mc.cores=B)
\end{Sinput}
\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}_h \mbox{\ as row vectors} \\
\mbox{\code{pred}, of type \code{pbart}, which is essentially a list} & \\
\mbox{\code{pred\$yhat.test}:\ } &
\wrap{\begin{array}{ccc}
\hat{y}_{11}& \dots & \hat{y}_{Q1} \\
\vdots & \ddots & \vdots \\
\hat{y}_{1M}& \dots & \hat{y}_{QM} \\
\end{array}} \where \hat{y}_{hm}=\mu_0+f_m(\bm{x}_h)  \\
\mbox{\code{pred\$prob.test}:\ } &
\wrap{\begin{array}{ccc}
\hat{p}_{11}& \dots & \hat{p}_{Q1} \\
\vdots & \ddots & \vdots \\
\hat{p}_{1M}& \dots & \hat{p}_{QM} \\
\end{array}} \where \hat{p}_{hm}=\Phi(\hat{y}_{hm}) \\ %\Phi(f_m(\bm{x}_h))  \\
\mbox{\code{pred\$prob.test.mean}:\ } &
\wrap{\hat{p}_{1}, \dots, \hat{p}_{Q}} \where \hat{p}_{h}=M^{-1} \sum_{m=1}^M  \hat{p}_{hm}   
\end{align*}

\subsection{Probit BART and Friedman's partial dependence function}

For an overview of Friedman's partial dependence function (including
the notation adopted in this article and its meaning), please see 
Section~\ref{Friedman-partial} which discusses continuous outcomes.
% \begin{comment}
% BART does not directly provide a summary of the effect of a single
% covariate, or a subset of covariates, on the outcome.  This is also
% the case for black-box, or nonparametric regression, models in general
% which have had to deal with this issue.  We recommend utilizing Friedman's
% partial dependence function \citep{Frie01} with BART to summarize the
% marginal effect due to a subset of the covariates, $\bm{x}_S$, by
% aggregating over the complement covariates, $\bm{x}_C$, i.e.,
% $\bm{x} =\wrap{\bm{x}_S,\bm{x}_C}$. The marginal dependence function
% is defined by fixing $\bm{x}_S$ while aggregating over the observed
% settings of the complement covariates in the cohort:
% $f(\bm{x}_S)={N^{-1}}\sum_{i=1}^N f(\bm{x}_S,\bm{x}_{iC})$.
% \end{comment}
For probit BART, the $f$ function is not directly of interest; rather,
the probability of an event is more interpretable:
$p(\bm{x}_{hS}) = {N^{-1}} \sum_{i=1}^N
\Phi(\mu_0+f(\bm{x}_{hS},\bm{x}_{iC}))$.  
% \begin{comment}
% 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{p}(\bm{x}_{hS}) = {M^{-1}} {N^{-1}} \sum_{m=1}^M \sum_{i=1}^N
% \Phi(\mu_0+f_m(\bm{x}_{hS},\bm{x}_{iC}))$ where $m$ indexes posterior
% samples.  
% Friedman's partial dependence function is a concept that is
% very flexible.  So flexible that we are unable to provide abstract
% functional support in the \pkg{BART} package; rather, we provide
% examples of the many practical uses in the \code{demo} directory.
% \end{comment}

\subsubsection{Probit BART example: chronic pain and obesity}
\label{nhanes}

We want to explore the hypothesis that obesity is a risk factor for
chronic lower-back pain (which includes buttock pain in this
definition).  A corollary to this hypothesis is that obesity is not
considered to be a risk factor for chronic neck pain.  A good source
of data for this question is available in the National Health and
Nutrition Examination Survey (NHANES) 2009-2010 Arthritis
Questionnaire.  5106 subjects were surveyed.  We will use probit BART
to analyze the dichotomous outcomes of chronic lower-back pain and
chronic neck pain.  We restrict our attention to the following
covariates: age, gender and anthropometric measurements including
weight (kg), height (cm), body mass index (kg/m$^2$) and waist
circumference (cm).  Also, note that survey sampling weights are
available to extrapolate the rates from the survey to the US
population as a whole.  We will concentrate on body mass index (BMI)
and gender, $\bm{x}_{hS}$, while utilizing Friedman's partial
dependence function as defined above and incorporating the survey
weights, i.e.,
$p_{hS}(\bm{x}_{hS}) = {\sum_{i=1}^N w_i
  \Phi(\mu_0+f(\bm{x}_{hS},\bm{x}_{iC}))}/ {\sum_{i'=1}^N w_{i'}}$.

The \pkg{BART} package provides two examples for the relationship
between chronic pain and BMI: \code{demo("nhanes.pbart1",
  package="BART")} for the probabilities and
\code{demo("nhanes.pbart2", package="BART")} for differences in these
probabilities.  In Figure~\ref{chronic-pain1}, the left panel for
lower-back pain and the right panel for neck pain, the
unweighted relationship between chronic pain, BMI and gender are
displayed: males (females) are represented by blue (red) solid lines
with corresponding 95\% credible intervals in dashed lines.  Although
there is a generous amount of uncertainty, it does not appear that the
probability of chronic lower-back pain increases with BMI for either
gender.  Conversely, chronic neck pain does appear to be rising, yet
again, the intervals are wide.  In both cases, these findings are not
anticipated given the original hypotheses.  Based on survey weights
(not shown), the results are basically the same.  In
Figure~\ref{chronic-pain2}, the unweighted relationship for females
between BMI and the difference in probability of chronic pain from a
baseline BMI of 25 (which is the upper limit of normal) with
corresponding 95\% credible intervals in dashed lines: the left panel
for lower-back pain (blue solid lines) and the right panel for
neck pain (red solid lines). Again, we have roughly the same
impression, i.e., there is no increase of lower-back chronic
pain with BMI and it is possibly dropping while neck pain might be
increasing, but the intervals are wide for both.  The results are
basically the same for males (not shown).
\begin{figure}%[t!]
\begin{center}
\includegraphics[scale=0.425]{figures/chronic-pain1.pdf}
\end{center}
\caption{\label{chronic-pain1}NHANES, BMI and the probability of
  chronic pain: the left panel for lower-back pain and the
  right panel for neck pain.  The unweighted Friedman's partial
  dependence relationship between chronic pain, BMI and gender are
  displayed as ascertained from NHANES data: males (females) are
  represented by blue (red) lines with the corresponding 95\% credible
  intervals (dashed lines).  We want to explore the hypothesis that
  obesity is a risk factor for chronic lower-back pain (which includes
  buttock pain in this definition).  A corollary to this hypothesis is
  that obesity is not considered to be a risk factor for chronic neck
  pain.  Although there is a generous amount of uncertainty, it does
  not appear that the probability of chronic lower-back pain increases
  with BMI for either gender.  Conversely, chronic neck pain does
  appear to be rising, yet again, the intervals are wide.  In both
  cases, these findings are not anticipated.  }
\end{figure}

\begin{figure}%[t!]
\begin{center}
\includegraphics[scale=0.425]{figures/chronic-pain2.pdf}
\end{center}
\caption{\label{chronic-pain2}NHANES, BMI and the probability of
  chronic pain for females only: the left panel for lower-back
  pain and the right panel for neck pain. The unweighted Friedman's
  partial dependence relationship between chronic pain and BMI are
  displayed as ascertained from NHANES data for females only:
  lower-back (blue) and neck pain (red) are presented with the
  corresponding 95\% credible intervals (dashed lines).  The
  difference in probability of chronic pain from a baseline BMI of 25
  (which is the upper limit of normal) is presented, i.e.,
  $p(x)-p(25)$.  We want to explore the hypothesis that obesity is a
  risk factor for chronic lower-back pain (which includes buttock pain
  in this definition).  A corollary to this hypothesis is that obesity
  is not considered to be a risk factor for chronic neck pain.
  Although there is a generous amount of uncertainty, it does not
  appear that the probability of chronic lower-back pain increases
  with BMI.  Conversely, chronic neck pain does appear to be rising,
  yet again, the intervals are wide.  In both cases, these findings
  are not anticipated. }
\end{figure}


\subsection{Logit BART for binary outcomes}

Assuming a Normal distribution of the unobserved latent, $z_i$ where
$y_i=\I{z_i>0}$, provides some challenges when estimating very small
or very large probabilities, $p_i$, since the Normal distribution has
relatively thin tails.  This restriction can be relaxed by assuming
the latents follow the Logistic distribution which has heavier tails. 
% and, therefore, is a better
%choice if the $p_i$ can be very close to zero or one.  
For Logistic latents, we employ a variant of the \citet{HolmHeld06}
technique by \citet{GramPols12} to create what we call logit BART.
However, it is important to recognize that logit BART is more
computationally intensive than probit BART.
  
The outcome, \code{y.train}, is provided as an integer with values 0
or 1.  Logit BART is provided by the \code{lbart} and \code{gbart}
functions.  Unlike probit BART where the auxiliary latents, $z_i$,
have a unit variance $\sd^2=1$; with Logisitic BART, we sample
truncated Normal latents, $z_i$, with a variance $\sd_i^2$ by the
\citet{Robe95} technique.  If ${\sd^2_i} = 4\psi_i^2 \where \psi_i$ is
sampled from the Kolmogorov-Smirnov distribution, then $z_i$ follow
the Logistic distribution.  Sampling from the Kolmogorov-Smirnov
distribution is described by \cite{Devr86}.  So, the conditionally
Normal latents, $z_i|\sd_i^2$, are the outcomes for a continuous BART
with a given heteroskedastic variance, $\sd_i^2$.  % Since Logistic
% latents are more flexible, there is no centering parameter, i.e.,
% $\mu_0=0$.  Therefore, the probabilities are $p_i=F(f(\bm{x}_i))$
% where $F$ is the standard Logistic distribution function.

The $z_i$ are centered around a known constant, $\mu_0$, which is
analagous to quasi-centering the probabilities, $p_i$, around
$p_0=F(\mu_0)$ where $F$ is the standard Logistic distribution
function.  The default value of $\mu_0$ is $F^{-1}(\bar{y})$ (which
you can over-ride with the \code{binaryOffset} argument to
\code{lbart} or the \code{offset} argument to \code{gbart}).
Therefore, the probabilities are $p_i=F(\mu_0+f(\bm{x}_i))$.

The input and output for \code{lbart} is essentially identical to
\code{pbart}.  Also, the \code{predict} function for objects of type
\code{lbart} is analogous.  The \code{gbart} function performs logit
BART when passed the \code{type='lbart'} argument.  

N.B.\ for
\code{lbart}, the thinning argument, \code{keepevery} defaults to 1
while for \code{gbart} with \code{type='lbart'}, \code{keepevery}
defaults to 10.

\subsection{Multinomial BART for categorical outcomes}

Several strategies for analyzing categorical outcomes have been
proposed from the Bayesian perspective
\citep{AlbeChib93,McCuRoss94,McCuPols00,ImaiVanD05,FruhFruh10,Scot11}
including two BART implementations \citep{KindWang16,Murr17}: our BART
implementations differ from these; although, since we are working on
the same problem, there are some similarities.
Generally, the literature has taken a logit approach.  Due to
the relative computational efficiency, we prefer probit to logit
(although, logit is available as an option).  To extend BART to
categorical outcomes, we have created two approaches to what we call
Multinomial BART.  The first approach works well when they are
relatively few categories while the second is preferable otherwise.

\subsubsection{Multinomial BART and conditional probability: mbart}

In the first approach, we fit a novel sequence of binary BART models
that bears some resemblance to continuation-ratio logits
\citep{Agre03}.  Let's assume that we have $K$ categories where each
are represented by mutually exclusive binary indicators:
$y_{i1}, \dots, y_{iK}$ for subjects indexed by $i=1, \dots, N$.  We
denote the probability of these outcome indicators via conditional
probabilities, $p_{ij}$, where $j=1, \dots, K$ as follows.
\begin{align*}
p_{i1}&=\P{y_{i1}=1}\\
p_{i2}&=\P{y_{i2}=1|y_{i1}=0} \\
p_{i3}&=\P{y_{i3}=1|y_{i1}=y_{i2}=0} \\
\vdots & \\
p_{i,K-1}&=\P{y_{i,K-1}=1|y_{i1}=\cdots=y_{i,K-2}=0} \\
p_{iK}&=\P{y_{i,K-1}=0|y_{i1}=\cdots=y_{i,K-2}=0} 
\end{align*}
Notice that $p_{iK}=1-p_{i,K-1}$ so we can specify the $K$
conditional probabilities via $K-1$ parameters.  Furthermore, these
conditional probabilities are, by construction, defined for subsets of
subjects: let $S_1=\{1, \dots, N\}$ and
$S_j = \{ i : y_{i1}=\cdots=y_{i,j-1}=0 \} \where j=2, \dots, K-1$.
Now, the unconditional probability of these outcome indicators,
$\pi_{ij}$, can be defined in terms of the conditional probablities
and their complements, $q_{ij}=1-p_{ij}$, for all subjects.
\begin{align*}
\pi_{i1}&=\P{y_{i1}=1}=p_{i1} \\
\pi_{i2}&=\P{y_{i2}=1}=p_{i2}q_{i1} \\
\pi_{i3}&=\P{y_{i3}=1}=p_{i3}q_{i2}q_{i1} \\
\vdots & \\
\pi_{i,K-1}&=\P{y_{i,K-1}=1}=p_{i,K-1}q_{i,K-2}\cdots q_{i1} \\
\pi_{iK}&=\P{y_{iK}=1}=q_{i,K-1}q_{i,K-2}\cdots q_{i1} 
\end{align*}
N.B.\ the conditional probability construction of $\pi_{ij}$
ensures that $\sum_{j=1}^K\pi_{ij}=1$.

Our modelling of these conditional probabilities based on a vector of
covariates $\bm{x}_i$ is what we call Multinomial BART:
\begin{align*}
y_{ij}|p_{ij} &\; ~ \; \B{p_{ij}} \where i \in S_j \mbox{ and } j=1, \dots, K-1 \\
p_{ij} & = \Phi(\mu_j+f_j(\bm{x}_i)) \\
 f_j & \prior \mathrm{BART}
\end{align*}
with $i$ indexing subjects, $i=1, \dots, N$; and the default value of
$\mu_j=\Phi^{-1}\wrap{\frac{\sum_i y_{ij}}{\sum_i
    \mathrm{I}\,\wrap[()]{i \in S_j}}}$.
This formulation yields the Multinomial likelihood:
$\wrap{\bm{y}|f_1,\dots,f_{K-1}} = \prod_{i=1}^N \prod_{j=1}^K
\pi_{ij}^{y_{ij}}$.

This approach is provided by the \pkg{BART} package as the
\code{mbart} function.  The input for \code{mbart} is essentially
identical to \code{gbart}, but the output is slightly different.  For
example, due to the way the model is estimated, the prediction for
\code{x.train} is not available; therefore, to request it set the
argument \code{x.test=x.train}.  By default, probit BART is employed
for computational efficiency, but logit BART can be specified with the
argument \code{type='lbart'}.
%Also, the \code{predict} function for objects of type \code{lbart} is analogous.
Notation: $M$~for the number of posterior samples, $B$~for the number
of threads (generally, $B=1$ for Windows), $N$~for the number of
observations in the training set, and $Q$~for the number of
observations in the test set.
\begin{Sinput}
R> set.seed(99)
R> post <- mbart(x.train, y.train, x.test, ndpost=M) 
R> post <- mc.mbart(x.train, y.train, x.test, ndpost=M, mc.cores=B, seed=99) 
\end{Sinput}
\begin{align*}
\mbox{Input: \code{x.train} and \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{\code{post}, of type \code{mbart}} & \\
\mbox{\code{post\$prob.test}:\ } &
\wrap{\begin{array}{ccccccc}
\hat{\pi}_{111} & \dots & \hat{\pi}_{1K1} & \dots & \hat{\pi}_{Q11} & \dots & \hat{\pi}_{QK1}\\
\vdots       & \ddots & \vdots    & \ddots & \vdots   & \ddots & \vdots \\
\hat{\pi}_{11M} & \dots & \hat{\pi}_{1KM} & \dots & \hat{\pi}_{Q1M} & \dots & \hat{\pi}_{QKM} \\
\end{array}} 
%& \where  \hat{p}_{ijm}=\Phi(f_{jm}(\bm{x}_i)) \mbox{ for \code{type='pbart'}}\\
%&\mbox{ or } \hat{p}_{ijm}=F(f_{jm}(\bm{x}_i)) \mbox{ for \code{type='lbart'}}\\
%&\mbox{ and $F(.)$ is the CDF of the Logistic distribution.}\\
\end{align*}
The columns of %\code{post\$prob.train} and
\code{post\$prob.test} represent different covariate settings crossed
with the $K$ categories.  The \code{predict} function for objects of
type \code{mbart} is analogous.

\subsubsection{Multinomial BART and the logit transformation: mbart2}

The second approach is inspired by the logit transformation and is
provided by the \code{mbart2} function which has a similar calling
convention to \code{mbart} described above.  Furthermore, as we shall
see, the computationally friendly probit is even applicable in this
instance.  Here, $y_i$ is categorical, i.e., $y_i \in \{1, \dots, K\}$
(technically, the \code{mbart2} function does not require the
categories to be $1, \dots, K$; it only requires that there are $K$
distinct categories).  Now, we have the following framework motivated
by the logit transformation.
\begin{align*}
%y_i|\bm{\pi}_i & ~ \M{1}{\bm{\pi}_i} \\
\P{y_i=j} & = \frac{\exp(\mu_j + f_j(\bm{x}_i))}
{\sum_{j'=1}^K \exp(\mu_{j'}+ f_{j'}(\bm{x}_i))} 
= \pi_{ij} \\
& \where f_j \prior \mathrm{BART},\ j=1, \dots, K
\end{align*}
Suppose for the moment, the centering parameters, $\mu_j$, are defined
as in logit BART.

It would appear that this definition has an identifiability issue
since
$\pi_{ij}= \frac{\exp(\mu_j + f_j(\bm{x}_i))} {\sum_{j'=1}^K
  \exp(\mu_{j'}+ f_{j'}(\bm{x}_i))} = \frac{\exp(\mu_j + f_j(\bm{x}_i)
  + c)} {\sum_{j'=1}^K \exp(\mu_{j'}+ f_{j'}(\bm{x}_i) + c)}$.
Identifiability could be restored by setting a single BART function to
zero, i.e., $f_{j'}(\bm{x}_i)=0$.  However, this is really unnecessary
since $\pi_{ij}$ is identified regardless.  % Although, the likelihood is 
% unidentified, due to our proper prior, we generate a proper posterior.

Computationally, this inference can be performed via a series of 
binary BARTs.  This can be shown by following the work of
\citep{HolmHeld06}: define $\P{y_i=c} \propto \exp f_c(\bm{x}_i)$.
Consider two cases: $\P{y_i=c}$ and $\P{y_i=j} \where j\not=c$.
The first case gives us the following in terms of $f_c$.
\begin{align*}
\P{y_i=c} & = \frac{\exp f_c(\bm{x}_i)}{\exp f_c(\bm{x}_i) 
+ \sum_{k\not=c} \exp f_k(\bm{x}_i)} \\
& = \frac{\exp f_c(\bm{x}_i)}{\exp f_c(\bm{x}_i) + \exp S}
\where S=\log \sum_{k\not=c} \exp f_k(\bm{x}_i) \\
& = \frac{ \exp -S}{ \exp -S}\frac{\exp f_c(\bm{x}_i)}
{\exp f_c(\bm{x}_i) + \exp S} \\
& = \frac{\exp (f_c(\bm{x}_i)-S)} {\exp (f_c(\bm{x}_i) -S)  +1}
\end{align*} \\
And the second case, where $j\not=c$, is as follows in terms of $f_c$.\\
\begin{align*}
\P{y_i=j} & = \frac{\exp f_j(\bm{x}_i)}{\exp f_c(\bm{x}_i) 
+ \sum_{k\not=c} \exp f_k(\bm{x}_i)} \\
& \propto \frac{1}{\exp f_c(\bm{x}_i) + \exp S} \\
& \propto \frac{1}{ \exp -S}\frac{1} {\exp f_c(\bm{x}_i) + \exp S} \\
& = \frac{1} {\exp (f_c(\bm{x}_i) -S)  +1}
\end{align*} 
Thus, the conditional inference for $f_c$ is equivalent to a binary indicator
$\I{y=c}$.  Therefore, \code{mbart2}
computes a full series of all $K$ BART functions for binary indicators.

The \code{mbart2} function defaults to \code{type='lbart'}, i.e.,
Logistic latents are used to compute the $f_j$'s which fits nicely
with the logit development of this approach.  However, the Logistic
latent fitting method can be computationally demanding.  Therefore,
Normal latents can be specified by \code{type='pbart'}.  This latter
setting would appear to contradict the development of this approach;
but notice that $\pi_{ij}$ is still a probability in this case and,
in our experience, the results produced are often reasonable.

\subsubsection{Multinomial BART example:  alligator food preference}

We demonstrate the usage of these functions by the American alligator
food preference example \citep{DelaLind99,Agre03}. In 1985, American
alligators were harvested by hunters from August 26 to September 30 in
peninsular Florida from lakes Oklawaha (Putnam County), George (Putnam
and Volusia counties), Hancock (Polk County) and Trafford (Collier
County). Lake, length and sex were recorded for each
alligator. Stomachs from a sample of alligators 1.09:3.89m long were
frozen prior to analysis. After thawing, stomach contents were removed
and separated and food items were identified and tallied. Volumes were
determined by water displacement. The stomach contents of 219
alligators were classified into five categories of primary food
preference: bird, fish (the most common primary food choice),
invertebrate (snails, insects, crayfish, etc.), reptile (turtles,
alligators), bird, and other (amphibians, plants, household pets,
stones, and other debris).  The length of alligators was dichotomized
into small, $\le$2.3m, vs.\ large, $>$2.3m.  We estimate the
probability of each food preference category for the marginal effect
of size by resorting to Friedman's partial dependence function
\citep{Frie01}.  We have supplied Figure~\ref{alligator} which
summarizes the BART results generated by the example
\code{alligator.R}: you can find this demo with the command
\code{demo("alligator", package="BART")}.  The
\code{mbart} function was used since the number of categories is
small.  The 95\% credible intervals are wide, but it appears that
large alligators are more likely to rely on a diet of fish while small
alligators are more likely to rely on invertebrates.  Although the
true probabilities are obviously unknown, we compared \code{mbart} to
an analysis by a single hidden-layer/feed-forward Neural Network via
the \pkg{nnet} \proglang{R} package \citep{Ripl07,VenaRipl13} and the
results were essentially identical (see the \code{demo} for details).
\begin{figure}%[t!]
\begin{center}
\includegraphics[scale=0.425]{figures/alligator.pdf}
\end{center}
\caption{\label{alligator}In 1985, American alligators were harvested
  by hunters in peninsular Florida from four lakes. Lake, length and
  sex were recorded for each alligator.  The stomach contents of 219
  alligators were classified into five categories based on the primary
  food preference: bird, fish, invertebrate, reptile and other.  The
  length of alligators was dichotomized into small, $\le$2.3m, vs.\
  large, $>$2.3m.  We estimate the probability of each food preference
  category for the marginal effect of size by resorting to Friedman's
  partial dependence function \citep{Frie01}.  The 95\% credible
  intervals are wide, but it appears that large alligators are more
  likely to rely on a diet of fish while small alligators are more
  likely to rely on invertebrates.  }
\end{figure}

\subsection{Converegence diagnostics for binary and categorical outcomes}\label{geweke}

How do you perform convergence diagnostics for BART?  For continuous
outcomes, convegence can easily be determined from the trace plots of
the the error standard deviation, $\sigma$.  However, for probit and
Multinomial BART with Normal latents, the error variance is fixed at 1
so this is not an option.  Similarly, for logit BART, $\sigma_i$, are
auxiliary latent variables not suitable for convergence diagnostics.
Therefore, we adapt traditional MCMC diagnostic approaches to BART.
We perform graphical checks via auto-correlation, trace plots and an
approach due to \cite{Gewe92}.

Geweke diagnostics are based on earlier work which characterizes
MCMC as a time series \citep{Hast70}.  Once this transition is
made, auto-regressive, moving-average (ARMA) process theory is
employed \citep{Silv86}.  Generally, we define our Bayesian estimator
as $\hat\theta_M = M^{-1}\sum_{m=1}^M \theta_m$.  We represent the
asymptotic variance of the estimator by
$ \sd^2_{\hat\theta} =\lim_{M \rightarrow \infty} \V{\hat\theta_M}$.
If we suppose that $\theta_m$ is an $\mathrm{ARMA}(p, q)$ process, then the
spectral density of the estimator is defined as
$\gamma(w) = (2\pi)^{-1} \sum_{m=-\infty}^{\infty} \V{\theta_0,
  \theta_m} \e{\mathrm{i} m w}$
where $\e{\mathrm{i} t w}=\cos(t w)+\mathrm{i} \sin(t w)$.  This
leads us to an estimator of the asymptotic variance which is
$\hat\sd^2_{\hat\theta} = {\hat{\gamma}^2(0)}$.
We divide our chain into two segments, $A$ and $B$, as follows:
$m \in A=\{1, \dots, M_A\} \where M_A= a M$; and
$m \in B=\{M-M_B+1, \dots, M\} \where M_B= b M$.  Note that 
$a+b<1$. Geweke suggests $a=0.1,\ b=0.5$ and recommends the
following Normal test for convergence.
\begin{align*}
\hat\theta_A & = M_A^{-1}\sum_{m \in A} \theta_m  &
\hat\theta_B & = M_B^{-1}\sum_{m \in B} \theta_m  \\
& \\
\hat\sd^2_{\hat\theta_A} & = \hat{\gamma}_{m \in A}^2(0) &
\hat\sd^2_{\hat\theta_B} & = \hat{\gamma}_{m \in B}^2(0) \\
& \\
Z_{AB} & = \frac{\sqrt{M}(\hat\theta_A-\hat\theta_B)}
{\sqrt{a^{-1}\hat\sd^2_{\hat\theta_A}+b^{-1}\hat\sd^2_{\hat\theta_B}}}~{\N{0}{1}}
%                                                  &  & 
\end{align*}

In our \pkg{BART} package, we supply \proglang{R} functions adapted
from the \pkg{coda} \proglang{R} package \citep{PlumBest06} to perform
Geweke diagnostics: {\tt spectrum0ar} and {\tt gewekediag}.  But, how
do we apply Geweke's diagnostic to BART?  We can check convergence for
any estimator of the form ${\theta}=h(f(\bm{x}))$, but often setting
$h$ to the identify function will suffice, i.e., ${\theta}=f(\bm{x})$.
However, BART being a Bayesian nonparametric technique means that we
have many potential estimators to check, i.e., essentially one
estimator for every possible choice of $\bm{x}$.

We have supplied Figures~\ref{geweke-pbart2-200},
\ref{geweke-pbart2-1000} and \ref{geweke-pbart2-5000} generated by
the example \code{geweke.pbart2.R}:\\ \code{demo("geweke.pbart2",
  package="BART")}.  The data are simulated by Friedman's
five-dimensional test function \citep{Frie91} where 50 covariates are
generated as $x_{ij}~\U{0, 1}$ but only the first 5 covariates have an
impact on the outcome at sample sizes $N = 200, 1000, 5000$.
\begin{align*}
{f(\bm{x}_i)} & = -1.5+\sin(\pi x_{i1} x_{i2}) + 2 (x_{i3}-0.5)^2 +x_{i4}+0.5x_{i5} \\
z_i & ~\N{f(\bm{x}_i)}{1} \\
{y_i} & =\I{z_i>0} 
\end{align*}
\begin{figure}%[t!]
\begin{center}
\includegraphics[scale=0.3]{figures/geweke-pbart2-200.pdf}
\end{center}
\caption{\label{geweke-pbart2-200}Geweke convergence diagnostics for
  probit BART: $N=200$.  In the upper left quadrant, we have plotted
  Friedman's partial dependence function for $f(x_{i4})$ vs.\ $x_{i4}$
  for 10 values of $x_{i4}$.  This is a check that can't be performed
  for real data, but it is informative in this case.  Notice that
  $f(x_{i4})$ vs.\ $x_{i4}$ is mainly directly proportional expected.
  In the upper right quadrant, we plot the auto-correlations of
  $f(\bm{x}_i)$ for 10 randomly selected $\bm{x}_i$ where $i$ indexes
  subjects.  Notice that there is very little auto-correlation.  In
  the lower left quadrant, we display the corresponding trace plots
  for these same settings.  The traces demonstrate that samples of
  $f(\bm{x}_i)$ appear to adequately traverse the sample space.  In
  the lower right quadrant, we plot the Geweke $Z_{AB}$ statistics for
  each subject $i$.  Notice that the $Z_{AB}$ exceed the 95\% limits
  only a handful of times.  Based on this figure, we conclude that the
  chains have converged.  }
\end{figure}

\begin{figure}%[t!]
\begin{center}
\includegraphics[scale=0.425]{figures/geweke-pbart2-1000.pdf}
\end{center}
\caption{\label{geweke-pbart2-1000}Geweke convergence diagnostics for
  probit BART: $N=1000$.  In the upper left quadrant, we have plotted
  Friedman's partial dependence function for $f(x_{i4})$ vs.\ $x_{i4}$
  for 10 values of $x_{i4}$.  This is a check that can't be performed
  for real data, but it is informative in this case.  Notice that
  $f(x_{i4})$ vs.\ $x_{i4}$ is directly proportional as expected.  In
  the upper right quadrant, we plot the auto-correlations of
  $f(\bm{x}_i)$ for 10 randomly selected $\bm{x}_i$ where $i$ indexes
  subjects.  Notice that there is very little auto-correlation.  In
  the lower left quadrant, we display the corresponding trace plots
  for these same settings.  The traces demonstrate that samples of
  $f(\bm{x}_i)$ appear to adequately traverse the sample space.  In
  the lower right quadrant, we plot the Geweke $Z_{AB}$ statistics for
  each subject $i$.  Notice that there appear to be a considerable
  number exceeding the 95\% limits.  Based on this figure, we conclude
  that convergence is questionable.  We would suggest that more
  thinning be employed via the \code{keepevery} argument to
  \code{pbart}; perhaps, \code{keepevery=50}.  }
\end{figure}

\begin{figure}%[t!]
\begin{center}
\includegraphics[scale=0.425]{figures/geweke-pbart2-5000.pdf}
\end{center}
\caption{\label{geweke-pbart2-5000}Geweke convergence diagnostics for
  probit BART: $N=5000$.  In the upper left quadrant, we have plotted
  Friedman's partial dependence function for $f(x_{i4})$ vs.\ $x_{i4}$
  for 10 values of $x_{i4}$.  This is a check that can't be performed
  for real data, but it is informative in this case.  Notice that
  $f(x_{i4})$ vs.\ $x_{i4}$ is directly proportional as expected.  In
  the upper right quadrant, we plot the auto-correlations of
  $f(\bm{x}_i)$ for 10 randomly selected $\bm{x}_i$ where $i$ indexes
  subjects.  Notice that there is some auto-correlation.  In the lower
  left quadrant, we display the corresponding trace plots for these
  same settings.  The traces demonstrate that samples of $f(\bm{x}_i)$
  appear to traverse the sample space, but there are some slower
  oscillations.  In the lower right quadrant, we plot the Geweke
  $Z_{AB}$ statistics for each subject $i$.  Notice that there appear
  to be far too many exceeding the 95\% limits.  Based on these
  figures, we conclude that convergence has not been attained.  We
  would suggest that more thinning be employed via the
  \code{keepevery} argument to \code{pbart}; perhaps,
  \code{keepevery=250}.  }
\end{figure}

The convergence for each of these data sets is graphically displayed
in Figures~\ref{geweke-pbart2-200}, \ref{geweke-pbart2-1000} and
\ref{geweke-pbart2-5000} where each figure is broken into four
quadrants.  In the upper left quadrant, we have plotted Friedman's
partial dependence function for $f(x_{i4})$ vs.\ $x_{i4}$ for 10 values of
$x_{i4}$.  This is a check that can't be performed for real data, but it
is informative in this case.  Notice that $f(x_{i4})$ vs.\ $x_{i4}$ is
directly proportional in each figure as expected.  In the upper right
quadrant, we plot the auto-correlations of $f(\bm{x}_i)$ for 10
randomly selected $\bm{x}_i$ where $i$ indexes subjects.  Notice that
there is very little auto-correlation for $N=200, 1000$, but a more
notable amount for $N=5000$.  In the lower left quadrant, we display
the corresponding trace plots for these same settings.  The traces
demonstrate that samples of $f(\bm{x}_i)$ appear to adequately
traverse the sample space for $N=200, 1000$, but less notably for
$N=5000$.  In the lower right quadrant, we plot the Geweke $Z_{AB}$
statistics for each subject $i$.  Notice that for $N=200$, the
$Z_{AB}$ exceed the 95\% limits only a handful of times.  Although,
there are 10 times more comparisons, $N=1000$ has seemingly more than
10 times as many values exceeding the 95\% limits.  And, for
$N=5000$, there are dramatically more values exceeding the 95\%
limits.  Based on these figures, we conclude that the chains have
converged for $N=200$; for $N=1000$, convergence is questionable; and,
for $N=5000$, convergence has not been attained.  We would suggest
that more thinning be employed for $N=1000, 5000$ via the
\code{keepevery} argument to \code{pbart}; perhaps,
\code{keepevery=50} for $N=1000$ and \code{keepevery=250} for
$N=5000$.

\subsection{BART and variable selection}

Bayesian variable selection techniques applicable to BART have been
studied by
\cite{ChipGeor10,ChipGeor13,BleiKape14,HahnCarv15,McCuCarv15,Line16}.
The \pkg{BART} package supports the sparse prior of \citet{Line16} by
specifying \code{sparse=TRUE} (the default is \code{sparse=FALSE}).
Let's represent the variable selection probabilities by
$s_j \where j=1, \., P$.  Now, replace the uniform variable selection
prior in BART with a Dirichlet prior.  Also, place a Beta prior on the
$\theta$ parameter.
\begin{align*}
\wrap{s_1, \., s_P}|\theta & \prior \Dir{\theta/P, \., \theta/P} \\
\frac{\theta}{\theta+\rho} & \prior \Beta{a}{b} 
\end{align*}
Typical settings are $b=1$ and $\rho=P$ (the defaults) which you can
over-ride with the \code{b} and \code{rho} arguments respectively.
The value $a=0.5$ (the default) is a sparse setting whereas an
alternative setting $a=1$ is not sparse; you can specify this
parameter with argument \code{a}.  If additional sparsity is desired,
then you can set the argument \code{rho} to a value smaller than $P$:
for more details, see Appendix~\ref{trees}.
Furthermore, Linero discusses two assumptions: Assumption 2.1 and
Assumption 2.2 (see \citet{Line16} for more details).  Basically,
Assumption 2.2 (2.1) is more (less) friendly to binary/ordinal
covariates and is (is not) the default corresponding to
\code{augment=FALSE} (\code{augment=TRUE}).

Let's return to the simulated probit BART example explored above which
is in the \pkg{BART} package: \code{demo("sparse.pbart",
  package="BART")}.  For sample sizes of $N=200, 1000, 5000$, there
are $P=100$ covariates, but only the first 5 are active.  In
Figure~\ref{varsel}, the 5 (95) active (inactive) covariates are red
(black) and circles (dots) are $>$ ($\le$) $P^{-1}$ which is chance
association represented by a black line.  For $N=200$, all five active
variables are identified, but notice that there are 20 false
positives.  For $N=1000$, all five active covariates are identified,
but notice that there are still 14 false positives.  For $N=5000$, all
five active covariates are identified and notice that there is only
one false positive.
\begin{figure}%[t!]
\begin{center}
\includegraphics[scale=0.425]{figures/sparse-pbart.pdf}
\end{center}
\caption{\label{varsel}Probit BART and variable selection example.
  For sample sizes of $N=200, 1000, 5000$, there are $P=100$
  covariates, but only the first 5 are active.  The 5 (95) active
  (inactive) covariates are red (black) and circles (dots) are $>$
  ($\le$) $P^{-1}$ which is chance association represented by a black
  line.  For $N=200$, all five active variables are identified, but
  notice that there are 20 false positives.  For $N=1000$, all five
  active covariates are identified, but notice that there are still 14
  false positives.  For $N=5000$, all five active covariates are
  identified and notice that there is only one false positive.  }
\end{figure}

We are often interested in the inter-relationship between covariates
within our model.  We can assess these relationships by inspecting the
binary trees.  For example, we can ascertain how often $x_1$ is chosen
as a branch decision rule leading to a branch decision rule with $x_2$
further up the tree or vice versa.  In this case, we call $x_1$ and
$x_2$ a concordant pair and we denote by $x_1 \leftrightarrow x_2$
which is a symmetric relationship, i.e., $x_1 \leftrightarrow x_2$
implies $x_2 \leftrightarrow x_1$.  If $\mathcal{B}_h$ is the number
of branches in tree $\mathcal{T}_h$, then the concordant pair
probability is:
$\kappa_{ij} = \P{x_i \leftrightarrow x_j \in
  \mathcal{T}_h|\mathcal{B}_h>1}$
for $i=1, \dots, P-1$ and $j=i+1, \dots, P$.  See an example of
calculating these probabilities in
\code{demo("trees.pbart", package="BART")}.

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

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

% \begin{comment}
% \subsection{Survival analysis with the Cox proportional hazard model}

% The inspiration for survival analysis with BART is the de facto
% standard: the Cox proportional hazard model \citep{Cox72}.  The data
% is $(s_i, \delta_i, \bm{x}_{i})$ where $i$ indexes subjects,
% $i=1, \dots, N$; $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_{(K)}<\infty$,
% and we consider the following time intervals: $(0,
% t_{(1)}], (t_{(1)}, t_{(2)}], \. (t_{(K-1)},
% t_{(K)}]$.  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 the parameters $\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.
% \end{comment}

\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.  Survival analysis has been studied by many, however,
most take a proportional hazards approach
\citep{Cox72,KalbPren80,KleiMoes06}.
%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,Fahr14}.  Relying on the
capabilities of BART, we do not stipulate a linear relationship with
the covariates nor proportional hazards.  

The data
is $(s_i, \delta_i, \bm{x}_{i})$ where $i$ indexes subjects,
$i=1, \dots, N$; $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_{(K)}<\infty$,
and we consider the following time intervals: $(0,
t_{(1)}], (t_{(1)}, t_{(2)}], \. (t_{(K-1)},
t_{(K)}]$. 

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.

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} & =  \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 $i$ indexing subjects, $i=1, \dots, N$; and
$\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), 
This formulation creates the likelihood of
$\wrap{\bm{y}|f}\ = \prod_{i=1}^N\prod_{j=1}^{n_i} p_{ij}^{y_{ij}}
(1-p_{ij})^{1-y_{ij}}$.
% (in contrast to the Cox model, this is the
% actual likelihood rather than the partial likelihood).
%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 either case, the default value of $\mu_0$
is $\Phi^{-1}(\bar{y})$
(which you can over-ride with the \code{offset} 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"}).
For computational efficiency, probit \citep{AlbeChib93} is the
default, but logit \citep{HolmHeld06,GramPols12} can be specified as
an option via \code{type="lbart"}.

Based on the posterior samples, we construct quantities of interest
with BART for survival analysis.  In discrete-time survival analysis,
the instantaneous hazard from continuous-time survival is essentially
replaced with the probability of an event in an interval, i.e., 
$p(t_{(j)},\bm{x})=\Phi(\mu_0+f(t_{(j)},\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;
and \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{Sinput}
R> times <- c(2.5, 1.5, 3.0)
R> delta <- c(  1,   1,   0)
R> surv.pre.bart(times=times, delta=delta)
\end{Sinput}
\begin{Soutput}
$y.train  $tx.train  $tx.test  $times    $K     
[1]              t          t  [1]       [1] 3  
    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{Soutput}  
%$

Here is a diagram of the input and output for the
\code{surv.pre.bart} function.  \code{pre} is a list that is generated
to contain the matrix \code{pre\$tx.train} and
the vector \code{pre\$y.train}. \\
\begin{Sinput}
R> pre <- surv.pre.bart(times, delta, x.train, x.test=x.train)  
\end{Sinput}
\begin{align*}
\mbox{\code{tx.train}\quad} & \mbox{\quad\code{y.train}} \\
\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
which is based on \code{gbart},
the default is \code{keepevery=10} 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. 
To avoid unnecessarily enlarged data sets, it is often prudent to
coarsen the time axis appropriately.  Although this might seem
drastic, times are often collected orders of magnitude more precisely
than necessary for the problem under study.  For example, cancer
registries often collect survival times in days while time in months 
or quarters would suffice for many typical applications.  
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 diagram of the input and output for the \code{surv.bart}
function for serial computation and \code{mc.surv.bart} for parallel
computation. 

Serial call  
\begin{Sinput}
R> set.seed(99)
R> post=surv.bart(x.train, times=times, delta=delta, x.test=x.train, 
+    ndpost=M)
\end{Sinput} 
Parallel call
\begin{Sinput}
R> post=mc.surv.bart(x.train, times=times, delta=delta, x.test=x.train, 
+    ndpost=M, mc.cores=B, seed=99) 
\end{Sinput}

% \begin{comment}
% \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{\ with rows of\ } \bm{x}_i \\
% \end{align*}
% \end{comment}
Input vector \code{times} with \code{K} distinct values and
\begin{align*}
%\mbox{Input vector \code{times}} &\\
%\mbox{with \code{K} distinct values} &\\
\mbox{\code{x.train}:\ } &
\wrap{\begin{array}{c}
\bm{x}_{1} \\
\bm{x}_{2} \\
\vdots \\
\bm{x}_{N} \\
\end{array}} \mbox{\ with rows of\ } \bm{x}_i \\
%\begin{align*}
%\mbox{Output \code{post}, type \code{survbart},} & \\
%\mbox{which is essentially a list } & \\
\mbox{Output \code{post},} & \\
\mbox{of type \code{survbart},} & \\
\mbox{is essentially a list } & \\
\mbox{including the matrix} & \\
\mbox{\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 \dots & \scriptstyle \hat{S}_1(t_{(K)}, \bm{x}_1) & \scriptstyle \dots
& \scriptstyle \hat{S}_1(t_{(1)}, \bm{x}_N)& \scriptstyle \dots & \scriptstyle \hat{S}_1(t_{(K)}, \bm{x}_N) \\
\vdots & \ddots & \vdots & \ddots & \vdots & \ddots & \vdots \\
\scriptstyle \hat{S}_M(t_{(1)}, \bm{x}_1)& \scriptstyle \dots & \scriptstyle \hat{S}_M(t_{(K)}, \bm{x}_1) & \scriptstyle \dots
& \scriptstyle \hat{S}_M(t_{(1)}, \bm{x}_N)& \scriptstyle \dots & \scriptstyle \hat{S}_M(t_{(K)}, \bm{x}_N) \\
\end{array}} & 
\end{align*}

Here is a diagram of the input and output for the \code{predict.survbart} function.\\
\begin{Sinput}
R> pred <- predict(post, pre$tx.test, mc.cores=B)
\end{Sinput}
%$
\begin{align*}
\mbox{Input: \code{x.test}\ } &
\wrap{\begin{array}{c}
\bm{x}_{1} \\
\bm{x}_{2} \\
\vdots \\
\bm{x}_{Q} \\
\end{array}} \mbox{\ with rows of \ } \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 \dots & \scriptstyle \hat{S}_1(t_{(K)}, \bm{x}_1) & \scriptstyle \dots
& \scriptstyle \hat{S}_1(t_{(1)}, \bm{x}_Q)& \scriptstyle \dots & \scriptstyle \hat{S}_1(t_{(K)}, \bm{x}_Q) \\
\vdots & \ddots & \vdots & \ddots & \vdots & \ddots & \vdots \\
\scriptstyle \hat{S}_M(t_{(1)}, \bm{x}_1)& \scriptstyle \dots & \scriptstyle \hat{S}_M(t_{(K)}, \bm{x}_1) & \scriptstyle \dots
& \scriptstyle \hat{S}_M(t_{(1)}, \bm{x}_Q)& \scriptstyle \dots & \scriptstyle \hat{S}_M(t_{(K)}, \bm{x}_Q) \\
\end{array}} 
\end{align*}

For an overview of Friedman's partial dependence function (including
the notation adopted in this article and its meaning), please see 
Section~\ref{Friedman-partial} which discusses continuous outcomes.
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 settings which, naturally, includes time,
$(t_{(j)}, \bm{x}_{hS})$.
% This marginal dependence function is defined by fixing time at
% $t_{(j)}$ and the subset at test observation setting $\bm{x}_{hS}$
% while aggregating over the observed settings of the complement
% covariates training observations:
% $f(t_{(j)}, \bm{x}_{hS})={N^{-1}}\sum_{i=1}^N f(t_{(j)}, \bm{x}_{hS},
% \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_{(j)}, \bm{x}_{hS}) = {N^{-1}} \sum_{i=1}^N S(t_{(j)},
\bm{x}_{hS}, \bm{x}_{iC})$.
% \begin{comment}
% 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 \pkg{BART} package; rather, we
% provide examples of the many practical uses in the \code{demo}
% directory.
% \end{comment}

\subsubsection{Survival analysis with BART example:  
advanced lung cancer}

Here we present an example that is available in the \pkg{BART} package:\\
\code{demo("lung.surv.bart", 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 logit BART and the results are
quite similar (not shown): 
\code{demo("lung.surv.lbart", package="BART")}.
Furthermore, we perform convergence diagnostics on the chain:\\
\code{demo("geweke.lung.surv.bart", package="BART")}.
\begin{figure}%[t!]
\begin{center}
\includegraphics[scale=0.45]{figures/lung.pdf}
\end{center}
\caption{\label{lung}Advanced lung cancer example: Friedman's partial
  dependence function with 95\% credible intervals: males (blue) vs.\
  females (red). A cohort of advanced lung cancer patients was
  recruited from the North Central Cancer Treatment Group. For
  survival time, these patients were followed for nearly 3 years or
  until lost to follow-up.}
\end{figure}

\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.
Suppose that we have two event times, $t_1$ and $t_2$, (let's say each
based on a different subject profile), then the concordance probability is
defined as $\kappa_{t_1, t_2}=\P{t_1<t_2}$.  A simple analytic example
with the Exponential distribution is as
follows. % where $t_1$ and $t_2$ are
%conditionally independent, i.e., $t_1|\lambda_1 \perp t_2|\lambda_2$,
\begin{align*}
t_i|\lambda_i & \ind \Exp{\lambda_i} \where i \in \{1, 2\} \\
\P{t_1<t_2|\lambda_1, \lambda_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|\lambda_1, \lambda_2} & = 1- \frac{\lambda_2}{\lambda_1+\lambda_2}
  = \frac{\lambda_1}{\lambda_1+\lambda_2}\\
& = \P{t_1<t_2|\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.  Suppose that we have two event times, $s_1$ and $s_2$, which
are conditionally independent, i.e.,
$s_1|(f, \bm{x}_1) \perp s_2|(f, \bm{x}_2)$.  First, we calculate
$ \P{s_1 < s_2|f, \bm{x}_1, \bm{x}_2} $ (from here on, we suppress $f$
and $\bm{x}_i$ for notational convenience).
\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 in this form.
Yet, we can arrive at 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{demo("concord.surv.bart", package="BART")}.

\subsection{Competing risks with BART}

Competing risks survival analysis
\citep{KalbPren80,FineGray99,KleiMoes06,NicovanH10,IshwGerd14,SparLoga19}
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.
We take two approaches to support competing risks with BART: both
approaches are extensions of BART survival analysis.  We flexibly
model the cause-specific hazards and eschew precarious restrictive
assumptions like linearity of covariate effects, proportionality
and/or parametric distributions of the outcomes.

\subsubsection{Competing risks with crisk.bart}

The first approach is supported by the function \code{crisk.bart} for
serial computation and\\ \code{mc.crisk.bart} for parallel
computation.  % These functions along with a related approach are
% described in \cite{SparLoga18}.
% We adapt 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 no such precarious
% restrictive assumptions.
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_{(K)} < \infty$.  % To accommodate
% 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, we create event indicators by melding absorbing events survival
analysis with mutually exclusive Multinomial categories where $i$
indexes subjects: $i=1, \dots, N$.
% 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*}
The likelihood is: $\wrap{\bm{y}|f_1, f_2} = \prod_{i=1}^N 
\prod_{j=1}^{n_i} p_{1ij}^{y_{1ij}} (1-p_{1ij})^{1-y_{1ij}}
\prod_{j'=1}^{n_i-y_{1in_i}} p_{2ij'}^{y_{2ij'}} (1-p_{2ij'})^{1-y_{2ij'}}$.
Now, 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}.

\subsubsection{Competing risks with crisk2.bart}

The second approach is supported by the function \code{crisk2.bart}
for serial computation and\\ \code{mc.crisk2.bart} for parallel
computation.  We take a similar approach as \citet{NicovanH10}.  We
model the probability for an event of either kind,
$p_{ij}=p(t_{(j)}, \bm{x}_i)$ (this is standard survival analysis);
and, given an event has occurred, the probability of a kind 1 event,
$\pi_{i}=\pi(t_{i}, \bm{x}_i)$.  Now, we create the corresponding
event indicators $y_{ij}$ and $u_i$ where $i$ indexes subjects:
$i=1, \dots, N$.
\begin{align*}
y_{ij} & = \I{\delta_i\not=0} \I{j=n_i} \where j=1, \dots, n_i\\ %  \where y_{1i0}=0 \\
y_{ij}|p_{ij} & ~ \B{p_{ij}} \\ % \where p_{1i0}=0 \\
 p_{ij} & = \Phi(\mu_y+f_y(t_{(j)}, {\bm{x}}_{i})) \where f_y \prior \mathrm{BART} \\
u_{i} & = \I{\delta_i=1} \where i \in \{i' : \delta_{i'}\not=0 \} \\
u_{i}|\pi_{i} & ~ \B{\pi_{i}} \\ %\where p_{2i0}=0\\
%y_{2ij}|(y_{1i\,j-1}=y_{1ij}=0, p_{2ij}) & ~ \B{p_{2ij}} \where p_{2i0}=0\\
 \pi_{i} & = \Phi(\mu_u+f_u(t_{i}, {\bm{x}}_{i})) \where f_u \prior \mathrm{BART} 
\end{align*}
The likelihood is: $\wrap{\bm{y},\bm{u}|f_y, f_u} = \prod_{i=1}^N 
\prod_{j=1}^{n_i} p_{ij}^{y_{ij}} (1-p_{ij})^{1-y_{ij}}
\prod_{i' : \delta_{i'}\not=0} \pi_{i'}^{u_{i'}} (1-\pi_{i'})^{1-u_{i'}}$.
Now, we can estimate the survival function and the cumulative
incidence functions similar to the first approach.
The returned object of type \code{crisk2bart} from \code{crisk2.bart}
or \code{mc.crisk2.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}.

\subsubsection{Competing risks with BART example: liver transplants}

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 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 to liver transplants 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.  Nevertheless, there still remains an acute shortage
of donor livers today.  The \code{transplant} data set is provided by
the \pkg{BART} \proglang{R} package as is this example:
\code{demo("liver.crisk.bart", 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{figure}%[t!]
\begin{center}
\includegraphics[scale=0.45]{figures/liver-BART.pdf}
\end{center}
\caption{ \label{liver-transplant}Liver transplant competing risks for
  type O patients estimated by BART and Aalen-Johansen. This data is
  from the Mayo Clinic liver transplant waiting list from 1990-1999.
  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 all blood types; whereas a
  donor liver from the other types will only be transplanted to a
  matching recipient.  Therefore, type O subjects on the waiting list
  were at a disadvantage since the pool of competitors was larger.}
\end{figure}

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

The \pkg{BART} package supports recurrent events \citep{SparRein18}
with \code{recur.bart} for serial computation and \code{mc.recur.bart}
for parallel computation.  Survival analysis is generally concerned
with absorbing events that a subject can only experience once like
mortality.  Recurrent events analysis is concerned with non-absorbing
events that a subject can experience more than once like hospital
admissions \citep{AndeGill82,WeiLin89,KalbPren02,SparRein18}.
Recurrent events analysis with 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 $K$ collectively distinct event start and end times for
all subjects by $0<t_{(1)}< \dots< t_{(K)}<\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 number 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 in Equation~\ref{recur}.  N.B.\ we follow the
recurrent events literature's favored terminology by using the term
``intensity'' rather than ``hazard'', but they 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 often 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,Fahr14}.  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 \pkg{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{Sinput}
R> library("BART")
R> times <- matrix(c(3, 8, 9, 4, 12, 12), nrow=2, ncol=3, byrow=TRUE)
R> tstop <- matrix(c(7, 8, 0, 7,  0,  0), nrow=2, ncol=3, byrow=TRUE)
R> delta <- matrix(c(1, 1, 0, 1,  0,  0), nrow=2, ncol=3, byrow=TRUE)
R> recur.pre.bart(times=times, delta=delta, tstop=tstop)
\end{Sinput}
\begin{Soutput}
$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{Soutput} 
%$
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 such that %correspond to 
$j \in R_i$ which is the risk set for subject $i$.
%$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}) \} }
% }$.
\begin{align*} 
R_i = \wrap[\{\}]{j : \wrap{ j \in \{1, \dots, n_i\} \mbox{\
and\ } {\cap_{k=1}^{N_i} \{t_{(j)} \not\in (t_{ik}, u_{ik}) \} } } }
\end{align*}
% 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)})$.
i.e., the risk set contains $j$ if $t_{(j)}$ is during the observation
period for subject $i$ and $t_{(j)}$ is not contained within an
already ongoing event for this subject.

Putting it all together, we arrive at the following recurrent events
discrete-time model with $i$ indexing subjects; $i=1, \dots, N$.
\begin{align}\label{recur}
y_{ij}|p_{ij} & ~  \B{p_{ij}} \where j \in R_i \notag \\
p_{ij} & = \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}
This produces the following likelihood:
$\wrap{\bm{y}|f} = \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$.
We center the BART function, $f$, %and, therefore, the probabilities, $p_{ij}$,
by $\mu_0=\Phi^{-1}\!\wrap[()]{\bar{y}}$ where
$\bar{y}=\frac{\sum_i \sum_{j \in R_i} y_{ij}} {\sum_i
  \sum_{j=1}^{n_i} \I[\,]{j \in R_i}}$.

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 for logit with
Logistic latents \citep{HolmHeld06,GramPols12} 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}}$.

%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 Equation~\ref{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}$.

\subsubsection{Recurrent events with BART example: bladder tumors}

An interesting example of recurrent events involves a clinical trial
conducted by the Veterans Administration Cooperative Urological
Research Group \citep{Byar90}.  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{demo("bladder.recur.bart", package="BART")}.  In
Figure~\ref{RI-Th-Pl}, notice that the relative intensity calculated
by Friedman's partial dependence function finds thiotepa inferior to
placebo from roughly 6 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 finds thiotepa inferior to vitamin B6 from roughly
3 to 24 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 superior to placebo throughout, but the
95\% credible intervals are wide; see Figure~\ref{RI-B6-Pl}.
\begin{figure}[h]
\begin{center}
\includegraphics[scale=0.45]{figures/RI-Th-Pl.pdf}
\end{center}
\caption{\label{RI-Th-Pl}Relative Intensity: Thiotepa vs.\ Placebo.
  The relative intensity function is as follows:
  $RI(t,\tilde{\bm{x}}_T(t),\tilde{\bm{x}}_P(t))=\frac{p(t,\tilde{\bm{x}}_T(t))}{p(t,\tilde{\bm{x}}_P(t))}$
  where $T$ is for Thiotepa and $P$ is for Placebo.  The blue lines
  are the relative intensity functions themselves and the red lines
  are their 95\% credible intervals.  The relative intensity is
  calculated by Friedman's partial dependence function, i.e., 
  aggregated over all other covariates.}
\end{figure}

\begin{figure}[h]
\begin{center}
\includegraphics[scale=0.45]{figures/RI-Th-B6.pdf}
\end{center}
\caption{\label{RI-Th-B6}  
Relative Intensity: Thiotepa vs.\ Vitamin B6.
  The relative intensity function is as follows:
  $RI(t,\tilde{\bm{x}}_T(t),\tilde{\bm{x}}_B(t))=\frac{p(t,\tilde{\bm{x}}_T(t))}{p(t,\tilde{\bm{x}}_{B}(t))}$
  where $T$ is for Thiotepa and $B$ is for Vitamin B6.  The blue lines
  are the relative intensity functions themselves and the red lines
  are their 95\% credible intervals.  The relative intensity is
  calculated by Friedman's partial dependence function, i.e., 
  aggregated over all other covariates.}
\end{figure}

\begin{figure}[h]
\begin{center}
\includegraphics[scale=0.45]{figures/RI-B6-Pl.pdf}
\end{center}
\caption{\label{RI-B6-Pl}
Relative Intensity: Vitamin B6 vs.\ Placebo.  
  The relative intensity function is as follows:
  $RI(t,\tilde{\bm{x}}_B(t),\tilde{\bm{x}}_P(t))=\frac{p(t,\tilde{\bm{x}}_B(t))}{p(t,\tilde{\bm{x}}_{P}(t))}$
  where $B$ is for Vitamin B6 and $P$ is for Placebo.  The blue lines
  are the relative intensity functions themselves and the red lines
  are their 95\% credible intervals.  The relative intensity is
  calculated by Friedman's partial dependence function, i.e., 
  aggregated over all other covariates.}
\end{figure}

\section{Discussion}

The \pkg{BART} \proglang{R} package provides a user-friendly reference
implementation of Bayesian Additive Regression Trees (BART).  BART is
a Bayesian nonparametric, tree-based ensemble, machine learning
technique with best-of-breed properties.  In the spirit of machine
learning, BART {\it learns} the relationship between the covariates,
$\bm{x}$, and the response variable arriving at $f(\bm{x})$ while not
burdening the user to pre-specify the functional form of $f$ nor the
interaction terms among the covariates.  By specifying an optional
sparse Dirichlet prior, BART is capable of variable selection: a form
of {\it learning} which is especially useful in high-dimensional
settings.  In the class of ensemble predictive models, BART's
out-of-sample predictive performance is competitive with other leading
members of this class.  Due to its membership in the class of Bayesian
nonparametric models, BART not only provides an estimate of
$f(\bm{x})$, but naturally generates the uncertainty as well.

There are user-friendly features that are inherent to BART itself
which, of course, are available in this package as well.  BART was
designed to be very flexible via its prior arguments while providing
the user robust, low information, default settings that will likely
produce a good fit without resorting to computationally demanding
cross-validation.  BART itself is relatively computationally
efficient, but larger data sets will naturally take more time to
estimate.  Therefore, the \pkg{BART} package provides the user with
simple and easy to use multi-threading to keep elapsed time to a
minimum.  Another important time-saver, the \pkg{BART} package allows
the user to save the trees from a BART model fit so that prediction
via the \proglang{R} \code{predict} function can take place at a later
time without having to re-fit the model.  And these predictions can
also take advantage of multi-threading.

The \pkg{BART} package has been written in \proglang{C++} for portability,
maintainability and efficiency; this allows BART to be called either
from \proglang{R} or from other computer source code written in many
languages.  The package supports missing data handling of the
covariates and provides the user with access to BART implementations
for several types of responses.  The \pkg{BART} package supports the
following:
\begin{itemize}
\item continuous outcomes;
\item binary outcomes via probit or logit transformation;
\item categorical outcomes;
\item time-to-event outcomes with right censoring including
\begin{itemize}
\item absorbing events,
\item competing risks, and
\item recurrent events.
\end{itemize}
\end{itemize}

In this article, we have provided the user with an overview of much
that is described in this section including (but not limited to): details of
the BART prior and its arguments, sparse variable selection,
prediction, multi-threading, support for the outcomes listed above and
missing data handling.  In addition, this article has provided primers
on important BART topics such as posterior computation, Friedman's partial
dependence function and convergence diagnostics.  With a
computational method such as BART, the user needs a reliable,
well-documented software package with a diverse set of examples.  With
this article, and the \pkg{BART} package itself, we believe that
interested users now have the tools to successfully employ
BART for their rigorous data analysis needs.

\clearpage
\appendix
\section{Getting and installing the BART R package}\label{BART}

The \pkg{BART} package \citep{McCuSpar18} is GNU General Public
License (GPL) software available on the Comprehensive \proglang{R} Archive
Network (CRAN).  You can install it from CRAN as follows.
\begin{Sinput}
R> options(repos=c(CRAN="https://cran.r-project.org"))
R> install.packages("BART", dependencies=TRUE)
\end{Sinput}
The examples in this article are included in the package.  You can run
the first example (described in Section~\ref{boston}) as follows.
%    ## there are diminishing returns so often 8 cores is sufficient
\begin{Sinput}
R> options(figures=".")
R> if(.Platform$OS.type=="unix") {
R>   options(mc.cores=min(8, parallel::detectCores()))
R> } else {
R>   options(mc.cores=1)
R> }
R> demo("boston", package="BART"))
\end{Sinput}
%$
%source(system.file('demo/boston.R', package='BART'))
As we shall see, these examples produce \proglang{R} objects
containing BART model fits.  But, these fits are Bayesian
nonparametric samples from the posterior and require statistical
summarization before they are readily interpretable.  Therefore, we
often employ graphical summaries (such as the figures in this article)
to visualize the BART model fit.  Note that the \code{figures} option
(in the code snippet above) specifies a directory where the Portable
Document Format (PDF) graphics files will be produced; if it is not
specified, then the graphics will be generated by \proglang{R},
however, no PDF files will be created.  Furthermore, some of these
BART model fits can take a few minutes so it is wise to utilize
multi-threading when it is available (for a discussion of efficient
computation with BART including multi-threading, see Appendix
Section~\ref{efficient}).  Returning to the snippet above, the option
\code{mc.cores} specifies the number of cores to employ in
multi-threading, e.g., there are diminishing returns so often 8
cores is sufficient.  And, finally, to run all of the examples in this
article (with the options as specified above), then do the following.
\code{demo("replication", package="BART")}

%  This requires delving into the details of the BART
% prior itself and the corresponding arguments to the \code{wbart} and
% \code{gbart} functions.

% \begin{comment}
% For continuous outcomes, 
% Bayesian Additive Regression Trees (BART) \citep{ChipGeor10}
% fit the basic model:
% $$ y_i = f(x_i) + \epsilon_i, \;\; \epsilon_i \sim N(0,\sigma^2) $$
% where $(f,\sigma) \prior \mathrm{BART}$.
% We use Markov chain Monte Carlo (MCMC) to get draws from the posterior
% distribution of the parameter $(f,\sigma)$.  In this section, we
% describe the functionality of \code{BART::wbart} which is the basic
% function in the \pkg{BART} \proglang{R} package.  But first, we
% delve into the details of the BART prior itself.
% \end{comment}

\section{Binary trees and the BART prior}\label{trees}

BART relies on an ensemble of $H$ binary trees which are a type of a
directed acyclic graph.  We exploit the wooden tree metaphor to its
fullest.  Each of these trees grows from the ground up starting out as
a root node.  The root node is generally a branch decision rule, but
it doesn't have to be; occasionally there are trees in the ensemble
which are only a root terminal node consisting of a single leaf output
value.  If the root is a branch decision rule, then it spawns a left
and a right node which each can be either a branch decision rule or a
terminal leaf value and so on.  In binary tree, $\mathcal{T}$, there
are $C$ nodes which are made of $B$ branches and $L$ leaves: $C=B+L$.
There is an algebraic relationship between the number of branches and
leaves which we express as $B= L-1$.

The ensemble of trees is encoded in an ASCII string
which is returned in the \code{treedraws\$trees} list item.  This
string can be easily %exported and 
imported by \proglang{R} with the following:\\
\begin{Sinput}
R> write(post$treedraws$trees, "trees.txt")
R> tc <- textConnection(post$treedraws$tree)
R> trees <- read.table(file=tc, fill=TRUE, row.names=NULL, header=FALSE,
+    col.names=c("node", "var", "cut", "leaf"))
R> close(tc)
R> head(trees)
\end{Sinput}
\begin{minipage}{9cm}
\begin{Soutput}
  node var cut         leaf
1 1000 200   1           NA
2    3  NA  NA           NA
3    1   0  66 -0.001032108
4    2   0   0  0.004806880
5    3   0   0  0.035709372
6    3  NA  NA           NA
\end{Soutput}
\end{minipage}
\begin{minipage}{4cm}
\usetikzlibrary{shadows}
 \begin{tikzpicture}
 [level distance=20mm,sibling distance=25mm,
   int/.style={fill=white,draw=black,drop shadow,circle,anchor=north},
   ter/.style={fill=white,rectangle,draw=black,drop shadow}]
 \node[int]  {$x_1$} [grow=up]
 child {node[ter] {$0.036$}
 edge from parent node [right,pos=0.3] {$> c_{1,67}$}}
 child {node[ter] {$0.005$}
 edge from parent node [left,pos=0.3] {$\le c_{1,67}$}};
 \end{tikzpicture}
\end{minipage}

The string is encoded via the following binary tree notation.  The
first line is an exception which has the number of MCMC samples, $M$,
in the field \code{node}; the number of trees, $H$, in the field
\code{var}; and the number of variables, $P$, in the field \code{cut}.
For the rest of the file, the field \code{node} is used for the number
of nodes in the tree when all other fields are \code{NA}; or for a
specific node when the other fields are present.  The nodes are
numbered in relation to the tree's tier level, 
$t(n)=\lfloor \log_2 n \rfloor$ or
\code{t=floor(log2(node))}, as follows.
\begin{table}[!h]\label{tree-schematic}
\begin{center}
\begin{tabular}{r|ccccccc} \hline
Tier & \\ %\hline
$t$ & \multicolumn{3}{c}{$2^t$} & $\dots$ &
\multicolumn{3}{c}{$2^{t+1}\!-\!1$} \\ 
$\vdots$ & \\
2 & 4 &   & 5 &   & 6 &   & 7 \\
1 &   & 2 &   &   &   & 3 &   \\
0 &   &   &   & 1 &   &   &   \\ \hline
\end{tabular}
\end{center}
\end{table}

The \code{var} field is the variable in the branch decision rule which
is encoded $0, \dots, P-1$ as a \proglang{C}/\proglang{C++} array
index (rather than an \proglang{R} index).  Similarly, the \code{cut}
field is the cutpoint of the variable in the branch decision rule
which is encoded $0, \dots, c_j-1$ for variable $j$; note that the
cutpoints are returned in the \code{treedraws\$cutpoints} list item.
The terminal leaf output value is contained in the field \code{leaf}.
It is not immediately obvious which nodes are branches vs.\ leaves
since, at first, it would appear that the \code{leaf} field is given
for both branches and leaves.  Leaves are always associated with
\code{var=0} and \code{cut=0}; however, note that this is also a valid
branch variable/cutpoint since these are \proglang{C}/\proglang{C++}
indices.
% The key insight is that the 
% first $B$ rows of each node are branches and the rest are leaves.
The key to discriminating between branches and leaves is via the
algebraic relationship between a branch, $n$, at tree tier $t(n)$
leading to its left, $l=2n$, and right, $r=2n+1$, nodes at tier
$t(n)+1$, i.e., for each node, besides root,
%$l=2^{t+1}+2k$ and $r=2^{t+1}+2k+1$, i.e., for each node, besides root,
you can determine from which branch it arose and those nodes that are
not a branch (since they have no leaves) are necessarily leaves.

%\subsection{The BART prior}

Underlying this methodology is the BART prior.  The BART prior
specifies a flexible class of unknown functions, $f$, from which we
can gather randomly generated fits to the given data via the
posterior.  N.B.\ we define $f$ as returning a scalar value, but BART
extensions which return multivariate values are conceivable.  Let the
function $g(\bm{x}; \mathcal{T}, \mathcal{M})$ assign a value based on
the input $\bm{x}$.  The binary decision tree $\mathcal{T}$ is
represented by a set of ordered triples, $(n, j, k)$, representing
branch decision rules: $n \in \mathcal{B}$ for node $n$ in the set of
branches $\mathcal{B}$, $j$ for covariate $x_j$ and $k$ for the
cutpoint $c_{jk}$.  The branch decision rules are of the form
$x_j\le c_{jk}$ which means branch left and $x_j>c_{jk}$, branch
right; or terminal leaves where it stops.  $\mathcal{M}$ represents
leaves and is a set of ordered pairs, $(n, \mu_n)$:
$n \in \mathcal{L}$ where $\mathcal{L}$ is the set of leaves
($\mathcal{L}$ is the complement of $\mathcal{B}$) and $\mu_n$ for the
outcome value.

%\{\mu_1, \dots, \mu_L\}$, the parameter values
%associated with the $L$ leaves one of which will be the output
%of $g$
%$T$ representing the structure of a binary tree,
%including interior decision rules as branches and terminal nodes as
%leaves; and $M=.  
The function, $f(\bm{x})$, is a sum of $H$ trees:
\begin{align}\label{BARTfunction}
f(\bm{x})=\sum_{h=1}^H g(\bm{x}; \mathcal{T}_h, \mathcal{M}_h)
\end{align} 
where $H$ is ``large'', let's say, 50, 100 or 200.

For a continuous outcome, $y_i$, we have the following BART regression
on the vector of covariates, $\bm{x}_i$:
\begin{align*}
y_i=\mu_0+f(\bm{x}_i)+\epsilon_i \where \epsilon_i \iid \N{0}{ w_i^2
  \sd^2}
\end{align*}
with $i$ indexing subjects $i=1, \dots, N$.  The unknown random
function, $f$, and the error variance, $\sd^2$, follow the BART prior
expressed notationally as
\begin{align*}
(f,\sd^2)\prior\mathrm{BART}(H, \mu_{0}, \tau, k, \alpha, \gamma;
\nu, \lambda, q)
\end{align*}
where $H$ is the number of trees, $\mu_0$ is a known constant which
centers ${y}$ and the rest of the parameters will be explained later
in this section (for brevity, we will often use the simpler shorthand
$(f,\sd^2)\prior\mathrm{BART}$).  The $w_i$ are known standard
deviation weight multiples which you can supply with the argument
\code{w} that is only available for continuous outcomes, hence, the
weighted BART name; the unit weight vector is the default.  The
centering parameter, $\mu_0$, can be specified via the \code{fmean}
argument where the default is taken to be $\b{y}$.

BART is a Bayesian nonparametric prior.  Using the Gelfand-Smith
generic bracket notation for the specification of random variable
distributions \citep{GelfSmit90}, we represent the BART prior in terms
of the collection of all trees, $\mathcal{T}$; collection of all leaves,
$\mathcal{M}$; and the error variance, $\sd^2$, as the following
product: $\wrap{\mathcal{T}, \mathcal{M}, \sd^2}=
\wrap{\sd^2}\wrap{\mathcal{T}, \mathcal{M}}=
\wrap{\sd^2}\wrap{\mathcal{T}}\wrap{\mathcal{M}|\mathcal{T}}$.
Furthermore, the individual trees themselves are independent:
$\wrap{\mathcal{T}, \mathcal{M}}=\prod_h
\wrap{\mathcal{T}_h}\wrap{\mathcal{M}_h|\mathcal{T}_h}$.
where $\wrap{\mathcal{T}_h}$ is the prior for the $h$th tree and
$\wrap{\mathcal{M}_h|\mathcal{T}_h}$ is the collection of leaves for
the $h$th tree.  And, finally, the collection of leaves for the
$h$th tree are independent:
$\wrap{\mathcal{M}_h|\mathcal{T}_h}=
\prod_n\wrap{\mu_{hn}|\mathcal{T}_h}$
where $n$ indexes the leaf nodes.

The tree prior: $\wrap{\mathcal{T}_h}$.  There are three prior
components of $\mathcal{T}_h$ which govern whether the tree branches
grow or are pruned.
% As can be seen
% in the tree schematic, Table~\ref{tree-schematic}, the nodes are
% numbered as follows: $n=1, \dots$.  
The first tree prior regularizes the probability of a branch at leaf
node $n$ in tree tier $t(n)=\lfloor\log_2 n\rfloor$ as
\begin{align}\label{regularity}
\P{B_n=1}=\alpha (t(n)+1)^{-\gamma}
\end{align}
where $B_n=1$ represents
a branch while $B_n=0$ is a leaf, $0<\alpha<1$ and $\gamma\ge 0$.  You
can specify these prior parameters with arguments, but the following
defaults are recommended: $\alpha$ is set by the parameter
\code{base=0.95} and $\gamma$ by \code{power=2}; for a detailed
discussion of these parameter settings, see \citet{ChipGeor98}.  Note
that this prior penalizes branch growth, i.e., in prior probability,
the default number of branches will likely be 1 or 2.  Next, there is
a prior dictating the choice of a splitting variable $j$ conditional on a
branch event $B_n$ which defaults to uniform probability $s_j=P^{-1}$ where
$P$ is the number of covariates (however, you can specify a Dirichlet
prior which is more appropriate if the number of covariates is large
\citep{Line16}; see below).  Given a branch event, $B_n$, and a
variable chosen, $x_j$, the last tree prior selects a cut point,
$c_{jk}$, within the range of observed values for $x_j$; this prior is
uniform.

We can also represent the probability of variable selection via the
sparse Dirichlet prior as
$\wrap{s_1, \dots, s_P}|\theta \prior \Dir{\theta/P, \dots, \theta/P}$ which
is specified by the argument \code{sparse=TRUE} while the default is
\code{sparse=FALSE} for uniform $s_j=P^{-1}$.  The prior parameter
$\theta$ can be fixed or random: supplying a positive number will specify
$\theta$ fixed at that value while the default \code{theta=0} is random
and its value will be learned from the data.  The random $\theta$ prior
is induced via $\theta/(\theta+\rho) \prior \Bet{a}{b}$ where the
parameter $\rho$ can be specified by the argument \code{rho} (which
defaults to \code{NULL} representing the value $P$; provide a value to
over-ride), the parameter $b$ defaults to 1 (which can be over-ridden
by the argument \code{b}) and the parameter $a$ defaults to 0.5 (which
can be over-ridden by the argument \code{a}).  The distribution of
\code{theta} controls the sparsity of the model: \code{a=0.5} induces
a sparse posture while \code{a=1} is not sparse and similar to the
uniform prior with probability $s_j=P^{-1}$.  
If additional sparsity is desired,
then you can set the argument \code{rho} to a value smaller than $P$.

Here, we take the opportunity to provide some insight into how and why
the sparse prior works as desired. The key to understanding the
inducement of sparsity is the distribution of the arguments to
the Dirichlet prior: $\theta/P$.  We are unaware of this result 
appearing elsewhere in the literature.  But, it can be
shown that $\theta/P ~ F(a, b, \rho/P)$ where $F(.)$ is the Beta Prime
distribution scaled by $\rho/P$ \citep{JohnKotz95}.  The non-sparse
setting is $(a, b, \rho/P)=(1, 1, 1)$.  As you can see in the
Figure~\ref{sparse}, sparsity is increased by reducing $\rho$, reducing
$a$ or reducing both.
\begin{figure}%[t!]
\begin{center}
\includegraphics[scale=0.48]{figures/sparse-beta-prime.pdf}
\end{center}
\caption{\label{sparse}
The distribution of $\theta/P$ and the sparse Dirichlet prior.
The key to understanding the inducement of sparsity is the 
distribution of the arguments to the Dirichlet prior: 
$\theta/P ~ F(a, b, \rho/P)$ where $F(.)$ is the Beta Prime
distribution scaled by $\rho/P$.  Here we plot the natural
logarithm of the scaled Beta Prime density, $f(.)$, at a
non-sparse setting and three sparse settings.  The non-sparse
setting is $(a, b, \rho/P)=(1, 1, 1)$ (solid black line).  As 
you can see in the figure, sparsity is increased by reducing $\rho$
(long dashed red line), reducing $a$ (short dashed blue line) or 
reducing both (mixed dashed gray line).
}
\end{figure}

Unlike matrices, data frames can contain categorical factors.
Therefore, factors can be supplied when \code{x.train} is a data
frame. Factors with multiple levels are transformed into dummy
variables with each level as their own binary
indicator; %while each continuous variable is a group all by itself
factors with only two levels are a binary indicator with a
single dummy variable.
%TODO
% By specifying \code{sparse=TRUE}, these groups are handled
% with a sparse Dirichlet prior for grouped variables
% \citep{YuanLin06,LineYang17}.  
% Suppose we have $P$ groups each with
% $K_j$ variables, then the probability of selecting 
% a particular variable is $u_{jk} = s_j t_{jk}$ where
% $\wrap{s_1, \dots, s_P} \prior \Dir{\theta/P, \dots, \theta/P}$ and
% $\wrap{t_{j1}, \dots, t_{jK_j}} \prior \Dir{\omega/K_j, \dots, \omega/K_j}$.
% % \citep{YuanLin06,LineYang17}.  Suppose we have $P$ groups each with
% % $K_j$ variables, then $s_{jk} = u_j v_{jk}$ where
% % $\wrap{u_1, \dots, u_P} \prior \Dir{\alpha/P, \dots, \alpha/P}$ and
% % $\wrap{v_{j1}, \dots, v_{jK_j}} \prior \Dir{\omega/K_j, \dots, \omega/K_j}$.
% The specification of the $\theta$ prior is as above.  The prior
% parameter $\omega$ is fixed and the default specification is set by
% the argument \code{omega=1}.

The leaf prior: $\wrap{\mu_{hn}|\mathcal{T}_h}$.  Given a tree, $\mathcal{T}_h$, there is
a prior on its leaf values, $\mu_{hn}|\mathcal{T}_h$ and we denote the
collection of all leaves in $\mathcal{T}_h$ by
$\mathcal{M}_h=\{(n, \mu_{hn}): n \in \mathcal{L}_h \}$.
%\wrap{\mu_{h1}, \dots, \mu_{hL_h}}$.  
Suppose that $y_i \in [y_{\min}, y_{\max}]$ for all $i$ and denote 
$\wrap{\mu_{1(i)}, \dots, \mu_{H(i)}}$ as the leaf output values from each 
tree corresponding to the vector of covariates, $\bm{x}_i$.
If $\mu_{h(i)}|\mathcal{T}_h \iid \N{0}{\sd_{\mu}^2}$, then the model 
estimate for subject~$i$ is 
$\mu_i=\E{y_i|\bm{x}_i}=\mu_0+\sum_h\mu_{h(i)}$ where
$\mu_i ~ \N{\mu_0}{H \sd_{\mu}^2}$.  We
%$\mu_{h(i)}|T_h \iid \N{\mu_{\mu}}{\sd_{\mu}^2}$, then
%$\E{y_i|\bm{x}_i}=\mu_i ~ \N{\mu_0+H\mu_{\mu}}{H \sd_{\mu}^2}$.  We
choose a value for $\sd_{\mu}$ which is the solution to the equations
%the system of equations created by the following $1-\alpha/2$ symmetric intervals: 
$y_{\min}=\mu_0-k\sqrt{H}\sd_{\mu}$
%$y_{\min}=\mu_0-|z_{\alpha/2}|\sqrt{H}\sd_{\mu}$
and $y_{\max}=\mu_0+k\sqrt{H}\sd_{\mu}$, i.e.,
%and $y_{\max}=\mu_0+|z_{\alpha/2}|\sqrt{H}\sd_{\mu}$, i.e.,
%$y_{\min}=\mu_0+H\mu_{\mu}-|z_{\alpha/2}|\sqrt{H}\sd_{\mu}$
%and $y_{\max}=\mu_0+H\mu_{\mu}+|z_{\alpha/2}|\sqrt{H}\sd_{\mu}$, i.e.,
%$\mu_{\mu}=\frac{y_{\max}-\mu_0+y_{\min}-\mu_0}{2H}$ and
$\sd_{\mu}=\frac{y_{\max}-y_{\min}}{2 k \sqrt{H}}$.
%$\sd_{\mu}=\frac{y_{\max}-y_{\min}}{2 |z_{\alpha/2}| \sqrt{H}}$.
%Since $y$ is centered around $\mu_0$, the solution for $\mu_{\mu}$
%will generally be near zero so we set it to zero.  
Therefore, we arrive at
$\mu_{hn} \prior \N{0}{\wrap{\frac{\tau}{2k\sqrt{H}}}^2} \where
\tau={y_{\max}-y_{\min}}$.  So, the prior for $\mu_{hn}$ is informed
by the data, $y$, but only weakly via the extrema,
$y_{\min}\mbox{\ and\ }y_{\max}$.  The parameter $k$ calibrates this
prior as follows.
\begin{align*}
\mu_i ~\N{\mu_0}{\wrap{\frac{\tau}{2 k}}^2}& \\
\P{y_{\min} \le \mu_i \le y_{\max}} &  = \Phi(k) - \Phi(-k)\\
\mbox{Since\ }\P{\mu_i \le y_{\max}} &= \P{z \le 2k \frac{y_{\max}-\mu_0}{\tau}} \approx
 \P{z \le k} = \Phi(k) \\
\mbox{Similarly\ }\P{\mu_i \le y_{\min}} &= \Phi(-k)
\end{align*}
The default value, $k=2$, corresponds to
$\mu_i$ falling within the extrema with approximately 0.95
probability.  Alternative choices of
$k$ can be supplied via the \code{k} argument.  We have found that
values of $k \in [1, 3]$ generally yield good results.  Note that
$k$ is a potential candidate parameter for choice via
cross-validation.

The error variance prior: $\wrap{\sd^2}$. The prior for $\sd^2$ is the
conjugate scaled inverse Chi-square distribution, i.e.,
$\nu \lambda \IC{\nu}$.  We recommend that the degrees of freedom,
$\nu$, be from 3 to 10 and the default is 3 which can be over-ridden
by the argument \code{sigdf}.  The $\lambda$ parameter can be specified by
the \code{lambda} argument which defaults to \code{NA}.  If
\code{lambda} is unspecified, then we determine a reasonable value for
$\lambda$ based on an estimate, $\widehat\sd$, (which can be specified by
the argument \code{sigest} and defaults to \code{NA}).  If
\code{sigest} is unspecified, the default value of \code{sigest} is
determined via linear regression or the sample standard deviation: if
$P<N$, then $y_i ~\N{\bm{x}_i'\widehat{\bm{\beta}}}{\widehat{\sd}^2}$;
otherwise, $\widehat{\sd}=s_y$.  Now we solve for $\lambda$ such that
$\P{\sd^2\le \widehat{\sd}^2}=q$.  This quantity, $q$, can be specified by
the argument \code{sigquant} and the default is 0.9 whereas we also
recommend considering 0.75 and 0.99.  Note that the pair $(\nu, q)$
are potential candidate parameters for choice via cross-validation.

Other important arguments for the BART prior.  We fix the number of
trees at $H$ which corresponds to the argument \code{ntree}.  The
default number of trees is 200 for continuous outcomes; as shown by
\citet{BleiKape14}, 50 is also a reasonable choice which is the
default for all other outcomes: cross-validation could be considered.
The number of cutpoints is provided by the argument \code{numcut} and
the default is 100.  The default number of cutpoints is achieved for
continuous covariates.  For continuous covariates, the cutpoints are
uniformly distributed by default, or generated via uniform quantiles
if the argument \code{usequants=TRUE} is provided.  By default,
discrete covariates which have fewer than 100 values will necessarily
have fewer cutpoints.  However, if you want a single discrete
covariate to be represented by a group of binary dummy variables, one
for each category, then pass the variable as a factor within a data
frame.

\section{Posterior computation for BART}\label{post}

In order to generate samples from the posterior for $f$, we sample the
structure of all the trees $\mathcal{T}_h$, for $h=1,\dots,H$; the
values of all leaves $\mu_{hn}$ for $n \in \mathcal{L}_h$
%$l=1,\dots,L_h$ 
within tree $h$; and, when appropriate, the error variance $\sd^2$.
Additionally, with the sparsity prior, there are samples of the vector
of splitting variable selection probabilities $[s_1,\dots,s_P]$ and,
when the sparsity parameter is random, samples of $\theta$.

The leaf and variance parameters are sampled from the posterior using
Gibbs sampling \citep{GemaGema84,GelfSmit90}.  Since the priors on these
parameters are conjugate, the Gibbs conditionals are specified
analytically. For the leaves, each $\mu_{hn}$ is drawn from a Normal
conditional density. The error variance, $\sd^2$, is drawn from a
scaled inverse Chi-square conditional.

Drawing a tree from the posterior requires a Metropolis-within-Gibbs
sampling scheme \citep{Muel91,Muel93}, i.e., a Metropolis-Hastings
(MH) step \citep{MetrRose53,Hast70} within Gibbs sampling.  For
single-tree models, four different proposal mechanisms are defined
\citep{ChipGeor98} (N.B.\ other MCMC tree sampling strategies have been
proposed: \citet{DeniMall98,WuTjel07,Prat16}).  The complementary
BIRTH/DEATH proposals are essential (the two other proposals
 are CHANGE and SWAP \citep{ChipGeor98}).
% However, because the branch regularization
% typically generates trees with few nodes, all four proposals are not
% necessary to explore the ensemble sample space. 
 For programming simplicity, the \pkg{BART} package only implements
 the BIRTH and DEATH proposals each with equal probability.  BIRTH
 selects a leaf and turns it into a branch, i.e., selects a new
 variable and cutpoint with two leaves ``born'' as its
 descendants. DEATH selects a branch leading to two terminal leaves
 and ``kills'' the branch by replacing it with a single
 leaf. % These two proposals are reversible, i.e., the MH ratio for the
% DEATH step is the reciprocal of the MH ratio for the BIRTH step and
% vice versa.
 To illustrate this discussion, we present the acceptance probability
 for a BIRTH proposal. Note that a DEATH proposal is the reversible
 inverse of a BIRTH proposal.

The algorithm assumes a fixed discrete set of possible split values
for each $x_j$. %(which is specified by the \code{numcut}
%argument: defaulting to 100).  
Furthermore, the leaf values, $\mu_{hn}$, are integrated over so
that our search in tree space is over a large, but discrete, set of
possibilities.  At the $m$th MCMC step, let $\mathcal{T}^m$ denote the
current state for the $h$th tree and $\mathcal{T}^*$ denotes the
proposed $h$th tree (subscript $h$ is suppressed for convenience).
$\mathcal{T}^*$ are identical $\mathcal{T}^m$ except that one terminal
leaf of $\mathcal{T}^m$ is replaced by a branch of $\mathcal{T}^*$
with two terminal leaves.  The proposed tree is accepted with the
following probability:
\begin{align*}
\pi_{\mathrm{BIRTH}}=\min\wrap[()]{1, 
\frac{\P{\mathcal{T}^*\,}}{\P{\mathcal{T}^m}}
\frac{\P{\mathcal{T}^m|\mathcal{T}^*\,}}
{\P{\mathcal{T}^*\,|\mathcal{T}^m}}}
\end{align*} where $\P{\mathcal{T}^m}$ and $\P{\mathcal{T}^*}$
are the posterior probabilities of ${T}^m$ and ${T}^*$
respectively. These are the targets of this sampling,
each consisting of a likelihood contribution and prior
contribution. Additionally, $\P{\mathcal{T}^m|\mathcal{T}^*}$
is the probability of proposing $\mathcal{T}^m$ given current
state $\mathcal{T}^*$ (a DEATH) and
$\P{\mathcal{T}^*\,|\mathcal{T}^m}$
is the probability of proposing $\mathcal{T}^*$ given current
state $\mathcal{T}^m$ (a BIRTH).

First, we describe the likelihood contribution to the posterior.  Let
$\bm{y}_n$ denote the partition of $\bm{y}$ corresponding to the leaf
node $n$ given the tree $\mathcal{T}$.  Because the leaf values are a
priori conditionally independent, we have
$\wrap{\bm{y}|\mathcal{T}}=\prod_n\wrap{\bm{y}_n|\mathcal{T}}$.  So,
for the ratio $\frac{\P{\mathcal{T}^*\,}}{\P{\mathcal{T}^m}}$ after
cancellation of terms in the numerator and denominator, we have
the likelihood contribution:
\begin{align*}
\frac{\P{\bm{y}_{\mathrm{L}},\bm{y}_{\mathrm{R}}|\mathcal{T}^*}}
{\P{\bm{y}_{\mathrm{LR}}|\mathcal{T}^m}}
&=\frac{\P{\bm{y}_{\mathrm{L}}|\mathcal{T}^*}
\P{\bm{y}_{\mathrm{R}}|\mathcal{T}^*}}
{\P{\bm{y}_{\mathrm{LR}}|\mathcal{T}^m}}
\end{align*}
where $\bm{y}_{\mathrm{L}}$ is the partition corresponding to the
newborn left leaf node; $\bm{y}_{\mathrm{R}}$, the partition for the
newborn right leaf node; and
$\bm{y}_{\mathrm{LR}}=\wrap{\bm{y}_{\mathrm{L}} \atop \bm{y}_{\mathrm{R}}}$.
N.B.\ the terms in the ratio are the predictive densities of a Normal
mean with a known variance and a Normal prior for the mean.

Similarly, the terms that the prior contributes to the posterior ratio
often cancel since there is only one ``place'' where the trees differ
and the prior draws components independently at different ``places''
of the tree.  Therefore, the prior contribution to
$\frac{\P{\mathcal{T}^*\,}}{\P{\mathcal{T}^m}}$ is 
\begin{align*}
\frac{\P{B_n=1}\P{B_l=0} \P{B_r=0} s_j} {\P{B_n=0}} & =
\frac{\alpha(t(n)+1)^{-\gamma}\wrap{1-\alpha(t(n)+2)^{-\gamma}}^2 s_j}
{1-\alpha(t(n)+1)^{-\gamma}}
\end{align*}
where $\P{B_n}$ is the branch regularity prior (see
Equation~\ref{regularity}), $s_j$ is the splitting variable selection
probability, $n$ is the chosen leaf node in tree $\mathcal{T}^m$,
$l=2n$ is the newborn left leaf node in tree $\mathcal{T}^*$ and
$r=2n+1$ is the newborn right leaf node in tree $\mathcal{T}^*$.

Finally, the ratio $ \frac{\P{\mathcal{T}^m|\mathcal{T}^*\,}}
{\P{\mathcal{T}^*\,|\mathcal{T}^m}}$ is 
\begin{align*}
\frac{\P{\mathrm{DEATH}|\mathcal{T}^*}
\P{n|\mathcal{T}^*}}
{\P{\mathrm{BIRTH}|\mathcal{T}^m}\P{n|\mathcal{T}^m} s_j}
\end{align*}
where $\P{n|\mathcal{T}} $ is the probability of choosing node
$n$ given tree $\mathcal{T}$.

N.B.\ $s_j$ appears in both the numerator and denominator
of the acceptance probability $\pi_{\mathrm{BIRTH}}$,
therefore, canceling which is mathematically convenient.
% Although not readily apparent, the acceptance probability is invariant
% to the splitting probability, $s_j$, for either the BIRTH or DEATH
% proposals.
% %regardless of whether a uniform or Dirichlet prior is specified. 
% For example, consider the BIRTH proposal. The denominator term
% $\P{\mathcal{T}^*|\mathcal{T}^m}$ contains $s_j$, but it also appears
% in the numerator term $\P{\mathcal{T}^*}$. Because of this
% cancellation, $\pi_{\mathrm{BIRTH}}$ does not depend on the variable
% splitting probability which is mathematically convenient.
% % since it is invariant to the prior chosen for $\bm{s}$.

Now, let's briefly discuss the posterior computation related to the
Dirichlet sparse prior.  If a Dirichlet prior is placed on the
variable splitting probabilities, $\bm{s}$, then its posterior samples
are drawn via Gibbs sampling with conjugate Dirichlet draws.  The
Dirichlet parameter is updated by adding the total variable
branch count over the ensemble, $m_j$, to the prior setting,
$\frac{\theta}{P}$, i.e.,
$\wrap{\frac{\theta}{P}+m_1, \dots, \frac{\theta}{P}+m_P}$.
% This is equivalent to any similar Multinomial-Dirichlet system where
% the components with the most observed counts get a larger share of
% probability from the Dirichlet.
In this way, the Dirichlet prior
induces a ``rich get richer'' variable selection strategy.
The sparsity parameter, $\theta$, is drawn on a fine grid of values
for the analytic posterior \citep{Line16}. This draw only depends on
$[s_1,\dots,s_P]$.

\section{Efficient computing with BART}\label{efficient}

If you had the task of creating an efficient implementation for a
black-box model such as BART, which tools would you use?
Surprisingly, linear algebra routines which are a traditional building
block of scientific computing will be of little use for a tree-based
method such as BART.  So what is needed?  Restricting ourselves to
widely available off-the-shelf hardware and open-source software,
we believe there are four key technologies necessary for a successful
BART implementation.
\begin{itemize}
\item an object-oriented language to facilitate working with trees and matrices
\item a parallel (or distributed) CPU computing framework for faster processing
\item a high-quality parallel random number generator
\item an interpreted shell for high-level data processing and analysis
\end{itemize}
In our implementation of BART, we pair the objected-oriented languages
of \proglang{R} and \proglang{C++} to satisfy these requirements.  In
this Section, we give a brief introduction to the concepts and
technologies harnessed for efficient computing by our \pkg{BART}
package.

\subsection{A brief history of multi-threading}

Writing multi-threaded programs is a fairly routine practice today
with a high-level language like \proglang{R} and corresponding
user-friendly interfaces such as the \pkg{parallel} \proglang{R} package
\citep{RCor18}.  Modern off-the-shelf laptops typically have 4 or
8 CPU cores placing reasonably priced multi-threaded hardware at your
fingertips.  Although, BART is often computationally undemanding, we
find it very convenient, with the aid of multi-threading, to run in
seconds that which would otherwise take minutes.  To highlight the
point that multi-threading is a mature technology, we now present a
brief history of multi-threading.  This is not meant to be exhaustive;
rather, we only provide enough detail to explain the capability and
popularity of multi-threading today.

Multi-threading
emerged rather early in the digital computer age with pioneers laying the
research groundwork in the 1960s.  In 1961, Burroughs released
the B5000 which was the first commercial hardware capable of
multi-threading \citep{Lync65}.  The B5000 performed asymmetric
multiprocessing which is commonly employed in modern hardware like
numerical co-processors and/or graphical processors today.  In 1962,
Burroughs released the D825 which was the first commercial hardware
capable of symmetric multiprocessing (SMP) with CPUs
\citep{AndeHoff62}.  In 1967, Gene Amdahl derived the theoretical
limits for multi-threading which came to be known as Amdahl's law
\citep{Amda67}.  If $B$ is the number of CPUs and $b$ is the fraction
of work that can't be parallelized, then the gain due to
multi-threading is $((1-b)/B+b)^{-1}$.

Now, fast-forward to the modern era of multi-threading.  Hardware and
software architectures in current use both directly, and indirectly,
led to the wide availability of pervasive multi-threading today.  In
2000, Advanced Micro Devices (AMD) released the AMD64 specification
that created a new 64-bit x86 instruction set which was capable of
co-existing with 16-bit and 32-bit x86 legacy instructions.  This was
an important advance since 64-bit math is capable of addressing vastly
more memory than 16-bit or 32-bit ($2^{64}$ vs.\ $2^{16}$ or $2^{32}$)
and multi-threading inherently requires more memory resources.  In
2003, version 2.6 of the Linux kernel incorporated full SMP support;
prior Linux kernels had either no support or very limited/crippled
support.  From 2005 to 2011, AMD released a series of Opteron chips
with multiple cores for multi-threading: 2 cores in 2005, 4 cores in
2007, 6 cores in 2009, 12 cores in 2010 and 16 cores in 2011.  From
2008 to 2010, Intel brought to market Xeon chips with their
hyper-threading technology that allows each core to issue two
instructions per clock cycle: 4 cores (8 threads) in 2008 and 8 cores
(16 threads) in 2010.  In today's era, most off-the-shelf hardware
available features 1 to 4 CPUs each of which is capable of
multi-threading.  Therefore, in the span of only a few
years, multi-threading rapidly trickled down from higher-end servers
to mass-market products such as desktops and laptops.  For example, the
consumer laptop that \pkg{BART} is developed on, purchased in 2016,
is capable of 8 threads (and hence many of the examples default to 8
threads).

\subsection{Modern multi-threading software frameworks}

% begin: addressing distributed computing major comment part 1
Up to this point, we have introduced multi-threading with respect to
parallelizing a task on a single system.  Here we want to make a
distinction between simple multi-threading on a single system and more
complex multi-threading on multiple systems simultaneously which is
often denoted by the term distributed computing.  On a single system,
various programming techniques can be used to create multi-threaded
software.  Basic multi-threading can be provided by the \code{fork}
system call which is often termed forking.  More advanced
multi-threading is provided by software frameworks such as
\pkg{OpenMP} and the Message Passing Interface (MPI).  Please note
that MPI can be employed for both simple multi-threading and for
distributed computing, e.g., MPI software initially written for a
single system could be extended to operate on multiple systems as
computational needs expand.  In the following, BART computations with
multi-threading are explored where the term multi-threading is used
for a single system and the term distributed computing is used for
multiple systems.
% end

In the late 1990s, MPI
%In the late 1990s, the Message Passing Interface (MPI)
\citep{WalkDong96} was introduced which is the dominant distributed
computing framework in use today \citep{GabrFagg04}.
%, i.e., distributed meaning tasks which span multiple computers.  
MPI support in \proglang{R} is built upon a fairly consistent
interface provided by the \pkg{parallel} package \citep{RCor18} which
is extended by other CRAN packages such as \pkg{snow}
\citep{TierRoss16} and \pkg{Rmpi} \citep{Yu17}.  To support MPI, new BART
software was created with a \proglang{C++} object schema
that is simple to program and maintain for distributed computing: we
call this the MPI BART code-base \citep{PratChip14}.  The \pkg{BART}
package source code is a descendent of MPI BART and its
programmer-friendly objects, although, the multi-threading MPI support
is now provided by \proglang{R} packages, e.g., \pkg{parallel},
\pkg{snow} and \pkg{Rmpi}.

The \pkg{BART} package supports multi-threading in two ways: 1)~via
\pkg{parallel} and related packages (which is how MPI is provided);
and 2)~via the OpenMP standard \citep{DaguMeno98}.  OpenMP takes
advantage of modern hardware by performing multi-threading on single
machines which often have multiple CPUs each with multiple cores.
Currently, the \pkg{BART} package only uses OpenMP for parallelizing
\code{predict} function calculations.  The challenge with OpenMP
(besides the \proglang{C}/\proglang{C++} programming required to
support it) is that it is not available on all platforms.
Operating system support can be detected by the GNU autotools
\citep{Calc10} which define a \proglang{C} pre-processor macro called
\code{_OPENMP} if it is available.  There are numerous exceptions for
operating systems so it is difficult to make universal statements.
But, generally, Microsoft Windows lacks OpenMP detection since the GNU
autotools do not natively exist on this platform.  For Apple macOS,
the standard Xcode toolkit does not provide OpenMP; however, the macOS
compilers on CRAN do provide OpenMP (see
\url{https://cran.r-project.org/bin/macosx/tools}).  Most Linux and
UNIX distributions provide OpenMP by default.  We provide the function
\code{mc.cores.openmp} which returns 1 if the \code{predict} function
is capable of utilizing OpenMP; otherwise, returns 0.

The \pkg{parallel} package provides multi-threading via forking.
Forking is available on Unix platforms, but not Windows (we use the
term Unix to refer to UNIX, Linux and macOS since they are all in the
UNIX family tree).  The \pkg{BART} package uses forking for posterior
sampling of the $f$ function, and also for the \code{predict} function
when OpenMP is not available.  Except for \code{predict}, all
functions that use forking start with \code{mc}.  And, regardless of
whether OpenMP or forking is employed, these functions accept the
argument \code{mc.cores} which controls the number of threads to be
used.  The \pkg{parallel} package provides the function
\code{detectCores} which returns the number of threads that your
hardware can support and, therefore, the \pkg{BART} package can use.

\subsection{BART implementations on CRAN}

Currently, there are four BART implementations on the Comprehensive R
Archive Network (CRAN); see the Appendix for a tabulated comparative
summary of their features.

\pkg{BayesTree} was the first released in 2006 \citep{ChipMcCu16}.
Reported bugs will be fixed, but no future improvements are planned;
so, we suggest choosing one of the newer packages such as \pkg{BART}.
The basic interface and work-flow of \pkg{BayesTree} has strongly
influenced the other packages which followed.  However, the
\pkg{BayesTree} source code is difficult to maintain and, therefore,
improvements were limited leaving it with relatively fewer features
than the other entries.

The second entrant is \pkg{bartMachine} which is written in
\proglang{Java} and was first released in 2013 \citep{KapeBlei16}.  It
provides advanced features like multi-threading, variable selection
\citep{BleiKape14}, a \code{predict} function, convergence diagnostics
and missing data handling.  However, the \proglang{R} to
\proglang{Java} interface can be challenging to deal with.
\proglang{R} is written in \proglang{C} and \proglang{Fortran},
consequentially, functions written in \proglang{Java} do not have a
natural interface to \proglang{R}.  This interface is provided by the
\pkg{rJava} \citep{Urba17} package which requires the Java Development
Kit (JDK).  Therefore, we highly recommend \pkg{bartMachine} for  
\proglang{Java} users. 

The third entrant is \pkg{dbarts} which is written in \proglang{C++}
and was first released in 2014 \citep{DoriChip16}.  It is a clone of
the \pkg{BayesTree} interface, but it does not share the source code;
\pkg{dbarts} source has been re-written from scratch for efficiency
and maintainability.  \pkg{dbarts} is a drop-in replacement for
\pkg{BayesTree}.  Although, it lacks multi-threading, the \pkg{dbarts}
serial implementation is the fastest, therefore, it is preferable when
multi-threading is unavailable such as on Windows.

The \pkg{BART} package which is written in \proglang{C++} was first
released in 2017 \citep{McCuSpar18}.  It provides advanced features
like multi-threading, variable selection \citep{Line16}, a
\code{predict} function and convergence diagnostics.  The source code
is a descendent of MPI BART.  Although, \proglang{R} is mainly written
in \proglang{C} and \proglang{Fortran} (at the time of this writing,
39.2\% and 26.8\% lines of source code respectively), \proglang{C++}
is a natural choice for creating \proglang{R} functions since they are
both object-oriented languages.  The \proglang{C++} interface to
\proglang{R} has been seamlessly provided by the \pkg{Rcpp} package
\citep{EddeFran11} which efficiently passes object references from
\proglang{R} to \proglang{C++} (and vice versa) as well as providing
direct access to the \proglang{R} random number generator.  The
source code can also be called from \proglang{C++} alone without an
\proglang{R} instance where the random number generation is provided
by either the standalone Rmath library \citep{RCor17} or the
\proglang{C++} \code{random} Standard Template Library.  Furthermore,
it is the only BART package to support categorical; and time-to-event
outcomes \citep{SparLoga16,SparRein18,SparLoga19}.  For one or more
missing covariates, record-level hot-decking imputation
\citep{deWaPann11} is employed that is biased towards the null, i.e.,
non-missing values from another record are randomly selected regardless
of the outcome.  This simple missing data imputation method is
sufficient for data sets with relatively few missing values; for more
advanced needs, we recommend the \pkg{sbart} package which utilizes
the Sequential BART algorithm \citep{DaniSing17,XuDani16} (N.B.\
\pkg{sbart} is also a descendent of MPI BART).

\begin{sidewaystable}
\begin{center}%\label{BARTcompare}
\begin{tabular}{l|llll}
Category                 & \pkg{BayesTree}      & \pkg{bartMachine}  & \pkg{dbarts}     & \pkg{BART}           \\ \hline
First release            & 2006                 & 2013               & 2014             & 2017                 \\
Authors                  & Chipman              & Kapelner           & Dorie,           & McCulloch, Sparapani \\
                         & \& McCulloch         & \& Bleich          & Chipman          & Gramacy, Spanbauer   \\
                         &                      &                    & \& McCulloch     & \& Pratola           \\
Source code              & \proglang{C++}       & \proglang{Java}    & \proglang{C++}   & \proglang{C++}       \\
% R package dependencies   & \pkg{nnet}        & \pkg{rJava},
%                                                   \pkg{car},      &
%                                                                                         & \pkg{Rcpp}, \pkg{nlme},
%                                                                                           \pkg{nnet},       \\
%                          &                      & \pkg{randomForest},&               & \pkg{parallel}, 
%                                                                                           \pkg{survival},   \\
%                          &                      & \pkg{missForest}&                  & \pkg{tools}       \\
R package dependencies   & None                 & \pkg{rJava}, \pkg{car},& None         & \pkg{Rcpp}           \\
excluding Recommended    &                      & \pkg{randomForest},&                  &                      \\
                         &                      & \pkg{missForest}   &                  &                      \\
Tree transition proposals& 4                    & 3                  & 4                & 2                    \\
Multi-threaded           & No                   & Yes                & No               & Yes                  \\
\code{predict} function  & No                   & Yes                & No               & Yes                  \\
Variable selection       & No                   & Yes                & No               & Yes                  \\
Continuous outcomes      & Yes                  & Yes                & Yes              & Yes                  \\
Binary outcomes probit      
                         & Yes                  & Yes                & Yes              & Yes                  \\
Binary outcomes logit    
                         & No                   & No                 & No               & Yes                  \\
Categorical outcomes     & No                   & No                 & No               & Yes                  \\
Time-to-event outcomes   & No                   & No                 & No               & Yes                  \\
Convergence diagnostics  & No                   & Yes                & No               & Yes                  \\
Thinning                 & Yes                  & No                 & Yes              & Yes                  \\
Missing data handling    & No                   & Yes                & No               & Yes                  \\
Cross-validation         & No                   & Yes                & Yes              & No                   \\
Partial dependence plots & Yes                  & Yes                & Yes              & No                   \\ \hline
% Citations                & \multicolumn{4}{l}{\cite{ChipMcCu16}}                                               \\ 
%                          &                      & \multicolumn{3}{l}{\cite{KapeBlei16}}                        \\ 
%                          &                      &                    & \multicolumn{2}{l}{\cite{DoriChip16}}   \\ 
%                          &                      &                    &                  & \cite{McCuSpar18}    \\ \hline
\end{tabular}\\
\end{center}
\begin{quote}
\cite{ChipMcCu16}\\
\cite{KapeBlei16}\\
\cite{DoriChip16}\\
\cite{McCuSpar18}
\end{quote}
\end{sidewaystable}

\subsection{MCMC is embarrassingly parallel}

In general, Bayesian Markov chain Monte Carlo (MCMC) posterior
sampling is considered to be embarrassingly parallel
\citep{RossTier07}, i.e., since the chains only share the data and
don't have to communicate with each other, parallel implementations
are considered to be trivial.  BART MCMC also falls into this class.

% begin: addressing distributed computing major comment part 2
However, to clarify this point before proceeding, the embarrassingly
parallel designation is in the context of simple multi-threading on
single systems.  An adaptation of distributed computing to large data
sets exhaustively divides the data into mutually exclusive partitions,
called shards, such that each system only processes a single shard.
With sharded distributed computing, the embarrassingly parallel moniker
does not apply.  Recently, two advanced techniques have been developed
for BART computations with sharding: Monte Carlo consensus
\citep{PratChip14} and modified likelihood inflating sampling
algorithm, or modified LISA, \citep{EnteCrai18}.  From here on, simple
multi-threading is assumed.
%end

Typical practice for Bayesian MCMC is to start in some initial state,
perform a limited number of samples to generate a new random starting
position and throw away the preceding samples which we call burn-in.
The amount of burn-in in the \pkg{BART} package is controlled by the
argument \code{nskip}: defaults to 100 with the exception of
time-to-event outcomes which default to 250.  The total
length of the chain returned is controlled by the argument
\code{ndpost} which defaults to 1000.  The theoretical gain due to
multi-threading can be calculated by what we call the MCMC Corollary
to Amdahl's Law.  Let $b$ be the burn-in fraction and $B$ be the
number of threads, then the gain limit is $((1-b)/B+b)^{-1}$.  (As an
aside, note that we can derive Amdahl's Law as follows where the
amount of work done is in the numerator and elapsed time is in the
denominator: $\frac{1-b+b}{(1-b)/B+b}=\frac{1}{(1-b)/B+b}$).  For
example, see the diagram in Figure~\ref{MCMC} where the burn-in
fraction, $b=\frac{100}{1100}=0.09$, and the number of CPUs, $B=5$,
results in an elapsed time of only $((1-b)/B+b)=0.27$ or a
$((1-b)/B+b)^{-1}=3.67$ fold reduction which is the gain in
efficiency.  In Figure~\ref{Amdahl}, we plot theoretical gains
on the y-axis and the number of CPUs on the x-axis
for two settings: $b \in \{0.025, 0.1\}$.
\begin{figure}%[t!]
\begin{center}
\includegraphics[scale=0.48]{figures/parallel.pdf}
\end{center}
\caption{\label{MCMC}The theoretical gain due to multi-threading can
  be calculated by Amdahl's Law.  Let $b$ be the burn-in fraction and
  $B$ be the number of threads, then the theoretical gain limit is
  $((1-b)/B+b)^{-1}$.  In this diagram, the burn-in fraction,
  $b=\frac{100}{1100}=0.09$, and the number of CPUs, $B=5$, results in
  an elapsed time of only $((1-b)/B+b)=0.27$ or a
  $((1-b)/B+b)^{-1}=3.67$ fold reduction which is the gain in
  efficiency.  }
\end{figure}

\begin{figure}%[t!]
\begin{center}
\includegraphics[scale=0.42]{figures/amdahl.pdf}
\end{center}
\caption{\label{Amdahl}The theoretical gain due to multi-threading can
  be calculated by Amdahl's Law.  Let $b$ be the burn-in fraction and
  $B$ be the number of threads, then the theoretical gain limit is
  $((1-b)/B+b)^{-1}$.  In this figure, the theoretical gains are on
  the y-axis and the number of CPUs, the x-axis, for two settings:
  $b \in \{0.025, 0.1\}$.  }
\end{figure}

\subsection{Multi-threading and random access memory}

The IEEE standard 754-2008 \citep{IEEE08} specifies that every
double-precision number consumes 8 bytes (64 bits).  Therefore, it is
quite simple to estimate the amount of random access memory (RAM)
required to store a matrix.  If $A$ is $m \times n$, then the amount
of RAM needed is $8 \times m \times n$ bytes.  Large matrices held in
RAM can present a challenge to system performance.  If you consume all
of the physical RAM, the system will ``swap'' segments out to virtual
RAM which are disk files and this can degrade performance and possibly
even crash the system.  On Unix, you can monitor memory and swap usage
with the \code{top} command-line utility.  And, within \proglang{R},
you can determine the size of an object with the \code{object.size}
function.

Mathematically, a matrix is represented as follows.
\begin{align*}
A=\wrap{\begin{array}{cccc}
a_{11} & a_{12} & \cdots & a_{1n} \\
a_{21} & a_{22} & \cdots & a_{2n} \\
\vdots& \vdots & \ddots & \vdots \\
a_{m1} & a_{m2} & \cdots & a_{mn} \\
\end{array}}
\end{align*}
\proglang{R} is a column-major language, i.e., matrices are laid out
in consecutive memory locations by traversing the columns:
$\wrap{a_{11}, a_{21}, \., a_{12}, a_{22}, \.}$.  \proglang{R} is
written in \proglang{C} and \proglang{Fortran} where
\proglang{Fortran} is a column-major language as well.  However,
\proglang{C} and \proglang{C++} are row-major languages, i.e.,
matrices are laid out in consecutive memory locations by traversing
the rows: $\wrap{a_{11}, a_{12}, \., a_{21}, a_{22}, \.}$.  So, if you
have written an \proglang{R} function in \proglang{C}/\proglang{C++},
then you need to be cognizant of the clash in paradigms (also note
that \proglang{R/Fortran} array indexing goes from 1 to $m$ while
\proglang{C}/\proglang{C++} indexing goes from 0 to $m-1$).  As you
might surmise, this is easily addressed with a transpose, i.e.,
instead of passing $A$ from \proglang{R} to
\proglang{C}/\proglang{C++} pass $A^{\top}$.

\proglang{R} is very efficient in passing objects; rather, than
passing an object (along with all of its memory consumption) on the
stack, it passes objects merely by a pointer referencing the original
memory location.  However, \proglang{R} follows copy-on-write memory
allocation, i.e., all objects present in the parent thread can be read
by a child thread without a copy, but when an object is
altered/written by the child, then a new copy is created in memory.
Therefore, if we pass $A$ from \proglang{R} to
\proglang{C}/\proglang{C++}, and then transpose, we will create
multiple copies of $A$ consuming $8 \times m \times n \times B$ where
$B$ is the number of children.  If $A$ is a large matrix, then you may
stress the system's limits.  The simple solution is for the parent to
create the transpose before passing $A$ and avoiding the multiple
copies, i.e., \code{A <- t(A)}.  And this is the philosophy that the
\pkg{BART} package follows for the multi-threaded BART functions; see
the documentation for the \code{transposed} argument.

\subsection{Multi-threading: interactive and batch processing}

Interactive jobs must take precedence over batch jobs to prevent the
user experience from suffering high latency.  For example, have you
ever experienced a system slowdown while you are typing and the
display of your keystrokes can not keep up; this should never happen
and is the sign of something amiss.  With large multi-threaded jobs,
it is surprisingly easy to naively degrade system performance.  But,
this can easily be avoided by operating system support provided by
\proglang{R}.  In the \pkg{tools} package \citep{RCor18}, there is
the \code{psnice} function.  Paraphrased from the \code{?psnice} help
page.
\begin{quote}
  Unix has a concept of process priority.  Priority is assigned values from 0
  to 39 with 20 being the normal priority and (counter-intuitively)
  larger numeric values denoting lower priority.  Adding to the
  complexity, there is a ``nice'' value, the amount by which the
  priority exceeds 20.  Processes with higher nice values
  will receive less CPU time than those with normal priority.
  Generally, processes with nice 19 are only run when the system would
  otherwise be idle.
\end{quote}
Therefore, by default, the \pkg{BART} package children have their nice
value set to 19.

% \begin{comment}
% \subsection{Continuous BART: serial and parallel implementations}

% Here we present snippets of \proglang{R} code to run BART in serial
% and parallel for a continuous outcome which we call continuous BART.
% While we only demonstrate continuous outcomes, the other outcomes
% are as similarly handled by the \pkg{BART} package as possible to
% present a consistent interface.
% The serial function is \code{wbart} and the parallel, \code{mc.wbart}.
% The 'w' in the name stands for weighted since you can provide known
% weights (with the \code{w} argument) for the following model:
% $y_i ~ \N{f(\bm{x}_i}{w_i^2\sd^2} \where (f, \sd) \prior \mathrm{BART}$
% and $i=1, \., N$ indexes subjects.  Now,
% we can perform the calculations in serial,\\ 
% \code{set.seed(99); post <- wbart(x.train, y.train, ndpost=M)}\\ 
% % keepevery=1)}\\
%  or in parallel  (when said support is available),\\
% \code{post <- {mc.wbart}(x.train, y.train, ndpost=M, %keepevery=1, 
% mc.cores=8, seed=99)}.\\
%   Notice the difference in how the seed is set;
% we will return to this detail later on.  The \pkg{BART} package
% allows \code{x.train} (and \code{x.test}) to be provided as matrices
% or data frames, but for simplicity we present them as matrices.
% \begin{align*}
% \mbox{Input: \code{x.train} and, optionally, \code{x.test}:\ } &
% \wrap{\begin{array}{c}
% \bm{x}_{1} \\
% \bm{x}_{2} \\
% \vdots \\
% \bm{x}_{N} \\
% \end{array}} \where \bm{x}_{i} \mbox{\ is the $i^{th}$ row} \\
% \mbox{Output: \code{post\$yhat.train} and \code{post\$yhat.test}:\ } &
% \wrap{\begin{array}{ccc}
% \hat{y}_{11}& \dots & \hat{y}_{N1} \\
% \vdots & \ddots & \vdots \\
% \hat{y}_{1M}& \dots & \hat{y}_{NM} \\
% \end{array}} \where \hat{y}_{im}=\mu_0+f_m(\bm{x}_i) 
% \end{align*}
% The \code{post} object returned is of type \code{wbart} which is
% essentially a list.  There are other items returned in the list, but
% here we only focus on \code{post\$yhat.train} and
% \code{post\$yhat.test}; the latter only being returned if
% \code{x.test} is provided.  In the above display, $m=1, \., M$ are the
% MCMC samples which are the rows of \code{post\$yhat.train} and
% \code{post\$yhat.test}.  Note that each outcome has a different return
% type, i.e., \code{post} object of type \code{wbart} (continuous),
% \code{pbart} (binary probit), \code{lbart} (binary logit),
% \code{mbart} (multinomial with either probit or logit),
% \code{survbart} (survival analysis), \code{criskbart} (competing
% risks) or \code{recurbart} (recurrent events).

% \subsection{Continuous BART: using predict with a previous fit}

% Often when we are fitting a BART model, we have not specified an
% \code{x.test} matrix of hypothetical values for the evaluation of $f$.
% For fire-and-forget packages like \pkg{BayesTree} and \pkg{dbarts}, we
% would have to re-fit the model every time we want to evaluate
% \code{x.test} which can be very time-consuming.  Therefore, the
% \pkg{BART} package takes a unique approach: it returns the ensemble of
% trees in the \code{post} object for later use; specifically, they are
% encoded in an ASCII character string, \code{post\$treedraws\$trees}.
% This allows us to construct \code{x.test} after the fact which is
% often convenient when it is large since we can partition it into
% smaller chunks.  Then we can evaluate predictions via the S3 method
% \code{predict.wbart}.  The
% predictions are generated in serial by default,\\
% \code{pred <- predict({post}, {x.test})} \\
% but can be parallelized (when said support is available),\\
% \code{pred <- predict({post}, {x.test}, {mc.cores=B})}.\\
% Notation: $M$~for the number of posterior samples, $B$~for the number
% of threads (generally, $B=1$ for Windows), $N$~for the number of
% observations in the training set, and $Q$~for the number of
% observations in the test set.
% \begin{align*}
% \mbox{Input: \code{x.test}:\ } &
% \wrap{\begin{array}{c}
% x_{1} \\
% x_{2} \\
% \vdots \\
% x_{Q} \\
% \end{array}} \where \bm{x}_{i} \mbox{\ is the $i^{th}$ row and\ } i=1, \., Q \\
% \mbox{Output matrix:\ } &
% \wrap{\begin{array}{ccc}
% \hat{y}_{11}& \dots & \hat{y}_{Q1} \\
% \vdots & \ddots & \vdots \\
% \hat{y}_{1M}& \dots & \hat{y}_{QM} \\
% \end{array}} \where \hat{y}_{im}=f_m(\bm{x}_i) 
% \end{align*}
% In the above display, $m=1, \., M$ are the
% MCMC samples which are the rows of the output matrix.
% \end{comment}

  
\subsection{Creating a BART executable}

Occasionally, you may need to create a BART executable that you can
run without an \proglang{R} instance.  This is especially useful if
you need to include BART in another \proglang{C++} program.  Or, when
you need to debug the \pkg{BART} package \proglang{C++} source code
which is more difficult to do when you are calling the function from
\proglang{R}.  Several examples of these are provided with the
\pkg{BART} package.  With \proglang{R}, you can find the
\code{Makefile} and the weighted BART example with
\code{system.file("cxx-ex/Makefile", package="BART")} and\\
\code{system.file("cxx-ex/wmain.cpp", package="BART")} respectively.
Note that these examples require the installation of the standalone
Rmath library \citep{RCor17} which is contained in the \proglang{R}
source code distribution.  Rmath provides common \proglang{R}
functions and random number generation, e.g., \code{pnorm} and
\code{rnorm}.  You will likely need to copy the \code{cxx-ex}
directory to your workspace.  Once done, you can build and run the
weighted BART executable example from the command line shell as
follows.\\
\code{sh\% make wmain.out \#\# to build\\sh\% ./wmain.out \#\# to run}\\
By default, these examples are based on the Rmath random number
generator.  However, you can specify the \proglang{C++} Standard
Template Library random number generator (contained in the STL
\code{random} header file) by uncommenting the following line in
the \code{Makefile} (by removing the pound, \#, symbols):\\
\code{\#\# CPPFLAGS = -I. -I/usr/local/include -DMATHLIB_STANDALONE -DRNG_random}\\
(which still requires Rmath for other purposes).  These examples were
developed on Linux and macOS, but they should be readily adaptable to
UNIX and Windows as well.

\pagebreak

\bibliography{references}

\end{document}

