---
title: "Tree based models - 2a"
author: ''
output:
  html_document:
    toc: yes
  pdf_document: null
  pdf_presentation: null
  slidy_presentation:
    css: slidy.css
    widescreen: yes
  word_document: null
---

#Classification trees

   + Regression trees are used for continuous outputs and classification trees for discrete outputs
   
   + Regression trees split based on LSE, and Classification trees split based on GI

   + When fitting a classification tree, rather than choosing a split to minimize
   the total sum of squares, we choose a split to minimize some other measure of
   "node impurity".  If $\hat p_L$ and $\hat p_R$ are estimates of the probability
   that AHD="Yes" in the left and right hand branches of a split, the default is
   to choose the split which minimizes

   $$G= \hat p_L (1- \hat p_L) + \hat p_R (1- \hat p_R)$$

   + $G$ is called the "Gini index". You may recognize it as the estimated total variance,
   of the binary variable AHD in this case. Minimizing impurity maximizes purity.
   
   + think of $$\hat p_L$ as the probability that Yes output is going to the left branch.

   + Note that in this case some of the predictor variables are discrete.  For example,
   "sex" and "Thal".  There is only one possible split using "sex".  Ignoring missing
   values, there are three ways to split using "Thal".  

   + Generally the tree growing process generates a tree with too many splits.  This is
   an example  of **overfitting**. 

   + We would like to restrict the size of the tree in some manner.  We do this by cut back the
   initially grown tree.  This is known as **pruning**.

   + How far back should we prune the tree?  And along which branches?


# Review fitting process

```{r}
heart=read.csv("http://faculty.marshall.usc.edu/gareth-james/ISL/Heart.csv")
head(heart)
#summary(heart)
# the outcome is AHD (Yes or No)
nmiss=apply(is.na(heart),1,sum) #number of missing values, by row
heartm=heart[nmiss!=0,]  #cases with missing values
heart=heart[nmiss==0,]   #remove cases with missing values
#head(heartm)
# build a prediction model using the tree library
library(tree)
heart.tree=tree(AHD~., data=heart)
plot(heart.tree)
text(heart.tree)
print(heart.tree)
```


   
Build the cross-validation  
   


```{r}
# cross-validation  use the misclassification error rate to guide pruning
heartcv=cv.tree(heart.tree, FUN=prune.misclass)
```



Plot the result of cv

```{r}
print(heartcv)
plot(heartcv$size,heartcv$dev,type='b')
```



Estimate the best tree size (complexity)


```{r}
bestsize=heartcv$size[heartcv$dev==min(heartcv$dev)]
bestsize

```



Now use the bestsize to prune the tree just right


```{r}
heart.prune=prune.misclass(heart.tree,best=bestsize)

```

Have a look at the pruned tree


```{r}
heart.prune
plot(heart.prune,ylab="error rate")
text(heart.prune, pretty=0)
```


Use our final model to make predictions

   
```{r}
heart.pred=predict(heart.prune,data=heart, type="class")
```



Evaluate the prediction performances (confusion matrix): compare truth to predictions

   
```{r}
table(heart.pred,heart$AHD)
```


But problem: all data is used . Overly optimistic: try test-train

# Review evaluation process


Randomly split the data into a training and a test set of approximately the same size.
   
```{r}
nrec=dim(heart)[1]
nrec
 ntrain=ifelse(nrec%%2==0,nrec/2,(nrec+1)/2)
 ntrain
 train=sample(1:nrec,ntrain,replace=F)
 train
 htrain=heart[train,]
 htrain
 htest=heart[-train,]
 htest
```
 
   + Carry out the fitting on the training set
  
```{r}
attach(heart)
AHD.test=AHD[-train]
AHD.train=AHD[train]
heart.train=heart[train,]
heart.test=heart[-train,]
heart.tree=tree(AHD~., data=heart.train)
heartcv=cv.tree(heart.tree, FUN=prune.misclass)
plot(heartcv$size,heartcv$dev,type='b')
bestsize=heartcv$size[heartcv$dev==min(heartcv$dev)]
heart.prune=prune.misclass(heart.tree,best=bestsize)
print(heart.prune)
text(heart.prune)
```

   + Calculate the error rate for the training set.
```{r}
pred.train=predict(heart.prune,data=heart.train,type="class")
ctable=table(AHD.train,pred.train)
print(ctable)
trainrate= (sum(ctable)-sum(diag(ctable)))/sum(ctable)
```

   + Calculate the error rate for the test set.
   
```{r}
pred.test=predict(heart.prune,heart.test,type="class")
ctable=table(AHD.test,pred.test)
print(ctable)
testrate= (sum(ctable)-sum(diag(ctable)))/sum(ctable)
```
   
   + The misclassification  rate on the training set is
`r round( 100*trainrate,1)` %.  The misclassification rate on the test set is
`r round( 100*testrate,1)` %.  The test error rate is a more realistic measure
of the performance of the method.





