*This was adapted from James, Witten, Hastie, and Ibshirani (2021) An Introduction to Statistical Learning.* The `tree` library is used to construct classification and regression trees. ```{r chunk1} library(tree) ``` We first use classification trees to analyze the `Carseats` data set. In these data, `Sales` is a continuous variable, and so we begin by recoding it as a binary variable. We use the `ifelse()` function to create a variable, called `High`, which takes on a value of `Yes` if the `Sales` variable exceeds $8$, and takes on a value of `No` otherwise. ```{r chunk2} library(ISLR2) attach(Carseats) High <- factor(ifelse(Sales <= 8, "No", "Yes")) ``` Finally, we use the `data.frame()` function to merge `High` with the rest of the `Carseats` data. ```{r chunk3} Carseats <- data.frame(Carseats, High) ``` We now use the `tree()` function to fit a classification tree in order to predict `High` using all variables but `Sales`. The syntax of the `tree()` function is quite similar to that of the `lm()` function. ```{r chunk4} tree.carseats <- tree(High ~ . - Sales, Carseats) ``` The `summary()` function lists the variables that are used as internal nodes in the tree, the number of terminal nodes, and the (training) error rate. ```{r chunk5} summary(tree.carseats) ``` 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. This is closely related to the entropy, defined in (8.7). 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. The argument `pretty = 0` instructs `R` to include the category names for any qualitative predictors, rather than simply displaying a letter for each category. ```{r chunk6} plot(tree.carseats) text(tree.carseats, pretty = 0) ``` The most important indicator of `Sales` appears to be shelving location, since the first branch differentiates `Good` locations from `Bad` and `Medium` locations. 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`. Branches that lead to terminal nodes are indicated using asterisks. ```{r chunk7} tree.carseats ``` In order to properly evaluate the performance of a classification tree on these data, we must estimate the test error rather than simply computing the training error. We 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. This approach leads to correct predictions for around $77 \%$ of the locations in the test data set. ```{r chunk8} 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) (104 + 50) / 200 ``` (If you re-run the `predict()` function then you might get slightly different results, due to "ties": for instance, this can happen when the training observations corresponding to a terminal node are evenly split between `Yes` and `No` response values.) 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 the number of terminal nodes of each tree considered (`size`) as well as the corresponding error rate and the value of the cost-complexity parameter used (`k`, which corresponds to $\alpha$ in (8.4)). ```{r chunk9} set.seed(7) cv.carseats <- cv.tree(tree.carseats, FUN = prune.misclass) names(cv.carseats) cv.carseats ``` Despite its name, `dev` corresponds to the number of cross-validation errors. The tree with 9 terminal nodes results in only 74 cross-validation errors. We plot the error rate as a function of both `size` and `k`. ```{r chunk10} par(mfrow = c(1, 2)) plot(cv.carseats$size, cv.carseats$dev, type = "b") plot(cv.carseats$k, cv.carseats$dev, type = "b") ``` We now apply the `prune.misclass()` function in order to prune the tree to obtain the nine-node tree. ```{r chunk11} prune.carseats <- prune.misclass(tree.carseats, best = 9) plot(prune.carseats) text(prune.carseats, pretty = 0) ``` How well does this pruned tree perform on the test data set? Once again, we apply the `predict()` function. ```{r chunk12} tree.pred <- predict(prune.carseats, Carseats.test, type = "class") table(tree.pred, High.test) (97 + 58) / 200 ``` Now $77.5 \%$ of the test observations are correctly classified, so not only has the pruning process produced a more interpretable tree, but it has also slightly improved the classification accuracy. If we increase the value of `best`, we obtain a larger pruned tree with lower classification accuracy: ```{r chunk13} prune.carseats <- prune.misclass(tree.carseats, best = 14) plot(prune.carseats) text(prune.carseats, pretty = 0) tree.pred <- predict(prune.carseats, Carseats.test, type = "class") table(tree.pred, High.test) (102 + 52) / 200 ``` ## Fitting Regression Trees Here we fit a regression tree to the `Boston` data set. First, we create a training set, and fit the tree to the training data. ```{r chunk14} set.seed(1) train <- sample(1:nrow(Boston), nrow(Boston) / 2) tree.boston <- tree(medv ~ ., Boston, subset = train) summary(tree.boston) ``` Notice that the output of `summary()` indicates that only four of the variables have been used in constructing the tree. In the context of a regression tree, the deviance is simply the sum of squared errors for the tree. We now plot the tree. ```{r chunk15} plot(tree.boston) text(tree.boston, pretty = 0) ``` The variable `lstat` measures the percentage of individuals with {lower socioeconomic status}, while the variable `rm` corresponds to the average number of rooms. The tree indicates that larger values of `rm`, or lower values of `lstat`, correspond to more expensive houses. For example, the tree predicts a median house price of $$45{,}400$ for homes in census tracts in which `rm >= 7.553`. It is worth noting that we could have fit a much bigger tree, by passing `control = tree.control(nobs = length(train), mindev = 0)` into the `tree()` function. Now we use the `cv.tree()` function to see whether pruning the tree will improve performance. ```{r chunk16} cv.boston <- cv.tree(tree.boston) plot(cv.boston$size, cv.boston$dev, type = "b") ``` In this case, the most complex tree under consideration is selected by cross-validation. However, if we wish to prune the tree, we could do so as follows, using the `prune.tree()` function: ```{r chunk17} prune.boston <- prune.tree(tree.boston, best = 5) plot(prune.boston) text(prune.boston, pretty = 0) ``` In keeping with the cross-validation results, we use the unpruned tree to make predictions on the test set. ```{r chunk18} yhat <- predict(tree.boston, newdata = Boston[-train, ]) boston.test <- Boston[-train, "medv"] plot(yhat, boston.test) abline(0, 1) mean((yhat - boston.test)^2) ``` In other words, the test set MSE associated with the regression tree is $35.29$. The square root of the MSE is therefore around $5.941$, indicating that this model leads to test predictions that are (on average) within approximately $$5{,}941$ of the true median home value for the census tract. \clearpage # On Your Own 1. This problem involves the OJ data set which is part of the `ISLR2` package. ```{r} data(OJ, package = 'ISLR2') ``` (a) Create a training set containing a random sample of 800 observations, and a test set containing the remaining observations. (b) Fit a tree to the training data, with Purchase as the response and the other variables as predictors. Use the summary() function to produce summary statistics about the tree, and describe the results obtained. What is the training error rate? How many terminal nodes does the tree have? (c) Type in the name of the tree object in order to get a detailed text output. Pick one of the terminal nodes, and interpret the information displayed. (d) Create a plot of the tree, and interpret the results. (e) Predict the response on the test data, and produce a confusion matrix comparing the test labels to the predicted test labels. What is the test error rate? (f) Apply the cv.tree() function to the training set in order to determine the optimal tree size. (g) Produce a plot with tree size on the x-axis and cross-validated classification error rate on the y-axis. (h) Which tree size corresponds to the lowest cross-validated classification error rate? (i) Produce a pruned tree corresponding to the optimal tree size obtained using cross-validation. If cross-validation does not lead to selection of a pruned tree, then create a pruned tree with five terminal nodes. (j) Compare the training error rates between the pruned and unpruned trees. Which is higher? (k) Compare the test error rates between the pruned and unpruned trees. Which is higher?