\documentclass[nojss]{jss}
  %\VignetteEngine{knitr::knitr}
  %\VignetteIndexEntry{Continuous outcomes with BART: Part 1}
%\usepackage{thumbpdf,lmodern}
\usepackage{verbatim}
\usepackage{statex2}
\usepackage[authoryear,round]{natbib}
\usepackage{rotating}

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

%% - \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{Continuous outcomes with BART: Part 1}
\Plaintitle{Continuous outcomes with BART: Part 1}
\Shorttitle{Continuous outcomes with BART: Part 1}

%% - \Abstract{} almost as usual
\Abstract{
  This short article delves into the BART prior with respect to 
  continuous outcomes and the {\bf BART} \proglang{R} package
  implementation.
}

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

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

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

\begin{document}

\maketitle 

\section{Introduction to Bayesian Additive Regression Trees}

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
nice 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{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
yielding desirable properties.  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.

\section{Continuous outcomes with BART}

In this section, we document the analysis of continuous outcomes with
the {\bf BART} \proglang{R} package.  This requires delving into the
details of the BART prior itself and the corresponding arguments to
the \code{wbart} function.

\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)
$$

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 {\bf BART} \proglang{R} package.  But first, we
delve into the details of the BART prior itself.
\end{comment}

\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.  Let the function $g(\bm{x}; T, M)$ assign a value based on
the input $\bm{x}$; $T$ representing the structure of a binary tree,
including interior decision rules as branches and terminal nodes as
leaves; and $M=\{\mu_1, \dots, \mu_L\}$, the parameter values
associated with the $L$ leaves one of which will be the output.  $T$
is composed of decision rules of the form $x_j\le c$ which means
branch left and $x_j>c$, branch right; or terminal leaves where it
stops.  The function, $f(\bm{x})$, is a sum of $H$ trees, i.e.,
$f(\bm{x})=\sum_{h=1}^H g(\bm{x}; T_h, M_h)$ 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$:
$y_i=\mu_0+f(\bm{x}_i)+\epsilon_i \where \epsilon_i \iid \N{0}{ w_i^2 \sd^2}$; the
unknown random function, $f$, and the error variance, $\sd^2$, follow
the BART prior expressed notationally as $(f,\sd^2)\prior\mathrm{BART}$
and $\mu_0$ is a known constant which centers ${y}$.
The $w_i$ are known standard deviation weights which you can 
supply with the argument \code{w} that is only available for 
continuous outcomes; the unit 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, but it is best understood as a
bundle of priors.  Using the Gelfand-Smith shorthand bracket notation
for the abstract specification of random variable distributions \citep{GelfSmit90},
we represent the BART prior as $\wrap{T_h}\wrap{\mu_{hl}|T_h}\wrap{\sigma^2}$
where $\wrap{T_h}$ is the tree prior, $\wrap{\mu_{hl}|T_h}$ is the leaf
prior and $\wrap{\sigma^2}$ is the error variance prior.

The tree prior: $\wrap{T_h}$.  There are three priors for $T_h$ which govern
whether the tree branches grow or are pruned.  The first is a
regularity prior dictating the probability of an interior branch,
event $B$, given the branch tier, $t$, is:
$\P{B|t}=\mbox{base} (1+t)^{-\mbox{power}}$.  You can specify these
prior parameters with arguments, but the following defaults are highly
recommended: \code{base=0.95} and \code{power=2}; for a detailed
discussion of these parameters, 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 conditional on a
branch event $B$ which defaults to uniform probability $1/P$ 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} which we will return to below.  Given a branch
event, $B$, and a variable chosen, $x_j$, the last tree prior selects
a cut point, $c$, within the range of observed values for $x_j$; this
prior is uniform.

We represent the probability of variable selection via the sparse
Dirichlet prior as\\
$\wrap{s_1, \dots, s_P} \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 $1/P$.  The prior parameter $\theta$
can be fixed or random; the default \code{alpha=0} means random, but
fixed can be specified by supplying a positive number to the argument
\code{alpha}.  Random $\theta$ is specified via
$\frac{\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 $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}); \code{a=0.5} induces a sparse
posture and \code{a=1}, non-sparse or dense.

If covariates are supplied via a data frame (rather than a matrix)
with the argument \code{x.train}, then factors with more than two
levels can be specified along with continuous variables.  The factors
are transformed into dummy variables with each factor in their own
group while each continuous variable is a group all by itself
(similarly, factors with two levels are a group with a single dummy
variable).  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}$.
% 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_{hl}|T_h}$.  Given a tree, $T_h$, there is a prior
on its leaf values, $\mu_{hl}|T_h$ and we denote the collection of all
leaves in $T_h$ by $M_h=\wrap{\mu_{h1}, \dots, \mu_{hL_h}}$.  Suppose
that $y_i \in [y_{\min}, y_{\max}]$ for all $i$ and
$\wrap{\mu_{1(i)}, \dots, \mu_{H(i)}}$ are the leaf output values from
each tree corresponding to the vector of covariates, $\bm{x}_i$.  If
$\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 values for $\mu_{\mu}$ and $\sd_{\mu}$ which are solutions to
the system of equations created by the $1-\alpha/2$ confidence
interval: $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 |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_{hl} \prior \N{0}{\wrap{\frac{\tau}{2k\sqrt{H}}}^2} \where
\tau={y_{\max}-y_{\min}}$.  So, the prior for $\mu_{hl}$ 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 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
good 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 specifed 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, $\hat\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'\hat{\bm{\beta}}}{\hat{\sd}^2}$;
otherwise, $\hat{\sd}=s_y$.  Now we solve for $\lambda$ such that
$\P{\sd^2\le \hat{\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 good 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; as shown by \citet{BleiKape14}, 50 is
also a reasonable choice and 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; however, discrete covariates which have
fewer than 100 values will necessarily have fewer cutpoints.  For
continuous covariates, the cutpoints are uniformly distributed by
default, or generated via uniform quantiles if the argument
\code{usequants=TRUE} is provided.

\subsection{Posterior samples returned}

The number of MCMC samples discarded for
burnin 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.  You can also 
thin some returned values, but not others.  The following arguments
default to \code{ndpost}, but can be over-ridden as needed.
\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 ensembles to return
\end{itemize}

\subsection{Typical use case}

Typically, when calling the \code{wbart} function to analyze some
data, 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}) function
is for serial (parallel) computation.  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. \\
\code{set.seed(99)}\\
\code{post <- {wbart}(x.train, y.train, x.test, ..., ndpost=M)} or \\
\code{post <- {mc.pbart}(x.train, y.train, x.test, ..., ndpost=M, {mc.cores=2, seed=99})}  \\
\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{\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}& \. & \hat{y}_{N1} \\
\vdots & \dots & \vdots \\
\hat{y}_{1M}& \. & \hat{y}_{NM} \\
\end{array}} \where  \hat{y}_{im}=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.

Often it is impractical to provide \code{x.test} in the call to
\code{wbart} 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 {\bf 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{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.\\
\code{pred <- {predict}(post, x.test, {mc.cores=1}, ...)} \\
\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{\code{pred} is a matrix:} & 
%\mbox{\code{pred}:\ } &
\wrap{\begin{array}{ccc}
\hat{y}_{11}& \. & \hat{y}_{Q1} \\
\vdots & \vdots & \vdots \\
\hat{y}_{1M}& \. & \hat{y}_{QM} \\
\end{array}} \where \hat{y}_{hm}=f_m(\bm{x}_h)  \\
\end{align*}

\bibliography{references}

\end{document}

