Machine Learning (ML) is about training an algorithm that can perform certain tasks. The general steps to get the trained model include:
The goal of this course is that you can answer the following questions:
Typically we can define roughly three types of ML tasks:
Supervised learning, you train algorithms to predict something (classes = classification or values = regression) from some other data (= features), and you provide it with correct examples of the execution of the task (called training data). A linear regression is an example of supervised learning. Given \(y = f(x)\) with \(x\) our input feature (e.g. precipitation), \(y\) our response (growth), and \(f\) an unknown function that maps \(x \rightarrow y\) . The goal of supervised learning is to train a ML algorithm to approximate \(f\) given observed \((x_i, y_i)\) pairs.
Unsupervised learning, on the other hand, is when you provide the features, but no examples of the correct execution of the task. Clustering techniques are examples of unsupervised learning. (In the example above, \(y\) would be unknown).
Reinforcement learning is a technique that mimics a game-like situation. The algorithm finds a solution through trial and error, receiving either rewards or penalties for each action. As in games, the goal is to maximize the rewards. We will talk more about this technique on the last day of the course.
For now, we will focus on the first two tasks, supervised and unsupervised learning (here a YouTube video explaining again the difference).
In ML, predictors (or the explaining variables) are often called features:
In supervised learning the response (y) and the features (x) are known:
In unsupervised learning, only the features are known:
In reinforcement learning an agent (ML model) is trained by interacting with an environment:
Have a look at the two textbooks on ML (Elements of statistical learning and introduction to statistical learning) in our further readings at the end of the GRIPS course - which of the following statements is true?In unsupervised learning, we want to identify patterns in data without having any examples (supervision) about what the correct patterns / classes are. As an example, consider the iris data set. Here, we have 150 observations of 4 floral traits:
iris = datasets::iris
colors = hcl.colors(3)
traits = as.matrix(iris[,1:4])
species = iris$Species
image(y = 1:4, x = 1:length(species) , z = traits,
ylab = "Floral trait", xlab = "Individual")
segments(50.5, 0, 50.5, 5, col = "black", lwd = 2)
segments(100.5, 0, 100.5, 5, col = "black", lwd = 2)
The observations are from 3 species and indeed those species tend to have different traits, meaning that the observations form 3 clusters.
pairs(traits, pch = as.integer(species), col = colors[as.integer(species)])![](A2-MachineLearningTasks_files/figure-html/chunk-chapter3-2"-1.png)
However, imagine we don’t know what species are, what is basically the situation in which people in the antique have been. The people just noted that some plants have different flowers than others, and decided to give them different names. This kind of process is what unsupervised learning does.
An example for an unsupervised learning algorithm is k-means clustering, one of the simplest and most popular unsupervised machine learning algorithms (see more on this in section “distance based algorithms”).
To start with the algorithm, you first have to specify the number of clusters (for our example the number of species). Each cluster has a centroid, which is the assumed or real location representing the center of the cluster (for our example this would be how an average plant of a specific species would look like). The algorithm starts by randomly putting centroids somewhere. Afterwards each data point is assigned to the respective cluster that raises the overall in-cluster sum of squares (variance) related to the distance to the centroid least of all. After the algorithm has placed all data points into a cluster the centroids get updated. By iterating this procedure until the assignment doesn’t change any longer, the algorithm can find the (locally) optimal centroids and the data points belonging to this cluster. Note that results might differ according to the initial positions of the centroids. Thus several (locally) optimal solutions might be found.
The “k” in K-means refers to the number of clusters and the ‘means’ refers to averaging the data-points to find the centroids.
A typical pipeline for using k-means clustering looks the same as for other algorithms. After having visualized the data, we fit a model, visualize the results and have a look at the performance by use of the confusion matrix. By setting a fixed seed, we can ensure that results are reproducible.
set.seed(123)
#Reminder: traits = as.matrix(iris[,1:4]).
kc = kmeans(traits, 3)
print(kc)K-means clustering with 3 clusters of sizes 50, 62, 38
Cluster means:
Sepal.Length Sepal.Width Petal.Length Petal.Width
1 5.006000 3.428000 1.462000 0.246000
2 5.901613 2.748387 4.393548 1.433871
3 6.850000 3.073684 5.742105 2.071053
Clustering vector:
[1] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
[38] 1 1 1 1 1 1 1 1 1 1 1 1 1 2 2 3 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2
[75] 2 2 2 3 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 3 2 3 3 3 3 2 3 3 3 3
[112] 3 3 2 2 3 3 3 3 2 3 2 3 2 3 3 2 2 3 3 3 3 3 2 3 3 3 3 2 3 3 3 2 3 3 3 2 3
[149] 3 2
Within cluster sum of squares by cluster:
[1] 15.15100 39.82097 23.87947
(between_SS / total_SS = 88.4 %)
Available components:
[1] "cluster" "centers" "totss" "withinss" "tot.withinss"
[6] "betweenss" "size" "iter" "ifault" Visualizing the results. Color codes true species identity, symbol shows cluster result.
plot(iris[c("Sepal.Length", "Sepal.Width")],
col = colors[as.integer(species)], pch = kc$cluster)
points(kc$centers[, c("Sepal.Length", "Sepal.Width")],
col = colors, pch = 1:3, cex = 3)
We see that there are are some discrepancies. Confusion matrix:
table(iris$Species, kc$cluster)
1 2 3
setosa 50 0 0
versicolor 0 48 2
virginica 0 14 36If you want to animate the clustering process, you could run
library(animation)
saveGIF(kmeans.ani(x = traits[,1:2], col = colors),
interval = 1, ani.width = 800, ani.height = 800)Elbow technique to determine the probably best suited number of clusters:
set.seed(123)
getSumSq = function(k){ kmeans(traits, k, nstart = 25)$tot.withinss }
#Perform algorithm for different cluster sizes and retrieve variance.
iris.kmeans1to10 = sapply(1:10, getSumSq)
plot(1:10, iris.kmeans1to10, type = "b", pch = 19, frame = FALSE,
xlab = "Number of clusters K",
ylab = "Total within-clusters sum of squares",
col = c("black", "red", rep("black", 8)))
Often, one is interested in sparse models. Furthermore, higher k than necessary tends to overfitting. At the kink in the picture, the sum of squares dropped enough and k is still low enough. But keep in mind, this is only a rule of thumb and might be wrong in some special cases.
Information criteria such as AIC or BIC can be also used to select the number of clusters and control complexity.
The two most prominent branches of supervised learning are regression and classification. The basic distinction between the two is that classification is about predicting a categorical variable, and regression is about predicting a continuous variable.
The random forest (RF) algorithm is possibly the most widely used machine learning algorithm and can be used for regression and classification. We will talk more about the algorithm later.
For the moment, we want to go through a typical workflow for a supervised regression: First, we visualize the data. Next, we fit the model and lastly we visualize the results. We will again use the iris data set that we used before. The goal is now to predict Sepal.Length based on the information about the other variables (including species).
Fitting the model:
library(randomForest)
set.seed(123)Sepal.Length is a numerical variable:
str(iris)'data.frame': 150 obs. of 5 variables:
$ Sepal.Length: num 5.1 4.9 4.7 4.6 5 5.4 4.6 5 4.4 4.9 ...
$ Sepal.Width : num 3.5 3 3.2 3.1 3.6 3.9 3.4 3.4 2.9 3.1 ...
$ Petal.Length: num 1.4 1.4 1.3 1.5 1.4 1.7 1.4 1.5 1.4 1.5 ...
$ Petal.Width : num 0.2 0.2 0.2 0.2 0.2 0.4 0.3 0.2 0.2 0.1 ...
$ Species : Factor w/ 3 levels "setosa","versicolor",..: 1 1 1 1 1 1 1 1 1 1 ...hist(iris$Sepal.Length)
The randomForest can be used similar to a linear regression model, we can specify the features using the formula syntax (~. means that all other variables should be used as features):
m1 = randomForest(Sepal.Length ~ ., data = iris) # ~.: Against all others.
print(m1)
Call:
randomForest(formula = Sepal.Length ~ ., data = iris)
Type of random forest: regression
Number of trees: 500
No. of variables tried at each split: 1
Mean of squared residuals: 0.1364625
% Var explained: 79.97In statistics we would use a linear regression model:
mLM = lm(Sepal.Length~., data = iris)As many other ML algorithms, the RF is not interpretable, so we don’t get coefficients that connect the variables to the response. But, at least we get the variable importance which is similar to an anova, telling us which variables were the most important ones:
varImpPlot(m1)
Our liner model would report linear effects, however, the lm cannot keep up with the flexibility of a random forest!
coef(mLM) (Intercept) Sepal.Width Petal.Length Petal.Width
2.1712663 0.4958889 0.8292439 -0.3151552
Speciesversicolor Speciesvirginica
-0.7235620 -1.0234978 And the finally, we can use the model to make predictions using the predict method:
plot(predict(m1), iris$Sepal.Length, xlab = "Predicted", ylab = "Observed")
abline(0, 1)
To understand the structure of a random forest in more detail, we can use a package from GitHub.
reprtree:::plot.getTree(m1, iris)
Here, one of the regression trees is shown.
With the random forest, we can also do classification. The steps are the same as for regression tasks, but we can additionally see how well it performed by looking at the confusion matrix. Each row of this matrix contains the instances in a predicted class and each column represents the instances in the actual class. Thus the diagonals are the correctly predicted classes and the off-diagonal elements are the falsely classified elements.
Species is a factor with three levels:
str(iris)'data.frame': 150 obs. of 5 variables:
$ Sepal.Length: num 5.1 4.9 4.7 4.6 5 5.4 4.6 5 4.4 4.9 ...
$ Sepal.Width : num 3.5 3 3.2 3.1 3.6 3.9 3.4 3.4 2.9 3.1 ...
$ Petal.Length: num 1.4 1.4 1.3 1.5 1.4 1.7 1.4 1.5 1.4 1.5 ...
$ Petal.Width : num 0.2 0.2 0.2 0.2 0.2 0.4 0.3 0.2 0.2 0.1 ...
$ Species : Factor w/ 3 levels "setosa","versicolor",..: 1 1 1 1 1 1 1 1 1 1 ...Fitting the model (syntax is the same as for the regression task):
set.seed(123)
library(randomForest)
m1 = randomForest(Species ~ ., data = iris)
print(m1)
Call:
randomForest(formula = Species ~ ., data = iris)
Type of random forest: classification
Number of trees: 500
No. of variables tried at each split: 2
OOB estimate of error rate: 4.67%
Confusion matrix:
setosa versicolor virginica class.error
setosa 50 0 0 0.00
versicolor 0 47 3 0.06
virginica 0 4 46 0.08varImpPlot(m1)
Predictions:
head(predict(m1)) 1 2 3 4 5 6
setosa setosa setosa setosa setosa setosa
Levels: setosa versicolor virginicaConfusion matrix:
table(predict(m1), as.integer(iris$Species))
1 2 3
setosa 50 0 0
versicolor 0 47 4
virginica 0 3 46Our model made a few errors.
Visualizing results ecologically:
plot(iris$Petal.Width, iris$Petal.Length, col = iris$Species, main = "Observed")
plot(iris$Petal.Width, iris$Petal.Length, col = predict(m1), main = "Predicted")
Visualizing one of the fitted models:
reprtree:::plot.getTree(m1, iris)
Confusion matrix:
| setosa | versicolor | virginica | |
|---|---|---|---|
| setosa | 50 | 0 | 0 |
| versicolor | 0 | 47 | 4 |
| virginica | 0 | 3 | 46 |
A demonstration with the iris dataset:
library(randomForest)
# scale your features if possible (some ML algorithms converge faster with scaled features)
iris_scaled = iris
iris_scaled[,1:4] = scale(iris_scaled[,1:4])
model = randomForest(Species~., data = iris_scaled)RandomForest is not based on a specific data generating model and thus we will not get effect estimates that tell us how the input features affect the response:
# no summary method available
print(model)
Call:
randomForest(formula = Species ~ ., data = iris_scaled)
Type of random forest: classification
Number of trees: 500
No. of variables tried at each split: 2
OOB estimate of error rate: 4%
Confusion matrix:
setosa versicolor virginica class.error
setosa 50 0 0 0.00
versicolor 0 47 3 0.06
virginica 0 3 47 0.06The confusion matrix explains where (for which species) the model makes wrong predictions / classifications on the OOB splits (OOB = out of bag). Each tree in the random forest is trained on a bootstrap of the data (bootstrap = sample with replacement from the original data, on average, each bootstrap will have 66% of the original data). Observations not used in a specific bootstrap are then used to validate the specific tree, bootstrap errors are at the end averaged for the n trees in the random forest.
While we don’t get effect estimates as in a lm, we get the variable importance which reports how important the specific predictors are:
varImpPlot(model)
Predictions
head(predict(model)) 1 2 3 4 5 6
setosa setosa setosa setosa setosa setosa
Levels: setosa versicolor virginicaThe model predicts the species class for each observation
Performance:
table(predict(model), as.integer(iris$Species))
1 2 3
setosa 50 0 0
versicolor 0 47 3
virginica 0 3 47Task:
Sepal.Length instead of Species (classification -> regression)Regression:
Random Forest automatically infers the type of the task, so we don’t have to change much:
model = randomForest(Sepal.Length~., data = iris_scaled)The OOB error is now “% Var explained” which is very similar to a \(R^2\):
print(model)
Call:
randomForest(formula = Sepal.Length ~ ., data = iris_scaled)
Type of random forest: regression
Number of trees: 500
No. of variables tried at each split: 1
Mean of squared residuals: 0.2031826
% Var explained: 79.55Plot observed vs predicted:
plot(iris_scaled$Sepal.Length, predict(model), xlim = c(-3, 3), ylim = c(-3, 3))
abline(a = c(0, 1))
Calculate \(R^2\):
cor(iris_scaled$Sepal.Length, predict(model))**2[1] 0.7987894