This is the first article about tree based methods using R. Carseats data in the chapter 8 lab of ISLR is used to perform classification analysis. Unlike the lab example, the rpart package is used to fit the CART model on the data and the caret package is used for tuning the pruning parameter (cp).

• Visualizations
• Pre-Processing
• Data Splitting
• Miscellaneous Model Functions
• Model Training and Tuning
• Using Custom Models
• Variable Importance
• Feature Selection: RFE, Filters, GA, SA
• Other Functions
• Parallel Processing

The pruning parameter in the rpart package is scaled so that its values are from 0 to 1. Specifically the formula is

$$R_{cp}\left(T\right)\equiv R\left(T\right) + cp*|T|*R\left(T_{1}\right)$$

where $$T_{1}$$ is the tree with no splits, $$\mid T\mid$$ is the number of splits for a tree and R is the risk.

Due to the inclusion of $$R\left(T_{1}\right)$$, when cp=1, the tree will result in no splits while it is not pruned when cp=0. On the other hand, in the original setup without the term, the pruning parameter ($$\alpha$$) can range from 0 to infinity.

Let’s get started.

The following packages are used.

1library(dplyr) # data minipulation
2library(rpart) # fit tree
3library(rpart.plot) # plot tree
4library(caret) # tune model


Carseats data is created as following while the response (Sales) is converted into a binary variable.

1require(ISLR)
2data(Carseats)
3Carseats = Carseats %>%
4  mutate(High=factor(ifelse(Sales<=8,"No","High"),labels=c("High","No")))
5# structure of predictors
6str(subset(Carseats,select=c(-High,-Sales)))

## 'data.frame':	400 obs. of  10 variables:
##  $CompPrice : num 138 111 113 117 141 124 115 136 132 132 ... ##$ Income     : num  73 48 35 100 64 113 105 81 110 113 ...
##  $Advertising: num 11 16 10 4 3 13 0 15 0 0 ... ##$ Population : num  276 260 269 466 340 501 45 425 108 131 ...
##  $Price : num 120 83 80 97 128 72 108 120 124 124 ... ##$ ShelveLoc  : Factor w/ 3 levels "Bad","Good","Medium": 1 2 3 3 1 1 3 2 3 3 ...
##  $Age : num 42 65 59 55 38 78 71 67 76 76 ... ##$ Education  : num  17 10 12 14 13 16 15 10 10 17 ...
##  $Urban : Factor w/ 2 levels "No","Yes": 2 2 2 2 2 1 2 2 1 1 ... ##$ US         : Factor w/ 2 levels "No","Yes": 2 2 2 2 1 2 1 2 1 2 ...

1# classification response summary
2with(Carseats,table(High))

## High
## High   No
##  164  236

1with(Carseats,table(High) / length(High))

## High
## High   No
## 0.41 0.59


The train and test data sets are split using createDataPartition().

1## Data Splitting
2set.seed(1237)
3trainIndex.cl = createDataPartition(Carseats$High, p=.8, list=FALSE, times=1) 4trainData.cl = subset(Carseats, select=c(-Sales))[trainIndex.cl,] 5testData.cl = subset(Carseats, select=c(-Sales))[-trainIndex.cl,]  Stratify sampling is performed by default. 1# training response summary 2with(trainData.cl,table(High))  ## High ## High No ## 132 189  1with(trainData.cl,table(High) / length(High))  ## High ## High No ## 0.411215 0.588785  1# test response summary 2with(testData.cl,table(High))  ## High ## High No ## 32 47  1with(testData.cl,table(High) / length(High))  ## High ## High No ## 0.4050633 0.5949367  The following resampling strategies are considered: cross-validation, repeated cross-validation and bootstrap. 1## train control 2trControl.cv = trainControl(method="cv",number=10) 3trControl.recv = trainControl(method="repeatedcv",number=10,repeats=5) 4trControl.boot = trainControl(method="boot",number=50)  There are two methods in the caret package: rpart and repart2. The first method allows the pruning parameter to be tuned. The tune grid is not set up explicitly and it is adjusted by tuneLength - equally spaced cp values are created from 0 to 0.3 in the package.  1set.seed(12357) 2fit.cl.cv = train(High ~ . 3 ,data=trainData.cl 4 ,method="rpart" 5 ,tuneLength=20 6 ,trControl=trControl.cv) 7 8fit.cl.recv = train(High ~ . 9 ,data=trainData.cl 10 ,method="rpart" 11 ,tuneLength=20 12 ,trControl=trControl.recv) 13 14fit.cl.boot = train(High ~ . 15 ,data=trainData.cl 16 ,method="rpart" 17 ,tuneLength=20 18 ,trControl=trControl.boot)  Repeated cross-validation and bootstrap produce the same best tuned cp value while cross-validation returns a higher value. 1# results at best tuned cp 2subset(fit.cl.cv$results,subset=cp==fit.cl.cv$bestTune$cp)

##           cp  Accuracy     Kappa AccuracySD   KappaSD
## 3 0.03110048 0.7382148 0.4395888 0.06115425 0.1355713

1subset(fit.cl.recv$results,subset=cp==fit.cl.recv$bestTune$cp)  ## cp Accuracy Kappa AccuracySD KappaSD ## 2 0.01555024 0.7384537 0.456367 0.08602102 0.1821254  1subset(fit.cl.boot$results,subset=cp==fit.cl.boot$bestTune$cp)

##           cp  Accuracy     Kappa AccuracySD    KappaSD
## 2 0.01555024 0.7178692 0.4163337 0.04207806 0.08445416


The one from repeated cross-validation is taken to fit to the entire training data.

Updated on Feb 10, 2015

• As a value of cp is entered in rpart(), the function fits the model up to the value and takes the result. Therefore it produces a pruned tree.
• If it is not set or set to be a low value (eg, 0), pruning can be done using the prune() function.
1## refit the model to the entire training data
2cp.cl = fit.cl.recv$bestTune$cp
3fit.cl = rpart(High ~ ., data=trainData.cl, control=rpart.control(cp=cp.cl))
4
5# Updated on Feb 10, 2015
6# prune if cp not set or too low
7# fit.cl = prune(tree=fit.cl, cp=cp.cl)


The resulting tree is shown as following. The plot shows expected losses and node probabilities in the final nodes. For example, the leftmost node has

• expected loss of 0.13 (= 8/60 (0.13333))
• node probability of 19% (= 60/321)
 1# plot tree
2cols <- ifelse(fit.cl$frame$yval == 1,"green4","darkred") # green if high
3prp(fit.cl
4    ,main="CART model tree"
5    ,extra=106           # display prob of survival and percent of obs
6    ,nn=TRUE             # display the node numbers
7    ,fallen.leaves=TRUE  # put the leaves on the bottom of the page
8    ,branch=.5           # change angle of branch lines
9    ,faclen=0            # do not abbreviate factor levels
10    ,trace=1             # print the automatically calculated cex
12    ,branch.lty=3        # draw branches using dotted lines
13    ,split.cex=1.2       # make the split text larger than the node text
14    ,split.prefix="is "  # put "is " before split text
15    ,split.suffix="?"    # put "?" after split text
16    ,col=cols, border.col=cols   # green if survived
17    ,split.box.col="lightgray"   # lightgray split boxes (default is white)
18    ,split.border.col="darkgray" # darkgray border on split boxes
19    ,split.round=.5)

## cex 0.7   xlim c(0, 1)   ylim c(-0.1, 1.1) The fitted model is applied to the test data.

1## apply to the test data
2pre.cl = predict(fit.cl, newdata=testData.cl, type="class")
3
4# confusion matrix
5cMat.cl = table(data.frame(response=pre.cl,truth=testData.cl\$High))
6cMat.cl

##         truth
## response High No
##     High   21  4
##     No     11 43

1# mean misclassification error
2mmse.cl = 1 - sum(diag(cMat.cl))/sum(cMat.cl)
3mmse.cl

##  0.1898734


More models are going to be implemented/compared in the subsequent articles.