Intro to Trees

Overview

Tree-based models basically consist of one or more nested if-then statements for the predictors that partition the data. Within these partitions, a specific model is used to predict the outcome. This recursive partitioning technique provides for exploration of the structure of a set of data (outcome and predictors) and identification of easy to visualize decision rules for predicting a categorical (Classification Tree) or continuous (Regression Tree) outcome.

In this tutorial we briefly describe the process of growing, examining, and pruning regression trees.

Outline

In this session we cover …

Introduction to Data (Boston Data)
Multivariate Regression Baseline
Regression Tree (CART method): rpart (rpart package)
Regression Tree (Conditional Inference method): ctree (partykit package)
Conclusion

Loading libraries used in this script.

library(MASS)               #the Boston Data
library(psych)              #data descriptives
library(ggplot2)            #data visualization
library(caret)              #training and cross validation, calls other model libraries
library(rattle)             #fancy tree plot
library(rpart)              #trees
library(rpart.plot)         #enhanced tree plots
library(RColorBrewer)       #color pallets
library(party)              #alternative decision tree algorithm
library(partykit)           #updated party functions
library(dplyr)              #data manipulation

1. Introduction to the Data

For this example we use data that accompanies the MASS package. No special reason these data were selected, other than they were used in some other examples we were working on. The data can be considered “typical” social science data, with a mix of nominal, count, and continuous variables. Of note, there are no missing data.

Reading in the Boston Data Exploration Data Set

#loading the data
data("Boston")

Data Descriptives

Lets have a quick look at the data file and the descriptives.

#data structure
head(Boston, 10)

##       crim   zn indus chas   nox    rm   age    dis rad tax ptratio  black
## 1  0.00632 18.0  2.31    0 0.538 6.575  65.2 4.0900   1 296    15.3 396.90
## 2  0.02731  0.0  7.07    0 0.469 6.421  78.9 4.9671   2 242    17.8 396.90
## 3  0.02729  0.0  7.07    0 0.469 7.185  61.1 4.9671   2 242    17.8 392.83
## 4  0.03237  0.0  2.18    0 0.458 6.998  45.8 6.0622   3 222    18.7 394.63
## 5  0.06905  0.0  2.18    0 0.458 7.147  54.2 6.0622   3 222    18.7 396.90
## 6  0.02985  0.0  2.18    0 0.458 6.430  58.7 6.0622   3 222    18.7 394.12
## 7  0.08829 12.5  7.87    0 0.524 6.012  66.6 5.5605   5 311    15.2 395.60
## 8  0.14455 12.5  7.87    0 0.524 6.172  96.1 5.9505   5 311    15.2 396.90
## 9  0.21124 12.5  7.87    0 0.524 5.631 100.0 6.0821   5 311    15.2 386.63
## 10 0.17004 12.5  7.87    0 0.524 6.004  85.9 6.5921   5 311    15.2 386.71
##    lstat medv
## 1   4.98 24.0
## 2   9.14 21.6
## 3   4.03 34.7
## 4   2.94 33.4
## 5   5.33 36.2
## 6   5.21 28.7
## 7  12.43 22.9
## 8  19.15 27.1
## 9  29.93 16.5
## 10 17.10 18.9

Our outcome of interest is medv: median value of owner-occupied homes in $1000s.

Note that there is no id variable. This is convenient for some tasks.

Descriptives

#sample descriptives
describe(Boston)

##         vars   n   mean     sd median trimmed    mad    min    max  range  skew
## crim       1 506   3.61   8.60   0.26    1.68   0.33   0.01  88.98  88.97  5.19
## zn         2 506  11.36  23.32   0.00    5.08   0.00   0.00 100.00 100.00  2.21
## indus      3 506  11.14   6.86   9.69   10.93   9.37   0.46  27.74  27.28  0.29
## chas       4 506   0.07   0.25   0.00    0.00   0.00   0.00   1.00   1.00  3.39
## nox        5 506   0.55   0.12   0.54    0.55   0.13   0.38   0.87   0.49  0.72
## rm         6 506   6.28   0.70   6.21    6.25   0.51   3.56   8.78   5.22  0.40
## age        7 506  68.57  28.15  77.50   71.20  28.98   2.90 100.00  97.10 -0.60
## dis        8 506   3.80   2.11   3.21    3.54   1.91   1.13  12.13  11.00  1.01
## rad        9 506   9.55   8.71   5.00    8.73   2.97   1.00  24.00  23.00  1.00
## tax       10 506 408.24 168.54 330.00  400.04 108.23 187.00 711.00 524.00  0.67
## ptratio   11 506  18.46   2.16  19.05   18.66   1.70  12.60  22.00   9.40 -0.80
## black     12 506 356.67  91.29 391.44  383.17   8.09   0.32 396.90 396.58 -2.87
## lstat     13 506  12.65   7.14  11.36   11.90   7.11   1.73  37.97  36.24  0.90
## medv      14 506  22.53   9.20  21.20   21.56   5.93   5.00  50.00  45.00  1.10
##         kurtosis   se
## crim       36.60 0.38
## zn          3.95 1.04
## indus      -1.24 0.30
## chas        9.48 0.01
## nox        -0.09 0.01
## rm          1.84 0.03
## age        -0.98 1.25
## dis         0.46 0.09
## rad        -0.88 0.39
## tax        -1.15 7.49
## ptratio    -0.30 0.10
## black       7.10 4.06
## lstat       0.46 0.32
## medv        1.45 0.41

#plots
pairs.panels(Boston)

#histogram of outcome
Boston %>%
  ggplot(aes(x=medv)) +
  geom_histogram(binwidth=1, boundary=.5, fill="white", color="black") + 
  labs(x = "Median Home Value")

#correlation matrix
round(cor(Boston), 2)

##          crim    zn indus  chas   nox    rm   age   dis   rad   tax ptratio
## crim     1.00 -0.20  0.41 -0.06  0.42 -0.22  0.35 -0.38  0.63  0.58    0.29
## zn      -0.20  1.00 -0.53 -0.04 -0.52  0.31 -0.57  0.66 -0.31 -0.31   -0.39
## indus    0.41 -0.53  1.00  0.06  0.76 -0.39  0.64 -0.71  0.60  0.72    0.38
## chas    -0.06 -0.04  0.06  1.00  0.09  0.09  0.09 -0.10 -0.01 -0.04   -0.12
## nox      0.42 -0.52  0.76  0.09  1.00 -0.30  0.73 -0.77  0.61  0.67    0.19
## rm      -0.22  0.31 -0.39  0.09 -0.30  1.00 -0.24  0.21 -0.21 -0.29   -0.36
## age      0.35 -0.57  0.64  0.09  0.73 -0.24  1.00 -0.75  0.46  0.51    0.26
## dis     -0.38  0.66 -0.71 -0.10 -0.77  0.21 -0.75  1.00 -0.49 -0.53   -0.23
## rad      0.63 -0.31  0.60 -0.01  0.61 -0.21  0.46 -0.49  1.00  0.91    0.46
## tax      0.58 -0.31  0.72 -0.04  0.67 -0.29  0.51 -0.53  0.91  1.00    0.46
## ptratio  0.29 -0.39  0.38 -0.12  0.19 -0.36  0.26 -0.23  0.46  0.46    1.00
## black   -0.39  0.18 -0.36  0.05 -0.38  0.13 -0.27  0.29 -0.44 -0.44   -0.18
## lstat    0.46 -0.41  0.60 -0.05  0.59 -0.61  0.60 -0.50  0.49  0.54    0.37
## medv    -0.39  0.36 -0.48  0.18 -0.43  0.70 -0.38  0.25 -0.38 -0.47   -0.51
##         black lstat  medv
## crim    -0.39  0.46 -0.39
## zn       0.18 -0.41  0.36
## indus   -0.36  0.60 -0.48
## chas     0.05 -0.05  0.18
## nox     -0.38  0.59 -0.43
## rm       0.13 -0.61  0.70
## age     -0.27  0.60 -0.38
## dis      0.29 -0.50  0.25
## rad     -0.44  0.49 -0.38
## tax     -0.44  0.54 -0.47
## ptratio -0.18  0.37 -0.51
## black    1.00 -0.37  0.33
## lstat   -0.37  1.00 -0.74
## medv     0.33 -0.74  1.00

Split Training and Test Data

For independent comparison of model predictions, we partition the data into a Training Set and an independent Test Set

#Setting the random seed for replication
set.seed(1234)

#renaming data set 
dat <- Boston

#Spliting training set into two parts based on outcome: 75% and 25%
index <- sample(1:nrow(dat), size=0.75*nrow(dat))
trainData <- dat[index,]
testData <- dat[-index,]

# #Using caret package function  
index <- createDataPartition(dat$medv, times=1, p=0.75, list=FALSE)
trainData <- dat[index,]
testData <- dat[-index,]

There are some nuanced distinctions between indexes created using the base sample() function and the caret package’s createDataPartition() function.

From the documentation for caret: For bootstrap samples, simple random sampling is used. For other data splitting, the random sampling is done within the levels of y when y is a factor in an attempt to balance the class distributions within the splits. For numeric y, the sample is split into groups sections based on percentiles and sampling is done within these subgroups. For createDataPartition, the number of percentiles is set via the groups argument. Also, for createDataPartition, very small class sizes (<= 3) the classes may not show up in both the training and test data.

Here, we proceed with the createDataPartition() version.

2. Regression - As A Preliminary Prediction Model

For “baseline”, lets run a regression, predicting medv from all other variables. This is also a classification model.

#Running exploratory linear regression
lm.fit <- lm(medv ~., data=trainData)
summary(lm.fit)

## 
## Call:
## lm(formula = medv ~ ., data = trainData)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -10.3410  -2.7428  -0.5283   1.7940  26.7312 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  34.190305   5.836646   5.858 1.04e-08 ***
## crim         -0.100053   0.034240  -2.922 0.003692 ** 
## zn            0.047498   0.015994   2.970 0.003177 ** 
## indus         0.042237   0.070894   0.596 0.551692    
## chas          1.516263   0.992904   1.527 0.127598    
## nox         -16.833976   4.473360  -3.763 0.000195 ***
## rm            4.181097   0.477051   8.764  < 2e-16 ***
## age           0.010570   0.015158   0.697 0.486034    
## dis          -1.314441   0.220885  -5.951 6.23e-09 ***
## rad           0.305016   0.078772   3.872 0.000128 ***
## tax          -0.013688   0.004561  -3.001 0.002877 ** 
## ptratio      -0.993075   0.148989  -6.665 9.70e-11 ***
## black         0.007941   0.003005   2.643 0.008575 ** 
## lstat        -0.547640   0.058684  -9.332  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 4.709 on 367 degrees of freedom
## Multiple R-squared:  0.7533, Adjusted R-squared:  0.7445 
## F-statistic: 86.18 on 13 and 367 DF,  p-value: < 2.2e-16

Fit of the regression is pretty good. $R^2 = 0.74$

Unfortunately, there do not seem to be any really good ways for visualizing these models (besides when there are only two predictors we can obtain a prediction plane in 3-d space).

Test of Prediction

However, we would like to assess on the Test Data. We look at the squared correlation between predicted scores and actual scores in the Test Data.

cor(predict(lm.fit, newdata=testData), testData$medv)^2

## [1] 0.6976338

Also pretty good!

3. Regression Tree (CART method) - As An Alternative Prediction Model

Traditional Classification and Regression Trees, as described by Brieman, Freidman, Olshen, and Stone (2017), can be generated through the rpart package. In the terminology of tree models, the data are recursively split into terminal nodes or leaves of the tree. To obtain a prediction for a new sample, we would follow the if-then statements defined by the tree using values of the new sample’s predictors until reaching a terminal node. The model formula in the terminal node would then be used to generate the prediction. In simple (traditional) trees, the model is a simple numeric value (yes/no, or a given numeric value). In other cases, the terminal node may be defined by a more complex function of the predictors (terminal nodes have models within them).

Tree-based and rule-based models are popular modeling tools for a number of reasons:

They generate a set of conditions that are highly interpretable and are easy to implement.
They can effectively handle many types of predictors (sparse, skewed, continuous, categorical, etc.) without the need for pre-processing.
These models do not require the user to specify the form of the predictors’ relationship to the response (e.g., linear, quadratic).
these models can (in some forms) effectively handle missing data and implicitly conduct feature selection. They have been extremely useful in many scenarios.

Basic implementation is done by Growing, Examining, Pruning - as illustrated below.

3a. Grow a Tree

To grow a traditional tree, we can use the rpart() function in the rpart package.

tree.fit <- rpart(formula, data=, method=,control=) where
+formula is in the format outcome ~ predictor1+predictor2+predictor3+etc.
+data= specifies the data frame +method= “class” for a classification tree; “anova” for a regression tree +control= optional parameters for controlling tree growth.

For example, control=rpart.control(minsplit=30, cp=0.001) requires that the minimum number of observations in a node be 30 before attempting a split and that a split must decrease the overall lack of fit by a factor of 0.001 (cost complexity factor) before being attempted.

rtree.fit <- rpart(medv ~ .,
                   data=trainData,
                   method="anova", #for regression tree
                   control=rpart.control(minsplit=30, cp=0.001))

3b. Examine the Tree

A collection of functions help us evaluate and examine the model.

+printcp(tree.fit) displays table of fits across cp (complexity parameter) values +rsq.rpart(tree.fit) plots approximate R-squared and relative error for different splits (2 plots). Labels are only appropriate for the “anova” method. +plotcp(tree.fit) plots the cross-validation results across cp values +print(tree.fit) print results +summary(tree.fit) detailed results including surrogate splits +plot(tree.fit) plot decision tree +text(tree.fit) label the decision tree plot +post(tree.fit, file=) create postscript plot of decision tree (there may be better ways to get good looking tree plots)

First we look at what the error looks like across the range of complexity parameters (depth of tree)

printcp(rtree.fit) # display the results

## 
## Regression tree:
## rpart(formula = medv ~ ., data = trainData, method = "anova", 
##     control = rpart.control(minsplit = 30, cp = 0.001))
## 
## Variables actually used in tree construction:
## [1] black   crim    dis     lstat   nox     ptratio rm     
## 
## Root node error: 32986/381 = 86.579
## 
## n= 381 
## 
##           CP nsplit rel error  xerror     xstd
## 1  0.5006424      0   1.00000 1.00801 0.094884
## 2  0.1646937      1   0.49936 0.57253 0.062125
## 3  0.0815862      2   0.33466 0.38689 0.050904
## 4  0.0302051      3   0.25308 0.29129 0.040776
## 5  0.0263779      4   0.22287 0.26741 0.041113
## 6  0.0152967      5   0.19649 0.24659 0.042500
## 7  0.0101745      6   0.18120 0.24487 0.042233
## 8  0.0085808      7   0.17102 0.23779 0.040761
## 9  0.0056861      8   0.16244 0.22322 0.038003
## 10 0.0054798      9   0.15676 0.22804 0.036624
## 11 0.0048701     10   0.15128 0.22529 0.036536
## 12 0.0041288     11   0.14641 0.22093 0.036165
## 13 0.0017826     12   0.14228 0.21553 0.033452
## 14 0.0016574     13   0.14050 0.21586 0.033446
## 15 0.0015984     14   0.13884 0.21636 0.033438
## 16 0.0011805     15   0.13724 0.21633 0.033427
## 17 0.0011707     16   0.13606 0.21606 0.033294
## 18 0.0010000     17   0.13489 0.21566 0.033295

rsq.rpart(rtree.fit) #produces 2 plots

## 
## Regression tree:
## rpart(formula = medv ~ ., data = trainData, method = "anova", 
##     control = rpart.control(minsplit = 30, cp = 0.001))
## 
## Variables actually used in tree construction:
## [1] black   crim    dis     lstat   nox     ptratio rm     
## 
## Root node error: 32986/381 = 86.579
## 
## n= 381 
## 
##           CP nsplit rel error  xerror     xstd
## 1  0.5006424      0   1.00000 1.00801 0.094884
## 2  0.1646937      1   0.49936 0.57253 0.062125
## 3  0.0815862      2   0.33466 0.38689 0.050904
## 4  0.0302051      3   0.25308 0.29129 0.040776
## 5  0.0263779      4   0.22287 0.26741 0.041113
## 6  0.0152967      5   0.19649 0.24659 0.042500
## 7  0.0101745      6   0.18120 0.24487 0.042233
## 8  0.0085808      7   0.17102 0.23779 0.040761
## 9  0.0056861      8   0.16244 0.22322 0.038003
## 10 0.0054798      9   0.15676 0.22804 0.036624
## 11 0.0048701     10   0.15128 0.22529 0.036536
## 12 0.0041288     11   0.14641 0.22093 0.036165
## 13 0.0017826     12   0.14228 0.21553 0.033452
## 14 0.0016574     13   0.14050 0.21586 0.033446
## 15 0.0015984     14   0.13884 0.21636 0.033438
## 16 0.0011805     15   0.13724 0.21633 0.033427
## 17 0.0011707     16   0.13606 0.21606 0.033294
## 18 0.0010000     17   0.13489 0.21566 0.033295

plotcp(rtree.fit) # visualize cross-validation results

#A good choice of cp for pruning is often the leftmost value for which the mean lies below the horizontal line

The detailed summary of the tree.

summary(rtree.fit) # detailed summary of splits

## Call:
## rpart(formula = medv ~ ., data = trainData, method = "anova", 
##     control = rpart.control(minsplit = 30, cp = 0.001))
##   n= 381 
## 
##             CP nsplit rel error    xerror       xstd
## 1  0.500642386      0 1.0000000 1.0080115 0.09488424
## 2  0.164693655      1 0.4993576 0.5725309 0.06212545
## 3  0.081586230      2 0.3346640 0.3868926 0.05090420
## 4  0.030205062      3 0.2530777 0.2912859 0.04077632
## 5  0.026377853      4 0.2228727 0.2674150 0.04111330
## 6  0.015296677      5 0.1964948 0.2465941 0.04249958
## 7  0.010174509      6 0.1811981 0.2448725 0.04223288
## 8  0.008580823      7 0.1710236 0.2377950 0.04076078
## 9  0.005686078      8 0.1624428 0.2232208 0.03800273
## 10 0.005479803      9 0.1567567 0.2280385 0.03662362
## 11 0.004870068     10 0.1512769 0.2252898 0.03653612
## 12 0.004128835     11 0.1464069 0.2209296 0.03616539
## 13 0.001782610     12 0.1422780 0.2155303 0.03345180
## 14 0.001657362     13 0.1404954 0.2158570 0.03344575
## 15 0.001598403     14 0.1388380 0.2163636 0.03343820
## 16 0.001180458     15 0.1372396 0.2163319 0.03342721
## 17 0.001170673     16 0.1360592 0.2160629 0.03329441
## 18 0.001000000     17 0.1348885 0.2156567 0.03329454
## 
## Variable importance
##      rm   lstat   indus     nox     dis     age ptratio     tax    crim      zn 
##      32      23       8       7       7       6       6       5       3       2 
##     rad   black 
##       1       1 
## 
## Node number 1: 381 observations,    complexity param=0.5006424
##   mean=22.61522, MSE=86.57861 
##   left son=2 (321 obs) right son=3 (60 obs)
##   Primary splits:
##       rm      < 6.941    to the left,  improve=0.5006424, (0 missing)
##       lstat   < 8.13     to the right, improve=0.4490711, (0 missing)
##       indus   < 6.66     to the right, improve=0.2605818, (0 missing)
##       ptratio < 19.9     to the right, improve=0.2529090, (0 missing)
##       nox     < 0.6615   to the right, improve=0.2508830, (0 missing)
##   Surrogate splits:
##       lstat   < 4.65     to the right, agree=0.892, adj=0.317, (0 split)
##       ptratio < 13.85    to the right, agree=0.871, adj=0.183, (0 split)
##       indus   < 3.985    to the right, agree=0.856, adj=0.083, (0 split)
##       zn      < 87.5     to the left,  agree=0.853, adj=0.067, (0 split)
## 
## Node number 2: 321 observations,    complexity param=0.1646937
##   mean=19.76885, MSE=36.77791 
##   left son=4 (132 obs) right son=5 (189 obs)
##   Primary splits:
##       lstat < 14.215   to the right, improve=0.4601722, (0 missing)
##       nox   < 0.6615   to the right, improve=0.3362854, (0 missing)
##       crim  < 5.84803  to the right, improve=0.3133585, (0 missing)
##       dis   < 2.2595   to the left,  improve=0.2862421, (0 missing)
##       indus < 16.57    to the right, improve=0.2617936, (0 missing)
##   Surrogate splits:
##       age   < 87.75    to the right, agree=0.832, adj=0.591, (0 split)
##       dis   < 2.5977   to the left,  agree=0.832, adj=0.591, (0 split)
##       nox   < 0.577    to the right, agree=0.816, adj=0.553, (0 split)
##       indus < 16.57    to the right, agree=0.804, adj=0.523, (0 split)
##       tax   < 431      to the right, agree=0.785, adj=0.477, (0 split)
## 
## Node number 3: 60 observations,    complexity param=0.08158623
##   mean=37.84333, MSE=77.77212 
##   left son=6 (36 obs) right son=7 (24 obs)
##   Primary splits:
##       rm      < 7.437    to the left,  improve=0.57673630, (0 missing)
##       lstat   < 4.68     to the right, improve=0.30000230, (0 missing)
##       ptratio < 18.15    to the right, improve=0.16210820, (0 missing)
##       black   < 392.79   to the right, improve=0.13260650, (0 missing)
##       rad     < 6        to the right, improve=0.06068443, (0 missing)
##   Surrogate splits:
##       lstat   < 4.015    to the right, agree=0.767, adj=0.417, (0 split)
##       indus   < 18.84    to the left,  agree=0.650, adj=0.125, (0 split)
##       ptratio < 14.75    to the right, agree=0.650, adj=0.125, (0 split)
##       black   < 389.885  to the right, agree=0.650, adj=0.125, (0 split)
##       zn      < 81.25    to the left,  agree=0.633, adj=0.083, (0 split)
## 
## Node number 4: 132 observations,    complexity param=0.03020506
##   mean=14.84621, MSE=18.109 
##   left son=8 (57 obs) right son=9 (75 obs)
##   Primary splits:
##       crim < 6.99237  to the right, improve=0.4168184, (0 missing)
##       nox  < 0.646    to the right, improve=0.3344787, (0 missing)
##       tax  < 567.5    to the right, improve=0.3042238, (0 missing)
##       rad  < 16       to the right, improve=0.2984618, (0 missing)
##       dis  < 2.1038   to the left,  improve=0.2646445, (0 missing)
##   Surrogate splits:
##       rad   < 16       to the right, agree=0.879, adj=0.719, (0 split)
##       tax   < 567.5    to the right, agree=0.864, adj=0.684, (0 split)
##       lstat < 19.23    to the right, agree=0.758, adj=0.439, (0 split)
##       nox   < 0.646    to the right, agree=0.750, adj=0.421, (0 split)
##       dis   < 2.2045   to the left,  agree=0.742, adj=0.404, (0 split)
## 
## Node number 5: 189 observations,    complexity param=0.02637785
##   mean=23.20688, MSE=21.07228 
##   left son=10 (175 obs) right son=11 (14 obs)
##   Primary splits:
##       lstat < 4.91     to the right, improve=0.21847500, (0 missing)
##       rm    < 6.543    to the left,  improve=0.18051610, (0 missing)
##       tax   < 223.5    to the right, improve=0.08739034, (0 missing)
##       dis   < 2.0274   to the right, improve=0.06715994, (0 missing)
##       nox   < 0.5125   to the right, improve=0.05887722, (0 missing)
## 
## Node number 6: 36 observations,    complexity param=0.01017451
##   mean=32.375, MSE=42.85743 
##   left son=12 (17 obs) right son=13 (19 obs)
##   Primary splits:
##       lstat < 5.495    to the right, improve=0.21753060, (0 missing)
##       nox   < 0.4885   to the right, improve=0.10696000, (0 missing)
##       indus < 5.405    to the right, improve=0.10611250, (0 missing)
##       zn    < 8.75     to the left,  improve=0.09600593, (0 missing)
##       rad   < 6        to the right, improve=0.09038041, (0 missing)
##   Surrogate splits:
##       dis  < 3.5466   to the left,  agree=0.806, adj=0.588, (0 split)
##       nox  < 0.44495  to the right, agree=0.778, adj=0.529, (0 split)
##       age  < 37.75    to the right, agree=0.778, adj=0.529, (0 split)
##       crim < 0.102345 to the right, agree=0.722, adj=0.412, (0 split)
##       zn   < 33.5     to the left,  agree=0.722, adj=0.412, (0 split)
## 
## Node number 7: 24 observations
##   mean=46.04583, MSE=18.00915 
## 
## Node number 8: 57 observations,    complexity param=0.005479803
##   mean=11.69474, MSE=13.73243 
##   left son=16 (46 obs) right son=17 (11 obs)
##   Primary splits:
##       nox   < 0.6695   to the right, improve=0.23092890, (0 missing)
##       lstat < 20.315   to the right, improve=0.17336840, (0 missing)
##       crim  < 15.57415 to the right, improve=0.16618560, (0 missing)
##       rm    < 5.9965   to the left,  improve=0.14893880, (0 missing)
##       dis   < 2.0643   to the left,  improve=0.09113234, (0 missing)
##   Surrogate splits:
##       rm    < 6.8285   to the left,  agree=0.842, adj=0.182, (0 split)
##       age   < 73.55    to the right, agree=0.842, adj=0.182, (0 split)
##       black < 5.165    to the right, agree=0.842, adj=0.182, (0 split)
##       lstat < 34.195   to the left,  agree=0.825, adj=0.091, (0 split)
## 
## Node number 9: 75 observations,    complexity param=0.004128835
##   mean=17.24133, MSE=8.150425 
##   left son=18 (37 obs) right son=19 (38 obs)
##   Primary splits:
##       nox   < 0.607    to the right, improve=0.2228033, (0 missing)
##       age   < 82.55    to the right, improve=0.2046128, (0 missing)
##       crim  < 0.64739  to the right, improve=0.1928180, (0 missing)
##       dis   < 1.9864   to the left,  improve=0.1853864, (0 missing)
##       indus < 16.01    to the right, improve=0.1618525, (0 missing)
##   Surrogate splits:
##       indus   < 16.01    to the right, agree=0.907, adj=0.811, (0 split)
##       tax     < 397      to the right, agree=0.893, adj=0.784, (0 split)
##       dis     < 2.38405  to the left,  agree=0.840, adj=0.676, (0 split)
##       crim    < 1.310635 to the right, agree=0.813, adj=0.622, (0 split)
##       ptratio < 19.9     to the right, agree=0.733, adj=0.459, (0 split)
## 
## Node number 10: 175 observations,    complexity param=0.01529668
##   mean=22.6, MSE=14.48983 
##   left son=20 (78 obs) right son=21 (97 obs)
##   Primary splits:
##       lstat < 9.95     to the right, improve=0.19899010, (0 missing)
##       rm    < 6.142    to the left,  improve=0.17168350, (0 missing)
##       tax   < 223.5    to the right, improve=0.10665180, (0 missing)
##       nox   < 0.5125   to the right, improve=0.06624460, (0 missing)
##       indus < 4.1      to the right, improve=0.06274371, (0 missing)
##   Surrogate splits:
##       nox   < 0.519    to the right, agree=0.749, adj=0.436, (0 split)
##       indus < 7.625    to the right, agree=0.709, adj=0.346, (0 split)
##       crim  < 0.16951  to the right, agree=0.703, adj=0.333, (0 split)
##       age   < 58.75    to the right, agree=0.703, adj=0.333, (0 split)
##       rm    < 6.03     to the left,  agree=0.697, adj=0.321, (0 split)
## 
## Node number 11: 14 observations
##   mean=30.79286, MSE=41.20209 
## 
## Node number 12: 17 observations
##   mean=29.14706, MSE=48.12249 
## 
## Node number 13: 19 observations
##   mean=35.26316, MSE=20.48233 
## 
## Node number 16: 46 observations,    complexity param=0.004870068
##   mean=10.82391, MSE=8.35095 
##   left son=32 (29 obs) right son=33 (17 obs)
##   Primary splits:
##       crim  < 10.753   to the right, improve=0.4181931, (0 missing)
##       lstat < 19.645   to the right, improve=0.3635462, (0 missing)
##       dis   < 1.944    to the left,  improve=0.2809175, (0 missing)
##       rm    < 5.9965   to the left,  improve=0.2518830, (0 missing)
##       nox   < 0.7065   to the left,  improve=0.2129638, (0 missing)
##   Surrogate splits:
##       lstat < 20.315   to the right, agree=0.761, adj=0.353, (0 split)
##       dis   < 1.944    to the left,  agree=0.717, adj=0.235, (0 split)
##       nox   < 0.7065   to the left,  agree=0.674, adj=0.118, (0 split)
##       rm    < 6.392    to the left,  agree=0.652, adj=0.059, (0 split)
##       black < 374.515  to the left,  agree=0.652, adj=0.059, (0 split)
## 
## Node number 17: 11 observations
##   mean=15.33636, MSE=19.80413 
## 
## Node number 18: 37 observations,    complexity param=0.001180458
##   mean=15.87568, MSE=4.988868 
##   left son=36 (10 obs) right son=37 (27 obs)
##   Primary splits:
##       lstat < 18.885   to the right, improve=0.21095130, (0 missing)
##       dis   < 1.74745  to the left,  improve=0.08317967, (0 missing)
##       black < 319.515  to the left,  improve=0.05845274, (0 missing)
##       rm    < 5.701    to the left,  improve=0.05692733, (0 missing)
##       crim  < 4.78225  to the left,  improve=0.04675386, (0 missing)
##   Surrogate splits:
##       dis < 1.61755  to the left,  agree=0.892, adj=0.6, (0 split)
##       rm  < 5.229    to the left,  agree=0.811, adj=0.3, (0 split)
##       age < 99.4     to the right, agree=0.757, adj=0.1, (0 split)
## 
## Node number 19: 38 observations,    complexity param=0.001598403
##   mean=18.57105, MSE=7.644688 
##   left son=38 (10 obs) right son=39 (28 obs)
##   Primary splits:
##       black   < 378.73   to the left,  improve=0.1815008, (0 missing)
##       crim    < 0.17127  to the right, improve=0.1660778, (0 missing)
##       rm      < 5.7825   to the left,  improve=0.1288988, (0 missing)
##       ptratio < 18.85    to the right, improve=0.1206865, (0 missing)
##       age     < 76.6     to the right, improve=0.1194402, (0 missing)
##   Surrogate splits:
##       crim    < 0.7954   to the right, agree=0.842, adj=0.4, (0 split)
##       dis     < 1.96865  to the left,  agree=0.789, adj=0.2, (0 split)
##       indus   < 16.01    to the right, agree=0.763, adj=0.1, (0 split)
##       ptratio < 20.95    to the right, agree=0.763, adj=0.1, (0 split)
##       lstat   < 14.62    to the left,  agree=0.763, adj=0.1, (0 split)
## 
## Node number 20: 78 observations,    complexity param=0.00178261
##   mean=20.70641, MSE=6.031626 
##   left son=40 (58 obs) right son=41 (20 obs)
##   Primary splits:
##       ptratio < 17.85    to the right, improve=0.12498650, (0 missing)
##       dis     < 5.58775  to the right, improve=0.11024300, (0 missing)
##       black   < 376.74   to the left,  improve=0.07808010, (0 missing)
##       age     < 69.15    to the right, improve=0.07721863, (0 missing)
##       tax     < 281.5    to the right, improve=0.07232533, (0 missing)
##   Surrogate splits:
##       indus < 18.84    to the left,  agree=0.769, adj=0.1, (0 split)
##       tax   < 208      to the right, agree=0.769, adj=0.1, (0 split)
## 
## Node number 21: 97 observations,    complexity param=0.008580823
##   mean=24.12268, MSE=16.08938 
##   left son=42 (29 obs) right son=43 (68 obs)
##   Primary splits:
##       rm    < 6.1245   to the left,  improve=0.18136500, (0 missing)
##       dis   < 3.2948   to the right, improve=0.13881060, (0 missing)
##       age   < 85.7     to the left,  improve=0.10962750, (0 missing)
##       indus < 13.375   to the left,  improve=0.08972878, (0 missing)
##       black < 371.25   to the right, improve=0.06721395, (0 missing)
##   Surrogate splits:
##       dis     < 9.20395  to the right, agree=0.732, adj=0.103, (0 split)
##       ptratio < 19.95    to the right, agree=0.732, adj=0.103, (0 split)
##       lstat   < 9.545    to the right, agree=0.732, adj=0.103, (0 split)
##       crim    < 0.625475 to the right, agree=0.722, adj=0.069, (0 split)
##       black   < 356.34   to the left,  agree=0.722, adj=0.069, (0 split)
## 
## Node number 32: 29 observations
##   mean=9.393103, MSE=4.998573 
## 
## Node number 33: 17 observations
##   mean=13.26471, MSE=4.619931 
## 
## Node number 36: 10 observations
##   mean=14.19, MSE=2.1389 
## 
## Node number 37: 27 observations
##   mean=16.5, MSE=4.602222 
## 
## Node number 38: 10 observations
##   mean=16.6, MSE=4.866 
## 
## Node number 39: 28 observations
##   mean=19.275, MSE=6.754018 
## 
## Node number 40: 58 observations,    complexity param=0.001170673
##   mean=20.19655, MSE=4.282747 
##   left son=80 (14 obs) right son=81 (44 obs)
##   Primary splits:
##       black < 380.125  to the left,  improve=0.15546080, (0 missing)
##       age   < 68.35    to the right, improve=0.11553240, (0 missing)
##       rm    < 5.9015   to the left,  improve=0.08101572, (0 missing)
##       dis   < 4.13665  to the right, improve=0.07350646, (0 missing)
##       nox   < 0.6055   to the left,  improve=0.06544715, (0 missing)
##   Surrogate splits:
##       nox  < 0.434    to the left,  agree=0.810, adj=0.214, (0 split)
##       rm   < 5.665    to the left,  agree=0.810, adj=0.214, (0 split)
##       age  < 98.95    to the right, agree=0.810, adj=0.214, (0 split)
##       dis  < 7.59015  to the right, agree=0.810, adj=0.214, (0 split)
##       crim < 9.94376  to the right, agree=0.793, adj=0.143, (0 split)
## 
## Node number 41: 20 observations
##   mean=22.185, MSE=8.163275 
## 
## Node number 42: 29 observations
##   mean=21.5069, MSE=3.737194 
## 
## Node number 43: 68 observations,    complexity param=0.005686078
##   mean=25.23824, MSE=17.19471 
##   left son=86 (54 obs) right son=87 (14 obs)
##   Primary splits:
##       dis   < 3.2948   to the right, improve=0.16041480, (0 missing)
##       lstat < 9.33     to the left,  improve=0.11064070, (0 missing)
##       age   < 74.15    to the left,  improve=0.08072486, (0 missing)
##       indus < 13.375   to the left,  improve=0.05832395, (0 missing)
##       crim  < 0.04738  to the left,  improve=0.05557239, (0 missing)
##   Surrogate splits:
##       age   < 68.45    to the left,  agree=0.897, adj=0.500, (0 split)
##       nox   < 0.5585   to the left,  agree=0.882, adj=0.429, (0 split)
##       crim  < 0.57207  to the left,  agree=0.853, adj=0.286, (0 split)
##       indus < 16.57    to the left,  agree=0.838, adj=0.214, (0 split)
##       rm    < 6.778    to the left,  agree=0.824, adj=0.143, (0 split)
## 
## Node number 80: 14 observations
##   mean=18.75, MSE=5.5925 
## 
## Node number 81: 44 observations
##   mean=20.65682, MSE=2.988363 
## 
## Node number 86: 54 observations,    complexity param=0.001657362
##   mean=24.39259, MSE=4.685501 
##   left son=172 (42 obs) right son=173 (12 obs)
##   Primary splits:
##       rm    < 6.603    to the left,  improve=0.21607430, (0 missing)
##       tax   < 332.5    to the right, improve=0.14006240, (0 missing)
##       lstat < 6.725    to the right, improve=0.08868318, (0 missing)
##       nox   < 0.491    to the right, improve=0.08245957, (0 missing)
##       dis   < 6.16455  to the right, improve=0.08188036, (0 missing)
##   Surrogate splits:
##       crim < 0.44489  to the left,  agree=0.815, adj=0.167, (0 split)
##       dis  < 3.54875  to the right, agree=0.815, adj=0.167, (0 split)
##       age  < 69.1     to the left,  agree=0.796, adj=0.083, (0 split)
##       rad  < 7.5      to the left,  agree=0.796, adj=0.083, (0 split)
## 
## Node number 87: 14 observations
##   mean=28.5, MSE=52.04714 
## 
## Node number 172: 42 observations
##   mean=23.85476, MSE=3.887239 
## 
## Node number 173: 12 observations
##   mean=26.275, MSE=2.923542

That is a lot of output, but here we can also look at the predictors used in the tree and their relative importance in the prediction. We see specifically that rm (average number of rooms per dwelling) and lstat (lower status of the population, percent) are driving much of the prediction.

This particular tree methodology can also handle missing data. When building the tree, missing data are ignored. For each split, a variety of alternatives (called surrogate splits) are evaluated. A surrogate split is one whose results are similar to the original split actually used in the tree. If a surrogate split approximates the original split well, it can be used when the predictor data associated with the original split are not available. In practice, several surrogate splits may be saved for any particular split in the tree.

Plotting the tree.

# plot tree (old schol way)
plot(rtree.fit, uniform=TRUE,
     main="Regression Tree for Median Home Value")
text(rtree.fit, use.n=TRUE, all=TRUE, cex=.8)

#create more attractive plot of tree 
#using prp() in the rpart.plot package
prp(rtree.fit)

#using Rattle package
fancyRpartPlot(rtree.fit)

We see the intuitive value of the tree method in the plot.

3c. Prune the Tree

Prune back the tree to avoid overfitting the data. Hastie et al. (2009) suggest selecting the tree size associated with the numerically smallest error. That is, the size of the tree is selected by examining the error using cross-validation, specifically the minimum of the xerror column (cross-validation error) printed by printcp( ).

Pruning is easily done using the function prune(fit, cp= ) by examining the cross-validated error results from printcp(), selecting the complexity parameter associated with minimum error, and placing it into the prune( ) function. Alternatively, this can be automated using tree.fit$cptable[which.min(tree.fit$cptable[,"xerror"]),"CP"].

# prune the tree based on minimim xerror
pruned.rtree.fit<- prune(rtree.fit, cp= rtree.fit$cptable[which.min(rtree.fit$cptable[,"xerror"]),"CP"])

# plot the pruned tree using prp() in the rpart.plot package 
prp(pruned.rtree.fit, main="Pruned Regression Tree for Median Home Value")

In this case the pruned tree is not that much smaller than the original tree.

There are, of course other approaches for pruning. Breiman et al. (1984) suggest using the cross-validation approach and applying a one-standard-error rule on the optimization criteria for identifying the simplest tree. That is, find the smallest tree that is within one standard error of the tree with smallest absolute error, which is the leftmost cp value for which the mean lies below the horizontal line placed 1 SE above the minmum of the curve by the minline in the plotcp() function.

# prune the tree based on 1 SE error 
pruned2.rtree.fit<- prune(rtree.fit, cp=.01)

# plot the pruned tree using prp() in the rpart.plot package
prp(pruned2.rtree.fit, main="Pruned Regression Tree for Median Home Value")

Test of Prediction

Finally, for comparison with the regression model, we examine the $R^2$ of the original and pruned trees. Note: The predictive value of the model would typically be established through cross-validation and test samples. We do the below only for didactic illustration.

#original tree
cor(predict(rtree.fit, newdata=testData),testData$medv)^2

## [1] 0.6660661

#pruned tree #1
cor(predict(pruned.rtree.fit, newdata=testData),testData$medv)^2

## [1] 0.6517376

#pruned tree #2
cor(predict(pruned2.rtree.fit, newdata=testData),testData$medv)^2

## [1] 0.6344976

We see here the trade off between “overfit” to training data and potential generalizability to new data. More formal evaluations would be done using cross-validation. But the smaller pruned tree is still doing pretty well (almost as well as the multiple regression).

4. Regression Tree (Conditional Inference Method) - As An Alternative Prediction Model

Traditional CART-based trees recursively perform univariate splits of the dependent variable based on values on a set of covariates. An information measures (such as the Gini coefficient) is used to select the current covariate. There is, however, a variable selection bias in the algorithms used in the traditional (rpart and related methods) algorithms. These approaches tend to select variables that have many possible splits or many missing values.

To overcome that bias, conditional inference trees were introduced. Unlike the other approaches, Conditional Inference Trees use a significance test procedure to select variables at each split. The significance test, or better: the multiple significance tests computed at each start of the algorithm (select covariate - choose split - recurse) are permutation tests that are used to obtain the the distribution of the test statistic under the null hypothesis (by calculating all possible values of the test statistic under rearrangements of the labels on the observed data points (see wikipedia).

More details can be found here https://stats.stackexchange.com/questions/12140/conditional-inference-trees-vs-traditional-decision-trees, and in the original paper here http://statmath.wu-wien.ac.at/~zeileis/papers/Hothorn+Hornik+Zeileis-2006.pdf.

The steps for implementation are largely the same: Grow, examine (and maybe prune).

4a. Grow a Tree

To grow a tree using the conditional inference method, we can use the party (party: A Laboratory for Recursive Partitioning) package or the updated package partykit (partykit: A Toolkit for Recursive Partytioning). This package provides nonparametric regression trees for nominal, ordinal, numeric, censored, and multivariate responses.

Specifically, regression or classification trees are obtained using the function +ctree(formula, data=, control=) where +formula is in the format outcome ~ predictor1+predictor2+predictor3+etc.
+data= specifies the data frame +control= optional parameters for controlling tree growth. For example, control=ctree_control(maxdepth=3) requires that the maximum depth of the tree is 3. The default maxdepth = Inf means that no restrictions are applied to tree size.

ctree.fit <- ctree(medv ~ ., 
                   data=trainData,
                   control=ctree_control(maxdepth=Inf))

4b. Examine the Tree

A collection of functions help us evaluate and examine the model.

+print(tree.fit) displays the details of the tree
+plot(tree.fit) plot decision tree

For our example, this is

print(ctree.fit) # display the results

## 
## Model formula:
## medv ~ crim + zn + indus + chas + nox + rm + age + dis + rad + 
##     tax + ptratio + black + lstat
## 
## Fitted party:
## [1] root
## |   [2] lstat <= 8.1
## |   |   [3] rm <= 7.42
## |   |   |   [4] lstat <= 4.86
## |   |   |   |   [5] crim <= 0.12744
## |   |   |   |   |   [6] rm <= 6.957: 29.650 (n = 10, err = 77.4)
## |   |   |   |   |   [7] rm > 6.957: 34.510 (n = 10, err = 60.3)
## |   |   |   |   [8] crim > 0.12744: 35.875 (n = 8, err = 761.0)
## |   |   |   [9] lstat > 4.86
## |   |   |   |   [10] rm <= 7.061: 24.686 (n = 63, err = 530.8)
## |   |   |   |   [11] rm > 7.061: 34.391 (n = 11, err = 28.0)
## |   |   [12] rm > 7.42: 46.046 (n = 24, err = 432.2)
## |   [13] lstat > 8.1
## |   |   [14] lstat <= 14.98
## |   |   |   [15] rm <= 6.715
## |   |   |   |   [16] lstat <= 9.81
## |   |   |   |   |   [17] crim <= 0.12816: 22.265 (n = 26, err = 350.0)
## |   |   |   |   |   [18] crim > 0.12816: 26.789 (n = 9, err = 655.2)
## |   |   |   |   [19] lstat > 9.81
## |   |   |   |   |   [20] black <= 376.73: 18.589 (n = 18, err = 181.4)
## |   |   |   |   |   [21] black > 376.73: 20.797 (n = 71, err = 388.3)
## |   |   |   [22] rm > 6.715: 29.214 (n = 7, err = 53.8)
## |   |   [23] lstat > 14.98
## |   |   |   [24] tax <= 437
## |   |   |   |   [25] crim <= 0.54452: 18.668 (n = 34, err = 252.2)
## |   |   |   |   [26] crim > 0.54452: 15.141 (n = 17, err = 45.9)
## |   |   |   [27] tax > 437
## |   |   |   |   [28] lstat <= 21.32: 14.177 (n = 44, err = 613.5)
## |   |   |   |   [29] lstat > 21.32: 10.262 (n = 29, err = 240.0)
## 
## Number of inner nodes:    14
## Number of terminal nodes: 15

plot(ctree.fit,
     main="Regression CTree for Median Home Value")

For comparison with the regression model, we examine the $R^2$ of the conditional inference tree. (Note: The predictive value of the model would typically be established through cross-validation across many test samples.)

#R-square conditional inference tree
cor(predict(ctree.fit, newdata=testData),testData$medv)^2

## [1] 0.6968245

Although the statistical approach ensures that the right-sized tree is grown without additional (post-)pruning or cross-validation, the depth of the tree here is rather large (6 levels and 13 terminal nodes), which of course makes interpretation more difficult than with less deep trees.

4c. Prune the Tree

Prune back the tree to avoid overfitting the data. This time we might simply prune for simplicity of plotting and interpretation. Pruning is done by regrowing with a different control parameter.

# regrow the tree with small depth, maxdepth = 3
pruned.ctree.fit<- ctree(medv ~ ., 
                         data=trainData,
                         control=ctree_control(maxdepth=3))

#examine pruned tree
print(pruned.ctree.fit) # display the results

## 
## Model formula:
## medv ~ crim + zn + indus + chas + nox + rm + age + dis + rad + 
##     tax + ptratio + black + lstat
## 
## Fitted party:
## [1] root
## |   [2] lstat <= 8.1
## |   |   [3] rm <= 7.42
## |   |   |   [4] lstat <= 4.86: 33.164 (n = 28, err = 1099.1)
## |   |   |   [5] lstat > 4.86: 26.128 (n = 74, err = 1440.8)
## |   |   [6] rm > 7.42: 46.046 (n = 24, err = 432.2)
## |   [7] lstat > 8.1
## |   |   [8] lstat <= 14.98
## |   |   |   [9] rm <= 6.715: 21.219 (n = 124, err = 2019.7)
## |   |   |   [10] rm > 6.715: 29.214 (n = 7, err = 53.8)
## |   |   [11] lstat > 14.98
## |   |   |   [12] tax <= 437: 17.492 (n = 51, err = 439.1)
## |   |   |   [13] tax > 437: 12.622 (n = 73, err = 1121.4)
## 
## Number of inner nodes:    6
## Number of terminal nodes: 7

plot(pruned.ctree.fit,
     main="(Pruned) Regression CTree for Median Home Value")

#R-square conditional inference tree
cor(predict(pruned.ctree.fit, newdata=testData),testData$medv)^2

## [1] 0.6153179

In this case the pruned tree provides an easier set of rules, but gives up prediction accuracy (in the hope for better generalization to other data).

5. Conclusion

In this session we walked through some very basics of implementing regression tree models. Classification trees operate in much the same way, just that the outcome is a nominal variable. While individual trees are not often used in practice much anymore, they provide a foundation for the forthcoming ensemble methods - where many trees are combined together. So, next we take a walk into the forest.

As always, thank you for playing!

Citations

B_Miner. (2021, November 22). Conditional inference trees vs traditional decision trees [Forum post]. Cross Validated. https://stats.stackexchange.com/q/12140

Breiman, L., Friedman, J. H., Olshen, R. A., & Stone, C. J. (1984). Classification and Regression Trees (Wadsworth International Group, Belmont, California, 1984).

Breiman, L., Friedman, J., Olshen, R. A., & Stone, C. J. (2017). Classification and Regression Trees. Chapman and Hall/CRC. https://doi.org/10.1201/9781315139470

Brownlee, J. (2016, February 7). How to Build an Ensemble Of Machine Learning Algorithms in R. MachineLearningMastery.Com. https://www.machinelearningmastery.com/machine-learning-ensembles-with-r/

Hastie, T., Tibshirani, R., & Friedman, J. (2009). The Elements of Statistical Learning. Springer. https://doi.org/10.1007/978-0-387-84858-7

Hothorn, T., Bühlmann, P., Dudoit, S., Molinaro, A., & Van Der Laan, M. J. (2006). Survival ensembles. Biostatistics, 7(3), 355–373. https://doi.org/10.1093/biostatistics/kxj011

Hothorn, T., Hornik, K., & Zeileis, A. (2006). Unbiased Recursive Partitioning: A Conditional Inference Framework. Journal of Computational and Graphical Statistics, 15(3), 651–674. https://doi.org/10.1198/106186006X133933

Hothorn, T., & Zeileis, A. (2015). partykit: A Modular Toolkit for Recursive Partytioning in R. Journal of Machine Learning Research, 16(118), 3905–3909.

Kuhn, M. (2008). Building Predictive Models in R Using the caret Package. Journal of Statistical Software, 28, 1–26. https://doi.org/10.18637/jss.v028.i05

Milborrow, S. (2024). Rpart.plot: Plot “rpart” Models: An Enhanced Version of “plot.rpart” (Version 3.1.2). https://CRAN.R-project.org/package=rpart.plot

Neuwirth, E. (2022). RColorBrewer: ColorBrewer Palettes (Version 1.1-3). https://CRAN.R-project.org/package=RColorBrewer

R Core Team. (2024). R: A Language and Environment for Statistical Computing. Foundation for Statistical Computing. https://www.R-project.org/

Revelle, W. (2024). psych: Procedures for Psychological, Psychometric, and Personality Research. Northwestern University. https://CRAN.R-project.org/package=psych

Robin, X., Turck, N., Hainard, A., Tiberti, N., Lisacek, F., Sanchez, J.-C., & Müller, M. (2011). pROC: An open-source package for R and S+ to analyze and compare ROC curves. BMC Bioinformatics, 12(1), 77. https://doi.org/10.1186/1471-2105-12-77

Strobl, C., Boulesteix, A.-L., Zeileis, A., & Hothorn, T. (2007). Bias in random forest variable importance measures: Illustrations, sources and a solution. BMC Bioinformatics, 8(1), 25. https://doi.org/10.1186/1471-2105-8-25

Therneau, T., & Atkinson, B. (2025). rpart: Recursive Partitioning and Regression Trees (Version 4.1.24) . https://CRAN.R-project.org/package=rpart

Venables, W. N., & Ripley, B. D. (2002). Modern Applied Statistics with S. Springer. https://doi.org/10.1007/978-0-387-21706-2

Wickham, H. (2016). ggplot2: Elegant Graphics for Data Analysis. Springer-Verlag. https://ggplot2.tidyverse.org/

Wickham, H., François, R., Henry, L., Müller, K., & Vaughan, D. (2023). dplyr: A Grammar of Data Manipulation (Version 1.1.4). https://CRAN.R-project.org/package=dplyr

Williams, G. (2011). Data Mining with {Rattle} and {R}: The art of excavating data for knowledge discovery. Springer. ttps://rd.springer.com/book/10.1007/978-1-4419-9890-3

Zeileis, A., Hothorn, T., & Hornik, K. (2008). Model-Based Recursive Partitioning. Journal of Computational and Graphical Statistics, 17(2), 492–514. https://doi.org/10.1198/106186008X319331