Introduction

This vignette visualizes classification results from rpart (CART), using tools from the package. The displays in this vignette are discussed in section 4 of Raymaekers and Rousseeuw (2021) (for the training data), and in section A.4 of the Supplementary Material (for the test data).

library(rpart)
library(classmap)

Titanic training data

We first load and inspect the data.

data("data_titanic")

traindata <- data_titanic[which(data_titanic$dataType == "train"), -13]

dim(traindata) 

## [1] 891  12

colnames(traindata)

##  [1] "PassengerId" "Pclass"      "Name"        "Sex"         "Age"        
##  [6] "SibSp"       "Parch"       "Ticket"      "Fare"        "Cabin"      
## [11] "Embarked"    "y"

# SibSp: number of siblings + spouses aboard
# Parch: number of parents + children aboard

str(traindata)

## 'data.frame':    891 obs. of  12 variables:
##  $ PassengerId: int  1 2 3 4 5 6 7 8 9 10 ...
##  $ Pclass     : int  3 1 3 1 3 3 1 3 3 2 ...
##  $ Name       : chr  "Braund, Mr. Owen Harris" "Cumings, Mrs. John Bradley (Florence Briggs Thayer)" "Heikkinen, Miss. Laina" "Futrelle, Mrs. Jacques Heath (Lily May Peel)" ...
##  $ Sex        : chr  "male" "female" "female" "female" ...
##  $ Age        : num  22 38 26 35 35 NA 54 2 27 14 ...
##  $ SibSp      : int  1 1 0 1 0 0 0 3 0 1 ...
##  $ Parch      : int  0 0 0 0 0 0 0 1 2 0 ...
##  $ Ticket     : chr  "A/5 21171" "PC 17599" "STON/O2. 3101282" "113803" ...
##  $ Fare       : num  7.25 71.28 7.92 53.1 8.05 ...
##  $ Cabin      : chr  "" "C85" "" "C123" ...
##  $ Embarked   : chr  "S" "C" "S" "S" ...
##  $ y          : Factor w/ 2 levels "casualty","survived": 1 2 2 2 1 1 1 1 2 2 ...

table(traindata$y)

## 
## casualty survived 
##      549      342

Now we want to predict the response y by CART (rpart). We set the seed as rpart is not deterministic. We also inspect the resulting classification tree.

set.seed(123) 
rpart.out <- rpart(y ~ Pclass + Sex + SibSp + 
                    Parch + Fare + Embarked,
                  data = traindata, method = 'class',
                  model = TRUE)

## plot the tree:
# pdf(file = "titanic_train_tree.pdf", width = 6, height = 6)
rpart.plot::rpart.plot(rpart.out, box.palette = "RdBu")

# dev.off()

We now inspect the variable importance, which will be used to calculate the variable weights in the farness computation.

rpart.out$variable.importance

##        Sex       Fare     Pclass      Parch   Embarked      SibSp 
## 124.426330  43.978667  31.163132  16.994228   9.889646   8.349703

Now we declare the types of the variables. This is used in the calculation of the Daisy distance, which in turn is needed for the farness computation. There are 5 nominal columns and one ordinal. The variables not listed here are interval-scaled by default.

mytype <- list(nominal = c("Name", "Sex", "Ticket", "Cabin", "Embarked"), ordratio = c("Pclass"))

Now we prepare for visualization and inspect the elements of the output.

x_train <- traindata[, -12]
y_train <- traindata[,  12]
vcrtrain <- vcr.rpart.train(x_train, y_train, rpart.out, mytype)
names(vcrtrain)

##  [1] "X"         "yint"      "y"         "levels"    "predint"   "pred"     
##  [7] "altint"    "altlab"    "PAC"       "figparams" "fig"       "farness"  
## [13] "ofarness"  "trainfit"

vcrtrain$predint[1:10] # prediction as integer

##  [1] 1 2 1 2 1 1 1 1 2 2

vcrtrain$pred[1:10]    # prediction as label

##  [1] "casualty" "survived" "casualty" "survived" "casualty" "casualty"
##  [7] "casualty" "casualty" "survived" "survived"

vcrtrain$altint[1:10]  # alternative label as integer

##  [1] 2 1 1 1 2 2 2 2 1 1

vcrtrain$altlab[1:10]  # alternative label

##  [1] "survived" "casualty" "casualty" "casualty" "survived" "survived"
##  [7] "survived" "survived" "casualty" "casualty"

# Probability of Alternative Class (PAC) of each object:
vcrtrain$PAC[1:3] 

## [1] 0.18890815 0.05294118 0.59459459

#
summary(vcrtrain$PAC)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
## 0.05294 0.18891 0.18891 0.27905 0.29630 0.94706

# f(i, g) is the distance from case i to class g:
vcrtrain$fig[1:3, ] # for the first 3 objects:

##           [,1]      [,2]
## [1,] 0.2621228 0.5079752
## [2,] 0.9505125 0.3728041
## [3,] 0.2303989 0.1644634

# The farness of an object i is the f(i, g) to its own class: 
vcrtrain$farness[1:3]

## [1] 0.2621228 0.3728041 0.1644634

#
summary(vcrtrain$farness)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##  0.0000  0.0000  0.2437  0.3286  0.6234  0.9613

# The "overall farness" of an object is defined as the 
# lowest f(i, g) it has to any class g (including its own):
summary(vcrtrain$ofarness)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##  0.0000  0.0000  0.1919  0.2906  0.5342  0.9578

sum(vcrtrain$ofarness > 0.99, na.rm = TRUE) 

## [1] 0

# No farness is considered outlying in these data.

confmat.vcr(vcrtrain) 

## 
## Confusion matrix:
##           predicted
## given      casualty survived
##   casualty      521       28
##   survived      130      212
## 
## The accuracy is 82.27%.

cols <- c("firebrick", "blue")

We first make the stacked plot, visualizing the confusion matrix:

stackedplot(vcrtrain, classCols = cols,
            main = "Stacked plot of rpart on Titanic training data")

Now we consider the silhouette plot

silplot(vcrtrain, classCols = cols, 
        main = "Silhouettes of rpart on Titanic training data")

##  classNumber classLabel classSize classAveSi
##            1   casualty       549       0.55
##            2   survived       342       0.27

# silplot.out <- silplot(vcrtrain, classCols = cols)
# ggsave("titanic_train_silhouettes.pdf", silplot.out,
#        width = 5, height = 5)

Now we construct the quasi residual plot of PAC vs. age for males only. The age variable is not in the model, but including it did not improve the classification. We see that the very young male passengers are often misclassified.

hist(x_train$Age)

# Quasi residual plot versus age, for males only:

# pdf("titanic_qrp_versus_age_males.pdf", width = 4.8, height = 4.8)
PAC <- vcrtrain$PAC[which(x_train$Sex == "male")]
feat <- x_train$Age[which(x_train$Sex == "male")]
qresplot(PAC, feat, xlab = "Age (years)", opacity = 0.5,
         main = "quasi residual plot for male passengers",
         plotLoess = TRUE)
text(x = 14, y = 0.60, "loess curve", col = "red", cex = 1)

# dev.off()

Now we construct the class maps. First of the casualties.

classmap(vcrtrain, "casualty", classCols = cols)

# classmap(vcrtrain, "casualty", classCols = cols, identify = TRUE)

# blue points top right: cases "a" and "b" in the paper
x_train[which(y_train == "casualty")[119], ]

##     PassengerId Pclass                       Name    Sex Age SibSp Parch
## 178         178      1 Isham, Miss. Ann Elizabeth female  50     0     0
##       Ticket    Fare Cabin Embarked
## 178 PC 17595 28.7125   C49        C

# Woman in 1st class, should have survived.
#
x_train[which(y_train == "casualty")[192], ]

##     PassengerId Pclass                         Name    Sex Age SibSp Parch
## 298         298      1 Allison, Miss. Helen Loraine female   2     1     2
##     Ticket   Fare   Cabin Embarked
## 298 113781 151.55 C22 C26        S

# Similar, but child.

# red point most to the right: case "c" in the paper
x_train[which(y_train == "casualty")[268], ]

##     PassengerId Pclass              Name  Sex Age SibSp Parch Ticket Fare
## 439         439      1 Fortune, Mr. Mark male  64     1     4  19950  263
##           Cabin Embarked
## 439 C23 C25 C27        S

Now the class map of the survivors.

classmap(vcrtrain, "survived", classCols = cols)

# classmap(vcrtrain, "survived", classCols = cols, identify = TRUE)

# red point with highest farness among highest PAC: case "d" in the paper
x_train[which(y_train == "survived")[c(14)], ] 

##    PassengerId Pclass                                                      Name
## 26          26      3 Asplund, Mrs. Carl Oscar (Selma Augusta Emilia Johansson)
##       Sex Age SibSp Parch Ticket    Fare Cabin Embarked
## 26 female  38     1     5 347077 31.3875              S

# near-coinciding points with highest farness: cases "e" and "f" in the paper
#
x_train[which(y_train == "survived")[265], ]

##     PassengerId Pclass                               Name  Sex Age SibSp Parch
## 680         680      1 Cardeza, Mr. Thomas Drake Martinez male  36     0     1
##       Ticket     Fare       Cabin Embarked
## 680 PC 17755 512.3292 B51 B53 B55        C

x_train[which(y_train == "survived")[287], ] 

##     PassengerId Pclass                   Name  Sex Age SibSp Parch   Ticket
## 738         738      1 Lesurer, Mr. Gustave J male  35     0     0 PC 17755
##         Fare Cabin Embarked
## 738 512.3292  B101        C

# man --> predicted as not survived. also paid the highest 
# fare in the data and has same ticket number as passenger 
# 265. Also embarked at the same place. It turns out that
# Gustave Lesueur was the man servant of the banker Thomas Cardeza:
  
# https://www.encyclopedia-titanica.org/titanic-survivor/thomas-cardeza.html
# https://www.encyclopedia-titanica.org/titanic-survivor/gustave-lesueur.html

# blue point bottom right: case "g" in the paper
x_train[which(y_train == "survived")[90], ]

##     PassengerId Pclass             Name    Sex Age SibSp Parch   Ticket
## 259         259      1 Ward, Miss. Anna female  35     0     0 PC 17755
##         Fare Cabin Embarked
## 259 512.3292              C

# # woman  + first class so predicted in survivors.
# # paid highest fare in whole dataset 

Titanic test data

In addition to the training data, we consider the titanic test data. First we load the data and inspect.

testdata <- data_titanic[which(data_titanic$dataType == "test"), -13]

dim(testdata)

## [1] 418  12

x_test <- testdata[, -12]
y_test <- testdata[, 12]
table(y_test)

## y_test
## casualty survived 
##      260      158

Now we prepare for visualization:

vcrtest <- vcr.rpart.newdata(x_test, y_test, vcrtrain)

confmat.vcr(vcrtest)

## 
## Confusion matrix:
##           predicted
## given      casualty survived
##   casualty      230       30
##   survived       64       94
## 
## The accuracy is 77.51%.

cols <- c("firebrick", "blue")

Now we can visualize, starting with the stacked plot:

stackedplot(vcrtest, classCols = cols,
            main = "Stacked plot of Titanic test data")

Now we construct the silhouette plot:

silplot(vcrtest, classCols = cols, 
        main = "Silhouettes of rpart on Titanic test data")

##  classNumber classLabel classSize classAveSi
##            1   casualty       260       0.47
##            2   survived       158       0.25

Finally, we make the class maps. First of the casualties, in which the misclassifications are all female (since all males are predicted as casualty by this model).

classmap(vcrtest, "casualty", classCols = cols) 

Now the class map of survivors:

classmap(vcrtest, "survived", classCols = cols)