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

#Classification trees

A **regression** tree models a continuous outcome. It is an alternative to
carrying out a multiple regression.

A **classification** tree models a discrete outcome.  In many cases, the
outcome variable is binary (having two values), in which case the classification
tree is an alternative to fitting a **logistic regression model**.  Logistic
regression is a topic covered in Stat4620 (Data Analysis).

Consider the **Heart** data set.

```{r}
heart=read.csv("http://faculty.marshall.usc.edu/gareth-james/ISL/Heart.csv")
head(heart)
```

The outcome variable AHD is binary, taking values "Yes" or "No".  This is
equivalent to a variable which equals 1 when AHD="Yes", and 0 otherwise.

```{r}
summary(heart)
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)
```


Let's predict AHD using a classification tree, with all other variables in
the data set used as potential predictors.

```{r}
library(tree)
heart.tree=tree(AHD~., data=heart)
plot(heart.tree)
text(heart.tree)
print(heart.tree)
```
   + 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?


```{r}
heartcv=cv.tree(heart.tree, FUN=prune.misclass)
print(heartcv)
plot(heartcv$size,heartcv$dev,type='b')
bestsize=heartcv$size[heartcv$dev==min(heartcv$dev)]

```

   + The argument **FUN=prune.misclass** used in the call to cv.tree says to use
the misclassification error rate to guide the pruning process, rather
than the deviance.  This is a reasonable choice for classification trees,
where the goal is generally to classify the discrete outcome.

   + Prune the tree back to the best size, `r bestsize`,  using **prune.misclass**

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

   + make a table of true class vs predicted class.  In the machine learning literature, this type of table is known as a **confusion matrix**.
   
```{r}
heart.pred=predict(heart.prune,data=heart, type="class")
table(heart.pred,heart$AHD)
```

## Using a test set to assess the misclassification rate

   + the error rate in the above analysis is unrealistically low in that
   all data were used to fit the tree, and all data were used to assess
   the performance of the classification
   
   + a better strategy to is intially split the dataset into **training** and
   **test** sets, do all of the fitting on the training set, and assess the
   error rate on the test set.
   
   ### 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.


   
   
# Gini index or coefficient

$$E=\frac{1}{n^2} \sum_{i,j} |x_i - x_j| = \bar{d_{ij}}$$
$M=\frac{1}{n}\sum x_i$
$$G=E/2M$$
difference / mean

There is a some alternative measures. For instance, Generalized Entropy Indexes, Atkinson Index, Piesch index, Kakwani index. Each measure can exhibit differnt properties.


$$G=\sum_{i \neq j } p_i p_j$$

$$1-p_1^2 -p_2^2$$

$$1-F$$



<img src="https://upload.wikimedia.org/wikipedia/commons/thumb/1/11/World_Income_Gini_Map_%282013%29.svg/1024px-World_Income_Gini_Map_%282013%29.svg.png">




