# Fitting a simple regression tree


First set the right working directory (you need to change the directory `dataDir` to the local directory 
in your machine) 

```{r }
dataDir <- "../finalDataSets"
```

We bring in the bodyfat dataset that we used in the regression chapter.

```{r }
body = read.csv(file.path(dataDir, "bodyfat_short.csv"), header = T,stringsAsFactors=TRUE)
```

To fit a decision tree, we use the function `rpart`.
`rpart` stands for recursive partitioning. It is available in the R package `rpart`

```{r }
library(rpart)
rt = rpart(BODYFAT ~ AGE + WEIGHT + HEIGHT + CHEST + ABDOMEN + HIP + THIGH, data = body)
```

For now we will focus on just visualizing the results, and then go back and explore the output `rt` more. 

There is a built in plotting command for the output `plot`. This plots just the tree. If you want the labels plotted on top (and you do!), you have to then call `text` afterward.

```{r }
par(mar=c(.1,.1,.1,.1))
plot(rt)
text(rt)
```

Now we will run `rpart` on the college data, to demonstrate how `rpart` handles categorical explanatory variables. First we read in the data:

```{r }
scorecard = read.csv(file.path(dataDir, "college.csv"), stringsAsFactors=TRUE)
```

We will convert the variable `CONTROL` into a factor. Recall that `CONTROL` is a categorical variable describing what kind of college the observation is.

```{r }
scorecard$CONTROL<-factor(scorecard$CONTROL, labels=c("public","private","for-profit"))
```

We will now fit a very simple tree, with just two explanatory variables. This is not because this is the best way to fit a tree, but for illustration purposes we want to make sure that the resulting decision tree uses the categorical `CONTROL` variable.

```{r }
collegeTree = rpart(RET_FT4 ~ TUITIONFEE_OUT + CONTROL, data = scorecard)
par(mar=rep(.1,4))
plot(collegeTree)
text(collegeTree)
```

Remember that the plot of the tree doesn't actually show the names of the levels of the categorical variable, but instead uses "a","b","c",... correponding to the levels of the categorical variable. We can find the levels of `CONTROL` and interpret that "b" corresponds to "private", and "c" corresponds to "for-profit". 

```{r }
levels(scorecard$CONTROL)
```

The illustrations of the tree are nice, but rather messy and hard to read. Moreover, it can be hard to remember whether you should go left or right when you satisfy the condition. And dealing with the levels of a categorical variable is quite difficult. We can see the "tree" structure in text format by using the `print` command on our output. 

```{r }
print(rt)
```

This allows us to more precisely see how to traverse the tree, as well as the levels of the categorical variables.


```{r }
(4698.255 + 4358.48)/17578.99
```

We can look at the choice of where to prune the tree using the function `printcp`. This allows us to decide whether the default R value (0.01) was a good idea. 

```{r }
printcp(rt)
```

Recall that `printcp` only prints out the results *up to* the value of `cp` given to the original call to `rpart`, in this case the default of `cp=0.01`. If you want to go beyond that value, you need to set `cp` higher in the original call to `rpart`.

Here we demonstrate changing the `cp` value to be slightly higher based on the results of `printcp` (though the difference between `0.0122` and `0.01` were probably not really very meaningful)

```{r }
rtd = rpart(BODYFAT ~ AGE + WEIGHT + HEIGHT + CHEST + ABDOMEN + HIP + THIGH, data = body, cp = 0.0122)
plot(rtd)
text(rtd)
```

# Classification
We will not demonstrate classification trees using the the email spam dataset from the chapter on logistic regression. 

```{r }
library(DAAG)
data(spam7)
spam = spam7
```

We again call `rpart` only now we need the additional argument `method="class"`. This tells `rpart` to build a classification tree instead of a regression tree.

```{r }
sprt = rpart(yesno ~ crl.tot + dollar + bang + money + n000 + make, method = "class", data = spam)
plot(sprt)
text(sprt)
```

We can similarly use `printcp` to evaluate the pruning of the classification tree.

```{r }
printcp(sprt)
```

The results of `printcp` show that the prediction error is continuing to decrease up until 0.01. This indicates that we might want to try smaller values of `cp` in our definition of our classification tree, so that we can evaluate where a good stopping point is.

```{r }
sprt = rpart(yesno ~ crl.tot + dollar + bang + money + n000 + make, method = "class", cp = 0.001, data = spam)
```

Now we can see the results beyond `cp=0.01`, and we can see that a smaller value will perform better.

```{r }
printcp(sprt)
```

We can use a better value (`0.0028`) and refit our tree.

```{r }
sprt = rpart(yesno ~ crl.tot + dollar + bang + money + n000 + make, method = "class", cp = 0.0028, data = spam)
plot(sprt)
text(sprt)
```

we can also make predictions from our decision trees. We'll demonstrate this for our classification tree, but the main workhorse is the function `predict` which can be used for any output of `rpart`, whether classification or regression trees.

First we need to make a data.frame with our new observation(s) we want to get predictions for. We only predict for a single observations, but you can make a data.frame with as many observations as you want predictions.

```{r }
x0 = data.frame(crl.tot = 100, dollar = 3, bang = 0.33, money = 1.2, n000 = 0, make = .3)
```

Now we call `predict`, giving it our fitted tree along with the new data

```{r }
predict(sprt, newdata=x0)
```

The following code is copied over from the previous chapter (07), and calculates the confusion matrix for our data.

```{r }
confusion <- function (y, yhat, thres)
  {
    n <- length(thres)
    conf <- matrix(0,length(thres),ncol=4)
    colnames(conf) <- c("C0","W1","W0","C1")
    rownames(conf)<-round(thres,3)
    for ( i in 1:n)
      {
        C0 <- sum((!y) & (yhat<=thres[i]))
        W1 <- sum((!y) & (yhat>thres[i]))
        W0 <- sum((y) & (yhat<=thres[i])) 
        C1 <- sum((y) & (yhat>thres[i]))
        conf[i,] <- c(C0, W1, W0, C1)
      }
    return(conf)
}
makeConfMatrix<-function(v){
  v<-paste0("$",c("C_0","W_1","W_0","C_1"),"=",v,"$")
  knitr::kable(
  rbind(
    c("","$\\hat{y} = 0$","$\\hat{y} = 1$"),
    c("$y = 0$",v[1],v[2]),
    c("$y = 1$",v[3],v[4])
  )    ,
  align="r",
  escape = FALSE
)
}
```

We can use our predictions and our observed data to calculate the confusion matrix for our classification tree:

```{r }
y = as.numeric(spam$yesno == "y")
y.hat = predict(sprt, spam)[,2]
v<-seq(0.1, 0.9, by = 0.05) #a sequence of possible cutoff values
tree.conf = confusion(y, y.hat, thres=v)
```

We can plot the total errors ($W1+W0$) as a function of our cutoff threshold.

```{r }
plot(v,tree.conf[,"W1"]+tree.conf[,"W0"], xlab="threshold", ylab="Total error", type = "l")

makeConfMatrix(confusion(y, y.hat, thres=0.5))
```

# Random Forests
We use the package `randomForest` to fit our random forests. The syntax is very similar to that of `rpart`

```{r }
library(randomForest)
ft = randomForest(BODYFAT ~ AGE + WEIGHT + HEIGHT + CHEST + ABDOMEN + HIP + THIGH, data = body, importance = TRUE)
```

We chose `importance=TRUE` to have `randomForest` calculate importance measures for our variables for later use.

We can view a summary of our results:

```{r }
ft
```

We can see these importance measures for each variable in our model.

```{r }
importance(ft)
```

Again we can use the `predict` function to predict a new observations from our random forest

```{r }
x0 = data.frame(AGE = 40, WEIGHT = 170, HEIGHT = 76, CHEST = 120, ABDOMEN = 100, HIP = 101, THIGH = 60)
predict(ft, x0)
```

We can similarly use random forests for classification. We need to make sure that our response (`yesno`) is a factor. Unlike building the classification tree in `rpart`, we do not have to explicitly mention `method= class`; the fact the response is given as a factor variable tells `randomForest` to use a classification forest as opposed to a regression forest. 

```{r }
sprf = randomForest(as.factor(yesno) ~ crl.tot + dollar + bang + money + n000 + make, data = spam)
sprf
```

Again, prediction is exactly the same using the `predict` function:

```{r }
x0 = data.frame(crl.tot = 100, dollar = 3, bang = 0.33, money = 1.2, n000 = 0, make = 0.3)
predict(sprf, x0)
```

