Chapter 10 Supervised learning

10.1 Introduction

Copy link

Often, we faced with omics data, we have annotations for some of the samples (or features) we have acquired:

• Among the 100s of patient tumours that were assayed using RNA sequencing or microarrays, we know the grade of the tumour for about half. We want to predict the grade of the other half using the gene expression profiles.

• We have performed a spatial proteomics experiment such as in Christoforou et al. (⊕Christoforou et al. 2016Christoforou, A, C M Mulvey, L M Breckels, A Geladaki, T Hurrell, P C Hayward, T Naake, et al. 2016. “A Draft Map of the Mouse Pluripotent Stem Cell Spatial Proteome.” Nat Commun 7: 8992. https://doi.org/10.1038/ncomms9992.) (see section 8.7.2) and know the sub-cellular localisation of some proteins. We want to predict the sub-cellular localisation of the other proteins using the protein profiles.

In both of these examples, the quantitative data are the data that we want to use to predict properties about samples or features; these are called the predictors. The grade of the samples in the first example and the protein sub-cellular localisation in the second one are the labels that we have in some cases, and want to predict otherwise. We can thus split our data in two parts, depending whether we have labels, or whether we want to predict them. The former are called labelled, and the latter unlabelled.

Figure 10.1: In supervised learning, the data are split in labelled or unlabelled data. The same applies when some of the columns are labelled.

In the figure above, the labels represent categories that need to be inferred from the predictors. This class of supervised learning is called classification. Classification are also split into binary classification when there are only two classes, or multi-label problem when, like above, there are more than two. The latter is a generalisation of the binary task. When the annotations are continuous values, the situation is referred to as a regression problem.

library("rWSBIM1322")
data(giris2)
summary(giris2)

##      BRCA1           BRCA2            TP53             A1CF         GRADE    
##  Min.   :3.870   Min.   :1.888   Min.   :0.4919   Min.   :0.1000   A   : 50  
##  1st Qu.:5.115   1st Qu.:2.700   1st Qu.:1.6242   1st Qu.:0.3044   B   : 50  
##  Median :5.800   Median :3.000   Median :4.3276   Median :1.3327   C   : 50  
##  Mean   :5.853   Mean   :3.058   Mean   :3.7796   Mean   :1.2111   NA's:100  
##  3rd Qu.:6.400   3rd Qu.:3.400   3rd Qu.:5.1885   3rd Qu.:1.8000             
##  Max.   :8.116   Max.   :4.400   Max.   :7.1766   Max.   :2.5685

table(is.na(giris2$GRADE))

## 
## FALSE  TRUE 
##   150   100

pca <- prcomp(giris2[, -5], scale = TRUE)
cls <- as.character(giris2$GRADE)
cls[is.na(cls)] <- "unknown"
factoextra::fviz_pca_ind(pca, col.ind = cls)

10.2 K-nearest neighbours classifier

Copy link

In this chapter, we’ll use a simple, but useful classification algorithm, k-nearest neighbours (kNN) to classify the giris patients. We will use the knn function from the class package.

K-nearest neighbours works by directly measuring the (Euclidean) distance between observations and inferring the class of unlabelled data from the class of its nearest neighbours. In the figure below, the unlabelled instances 1 and 2 will be assigned classes A (blue) and B (red) as their closest neighbours are red and blue, respectively.

Figure 10.2: Schematic illustrating the k nearest neighbours algorithm.

Typically in machine learning, there are two clear steps, where one first trains a model and then uses the model to predict new outputs (class labels in this case). In kNN, these two steps are combined into a single function call to knn.

giris2_labelled <- giris2[!is.na(giris2$GRADE), ]
giris2_unlabelled <- giris2[is.na(giris2$GRADE), ]

The knn function takes, among others¹³, the following arguments:

1. the labelled predictors, that will be used to train the model,
2. the unlabelled predictors, on which the model will be applied (see below why this is called test),
3. the labels (the length of this vector must match the number of rows of the labelled predictors).

library("class")
knn_res <- knn(giris2_labelled[, -5],
               giris2_unlabelled[, -5],
               giris2_labelled[, 5])
class(knn_res)

## [1] "factor"

table(knn_res)

## knn_res
##  A  B  C 
## 31 38 31

► Solution

giris2_res <- giris2
giris2_res[is.na(giris2$GRADE), 5] <- knn_res
p1 <- factoextra::fviz_pca_ind(pca, col.ind = cls)
p2 <- factoextra::fviz_pca_ind(pca, col.ind = giris2_res$GRADE)
library("patchwork")
p1 / p2

Figure 10.3: Visualisation of the labelled and unlabelled data before (top) and after (bottom) classification.

## all pairwise distances
dd <- as.matrix(dist(giris2[, -5]))
## extract only those for sample 167
giris2$dist_167 <- dd[167, ]

## look at the 15 nearest neighbours
library("tidyverse")
vote <- giris2 %>%
    head(150) %>%
    arrange(dist_167) %>%
    head(15) %>%
    mutate(vote = cumsum(count = ifelse(GRADE == "C", 1, -1))) %>%
    mutate(grade = ifelse(vote <= 0, "B", "C"))
vote

##    BRCA1 BRCA2 TP53 A1CF GRADE  dist_167 vote grade
## 1    6.8   2.8  4.8  1.4     B 0.7239656   -1     B
## 2    6.7   3.0  5.0  1.7     B 0.7922045   -2     B
## 3    6.9   3.1  4.9  1.5     B 0.8065525   -3     B
## 4    6.7   2.5  5.8  1.8     C 0.8444024   -2     B
## 5    7.2   3.0  5.8  1.6     C 0.8592955   -1     B
## 6    6.9   3.1  5.4  2.1     C 0.8892644    0     B
## 7    6.8   3.0  5.5  2.1     C 0.8956019    1     C
## 8    6.4   2.7  5.3  1.9     C 0.9116501    2     C
## 9    7.0   3.2  4.7  1.4     B 0.9623517    1     C
## 10   6.9   3.1  5.1  2.3     C 0.9723507    2     C
## 11   6.3   2.5  5.0  1.9     C 0.9726503    3     C
## 12   6.5   3.0  5.2  2.0     C 0.9802555    4     C
##  [ reached 'max' / getOption("max.print") -- omitted 3 rows ]

vote %>%
    ggplot(aes(x = 1:15, y = vote)) +
    xlab("Number of nearest neighbours") +
    geom_point(aes(colour = grade), size = 3) + geom_line() +
    geom_hline(yintercept = 0)

As we have seen, the number of nearest neighbours k has an important influence on the classification results. We can now refine our understanding of the knn function; it has the following arguments:

1. the labelled predictors, that will be used to train the model,
2. the unlabelled predictors, on which the model will be applied (see below why this is called test),
3. the labels (the length of this vector must match the number of rows of the labelled predictors),
4. the number of neighbours to use,
5. when set to TRUE, the prob argument also return the proportion of votes in favour of the assigned class.

knn_1 <- knn(giris2_labelled[, -5], giris2_unlabelled[, -5], giris2_labelled[, 5], k = 1)
knn_11 <- knn(giris2_labelled[, -5], giris2_unlabelled[, -5], giris2_labelled[, 5], k = 11)
table(knn_1, knn_11)

##      knn_11
## knn_1  A  B  C
##     A 31  0  0
##     B  0 34  4
##     C  0  3 28

10.3 Model selection

Copy link

There is no way to decide which set of results above, using k=1, 2, … 11 or any other value for k is better. At least not while proceeding as above. To be able to make an informed decision about the model parameter k, we need to to measure the performance of the kNN classifier for different values of k. To do so, we are going to create training and testing sets using the labelled data. Each of these will be composed by a certain proportion of the original labelled data.

Below, we denote the training predictors \(X_{tr}\) and labels \(Y_{tr}\) and the testing predictors \(X_{te}\) and labels \(Y_{te}\).

Figure 10.4: Training and testing sets.

We are now going to do the following procedure

1. hide the testing labels \(Y_{te}\),
2. train a classifier model using \(X_{tr}\) and \(Y_{tr}\),
3. apply it on \(X_{te}\) to obtain a new \(Y_{te}^{predicted}\), and
4. compare \(Y_{te}\) to \(Y_{te}^{predicted}\).

There are numerous different ways to measure the performance of a classifier using \(Y_{te}\) and \(Y_{te}^{predicted}\). We are going to focus on the classification accuracy, counting the proportion of correct results.

set.seed(1)
i <- sample(150, 100)
giris2_train <- giris2_labelled[i, ]
giris2_test <- giris2_labelled[-i, ]

test_1 <- knn(giris2_train[, -5], giris2_test[, -5], giris2_train[, 5], k = 1)
test_11 <- knn(giris2_train[, -5], giris2_test[, -5], giris2_train[, 5], k = 11)
table(giris2_test[, 5], test_1)

##    test_1
##      A  B  C
##   A 16  0  0
##   B  0 18  1
##   C  0  1 14

table(giris2_test[, 5], test_11)

##    test_11
##      A  B  C
##   A 16  0  0
##   B  0 19  0
##   C  0  1 14

sum(test_1 == giris2_test[, 5])/length(test_1)

## [1] 0.96

sum(test_11 == giris2_test[, 5])/length(test_11)

## [1] 0.98

There are now two adjustments that we can do improve on the procedure above:

1. Test all odd value of k from 1 to 11 (or higher, if deemed useful), to have a better granularity of the model parameter we test.

2. Testing more than on training/testing split. Indeed, relying on a single random split, we rely on a random configuration, that could affect that results either overly optimistically, or negatively. We prefer to repeat the split a certain number of time and calculate an average performance over all splits.

10.4 k-fold cross-validation

Copy link

There exist several principled ways to split data into training and testing partitions. Here, we are going to focus on k-fold cross-validation, where data of size \(n\) a repeatedly split into a training set of size around \(\frac{n (k - 1)}{k}\) and a testing set of size around \(\frac{n}{k}\).

In practice, the data are split into k partitions of size \(\frac{n}{k}\), and the training/testing procedure is repeated \(k\) times using \(k - 1\) partition as training data and 1 partition as testing data.

The figure below illustrates the cross validation procedure, creating 3 folds. One would typically do a 10-fold cross validation (if the size of the data permits it). We split the data into 3 random and complementary folds, so that each data point appears exactly once in each fold. This leads to a total test set size that is identical to the size of the full annotated dataset but is composed of out-of-sample predictions (i.e. a sample is never used for training and testing).

Figure 10.5: The data is split into 3 folds. At each training/testing iteration, a different fold is used as test partition (white) and the two other ones (blue) are used to train the model.

After cross-validation, all models used within each fold are discarded, and a new model is built using the whole labelled dataset, with the best model parameter(s), i.e those that generalised over all folds. This makes cross-validation quite time consuming, as it takes x+1 (where x in the number of cross-validation folds) times as long as fitting a single model, but is essential.

folds <- caret::createFolds(giris2_labelled[[5]], k = 10)
folds

## $Fold01
##  [1]   2   9  11  15  40  51  54  62  76  95 103 107 115 140 143
## 
## $Fold02
##  [1]   1  20  21  41  48  63  66  71  78  80 104 111 135 136 142
## 
## $Fold03
##  [1]   3  31  34  39  46  74  89  90  94  96 105 117 120 131 139
## 
## $Fold04
##  [1]  17  18  27  28  49  59  64  75  97  99 112 121 134 141 144
## 
## $Fold05
##  [1]   8  29  35  37  42  53  69  79  81  82 108 110 119 123 125
## 
## $Fold06
##  [1]  13  19  22  24  50  55  68  72  91 100 116 122 127 132 137
## 
## $Fold07
##  [1]   6  10  33  38  43  52  56  65  73  87 101 114 124 126 130
## 
## $Fold08
##  [1]   5  12  23  45  47  57  77  84  85  88 106 109 146 149 150
## 
## $Fold09
##  [1]   4  14  16  26  44  58  61  92  93  98 102 118 129 138 147
## 
## $Fold10
##  [1]   7  25  30  32  36  60  67  70  83  86 113 128 133 145 148

all(sort(unlist(folds)) == 1:150)

## [1] TRUE

knn_accuracy <- function(fld, k) {
    ## Uses global variables!
    train <- giris2_labelled[-fld, ]
    test <- giris2_labelled[fld, ]
    test_k <- knn(train[, -5],
                  test[, -5],
                  train[, 5],
                  k)
    sum(test_k == test[, 5])/length(test_k)
}

accs <- sapply(folds, knn_accuracy, k = 11)
accs

##    Fold01    Fold02    Fold03    Fold04    Fold05    Fold06    Fold07    Fold08 
## 0.9333333 0.9333333 0.9333333 1.0000000 1.0000000 1.0000000 1.0000000 0.9333333 
##    Fold09    Fold10 
## 1.0000000 1.0000000

mean(accs)

## [1] 0.9733333

► Solution

accs <- tibble(k = seq(1, 30, 2)) %>%
    rowwise() %>%
    mutate(acc = mean(sapply(folds, knn_accuracy, k))) %>%
    ungroup()
accs

## # A tibble: 15 × 2
##        k   acc
##    <dbl> <dbl>
##  1     1 0.953
##  2     3 0.96 
##  3     5 0.967
##  4     7 0.98 
##  5     9 0.96 
##  6    11 0.98 
##  7    13 0.98 
##  8    15 0.973
##  9    17 0.973
## 10    19 0.98 
## 11    21 0.98 
## 12    23 0.98 
## 13    25 0.96 
## 14    27 0.967
## 15    29 0.953

accs %>%
    ggplot(aes(x = k, y = acc)) +
    geom_point() +
    geom_line()

Figure 10.6: Average classification accuracy as a function of the parameter k.

As we can see from the last exercise, most values of k give a good classification accuracy, ranging from 96% to 98%, and it would be difficult to identify a single best value for k. Indeed, the three groups we are typing to assign new data to are relatively well separated and the little uncertainty that we observe is most likely due to some instances that lie at the boundary between classes B and C.

There exist a variety of packages (for example caret, tidymodels or mlr) that automate model and parameter selection using cross-validation procedures as illustrated above, so that in practice, these tasks can be automated. This becomes particularly useful when more than one model parameter (called hyper-parameters) need to tuned, different classifiers need to be assessed, and/or several model performance metrics are to be computed.

10.5 Variance and bias trade-off

Copy link

A recurrent goal in machine learning is to find a trade-off between variance and bias when fitting models. It is important to build a good model using the data at hand, i.e. that learns from the data, but not too much. If the model is too specific for the data that was used to build it (the model would be said to be over fitting the data) and would not be applicable to new data (the model would have high bias). On the other hand, a model that is too general, that wouldn’t be applicable enough to the data, or type of data at hand¹⁴, it would perform poorly, and hence lead to high variance.

This variance and bias trade-off can be illustrated in different ways. The figure below (reproduced with permission from (⊕James et al. 2014James, Gareth, Daniela Witten, Trevor Hastie, and Robert Tibshirani. 2014. An Introduction to Statistical Learning: With Applications in r. Springer Publishing Company, Incorporated.)) shows the prediction error as a function of the model complexity. In our example above, we have used accuracy, a metric that we want to optimise, while in the prediction error is a measure to be minimised.

Figure 10.7: Effect of model complexity on the prediction error (lower is better) on the training and testing data. The former decreases for lower values of \(k\) while the test error reaches a minimum around \(k = 10\) before increasing again. Reproduced with permission from James et al. 2014)

The model complexity represent the ability of the model (or model parameter) to learn/adapt to the training data. In our case, the complexity would be illustrated by\(\frac{1}{k}\): when k=1, the model becomes very complex as it adapts to every single neighbour. The other extreme, when k becomes large, the classification will tend to using an average of so many neighbours that doesn’t reflect any specificity of our data.

Figure 10.8: The kNN classifier using k=1 (left, solid classification boundaries) and k=100 (right, solid classification boundaries) compared the Bayes decision boundaries (see original material for details). Reproduced with permission from James et al. 2014).

This last figure, reproduced from (⊕Holmes and Huber 2019Holmes, Susan, and Wolfgang Huber. 2019. Modern Statistics for Modern Biology. Cambridge Univeristy Press.), frames the variance-bias trade-off as one between accuracy and precision. An over-fit model with high bias and low variance is one that is precise but not accurate: it could work extremely well on the training data but miss on new data, thus lacking in generalisation power. Conversely, an under-fit model could be accurate but with low precision: on average, it works well, but is unable to provide a precise answer. Ideally, we want models that achieve good precision while still being applicable and accurate with new data.

Figure 10.9: Precision and accuracy when shooting.

10.6 Additional exercises

Copy link

Load the hyperLOPIT2015_se data from the pRolocdata package. See section 8.7.2) for details. The features variable column markers defines the proteins for which the sub-cellular localisation is known - these are the labelled data. Those that are marked unknown are of unknown location - these are the unlabelled data.

cls <- c("#E41A1C", "#377EB8", "#238B45", "#FF7F00", "#FFD700",
         "#333333", "#00CED1", "#A65628", "#F781BF", "#984EA3",
         "#9ACD32", "#B0C4DE", "#00008A", "#8B795E", "#E0E0E060")