---
title: "BART::wbart: BART for Numeric Outcomes"
author: "Robert McCulloch and Rodney Sparapani"
output:
  pdf_document :
    toc: true
    toc_depth: 3
    number_sections: yes
vignette: >
  %\VignetteIndexEntry{wbart, BART for Numeric Outcomes}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
```
```{r colfun,echo=FALSE}
library(knitr)
#Color Format
colFmt = function(x,color){
   outputFormat = opts_knit$get("rmarkdown.pandoc.to")
   if(outputFormat == 'latex')
     paste("\\textcolor{",color,"}{",x,"}",sep="")
   else if(outputFormat == 'html')
     paste("<font color='",color,"'>",x,"</font>",sep="")
   else
    x
}
```


\  

\  


<!-- section: %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%  -->
# BART

\  

BART is Bayesian Additive Regression Trees (see Chipman, George, and McCulloch).


We fit the basic model:

$$
Y_i = f(x_i) + \epsilon_i, \;\; \epsilon_i \sim N(0,\sigma^2)
$$

\  

We use Markov Chain Monte Carlo to get draws from the posterior distribution
of the parameter $(f,\sigma)$.

In this vignette we look at BART::wbart which is the basic function in the R package **BART**.

\  

<!-- slide: %%%%%%%%%%%%%%%%%%%%%%%%%% -->
***
## Boston Housing Data

\ 

Let's just use the good old Boston housing data.

We'll predict the median house value, y=mdev, from  x1 = rm (number of rooms) and x2=lsat (% lower status).

\  

```{r s1-1, include=TRUE, echo=TRUE,cache=TRUE}
library(MASS)
x = Boston[,c(6,13)] #rm=number of rooms and lstat= percent lower status
y = Boston$medv # median value
head(cbind(x,y))
```

\  

<!-- slide: %%%%%%%%%%%%%%%%%%%%%%%%%% -->
***
## A Quick Look at the Data  

\  

```{r pl-dat, include=TRUE, echo=TRUE,out.width='80%',fig.align='center', dependson="s1-1"}
par(mfrow=c(2,2))
par(mai=c(.8,.8,.2,.2))
plot(x[,1],y,xlab="x1=rm",ylab="y=mdev",cex.axis=1.3,cex.lab=1.2)
plot(x[,2],y,xlab="x2=lstat",ylab="y=mdev",cex.axis=1.3,cex.lab=1.2)
plot(x[,1],x[,2],xlab="x1=rm",ylab="x2=lstat",cex.axis=1.3,cex.lab=1.2)
```

<!-- slide: %%%%%%%%%%%%%%%%%%%%%%%%%% -->
***
## Run wbart


```{r s1-2, include=TRUE, echo=TRUE,cache=TRUE,dependson="s1-1",message=FALSE}
library(BART) #BART package
set.seed(99) #MCMC, so set the seed
nd=200 # number of kept draws
burn=50 # number of burn in draws
bf = wbart(x,y,nskip=burn,ndpost=nd)
```

<!-- slide: %%%%%%%%%%%%%%%%%%%%%%%%%% -->
***
## Results returned with a list

\  
 
We stored the results of running wbart in the list **bf**.

\  


```{r, include=TRUE, echo=TRUE,dependson=paste0("s1-",1:2),collapse=TRUE}
names(bf)
length(bf$sigma)
length(bf$yhat.train.mean)
dim(bf$yhat.train)
```

\  

*`r colFmt("Remember","red")`*, the training data has $n=$ `r length(y)` observations, we had
burn=`r burn` burn-in draws and nd=`r nd` kept draws.

Let's look at a couple key list components:

* **sigma**: burnin + kept (burn+nd) draws of $\sigma$.
* **yhat.train.mean**: $j^{th}$ value is posterior mean of $f(x_j)$, \ $f$ evaluated at the $j^{th}$
           training observation.
* **yhat.train**: $i,j$ value  is the $i^{th}$ kept MCMC draw of \ $f(x_j)$.

\  

<!-- slide: %%%%%%%%%%%%%%%%%%%%%%%%%% -->
***
## Assess Convegence with $\sigma$ Draws

\  

As with any high-dimensional MCMC, assessing convergence may be tricky.  

A nice simple thing to look at is the draws of $\sigma$.
The parameter $\sigma$ is the only identified parameter in the model and it also
gives us a sense of the size of the errors.  

\  


```{r, include=TRUE, echo=TRUE,dependson="s1-2", out.width='50%',fig.align='center'}
plot(bf$sigma)
abline(v=burn,lwd=2,col="red")
```

\  

Look's like it burned in almost right away.  
Just one initial draw looking a bit bigger than the rest.
Hopefully, subsequent variation is legitimate posterior variation.  

In a more difficult problem you may see the $\sigma$ draws initially declining as the MCMC search for fit.

\  

<!-- slide: %%%%%%%%%%%%%%%%%%%%%%%%%% -->
***
## Look at in-sample Fit and Compare to a Linear Fit

\  

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


\  

```{r , include=TRUE, echo=TRUE,out.width='70%',fig.align='center',dependson=paste0("s1-",1:2)}
lmf = lm(y~.,data.frame(x,y))
fitmat = cbind(y,bf$yhat.train.mean,lmf$fitted.values)
colnames(fitmat)=c("y","BART","Linear")
cor(fitmat)
pairs(fitmat)
```

\  

The BART fit is noticeably different from the linear fit.

\  

<!-- slide: %%%%%%%%%%%%%%%%%%%%%%%%%% -->
***
## A Quick Look at the Uncertainty

We order the observations by the fitted house value (yhat.train.mean)  and then use boxplots to display the draws of
$f(x)$ 
in each column of yhat.train.

```{r , include=TRUE, echo=TRUE,out.width='70%',fig.align='center',dependson="s1-2"}
ii = order(bf$yhat.train.mean) #order observations by predicted value
boxplot(bf$yhat.train[,ii]) #boxplots of f(x) draws
```

Substantial predictive uncertainty, but you are still pretty
sure some houses should cost more that others!!

\  


<!-- section: %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%  -->
# Using predict.wbart

We can get out of sample predictions two ways. 

First, we can can just ask for them when we call wbart by supply a matrix of test $x$ value. 

Second, we can call a predict method.


<!-- slide: %%%%%%%%%%%%%%%%%%%%%%%%%% -->
***
##  Train and Test Data Sets

\ 

Let's split our data into train and test subsets.  

\  

```{r ttsplit, include=TRUE, echo=TRUE,cache=TRUE,dependson="s1-1"}
n=length(y) #total sample size
set.seed(14)  # Dave Keon, greatest Leaf of all time!
ii = sample(1:n,floor(.75*n)) # indices for train data, 75% of data
xtrain=x[ii,]; ytrain=y[ii] # training data
xtest=x[-ii,]; ytest=y[-ii] # test data
cat("train sample size is ",length(ytrain)," and test sample size is ",length(ytest),"\n")
```


And now we can run `wbart` using the train to learn and predict at `xtest`.

First, we'll just give `xtest` to the `wbart` call.

```{r pred1,include=TRUE, echo=TRUE,cache=TRUE,message=FALSE,results="hide",dependson="ttsplit"}
set.seed(99)
bfp1 = wbart(xtrain,ytrain,xtest) #predict.wbart wants a matrix
```

<!-- slide: %%%%%%%%%%%%%%%%%%%%%%%%%% -->
***

```{r , include=TRUE, echo=TRUE,dependson="pred1",collapse=TRUE}
dim(bfp1$yhat.test)
length(bfp1$yhat.test.mean)
```

\  

Now

* **yhat.test**: $i,j$ value  is the $i^{th}$ kept MCMC draw of $f(x_j)$ where $x_j$ is the
            $j^{th}$ row of `xtest`.
* **yhat.test.mean**: $j^{th}$ value is posterior mean of $f(x_j)$, $f$ evaluated at the $j^{th}$
           row of `xtest`.
           

\ 

<!-- slide: %%%%%%%%%%%%%%%%%%%%%%%%%% -->
***

\  

Alternatively, we could run wbart saving all the MCMC results and then call predict.wbart.  

\  


```{r pred2,include=TRUE, echo=TRUE,cache=TRUE,message=FALSE,results="hide",dependson="ttsplit"}
set.seed(99)
bfp2 = wbart(xtrain,ytrain)
yhat = predict(bfp2,as.matrix(xtest)) #predict wants a matrix
```

Then `yhat` and `bfp1$yhat.test` are the same.

```{r , include=TRUE, echo=TRUE,dependson="pred2",collapse=TRUE}
dim(yhat)
summary(as.double(yhat-bfp1$yhat.test))
```

\  
\  


<!-- section: %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%  -->
#Thining
***


In our simple Boston housing data set `wbart` runs pretty fast.

But with more data and longer runs you may want to speed things up
by saving less and then using predict.

Let's just keep a thinned subset of 200 tree ensembles.


```{r bfthin,include=TRUE, echo=TRUE,cache=TRUE,message=FALSE,results="hide",dependson="ttsplit"}
set.seed(4) #Bobby Orr, let's change the seed
bfthin = wbart(xtrain,ytrain,nskip=1000,ndpost=10000,
                     nkeeptrain=0,nkeeptest=0,nkeeptestmean=0,nkeeptreedraws=200)
yhatthin = predict(bfthin,as.matrix(xtest)) #predict wants a matrix
```



```{r, include=TRUE, echo=TRUE,collapse=TRUE}
dim(bfthin$yhat.train)
dim(yhatthin)
```

\  

Now there are no kept draws of $f(x)$ for training $x$, and we have 200 tree ensembles
to use in predict.wbart.  

\  

<!-- slide: %%%%%%%%%%%%%%%%%%%%%%%%%% -->
***
##The thinning arguments:

\  

* **nkeeptrain** : number of $f(x)$ draws to save for training $x$.
* **nkeeptest** : number of $f(x)$ draws to save for test $x$.
* **nkeeptestmeam** : number of draws to use in computing yhat.test.mean.
* **nkeeptreedraws** : number of tree ensembles to keep.

\  

The default values are to keep all the draws (e.g. nkeeptrain=ndpost).

Of course, if you keep 100 out of 100,000, you keep every 1,000th draw.  

\  

<!-- slide: %%%%%%%%%%%%%%%%%%%%%%%%%% -->
***
## Let's have a look at the predictions

\  




```{r, include=TRUE, echo=TRUE,out.width='60%',fig.align='center',dependson=c("ttsplit","bfthin","pred1")}
fmat=cbind(ytest,bfp1$yhat.test.mean,apply(yhatthin,2,mean))
colnames(fmat) = c("y","BARTpred","BARTpredThin")
pairs(fmat)
```

\  

*`r colFmt("Recall","red")`*, the **BARTpred** predictions are from a run are from a BART run with seed=99 and all default values.  

\  

The **BARTpredThin** are from 200 kept trees out of a long run with 1,000 burn-in and 10,000 kept draws and seed=4.  


\ 

*`r colFmt("Interesting how similar they are !!!!","red")`*

\  



