Regression Tree Introduction

We will use the regression tree (= decision tree) for predictions. The quality of models will be compared using the root-mean-squared-error (RMSE) and the mean-absolute-percentage-error (MAPE).

In the Statistical Learning tutorial we introduced the faithful dataset and predicted eruption times using linear regression. Now, we will do the prediction using the regression tree method.

We will begin by developing some basic understanding of the regression tree algorithm.

First we extract 10 records from the original dataset.

ns = 10  # intended number of sample records
set.seed(1) # fix random generator (for reproducibility)
I = sample(1:nrow(faithful), ns) # get the sample indices
D = faithful[I,] # extract the sample
D

x = D$waiting; y = D$eruptions; # for easy access
cat('average of y is ', mean(y))

Let us visualise the data.

# library(ggplot2)
g <- ggplot(aes(x=waiting,y=eruptions), data = D)+
  geom_point(size=4, color="brown")+
  geom_text(label=1:ns, nudge_x = .5, nudge_y = .05)+
  xlab("x = waiting time [min]")+ ylab("y=eruption time[min]")+
  theme(text = element_text(size=12))
g

How would you split the data to minimise MSE? Recall the mean squared error is: \[\text{MSE} = \frac{1}{n}\sum_{k=1}^n \left(y_k -\hat{y_k}\right)^2, \] where \(n\) are the number of observations, \(y_k\) is the observed value (response, output), \(\hat{y}_k = \hat{f}(x_k)\) is the modelled value (prediction), and \(x_k\) is the input.

Before using the MSE, try it intuitively.

Probably, you did the split in a way similar to:
\[x_s = \frac{x_1+x_5}{2} = \frac{63+73}{2}=68,\] i.e. xs = (x[1]+x[5])/2.

Now, we can use this to predict (approximate) the \(y\)-values on the left \(c_1\) and right side \(c_2\): \[c_1 = \frac{1}{5}(y_1+y_2+y_6+y_8+y_9) = 2.04\] That means the MSE on the left side is: \[ \text{MSE}_\text{left} = \frac{1}{5}( (y_1-c_1)^2+(y_2-c_1)^2+(y_6-c_1)^2+(y_8-c_1)^2+(y_9-c_1)^2) \]

\[c_2 = \frac{1}{5}(y_5+y_7+y_{10}+y_4+y_3) = 4.1\] Similarly, to before we compute the \(\text{MSE}_\text{right}\). Other splits will give a higher MSE.

Let us visualise this additional information.

c1 = mean(y[x<68]) # left region
c2 = mean(y[x>=68]) # right region
g2 <- g + geom_vline(xintercept = 68, size=2, color = "blue")+ 
  geom_text(x=68,y=3, label="68", color="red")+
  geom_segment(aes(x=c(45),xend=c(68), y=c(c1), yend=c(c1)), 
               size=2, color="red", linetype=4)+
  geom_segment(aes(x=c(68),xend=c(95), y=c(c2), yend=c(c2)), 
               size=2, color="red", linetype=4)+
  geom_hline(yintercept = 3.1, size=2, color = "red")
g2

#library(pacman); # package manager
#p_load(rpart, rattle, rpart.plot, RColorBrewer) # load packages
dt <- rpart( eruptions ~ waiting, D, minsplit=2, maxdepth=1)
fancyRpartPlot(dt, sub="decision stump")

dt <- rpart( eruptions ~ waiting, D, minsplit=2, maxdepth=2)
fancyRpartPlot(dt, sub="Three split nodes - a level 2 regression tree.")

c3 = mean(y[c(9,8,2,6)])
c4 = y[1]
c5 = mean(y[c(5,7,10,4)])
c6 = y[3]

g2 + 
  geom_segment(aes(x=c(45),xend=c(61), y=c(c3), yend=c(c3)), 
               size=1, color="green", linetype=4)+
  geom_segment(aes(x=c(61),xend=c(68), y=c(c4), yend=c(c4)), 
               size=1, color="green", linetype=4)+
  geom_vline(xintercept = 61, size=1, color = "blue", linetype=4)+
  geom_vline(xintercept = 87, size=1, color = "blue", linetype=4)+
  geom_segment(aes(x=c(68),xend=c(87), y=c(c5), yend=c(c5)), 
               size=1, color="green", linetype=4)+
  geom_segment(aes(x=c(87),xend=c(95), y=c(c6), yend=c(c6)), 
               size=1, color="green", linetype=4)+
  geom_hline(yintercept = 3.1, size=1, color = "red", linetype=4)

predict(dt, data.frame(waiting=c(50, 60, 80, 65, 100)))

# training & test data split
n  = nrow(faithful); ns = floor(.7*n);
set.seed(1); I = sample(1:n,ns);
train = faithful[I,]; test = faithful[-I,];

# regression tree
dt <- rpart(eruptions ~ waiting, train, minsplit=3, maxdepth=7)
# fancyRpartPlot(dt)

# Quality measures
mse  <- function(y,yh) {mean((y-yh)^2)}
rmse <- function(y,yh) {sqrt(mse(y,yh))}
mape <- function(y,yh) {mean(abs( (y-yh)/y ))}
derr <- function(y,yh) {sprintf('MSE = %.3f, RMSE = %.3f, MAPE = %.2f%%\n',
                                    mse(y,yh),rmse(y,yh),100*mape(y,yh))}
y  = test$eruptions
yh = predict(dt, newdata = test)
cat("Regression tree: ", derr(y,yh))

# compare with linear regression
lr <- lm(eruptions ~ waiting, train)
yh = predict(lr, newdata = test)
cat("linear regression: ", derr(y,yh))

Regression Tree for Wage Data

We will use the Wage dataset from the ISLR package. Please see the related book James et al. (2013) for more details.

# library(ISLR) #loads the Wage data (pre-loaded)
data(Wage); head(Wage, n=3L);
str(Wage)

There are 3000 observations and 11 variables. Most variables are factors. We observe that the data is given per year. The salary (wage) is thousands of dollars.

Does the wage depend on the person’s age?

a = ggplot(Wage, aes(age, wage))
a+geom_point(colour="gray")+geom_smooth(method="lm")+geom_smooth()
# plot(Wage$age,Wage$wage) # built-in plot (avoids ggplot)

Visually, it appears that there is a small correlation. We not two distinct groups one with continuously high salary and one that increases initially before dropping. The group distinction might be well modeled with decision trees. The nonlinearity in the second group may require a new polynomial feature to improve the model quality.

What is the correlation coefficient?

cor(Wage$age,Wage$wage)

The correlation value indicates that there is only a very week correlation. That means, to infer the salary from the age variable is risky.

We will now derive a few models and compare their test-accuracy.

Split the data into 70% training samples and 30% test samples. We will also remove the logwage (column 10) since this is the logarithmic wage (see also cor(Wage$wage, Wage$logwage)). Hence, defeating the purpose of predicting the wage.

n = nrow(Wage); ns = floor(.7*n);
set.seed(1) # random generator seed
I = sample(1:n, ns)
train = Wage[I,-10]
test = Wage[-I,-10]
cat("train: ",ns, " observations")

Note, once the best model is derived the entire data-set can be used as training data.

pm <- lm(wage ~ age, train)
predict(pm, data.frame(age=c(20,25,30)))

rmse <- function(y,yh) {sqrt(mean((y-yh)^2))}

pm <- lm(wage ~ age, train)
yh = predict(pm, newdata=test)
y  = test$wage
rmse(y,yh)

mape     <- function(y,yh) {mean(abs((y-yh)/y))} 

pm <- lm(wage ~ age, train)
yh = predict(pm, newdata=test)
y  = test$wage
mape(y,yh)

cor(Wage$age^2,Wage$wage)

pm <- lm(wage ~ age + I(age^2), train )
mape ( test$wage, predict(pm, newdata=test) )

dt <- rpart( wage ~ age, train)
yh = predict(dt, newdata = test)
y  = test$wage
mape(y, yh); # rmse(y, yh)
fancyRpartPlot(dt)

dt <- rpart( wage ~ ., train)
mape(train$wage, predict(dt)); # training error
mape(test$wage, predict(dt, newdata=test)); # test error 

printcp(dt)
plotcp(dt)

dtp <- prune(dt, cp=0.021)
mape(test$wage, predict(dtp, newdata=test)); # test error 
plot(dtp,margin=0.2); text(dtp, pretty=0, use.n=TRUE, all=TRUE, cex=.8)

dt.pp <- rpart( wage ~ ., train, control = rpart.control(minbucket=100))
mape(train$wage, predict(dt.pp)); # training error
mape(test$wage, predict(dt.pp, newdata=test)); # test error 
fancyRpartPlot(dt.pp)

Warm-up and Kaggle

What is the proportion of survivors?

T = table(titanic_train$Survived)
prop.table(T)

Answer: 38.4% survived.

Create a file that you can submit at Kaggle’s competition. Predict that there are no survivors.

# initialise such that there are no survivors
ns = rep(0,nrow(titanic_test)) #
# create submission file for Kaggle competition
submit <- data.frame(PassengerId = titanic_test$PassengerId, Survived = ns)
# write.csv(submit,file="no-survivors.csv", row.names = FALSE)
head(submit, n=3L)

The file can be submitted using this link on Kaggle.

Decision trees

Plot a simple decision tree using Survived ~ Sex.

# library(rpart)
dtp <- rpart(Survived ~ Sex, data=titanic_train, method="class" )
plot(dtp,margin=0.2); text(dtp, pretty=0, use.n=TRUE, all=TRUE, cex=.8)

This is pretty ugly. We want to have a better decision tree visualisation and will use another package.

# library(pacman); p_load(rattle, rpart.plot, RColorBrewer)
dt <- rpart(Survived ~ Sex, data=titanic_train, method="class" )
fancyRpartPlot(dt)

This is much better. What is the meaning? The green colours mean passengers have died and blue means they survived. The root node (1) states that 62% of the passengers died. Hence, 38% survived. This agrees with the table we created earlier. Now, there is the questions Sex = male (is the passenger male)? If the answer is yes then we follow the left branch, otherwise the right one. Going to the left we end up at the green node 2. That means, males have a 81% chance of dying; and they constitute 65% of the data. So, a tree that is very easy to interpret.

What is the training-accuracy of this tree?

yh = predict(dt,type="class")
y  = titanic_train$Survived
mean(y == yh)

Answer: 78.7%

Let us now get more factors into the tree. First we will convert several fields to factors. Furthermore, we will set the age na fields to its average. The Name column will be dropped.

train = titanic_train; test = titanic_test;
cols <- c("Pclass", "Sex", "Ticket", "Cabin", "Embarked")
train$Survived = as.factor(train$Survived)
train[cols] <- lapply(train[cols], factor)
test[cols]  <- lapply(test[cols], factor)
mv = mean(c(train$Age,test$Age), na.rm =TRUE)
train$Age[is.na(train$Age)] <- mv
test$Age [is.na(test$Age) ] <- mv

# remove columns , e.g. train$Name = NULL
train = subset( train, select = -c(PassengerId, Name, Ticket, Cabin))
test  = subset( test, select = -c(PassengerId, Name, Ticket, Cabin))
cat("train: ",nrow(train)," rows"); head(train)
cat("test: " ,nrow(test)," rows"); head(test)

Please note that this test set is for Kaggle submissions. That means it does not include the Survived field.

Is the new train-accuracy better using all factors?

dtc <- rpart(Survived ~ ., data=train, method="class" )
yh = predict(dtc,type="class"); y  = train$Survived
mean(y == yh)

The accuracy has improved.

Let’s visualise the new model.

fancyRpartPlot(dtc)

Unfortunately, over-fitting can lead to worse test-accuracy.

yh <- predict(dtc, test, type="class")
head(yh, n=3L)

printcp(dtc)
plotcp(dtc)

dtp <- prune(dtc, cp=0.027)
#plot(dtp,margin=0.2); text(dtp, pretty=0, use.n=TRUE, all=TRUE, cex=.8)
fancyRpartPlot(dtp)

Random Forest

Let us now use the random forest method. This method creates several decision trees by varying the number of features and amount of data used. These trees then have a majority vote of which group to choose for a given observation.

The random forest function cannot deal with factors that have more than 53 categories. So we had to remove Ticket and Cabine from the input data.

#library(randomForest)
set.seed(7)
rf <- randomForest(Survived ~ ., data = train)
summary(rf)

What is the training-accuracy of this model?

yh = predict(rf) 
mean( yh == train$Survived)

Decision Tree

Create a classification tree and determine it’s test-accuracy.

#library(rpart); library(pacman); p_load(rattle, rpart.plot, RColorBrewer)
dt <- rpart(Species ~ ., data=train, method="class" )
# fancyRpartPlot(dt) # not much use too complex
yh <- predict(dt, type="class", newdata = test)
mean(yh == test$Species)

Random Forest

Create a random forest and determine it’s test-accuracy.

# library(randomForest)
set.seed(1)
rf <- randomForest(Species ~ ., data=train, importance=TRUE)
yh <- predict(rf,type="class", newdata = test)
mean(yh == test$Species)

kNN

Create a kNN classifier and determine it’s test-accuracy.

# library (class)
knn.cf=knn (train[,1:4], test[,1:4], train$Species, k=3, prob=TRUE)
table(knn.cf , test$Species)
mean(knn.cf == test$Species)

k-means

The \(k\)-means algorithm is a popular heuristic, which finds \(k\) clusters by minimising the distance between observation and cluster centre.

data(iris); head(iris); summary(iris); str(iris)

ggplot(iris, aes(Petal.Length, Petal.Width, color = Species)) +   geom_point()+ theme(legend.position = c(.18,.73))
  unique(iris$Species)

set.seed(7)
cm <- kmeans(iris[, 3:4], 3, nstart = 30)
cm
table(cm$cluster, iris$Species)

cm$cluster <- as.factor(cm$cluster)
ggplot(iris, aes(Petal.Length, Petal.Width, color = cm$cluster)) + geom_point()+ theme(legend.position = "none")

Simple Call

Hierarchically cluster the iris petal data using complete linkage. First scale the data by mean and standard-deviation. Then determine all distances between all observations.

What distance measure was used by default?

# p_load(ggplot2, latex2exp) # for plots
# p_load(ggdendro, dendextend) # for dendograms
md = "complete"
hc = iris[, 3:4] %>% scale %>% dist %>% 
  hclust(method = md)
hc

Determine the top three clusters and compare them to the actual species.

I = cutree (hc , 3)
table(I, iris$Species)

We see that the setosa and virginica species were modelled perfectly. However, only 50% of the versicolor irises were identified correctly. The problem is that we used the complete linkage, which was inappropriate in this case.

Visualise the modelled clusters.

X = iris[, 3:4];
qplot(X[,1],X[,2],geom = "point",
        color=as.factor(I))+
    xlab(TeX("x_1"))+ylab(TeX("x_2"))+
    theme(legend.position = "none")

Dendogram

Let’s just use only a small sample (30 observations) to gain more insights about dendograms and their linkages. We begin by defining a function that prepares ourdendogram for plotting.

ppdend <- function(D, method="average", nb_clusters=3){
  return(D %>% scale %>% dist %>% 
    hclust(method = method) %>% as.dendrogram %>%
    set("branches_k_color", k=nb_clusters) %>% set("branches_lwd", .6) %>%
    set("labels_colors") %>% set("labels_cex", c(.6,.6)) %>% 
    set("leaves_pch", 19) %>% set("leaves_col", c("blue", "purple")))
}

set.seed(1); I = sample(1:150, 30);
X = iris[I, 3:4]; 

par(mfrow =c(2,3))
plot( ppdend(X,"single"),main="single")
plot( ppdend(X,"complete"),main="complete")
plot( ppdend(X,"average"),main="average")
plot( ppdend(X,"centroid"),main="centroid")
plot( ppdend(X,"ward.D"),main="ward.D")
plot( ppdend(X,"ward.D2"),main="ward.D2")

Let us have a look at a few more dendogram visualisations.

dend <- ppdend(X)
plot( dend ) # as before simple plot
ggplot(as.ggdend(dend)) # as ggplot
ggplot(as.ggdend(dend), labels = FALSE) + 
  scale_y_reverse(expand = c(0.2, 0)) +
  coord_polar(theta="x") # as ploar coordinate plot

Acknowledgment

This tutorial was created using RStudio, R, rmarkdown, and many other tools and libraries. The packages learnr and gradethis were particularly useful. I’m very grateful to Prof. Andy Field for sharing his disovr package, which allowed me to improve the style of this tutorial and get more familiar with learnr. Allison Horst wrote a very instructive blog “Teach R with learnr: a powerful tool for remote teaching”, which encouraged me to continue with learnr. By the way, I find her statistic illustrations amazing.