While the last three articles illustrated the CART model for both classification (with equal/unequal costs) and regression tasks, this article is rather technical as it compares three packages: rpart, caret and mlr. For those who are not familiar with the last two packages, they are wrappers (or frameworks) that implement a range of models (or algorithms) in a unified way. They are useful because inconsistent API could be a drawback of R (like other open source tools) and it would be quite beneficial if there is a way to implement different models in a standardized way. In line with the earlier articles, the Carseats data is used for a classification task.

Before getting started, I should admit the names are not defined effectively. I hope the below list may be helpful to follow the script.

  • package: rpt - rpart, crt - caret, mlr - mlr
  • model fitting: ftd - fit on training data, ptd - fit on test data
  • parameter selection
    • bst - best cp by 1-SE rule (recommended by rpart)
    • lst - best cp by highest accuracy (caret) or lowest mmce (mlr)
  • cost: eq - equal cost, uq - unequal cost (uq case not added in this article)
  • etc: cp - complexity parameter, mmce - mean misclassification error, acc - Accuracy (caret), cm - confusion matrix

Also the data is randomly split into trainData and testData. In practice, the latter is not observed and it is used here for evaludation.

Let’s get started.

The following packages are used.

1library(reshape2)
2library(plyr)
3library(dplyr)
4library(ggplot2)
5library(rpart)
6library(caret)
7library(mlr)

The Sales column is converted into a binary variables.

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

Balanced splitting of data can be performed in either of the packages as shown below.

 1## split data
 2# caret
 3set.seed(1237)
 4trainIndex = createDataPartition(Carseats$High, p=0.8, list=FALSE, times=1)
 5trainData.caret = data[trainIndex,]
 6testData.caret = data[-trainIndex,]
 7
 8# mlr
 9set.seed(1237)
10split.desc = makeResampleDesc(method="Holdout", stratify=TRUE, split=0.8)
11split.task = makeClassifTask(data=data,target="High")
12resampleIns = makeResampleInstance(desc=split.desc, task=split.task)
13trainData.mlr = data[resampleIns$train.inds[[1]],]
14testData.mlr = data[-resampleIns$train.inds[[1]],]
15
16# response summaries
17train.res.caret = with(trainData.caret,rbind(table(High),table(High)/length(High)))
18train.res.caret
##            High         No
## [1,] 132.000000 189.000000
## [2,]   0.411215   0.588785
1train.res.mlr = with(trainData.mlr,rbind(table(High),table(High)/length(High)))
2train.res.mlr
##             High          No
## [1,] 131.0000000 188.0000000
## [2,]   0.4106583   0.5893417
1test.res.caret = with(testData.caret,rbind(table(High),table(High)/length(High)))
2test.res.caret
##            High         No
## [1,] 32.0000000 47.0000000
## [2,]  0.4050633  0.5949367
1test.res.mlr = with(testData.mlr,rbind(table(High),table(High)/length(High)))
2test.res.mlr
##            High         No
## [1,] 33.0000000 48.0000000
## [2,]  0.4074074  0.5925926

Same to the previous articles, the split by the caret package is taken.

1# data from caret is taken in line with previous articles
2trainData = trainData.caret
3testData = testData.caret

Note that two custom functions are used: bestParam() and updateCM(). The former searches the cp values by the 1-SE rule (bst) and at the lowest xerror (lst) from the cp table of a rpart object. The latter produces a confusion matrix with model and use error, added to the last column and row respectively. Their sources can be seen here.

1source("src/mlUtils.R")

At first, the model is fit using the rpart package and bst and lst cp values are obtained.

1### rpart
2## train on training data
3set.seed(12357)
4mod.rpt.eq = rpart(High ~ ., data=trainData, control=rpart.control(cp=0))
5mod.rpt.eq.par = bestParam(mod.rpt.eq$cptable,"CP","xerror","xstd")
6mod.rpt.eq.par
##            lowest       best
## param  0.01136364 0.10606061
## error  0.65909091 0.70454545
## errStd 0.06033108 0.06157189

The selected cp values can be check graphically below.

 1# plot xerror vs cp
 2df = as.data.frame(mod.rpt.eq$cptable)
 3best = mod.rpt.eq.par
 4ubound = ifelse(best[2,1]+best[3,1]>max(df$xerror),max(df$xerror),best[2,1]+best[3,1])
 5lbound = ifelse(best[2,1]-best[3,1]<min(df$xerror),min(df$xerror),best[2,1]-best[3,1])
 6
 7ggplot(data=df[1:nrow(df),], aes(x=CP,y=xerror)) + 
 8  geom_line() + geom_point() +   
 9  geom_abline(intercept=ubound,slope=0, color="purple") + 
10  geom_abline(intercept=lbound,slope=0, color="purple") + 
11  geom_point(aes(x=best[1,2],y=best[2,2]),color="red",size=3) + 
12  geom_point(aes(x=best[1,1],y=best[2,1]),color="blue",size=3)

The original tree is pruned with the 2 cp values, resulting in 2 separate trees, and they are fit on the training data.

 1## performance on train data
 2mod.rpt.eq.lst.cp = mod.rpt.eq.par[1,1]
 3mod.rpt.eq.bst.cp = mod.rpt.eq.par[1,2]
 4# prune
 5mod.rpt.eq.lst = prune(mod.rpt.eq, cp=mod.rpt.eq.lst.cp)
 6mod.rpt.eq.bst = prune(mod.rpt.eq, cp=mod.rpt.eq.bst.cp)
 7
 8# fit to train data
 9mod.rpt.eq.lst.ftd = predict(mod.rpt.eq.lst, type="class")
10mod.rpt.eq.bst.ftd = predict(mod.rpt.eq.bst, type="class")

Details of the fitting is kept in a list (mmce).

  • pkg: package name
  • isTest: fit on test data?
  • isBest: cp by 1-SE rule?
  • isEq: equal cost?
  • cp: cp value used
  • mmce: mean misclassification error
 1# fit to train data
 2mod.rpt.eq.lst.ftd = predict(mod.rpt.eq.lst, type="class")
 3mod.rpt.eq.bst.ftd = predict(mod.rpt.eq.bst, type="class")
 4
 5# confusion matrix
 6mod.rpt.eq.lst.ftd.cm = table(actual=trainData$High,fitted=mod.rpt.eq.lst.ftd)
 7mod.rpt.eq.lst.ftd.cm = updateCM(mod.rpt.eq.lst.ftd.cm, type="Fitted")
 8
 9mod.rpt.eq.bst.ftd.cm = table(actual=trainData$High,fitted=mod.rpt.eq.bst.ftd)
10mod.rpt.eq.bst.ftd.cm = updateCM(mod.rpt.eq.bst.ftd.cm, type="Fitted")
11
12# misclassification error
13mod.rpt.eq.lst.ftd.mmce = list(pkg="rpart",isTest=FALSE,isBest=FALSE,isEq=TRUE
14                               ,cp=round(mod.rpt.eq.lst.cp,4)
15                               ,mmce=mod.rpt.eq.lst.ftd.cm[3,3])
16mmce = list(unlist(mod.rpt.eq.lst.ftd.mmce))
17
18mod.rpt.eq.bst.ftd.mmce = list(pkg="rpart",isTest=FALSE,isBest=TRUE,isEq=TRUE
19                               ,cp=round(mod.rpt.eq.bst.cp,4)
20                              ,mmce=mod.rpt.eq.bst.ftd.cm[3,3])
21mmce[[length(mmce)+1]] = unlist(mod.rpt.eq.bst.ftd.mmce)
22
23ldply(mmce)
##     pkg isTest isBest isEq     cp mmce
## 1 rpart  FALSE  FALSE TRUE 0.0114 0.16
## 2 rpart  FALSE   TRUE TRUE 0.1061 0.29

The pruned trees are fit into the test data and the same details are added to the list (mmce).

 1## performance of test data
 2# fit to test data
 3mod.rpt.eq.lst.ptd = predict(mod.rpt.eq.lst, newdata=testData, type="class")
 4mod.rpt.eq.bst.ptd = predict(mod.rpt.eq.bst, newdata=testData, type="class")
 5
 6# confusion matrix
 7mod.rpt.eq.lst.ptd.cm = table(actual=testData$High, fitted=mod.rpt.eq.lst.ptd)
 8mod.rpt.eq.lst.ptd.cm = updateCM(mod.rpt.eq.lst.ptd.cm)
 9
10mod.rpt.eq.bst.ptd.cm = table(actual=testData$High, fitted=mod.rpt.eq.bst.ptd)
11mod.rpt.eq.bst.ptd.cm = updateCM(mod.rpt.eq.bst.ptd.cm)
12
13# misclassification error
14mod.rpt.eq.lst.ptd.mmce = list(pkg="rpart",isTest=TRUE,isBest=FALSE,isEq=TRUE
15                               ,cp=round(mod.rpt.eq.lst.cp,4)
16                               ,mmce=mod.rpt.eq.lst.ptd.cm[3,3])
17mmce[[length(mmce)+1]] = unlist(mod.rpt.eq.lst.ptd.mmce)
18
19mod.rpt.eq.bst.ptd.mmce = list(pkg="rpart",isTest=TRUE,isBest=TRUE,isEq=TRUE
20                               ,cp=round(mod.rpt.eq.bst.cp,4)
21                              ,mmce=mod.rpt.eq.bst.ptd.cm[3,3])
22mmce[[length(mmce)+1]] = unlist(mod.rpt.eq.bst.ptd.mmce)
23
24ldply(mmce)
##     pkg isTest isBest isEq     cp mmce
## 1 rpart  FALSE  FALSE TRUE 0.0114 0.16
## 2 rpart  FALSE   TRUE TRUE 0.1061 0.29
## 3 rpart   TRUE  FALSE TRUE 0.0114 0.19
## 4 rpart   TRUE   TRUE TRUE 0.1061  0.3

Secondly the caret package is employed to implement the CART model.

1### caret
2## train on training data
3trControl = trainControl(method="repeatedcv", number=10, repeats=5)
4set.seed(12357)
5mod.crt.eq = caret::train(High ~ .
6                          ,data=trainData
7                          ,method="rpart"
8                          ,tuneLength=20
9                          ,trControl=trControl)

Note that the caret package select the best cp value that corresponds to the lowest Accuracy. Therefore the best cp by this package is labeled as lst to be consistent with the rpart package. And the bst cp is selected by the 1-SE rule. Note that, as the standard error of Accuracy is relatively wide, an adjustment is maded to select the best cp value and it can be checked in the graph below.

 1# lowest CP
 2# caret: maximum accuracy, rpart: 1-SE
 3# according to rpart, caret's best is lowest
 4df = mod.crt.eq$results
 5mod.crt.eq.lst.cp = mod.crt.eq$bestTune$cp
 6mod.crt.eq.lst.acc = subset(df,cp==mod.crt.eq.lst.cp)[[2]]
 7
 8# best cp by 1-SE rule - values are adjusted from graph 
 9mod.crt.eq.bst.cp = df[17,1]
10mod.crt.eq.bst.acc = df[17,2]
 1# CP by 1-SE rule
 2maxAcc = subset(df,Accuracy==max(Accuracy))[[2]]
 3stdAtMaxAcc = subset(df,subset=Accuracy==max(Accuracy))[[4]]
 4# max cp within 1 SE
 5maxCP = subset(df,Accuracy>=maxAcc-stdAtMaxAcc)[nrow(subset(df,Accuracy>=maxAcc-stdAtMaxAcc)),][[1]]
 6accAtMaxCP = subset(df,cp==maxCP)[[2]]
 7
 8# plot Accuracy vs cp
 9ubound = ifelse(maxAcc+stdAtMaxAcc>max(df$Accuracy),max(df$Accuracy),maxAcc+stdAtMaxAcc)
10lbound = ifelse(maxAcc-stdAtMaxAcc<min(df$Accuracy),min(df$Accuracy),maxAcc-stdAtMaxAcc)
11
12ggplot(data=df[1:nrow(df),], aes(x=cp,y=Accuracy)) + 
13  geom_line() + geom_point() +   
14  geom_abline(intercept=ubound,slope=0, color="purple") + 
15  geom_abline(intercept=lbound,slope=0, color="purple") + 
16  geom_point(aes(x=mod.crt.eq.bst.cp,y=mod.crt.eq.bst.acc),color="red",size=3) + 
17  geom_point(aes(x=mod.crt.eq.lst.cp,y=mod.crt.eq.lst.acc),color="blue",size=3)

Similar to above, 2 trees with the respective cp values are fit into the train and test data and the details are kept in mmce. Below is the update by fitting from the train data.

 1## performance on train data
 2# refit from rpart for best cp - cp by 1-SE not fitted by caret
 3# note no cross-validation necessary - xval=0 (default: xval=10)
 4set.seed(12357)
 5mod.crt.eq.bst = rpart(High ~ ., data=trainData, control=rpart.control(xval=0,cp=mod.crt.eq.bst.cp))
 6
 7# fit to train data - lowest from caret, best (1-SE) from rpart
 8mod.crt.eq.lst.ftd = predict(mod.crt.eq,newdata=trainData)
 9mod.crt.eq.bst.ftd = predict(mod.crt.eq.bst, type="class")
10
11# confusion matrix
12mod.crt.eq.lst.ftd.cm = table(actual=trainData$High,fitted=mod.crt.eq.lst.ftd)
13mod.crt.eq.lst.ftd.cm = updateCM(mod.crt.eq.lst.ftd.cm, type="Fitted")
14
15mod.crt.eq.bst.ftd.cm = table(actual=trainData$High,fitted=mod.crt.eq.bst.ftd)
16mod.crt.eq.bst.ftd.cm = updateCM(mod.crt.eq.bst.ftd.cm, type="Fitted")
17
18# misclassification error
19mod.crt.eq.lst.ftd.mmce = list(pkg="caret",isTest=FALSE,isBest=FALSE,isEq=TRUE
20                               ,cp=round(mod.crt.eq.lst.cp,4)
21                               ,mmce=mod.crt.eq.lst.ftd.cm[3,3])
22mmce[[length(mmce)+1]] = unlist(mod.crt.eq.lst.ftd.mmce)
23
24mod.crt.eq.bst.ftd.mmce = list(pkg="caret",isTest=FALSE,isBest=TRUE,isEq=TRUE
25                               ,cp=round(mod.crt.eq.bst.cp,4)
26                               ,mmce=mod.crt.eq.bst.ftd.cm[3,3])
27mmce[[length(mmce)+1]] = unlist(mod.crt.eq.bst.ftd.mmce)

Below is the update by fitting from the test data. The updated fitting details can be checked.

 1## performance of test data
 2# fit to test data - lowest from caret, best (1-SE) from rpart
 3mod.crt.eq.lst.ptd = predict(mod.crt.eq,newdata=testData)
 4mod.crt.eq.bst.ptd = predict(mod.crt.eq.bst, newdata=testData)
 5
 6# confusion matrix
 7mod.crt.eq.lst.ptd.cm = table(actual=testData$High, fitted=mod.rpt.eq.lst.ptd)
 8mod.crt.eq.lst.ptd.cm = updateCM(mod.crt.eq.lst.ptd.cm)
 9
10mod.crt.eq.bst.ptd.cm = table(actual=testData$High, fitted=mod.rpt.eq.bst.ptd)
11mod.crt.eq.bst.ptd.cm = updateCM(mod.crt.eq.bst.ptd.cm)
12
13# misclassification error
14mod.crt.eq.lst.ptd.mmce = list(pkg="caret",isTest=TRUE,isBest=FALSE,isEq=TRUE
15                               ,cp=round(mod.crt.eq.lst.cp,4)
16                               ,mmce=mod.crt.eq.lst.ptd.cm[3,3])
17mmce[[length(mmce)+1]] = unlist(mod.crt.eq.lst.ptd.mmce)
18
19mod.crt.eq.bst.ptd.mmce = list(pkg="caret",isTest=TRUE,isBest=TRUE,isEq=TRUE
20                               ,cp=round(mod.crt.eq.bst.cp,4)
21                               ,mmce=mod.crt.eq.bst.ptd.cm[3,3])
22mmce[[length(mmce)+1]] = unlist(mod.crt.eq.bst.ptd.mmce)
23
24ldply(mmce)
##     pkg isTest isBest isEq     cp mmce
## 1 rpart  FALSE  FALSE TRUE 0.0114 0.16
## 2 rpart  FALSE   TRUE TRUE 0.1061 0.29
## 3 rpart   TRUE  FALSE TRUE 0.0114 0.19
## 4 rpart   TRUE   TRUE TRUE 0.1061  0.3
## 5 caret  FALSE  FALSE TRUE 0.0156 0.16
## 6 caret  FALSE   TRUE TRUE 0.2488 0.29
## 7 caret   TRUE  FALSE TRUE 0.0156 0.19
## 8 caret   TRUE   TRUE TRUE 0.2488  0.3

Finally the mlr package is employed.

At first, a taks and learner are set up.

1### mlr
2## task and learner
3tsk.mlr.eq = makeClassifTask(data=trainData,target="High")
4lrn.mlr.eq = makeLearner("classif.rpart",par.vals=list(cp=0))

Then a grid of cp values is generated followed by tuning the parameter. Note that, as the tuning optimization path does not include a standard-error-like variable, only the best cp values are taken into consideration.

 1## tune parameter
 2cpGrid = function(index) {
 3  start=0
 4  end=0.3
 5  len=20
 6  inc = (end-start)/len
 7  grid = c(start)
 8  while(start < end) {
 9    start = start + inc
10    grid = c(grid,start)
11  } 
12  grid[index]
13}
14# create tune control grid
15ps.mlr.eq = makeParamSet(
16    makeNumericParam("cp",lower=1,upper=20,trafo=cpGrid)
17  )
18ctrl.mlr.eq = makeTuneControlGrid(resolution=c(cp=(20)))
19
20# tune cp
21rdesc.mlr.eq = makeResampleDesc("RepCV",reps=5,folds=10,stratify=TRUE)
22set.seed(12357)
23tune.mlr.eq = tuneParams(learner=lrn.mlr.eq
24                         ,task=tsk.mlr.eq
25                         ,resampling=rdesc.mlr.eq
26                         ,par.set=ps.mlr.eq
27                         ,control=ctrl.mlr.eq)
28# tuned cp
29tune.mlr.eq$x
## $cp
## [1] 0.015
1# optimization path
2# no standard error like values, only best cp is considered
3path.mlr.eq = as.data.frame(tune.mlr.eq$opt.path)
4path.mlr.eq = transform(path.mlr.eq,cp=cpGrid(cp))
5head(path.mlr.eq,3)
##      cp mmce.test.mean dob eol error.message exec.time
## 1 0.000      0.2623790   1  NA          <NA>     0.948
## 2 0.015      0.2555407   2  NA          <NA>     0.944
## 3 0.030      0.2629069   3  NA          <NA>     0.945

Using the best cp value, the learner is updated followed by training the model.

1# obtain fitted and predicted responses
2# update cp
3lrn.mlr.eq.bst = setHyperPars(lrn.mlr.eq, par.vals=tune.mlr.eq$x)
4
5# train model
6trn.mlr.eq.bst = train(lrn.mlr.eq.bst, tsk.mlr.eq)

Then the model is fit into the train and test data and the fitting details are updated in mmce. The overall fitting results can be checked below.

 1# fitted responses
 2mod.mlr.eq.bst.ftd = predict(trn.mlr.eq.bst, tsk.mlr.eq)$data
 3mod.mlr.eq.bst.ptd = predict(trn.mlr.eq.bst, newdata=testData)$data
 4
 5# confusion matrix
 6mod.mlr.eq.bst.ftd.cm = table(actual=mod.mlr.eq.bst.ftd$truth, fitted=mod.mlr.eq.bst.ftd$response)
 7mod.mlr.eq.bst.ftd.cm = updateCM(mod.mlr.eq.bst.ftd.cm)
 8
 9mod.mlr.eq.bst.ptd.cm = table(actual=mod.mlr.eq.bst.ptd$truth, fitted=mod.mlr.eq.bst.ptd$response)
10mod.mlr.eq.bst.ptd.cm = updateCM(mod.mlr.eq.bst.ptd.cm)
11
12# misclassification error
13mod.mlr.eq.bst.ftd.mmce = list(pkg="mlr",isTest=FALSE,isBest=FALSE,isEq=TRUE
14                               ,cp=round(tune.mlr.eq$x[[1]],4)
15                               ,mmce=mod.mlr.eq.bst.ftd.cm[3,3])
16mmce[[length(mmce)+1]] = unlist(mod.mlr.eq.bst.ftd.mmce)
17
18mod.mlr.eq.bst.ptd.mmce = list(pkg="mlr",isTest=TRUE,isBest=FALSE,isEq=TRUE
19                               ,cp=round(tune.mlr.eq$x[[1]],4)
20                               ,mmce=mod.mlr.eq.bst.ptd.cm[3,3])
21mmce[[length(mmce)+1]] = unlist(mod.mlr.eq.bst.ptd.mmce)
22
23ldply(mmce)
##      pkg isTest isBest isEq     cp mmce
## 1  rpart  FALSE  FALSE TRUE 0.0114 0.16
## 2  rpart  FALSE   TRUE TRUE 0.1061 0.29
## 3  rpart   TRUE  FALSE TRUE 0.0114 0.19
## 4  rpart   TRUE   TRUE TRUE 0.1061  0.3
## 5  caret  FALSE  FALSE TRUE 0.0156 0.16
## 6  caret  FALSE   TRUE TRUE 0.2488 0.29
## 7  caret   TRUE  FALSE TRUE 0.0156 0.19
## 8  caret   TRUE   TRUE TRUE 0.2488  0.3
## 9    mlr  FALSE  FALSE TRUE  0.015 0.16
## 10   mlr   TRUE  FALSE TRUE  0.015 0.19

It is shown that the mmce values are identical and it seems to be because the model is quite stable with respect to cp. It can be checked in the following graph.

 1# mmce vs cp
 2ftd.mmce = c()
 3ptd.mmce = c()
 4cps = mod.crt.eq$results[[1]]
 5for(i in 1:length(cps)) {
 6  if(i %% 2 == 0) {
 7    set.seed(12357)
 8    mod = rpart(High ~ ., data=trainData, control=rpart.control(cp=cps[i]))
 9    ftd = predict(mod,type="class")
10    ptd = predict(mod,newdata=testData,type="class")
11    ftd.mmce = c(ftd.mmce,updateCM(table(trainData$High,ftd))[3,3])
12    ptd.mmce = c(ptd.mmce,updateCM(table(testData$High,ptd))[3,3])    
13  }
14}
15mmce.crt = data.frame(cp=as.factor(round(cps[seq(2,length(cps),2)],3))
16                      ,fitted=ftd.mmce
17                      ,predicted=ptd.mmce)
18mmce.crt = melt(mmce.crt,id=c("cp"),variable.name="data",value.name="mmce")
19ggplot(data=mmce.crt,aes(x=cp,y=mmce,fill=data)) + 
20  geom_bar(stat="identity", position=position_dodge())

It may not be convicing to use a wrapper by this article about a single model. For example, however, if there are multiple models with a variety of tuning parameters to compare, the benefit of having one can be considerable. In the following articles, a similar approach would be taken, which is comparing individual packages to the wrappers.