\documentclass[nojss]{jss}
%\VignetteEngine{knitr::knitr}
%\VignetteIndexEntry{Binary and categorical outcomes with BART}
%\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{Rodney Sparapani\\Medical College of Wisconsin
\And Robert McCulloch\\Arizona State University}
\Plainauthor{Rodney Sparapani, Robert McCulloch}

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

\title{Binary and categorical outcomes with BART}
\Plaintitle{Binary and categorical outcomes with BART}
\Shorttitle{Binary/categorical outcomes with BART}

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

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

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

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

\begin{document}

\maketitle

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

\section{Binary and categorical outcomes with BART}

The {\bf BART} package supports binary outcomes via probit BART with
Normal latents and logistic BART with Logistic latents.  Categorical
outcomes are supported with multinomial BART with Logistic latents.
Convergence diagnostics are provided and variable selection as well.

\subsection{Probit BART for binary outcomes}

To extend BART to binary outcomes, we employ the technique of
\cite{AlbeChib93} to create what we call probit BART.  Probit BART is
provided by the {\bf BART} package as the \code{pbart} function.  In
this case, the outcome, \code{y.train}, is provided as an integer with
values 0 or 1.  Given $y_i$, we introduce 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 where $i$ indexes subject and $\Phi$ is the
standard Normal cumulative distribution function.
\begin{align*}
y_{i}|p_{i} & \ind \B{p_{i}} \\
p_{i}|f & =  \Phi({\mu}_0+f(\bm{x}_i)) \where f \prior \mathrm{BART} \\
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*}
The $z_i$ are centered around a known constant, $\mu_0$, which is
tantamount to centering the probabilities, $p_i$, around
$p_0=\Phi(\mu_0)$.  If $\mu_0=0$, which is the default, then the $p_i$
are centered around 0.5; to specify a different value, say -1, pass
the argument \code{binaryOffset=-1} 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.  This setup leads to the following Bernoulli likelihood:
$\wrap{\bm{y}|f} = \prod_{i=1}^N p_{i}^{y_i}(1-p_{i})^{1-y_i}$. 

In the following, we assume that \code{binaryOffset=0} for convenience
(which is the default).  The \code{pbart} (\code{mc.pbart}) function
is 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. \\
\code{set.seed(99)}\\
\code{post <- {pbart}(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{pbart}, 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.
Although, \code{post\$yhat.train} and \code{post\$yhat.test}, when requested,
are returned, generally, \code{post\$prob.train} and \code{post\$prob.test}
are 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}& \. & \hat{p}_{N1} \\
\vdots & \dots & \vdots \\
\hat{p}_{1M}& \. & \hat{p}_{NM} \\
\end{array}} \where  \hat{p}_{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 {\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{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.\\
\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}, of type \code{pbart}, which is essentially a list} & \\
\mbox{\code{pred\$yhat.test}:\ } &
\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)  \\
\mbox{\code{pred\$prob.test}:\ } &
\wrap{\begin{array}{ccc}
\hat{p}_{11}& \. & \hat{p}_{Q1} \\
\vdots & \vdots & \vdots \\
\hat{p}_{1M}& \. & \hat{p}_{QM} \\
\end{array}} \where \hat{p}_{hm}=\Phi(f_m(\bm{x}_h))  \\
\mbox{\code{pred\$prob.test.mean}:\ } &
\wrap{\hat{p}_{1}, \., \hat{p}_{Q}} \where \hat{p}_{h}=M^{-1} \sum_{m=1}^M  \hat{p}_{hm}   \\
\end{align*}

\subsection{Friedman's partial dependence function}

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 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})$.  For
probit BART, the $f$ function is not directly of interest; rather, the
probability of an event is more interpretable:
$p(\bm{x}_S) = {N^{-1}} \sum_{i=1}^N
\Phi(\mu_0+f(\bm{x}_S,\bm{x}_{iC}))$.
Other marginal functions can be obtained in a similar fashion.
Estimates can be derived via functions of the posterior samples such
as means, quantiles, e.g.,
$\hat{p}(\bm{x}_S) = {M^{-1}} {N^{-1}} \sum_{m=1}^M \sum_{i=1}^N
\Phi(\mu_0+f_m(\bm{x}_S,\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 {\bf BART}
package; rather, we provide examples of the many practical uses in the
\code{demo} directory.

\subsection{Logistic BART for binary outcomes}

Note that the distribution of the latent $z_i$ is not identifiable
from the data so it is essentially a parametric assumption.  This
assumption 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 the technique of \cite{HolmHeld06} to
create what we call logistic BART.  However, it is important to
recognize that logistic BART is more computationally intensive 
than probit BART.
  
The outcome, \code{y.train}, is provided as an integer with values 0
or 1.  Logistic BART is provided by the \code{lbart} function.  Unlike
probit BART where the auxiliary latents, $z_i$, have a fixed variance
$\sd^2=1$; with Logisitic BART, we sample truncated Normal latents,
$z_i$, with a random variance $\sd_i^2$ \citep{Robe95}.  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 known 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 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.

\subsection{Multinomial BART for categorical outcomes}

To extend BART to categorical outcomes, we employ as many logistic
BARTs as there are categories to create what we call multinomial
BART.  Multinomial BART is provided by the {\bf BART} package as the
\code{mbart} function.  In this case, the outcome, \code{y.train}, is
provided as an integer with values $1, \dots, C$ which generate the
corresponding latents $z_{i1}, \dots, z_{iC}$ and probabilities 
$p_{i1}, \dots, p_{iC}$ which are constrained to sum to one.

The input for \code{mbart} is essentially identical to \code{pbart}.
The output is slightly different.
%Also, the \code{predict} function for objects of type \code{lbart} is analogous.
\code{set.seed(99)}\\
\code{post <- {mbart}(x.train, y.train, x.test, ..., ndpost=M)} or \\
\code{post <- {mc.mbart}(x.train, y.train, x.test, ..., ndpost=M, {mc.cores=2, seed=99})}  \\
\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}} \mbox{\ or\ } \bm{x}_{i} \\
\mbox{\code{post}, of type \code{mbart}} & \\
\mbox{\code{post\$prob.train} and \code{post\$prob.test}:\ } &
\wrap{\begin{array}{ccccccc}
\hat{p}_{111} & \. & \hat{p}_{1C1} & \. & \hat{p}_{N11} & \. & \hat{p}_{NC1}\\
\vdots       & \vdots & \vdots    & \vdots & \vdots   & \vdots & \vdots \\
\hat{p}_{11M} & \. & \hat{p}_{1CM} & \. & \hat{p}_{N1M} & \. & \hat{p}_{NCM} \\
\end{array}} \\
& \where  \hat{p}_{icm}=F(f_{cm}(\bm{x}_i)) \\
\end{align*}
The columns of \code{post\$prob.train} and \code{post\$prob.test}
represent different covariate settings crossed with the $C$
categories.  Also, the \code{predict} function for objects of type
\code{mbart} is analogous.

\subsection{Converegence diagnostics for dichotomous and categorical outcomes}

How do you perform convergence diagnostics for BART?  For continuous
outcomes, convegence can easily be determined from the trace plots of
the the error variance, $\sigma^2$.  However, for probit BART with
Normal latents, the error variance is fixed at 1 so this is not an
option.  Similarly, for logistic and multinomial BART, $\sigma_i^2$,
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 {\bf BART} package, we supply R functions adapted from the {\bf
  coda} 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-100},
\ref{geweke-pbart2-1000} and \ref{geweke-pbart2-10000} generated by the
example \code{geweke.pbart2.R}:\\  
\code{system.file('demo/geweke.pbart2.R', 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 = 100, 1000, 10000$.
\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*}

The convergence for each of these data sets is graphically displayed
in Figures~\ref{geweke-pbart2-100}, \ref{geweke-pbart2-1000} and
\ref{geweke-pbart2-10000} 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=100, 1000$, but a more
notable amount for $N=10000$.  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=100, 1000$, but less notably for
$N=10000$.  In the lower right quadrant, we plot the Geweke $Z_{AB}$
statistics for each subject $i$.  Notice that for $N=100$, 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=10000$, there are dramatically more values exceeding the 95\%
limits.  Based on these figures, we conclude that the chains have
converged for $N=100$; for $N=1000$, convergence is questionable; and,
for $N=10000$, convergence has not been attained.  We would suggest
that more thinning be employed for $N=1000, 10000$ via the
\code{keepevery} argument to \code{pbart}; perhaps,
\code{keepevery=50} for $N=1000$ and \code{keepevery=250} for
$N=10000$.

\subsection{BART and variable selection}

Several methods have been proposed for variable selection with BART
\citep{ChipGeor10,BleiKape14,HahnCarv15,McCuCarv15,Line16}.  The {\bf
  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} & \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 or dense; you can specify this
parameter with argument \code{a}.  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 (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 {\bf BART} package:
\code{system.file('demo/sparse.pbart.R', package='BART')}.
For sample sizes of $N=100, 1000, 10000$, 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=100$, only $s_2\le P^{-1}$, but
notice that there are 34 false positives.  For $N=1000$, all five active
covariates are identified, but notice that there are 18 false positives.
For $N=10000$, all five active covariates are identified and notice
that there are only two false positives.

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 $B_h$
is the number of branches in tree ${T}_h$, then the concordant
pair probability is:
$\kappa_{ij} = \P{x_i \leftrightarrow x_j \in {T}_h|B_h>1}$
for $i=1, \dots, P-1$ and $j=i+1, \dots, P$.  See an example of 
calculating these probabilities in \code{system.file('demo/trees.pbart.R',
package='BART')}.

\subsection{Motivating example: chronic pain and obesity}

We want to test 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 sampling weights are available to
extrapolate the rates from the survey to the US as a whole.  We will
concentrate on body mass index (BMI) and gender, $\bm{x}_S$, while
utilizing Friedman's partial dependence function as defined above and 
also incorporating the sampling weights, i.e., 
$p_S(\bm{x}_S) = {\sum_{i=1}^N w_i \Phi(\mu_0+f(\bm{x}_S,\bm{x}_{iC}))}/
{\sum_{i'=1}^N w_{i'}}$.

The {\bf BART} package provides two examples:\\
\code{system.file('demo/nhanes.pbart1.R', package='BART')} for chronic lower back pain and\\
\code{system.file('demo/nhanes.pbart2.R', package='BART')} for chronic
neck pain.  In Figure~\ref{clbp}, the unweighted relationship between
chronic lower back pain, BMI and gender are displayed: males (females)
are represented by blue (red) lines.  As you can see, there is a
non-linear relationship between the probability of chronic lower back
pain and BMI for both genders where females have a parallel higher
probability than males.  For frail and underweight, the probability is
high and drops as BMI increases until about 35 kg/m$^2$ and afterwards
increases until about 65 kg/m$^2$ and then is flat.  Based on sampling
weights, the results are basically the same (not shown).  In
Figure~\ref{neck}, the unweighted relationship between chronic neck
pain, BMI and gender are displayed: males (females) are represented by
blue (red) lines.  As you can see, there appears to be no relationship
between the probability of chronic neck pain and BMI for both
genders where females have a nearly parallel higher probability than
males.  Based on sampling weights (not shown), the results are
basically the same.

\bibliography{references}

\begin{figure}[t!]
\begin{center}
\includegraphics[scale=0.425]{figures/geweke-pbart2-100.pdf}
\end{center}
\caption{Geweke convergence diagnostics for probit BART: $N=100$.
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 the
$Z_{AB}$ exceed the 95\% limits only a handful of times.  Based on this figure, 
we conclude that the chains have converged.
\label{geweke-pbart2-100}}
\end{figure}

\begin{figure}[t!]
\begin{center}
\includegraphics[scale=0.425]{figures/geweke-pbart2-1000.pdf}
\end{center}
\caption{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}.
\label{geweke-pbart2-1000}}
\end{figure}

\begin{figure}[t!]
\begin{center}
\includegraphics[scale=0.425]{figures/geweke-pbart2-10000.pdf}
\end{center}
\caption{Geweke convergence diagnostics for probit BART: $N=10000$.
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}.
\label{geweke-pbart2-10000}}
\end{figure}

\begin{figure}[t!]
\begin{center}
\includegraphics[scale=0.425]{figures/sparse-pbart.pdf}
\end{center}
\caption{
  Probit BART and variable selection example.
%: \code{system.file('demo/sparse.pbart.R', package='BART')}.  
  For sample sizes of $N=100, 1000, 10000$, 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=100$, only $s_2\le P^{-1}$, but notice that there are
  34 false positives.  For $N=1000$, all five active covariates are
  identified, but notice that there are 18 false positives.  For
  $N=10000$, all five active covariates are identified and notice that
  there are only two false positives.  \label{varsel}}
\end{figure}

\begin{figure}[t!]
\begin{center}
\includegraphics[scale=0.425]{figures/clbp.pdf}
\end{center}
\caption{Friedman's partial dependence function:  BMI and
probability of chronic lower back pain.  The unweighted relationship between
chronic lower back pain, BMI and gender are displayed: males (females)
are represented by blue (red) lines.  As you can see, there is a
non-linear relationship between the probability of chronic lower back
pain and BMI for both genders where females have a parallel higher
probability than males.  For frail and underweight, the probability is
high and drops as BMI increases until about 35 kg/m$^2$ and afterwards
increases until about 65 kg/m$^2$ and then is flat.  Based on sampling
weights (not shown), the results are basically the same 
\label{clbp}}
\end{figure}

\begin{figure}[t!]
\begin{center}
\includegraphics[scale=0.425]{figures/neck.pdf}
\end{center}
\caption{Friedman's partial dependence function: BMI and 
probability of chronic neck pain.  The unweighted relationship between chronic neck
pain, BMI and gender are displayed: males (females) are represented by
blue (red) lines.  As you can see, there appears to be no relationship
between the probability of chronic neck pain and BMI for both
genders where females have a nearly parallel higher probability than
males.  Based on sampling weights (not shown), the results are
basically the same.
\label{neck}}
\end{figure}

\end{document}

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