We will show how to grow a tree for classification. Regression trees are developed in a very similar way.

As before, these codes are taken from the Introduction to Statistical Learning book.

In these data, Sales is a continuous variable, and so we begin by recoding it as a binary variable.

library(tree)
library(ISLR)
attach(Carseats)
High=ifelse(Sales<=8,"No","Yes")
Carseats=data.frame(Carseats,High)

tree.carseats=tree(High~.-Sales,Carseats)
summary(tree.carseats)

## 
## Classification tree:
## tree(formula = High ~ . - Sales, data = Carseats)
## Variables actually used in tree construction:
## [1] "ShelveLoc"   "Price"       "Income"      "CompPrice"   "Population" 
## [6] "Advertising" "Age"         "US"         
## Number of terminal nodes:  27 
## Residual mean deviance:  0.4575 = 170.7 / 373 
## Misclassification error rate: 0.09 = 36 / 400

We see that the training error rate is 9%.

For classification trees, the deviance reported in the output of summary() is given by \[ -2 \sum_{m} \sum_k n_{mk} \log(\hat{p}_{mk}) \] where \(n_{mk}\) is the number of observations in the \(m\)th terminal node that belong to the \(k\)th class.

A small deviance indicates a tree that provides a good fit to the (training) data.

The residual mean deviance reported is simply the deviance divided by \(n - |T_0|\), which in this case is 400 - 27 = 373.

• One of the most attractive properties of trees is that they can be graphically displayed.
• We use the plot() function to display the tree structure, and the text() function to display the node labels.

plot(tree.carseats)

plot(tree.carseats)
text(tree.carseats,pretty=0)

• If we just type the name of the tree object, R prints output corresponding to each branch of the tree.
• R displays the split criterion (e.g. Price<92.5), the number of observations in that branch, the deviance, the overall prediction for the branch (Yes or No), and the fraction of observations in that branch that take on values of Yes and No.

• Need test error for proper evaluation.
• Split the observations into a training set and a test set, build the tree using the training set, and evaluate its performance on the test data.
• The predict() function can be used for this purpose.
• In the case of a classification tree, the argument type="class" instructs R to return the actual class prediction.

set.seed(2)
train=sample(1:nrow(Carseats), 200)
Carseats.test=Carseats[-train,]
High.test=High[-train]
tree.carseats=tree(High~.-Sales,Carseats,subset=train)
tree.pred=predict(tree.carseats,Carseats.test,type="class")
table(tree.pred,High.test)

##          High.test
## tree.pred No Yes
##       No  86  27
##       Yes 30  57

This approach leads to correct predictions for 0.715 cases.

• What you just saw was the full tree.
• Next, we consider whether pruning the tree might lead to improved results.
• The function cv.tree() performs cross-validation in order to determine the optimal level of tree complexity; cost complexity pruning is used in order to select a sequence of trees for consideration.

• We use the argument FUN=prune.misclass in order to indicate that we want the classification error rate to guide the cross-validation and pruning process, rather than the default for the cv.tree() function, which is deviance.

• The cv.tree() function reports

1. the number of terminal nodes of each tree considered (size),
2. the corresponding error rate, and
3. the value of the cost-complexity parameter used (k corresponds to our \(\alpha\))

• Cost-complexity pruning is similar to regularization idea. You put a penalty for the size of the tree - larger trees are penalized more.

\[ \text{minimize} \sum_{m = 1}^{|T|} \sum_{x_i \in R_m} 1(y_i \neq y_{R_m}) + \alpha |T| \]

• Just like regularized methods, \(\alpha\) is a tuning parameter - we can choose a suitable value by cross-validation.

• Note that, despite the name, dev corresponds to the cross-validation error rate in this instance.
• The tree with 9 terminal nodes results in the lowest cross-validation error rate, with 50 cross-validation errors.

par(mfrow=c(1,2))
plot(cv.carseats$size,cv.carseats$dev,type="b")
plot(cv.carseats$k,cv.carseats$dev,type="b")

• How does the pruned tree do in terms of prediction?

tree.pred=predict(prune.carseats,Carseats.test,type="class")
table(tree.pred,High.test)

##          High.test
## tree.pred No Yes
##       No  94  24
##       Yes 22  60

Now precentage of correct prediction is 0.77.

• You can choose a different k, e.g. \(k = 15\).

prune.carseats=prune.misclass(tree.carseats,best=15)
tree.pred=predict(prune.carseats,Carseats.test,type="class")
table(tree.pred,High.test)

##          High.test
## tree.pred No Yes
##       No  86  22
##       Yes 30  62

• Percentage of correct predictions will be lower: 0.74

• Compare to the best model by CV: 0.77

par(mfrow=c(1,2))
plot(Carseats.test$ShelveLoc,glm.probs); abline(h=0.5)
plot(Carseats.test$Price,glm.probs);abline(h=0.5)

1. Plot the ROC curve for both logistic and decision tree (with pruning).
2. Also apply k-NN with CV and LDA on the same data-set and compare accuracies.
3. Can you create an artificial data-set where decision tree will do better than logistic or LDA?

• Linear Regression: \[ f(X) = \beta_0 + \sum_{j=1}^{p}X_j \beta_j \]

• Regression Tree Method: \[ f(X) = \sum_{m=1}^{M} c_m 1_{(X \in R_m)} \] where \(R_1, \ldots, R_m\) is a partition of the \(X\)-space.

• Which model is better? It depends on the problem at hand.

If true decision boundary is linear, classical methods work well, if it's non-linear tree-based methods might work better.

tree.car=tree(Sales~.,Carseats,subset=train)
summary(tree.car)

## 
## Regression tree:
## tree(formula = Sales ~ ., data = Carseats, subset = train)
## Variables actually used in tree construction:
## [1] "High"      "ShelveLoc" "CompPrice" "Price"    
## Number of terminal nodes:  10 
## Residual mean deviance:  1.798 = 341.6 / 190 
## Distribution of residuals:
##     Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
## -3.64200 -0.93070 -0.08071  0.00000  1.10800  3.09100

plot(tree.car, type = "unif")
text(tree.car,pretty=0,cex=0.7)

plot(cv.car)

plot(cv.car$size,cv.car$dev,type='b')

prune.car=prune.tree(tree.car,best=9)
plot(prune.car,type = "unif")
text(prune.car,pretty=0,cex=0.7)

Here we apply bagging and random forests to the Boston data, using the randomForest package in R.

• Ensemble Learning: Generae many predictors and aggregate their results.
• Two well-known methods: Boosting (Shapire et al. , 1998) and Bagging (Breiman, 1996)
• In boosting, successive trees give extra weight to points incorrectly predicted by earlier predictors. In the end, a weighted vote is taken for prediction.
• In bagging, successive trees do not depend on earlier trees - each is independently constructed using a bootstrap sample of the data set.

library(MASS);data(Boston)
str(Boston)

## 'data.frame':    506 obs. of  14 variables:
##  $ crim   : num  0.00632 0.02731 0.02729 0.03237 0.06905 ...
##  $ zn     : num  18 0 0 0 0 0 12.5 12.5 12.5 12.5 ...
##  $ indus  : num  2.31 7.07 7.07 2.18 2.18 2.18 7.87 7.87 7.87 7.87 ...
##  $ chas   : int  0 0 0 0 0 0 0 0 0 0 ...
##  $ nox    : num  0.538 0.469 0.469 0.458 0.458 0.458 0.524 0.524 0.524 0.524 ...
##  $ rm     : num  6.58 6.42 7.18 7 7.15 ...
##  $ age    : num  65.2 78.9 61.1 45.8 54.2 58.7 66.6 96.1 100 85.9 ...
##  $ dis    : num  4.09 4.97 4.97 6.06 6.06 ...
##  $ rad    : int  1 2 2 3 3 3 5 5 5 5 ...
##  $ tax    : num  296 242 242 222 222 222 311 311 311 311 ...
##  $ ptratio: num  15.3 17.8 17.8 18.7 18.7 18.7 15.2 15.2 15.2 15.2 ...
##  $ black  : num  397 397 393 395 397 ...
##  $ lstat  : num  4.98 9.14 4.03 2.94 5.33 ...
##  $ medv   : num  24 21.6 34.7 33.4 36.2 28.7 22.9 27.1 16.5 18.9 ...

• Bagging is a special case of a random forest with \(m = p\).
• randomForest() can perform both random forests and bagging.

library(randomForest)
train = sample(1:nrow(Boston), nrow(Boston)/2)
boston.test=Boston[-train,"medv"]
bag.boston=randomForest(medv~.,data=Boston,subset=train,mtry=13,
                        importance=TRUE)

• The argument mtry=13 indicates that all 13 predictors should be considered for each split of the tree.
• In other words, bagging should be done.
• importance=TRUE: should importance of the predictors be assessed?

bag.boston

## 
## Call:
##  randomForest(formula = medv ~ ., data = Boston, mtry = 13, importance = TRUE,      subset = train) 
##                Type of random forest: regression
##                      Number of trees: 500
## No. of variables tried at each split: 13
## 
##           Mean of squared residuals: 12.35287
##                     % Var explained: 85.91

yhat.bag = predict(bag.boston,newdata=Boston[-train,])
mean((yhat.bag-boston.test)^2)

## [1] 13.67243

plot(yhat.bag, boston.test); abline(0,1)

tree.boston=tree(medv~.,Boston,subset=train)
yhat=predict(tree.boston,newdata=Boston[-train,])
boston.test=Boston[-train,"medv"]
mean((yhat-boston.test)^2)

## [1] 27.07243

plot(yhat,boston.test);abline(0,1)

• By default ntree = 500.
• We could change the number of trees grown by randomForest() using the ntree argument:

bag.boston=randomForest(medv~.,data=Boston,subset=train,mtry=13,
                        ntree=25)
yhat.bag = predict(bag.boston,newdata=Boston[-train,])
mean((yhat.bag-boston.test)^2)

## [1] 14.35049

• This is still bagging as mtry = 13, but the test error has decreased.

• Growing a random forest proceeds in exactly the same way, except that we use a smaller value of the mtry argument.
• By default, randomForest() uses \(p/3\) variables when building a random forest of regression trees, and \([\sqrt{p}]\) variables when building a random forest of classification trees.
• Here we use mtry = 6.

set.seed(1)
rf.boston=randomForest(medv~.,data=Boston,subset=train,mtry=6,importance=TRUE)
yhat.rf = predict(rf.boston,newdata=Boston[-train,])
(mse.rf <- mean((yhat.rf-boston.test)^2))

## [1] 12.24709

The test set MSE is 12.2470904: this indicates that random forests yielded an improvement over bagging in this case.

varImpPlot(rf.boston)

ntreeset = seq(10,300,by = 10)
mse.bag = rep(0,length(ntreeset));mse.rf = rep(0,length(ntreeset))

for(i in 1:length(ntreeset)){
  nt = ntreeset[i]
  bag.boston=randomForest(medv~.,data=Boston,subset=train,mtry=13,ntree=nt)
  yhat.bag = predict(bag.boston,newdata=Boston[-train,])
  mse.bag[i] = mean((yhat.bag-boston.test)^2)
  
  rf.boston=randomForest(medv~.,data=Boston,subset=train,ntree=nt)
  yhat.bag = predict(rf.boston,newdata=Boston[-train,])
  mse.rf[i] = mean((yhat.bag-boston.test)^2)
}

plot(ntreeset,mse.bag,type="l",col=2,ylim=c(5,20))
lines(ntreeset,mse.rf,col=3)
legend("bottomright",c("Bagging","Random Forest"),col=c(2,3),lty=c(1,1))

ntreeset = seq(10,300,by = 10)
mse.bag = rep(0,length(ntreeset));mse.rf1 = rep(0,length(ntreeset));
mse.rf2 = rep(0,length(ntreeset));

for(i in 1:length(ntreeset)){
  nt = ntreeset[i]
  bag.boston=randomForest(medv~.,data=Boston,subset=train,mtry=13,ntree=nt)
  yhat.bag = predict(bag.boston,newdata=Boston[-train,])
  mse.bag[i] = mean((yhat.bag-boston.test)^2)
  
  rf.boston=randomForest(medv~.,data=Boston,mtry = 6, subset=train,ntree=nt)
  yhat.bag = predict(rf.boston,newdata=Boston[-train,])
  mse.rf1[i] = mean((yhat.bag-boston.test)^2)
  
  rf.boston=randomForest(medv~.,data=Boston,mtry = 4, subset=train,ntree=nt)
  yhat.bag = predict(rf.boston,newdata=Boston[-train,])
  mse.rf2[i] = mean((yhat.bag-boston.test)^2)
}

plot(ntreeset,mse.bag,type="l",col=2,ylim=c(5,20))
lines(ntreeset,mse.rf1,col=3); lines(ntreeset,mse.rf2,col=4)
legend("bottomright",c("Bagging","Random Forest (m = 6)",
                       "Random Forest (m = 4)"),col=c(2,3,4), lty=c(1,1,1))