List of Figures

Chapter 1. Introduction to machine learning

Figure 1.1. An algorithm for making and serving a cake. We start at the top and, after performing each operation, follow the next arrow. Diamonds are decision points, where the arrow we follow next depends on the state of our cake. Dotted arrows show routes that iterate back to previous operations. This algorithm takes ingredients as its input and outputs cake with either ice cream or custard!

Figure 1.2. Traditional AI vs. machine learning AI. In traditional AI applications, we provide the computer with a complete set of rules. When it’s given data, it outputs the relevant answers. In machine learning, we provide the computer with data and the answers, and it learns the rules for itself. When we pass new data through these rules, we get answers for this new data.

Figure 1.3. Any straight line can be described by its slope (the change in y divided by the change in x) and its intercept (where it crosses the y-axis when x = 0). The equation y = intercept + slope * x can be used to predict the value of y given a value of x.

Figure 1.4. A hypothetical algorithm for learning the parameters of a straight line. This algorithm takes two continuous variables as inputs and fits a straight line through the mean. It iteratively rotates the line until it finds a solution that minimizes the sum of squares. The parameters of the line are output as the learned model.

Figure 1.5. Supervised vs. unsupervised machine learning. Supervised algorithms take data that is already labeled with a ground truth and build a model that can predict the labels of unlabeled, new data. Unsupervised algorithms take unlabeled data and learn patterns within it, such that new data can be mapped onto these patterns.

Figure 1.6. Classification, regression, dimension reduction, and clustering. Classification and regression algorithms build models that predict categorical and continuous variables of unlabeled, new data, respectively. Dimension-reduction algorithms create a new representation of the original data in fewer dimensions and map new data onto this representation. Clustering algorithms identify clusters within the data and map new data onto these clusters.

Figure 1.7. The relationship between artificial intelligence (AI), machine learning, and deep learning. Deep learning comprises a collection of techniques that form a subset of machine learning techniques, which themselves are a subfield of AI.

Chapter 2. Tidying, manipulating, and plotting data with the tidyverse

Figure 2.1. A scatter plot created with ggplot2. The Sepal.Length variable is mapped to the x aesthetic, and the Sepal.Width variable is mapped to the y aesthetic. A black-and-white theme was applied by adding the theme_bw() layer.

Figure 2.2. The same scatter plot as in figure 2.1, with 2D density contours and a smoothed line added as layers using the geom_density_2d() and geom_smooth functions, respectively.

Figure 2.3. The same scatter plot as in figure 2.1, with the Species variable mapped to the shape and col aesthetics

Figure 2.4. The same data is shown, but with different iris species plotted on separate subplots or facets.

Figure 2.5. The map() function takes a vector or list as input, applies a function to each element individually, and returns a list of the returned values.

Figure 2.6. The result of “walking” the hist() function over each element of our list using walk()

Figure 2.7. The result of “walking” the hist() function over each element of our list using iwalk()

Figure 2.8. The pmap() function was used to iterate the rnorm() function over three vectors of arguments. The output from pmap() was piped into iwalk() to iterate the hist() function over each random sample.

Chapter 3. Classifying based on similarities with k-nearest neighbors

Figure 3.1. Body length and aggression of reptiles. Labeled cases for adders, grass snakes, and slow worms are indicated by their shape. New, unlabeled data are shown by black crosses.

Figure 3.2. The first step of the kNN algorithm: calculating distance. The lines represent the distance between one of the unlabeled cases (the cross) and each of the labeled cases.

Figure 3.3. The second step of the kNN algorithm: ranking the neighbors. The lines represent the distance between one of the unlabeled cases (the cross) and each of the labeled cases. The numbers represent the ranked distance between the unlabeled case (the cross) and each labeled case (1 = closest).

Figure 3.4. The final step of the kNN algorithm: identifying the k-nearest neighbors and taking the majority vote. Lines connect the unlabeled data with their one, three, and five nearest neighbors. The majority vote in each scenario is indicated by the shape drawn under each cross.

Figure 3.5. Plotting the relationships between variables in diabetesTib. All three combinations of the continuous variables are shown, shaded by class.

Figure 3.6. Defining a task in mlr. A task definition consists of the data containing the predictor variables and, for classification and regression problems, a target variable we want to predict. For unsupervised learning, the target is omitted.

Figure 3.7. Defining a learner in mlr. A learner definition consists of the class of algorithm you want to use, the name of the individual algorithm, and, optionally, any additional arguments to control the algorithm’s behavior.

Figure 3.8. Training a model in mlr. Training a model simply consists of combining a learner with a task.

Figure 3.9. A summary of the predict() and performance() functions of mlr. predict() passes observations into a model and outputs the predicted values. performance() compares these predicted values to the cases’ true values and outputs one or more performance metrics summarizing the similarity between the two.

Figure 3.10. The bias-variance trade-off. Generalization error is the proportion of erroneous predictions a model makes and is a result of overfitting and underfitting. The error associated with overfitting (too complex a model) is variance. The error associated with underfitting (too simple a model) is bias. The error associated with overfitting (too complex a model) is variance. An optimal model balances this trade-off.

Figure 3.11. Examples of underfitting, optimal fitting, and overfitting for a two-class classification problem. The dotted line represents a decision boundary.

Figure 3.12. Holdout CV. The data is randomly split into a training set and test set. The training set is used to train the model, which is then used to make predictions on the test set. The similarity of the predictions to the true values of the test set is used to evaluate model performance.

Figure 3.13. K-fold CV. The data is randomly split into near equally sized folds. Each fold is used as the test set once, with the rest of the data used as the training set. The similarity of the predictions to the true values of the test set is used to evaluate model performance.

Figure 3.14. Leave-one-out CV is the extreme of k-fold, where we reserve a single case as the test set and train the model on the remaining data. The similarity of the predictions to the true values of the test set is used to evaluate model performance.

Figure 3.15. The MMCE values from fitting the kNN model with different values of k during a grid search

Figure 3.16. Nested CV. The dataset is split into folds. For each fold, the training set is used to create sets of inner k-fold CV. Each of these inner sets cross-validates a single hyperparameter value by splitting the data into training and test sets. For each fold in these inner sets, a model is trained using the training set and evaluated on the test set, using that set’s hyperparameter value. The hyperparameter from each inner CV loop that gives the best-performing model is used to train the models on the outer loop.

Chapter 4. Classifying based on odds with logistic regression

Figure 4.1. Logistic regression learns models that output the probability (p) of new data belonging to each of the classes. Typically, new data is assigned the class to which it has the highest probability of belonging. The dotted arrow indicates that there are additional steps in calculating the probabilities, which we’ll discuss in section 4.1.1.

Figure 4.2. Plotting copper content against class. The y-axis displays the categorical class membership as if it were a continuous variable, with forgeries and originals taking the values of 0 and 1, respectively. The solid line represents a poor attempt to model a linear relationship between copper content and class. The dashed line at y = 0.5 indicates the threshold of classification.

Figure 4.3. Modeling the data with the logistic function. The S-shaped curve represents the logistic function fitted to the data. The center of the curve passes through the mean of copper content and maps it between 0 and 1.

Figure 4.4. Plotting the odds of being an original against copper content. The probabilities derived from the logistic function were converted into odds and plotted against copper content. Odds can take any positive value. The straight line represents a poor attempt to model a linear relationship between copper content and odds.

Figure 4.5. Plotting the log odds of being an original against copper content. The odds were converted into log odds using the logit function and plotted against copper content. Log odds are unbounded and can take any value. The straight line represents the linear relationship between copper content and log odds.

Figure 4.6. Summary of how logistic regression models predict class membership. Data is converted into log odds (logits), which are converted into odds and then into the probability of belonging to the “positive” class. Cases are assigned to the positive class if their probability exceeds a threshold probability (0.5 by default).

Figure 4.7. Visualizing a logistic regression model with two predictors. Copper content and lead content are plotted on the x- and z-axes, respectively. Log odds are plotted on the y-axis. The plane shown inside the plot represents the linear model that combines the intercept and the slopes of copper content and lead content to predict log odds.

Figure 4.8. The process of classifying new paintings. The predictor variable values of three paintings are converted into log odds based on the learned model parameters (intercept and slopes). The log odds are converted into probabilities (p), and if p > 0.5, the case is classified as the “positive” class.

Figure 4.9. Summary of the softmax function. In the binomial case, only one logit is needed per case (the logit for the positive class). Where there are multiple classes (a, b, and c in this example), the model estimates one logit per class for each case. The softmax function maps these logits to probabilities that sum to one. The case is assigned to the class with the largest probability.

Figure 4.10. Faceted plot of Survived against FamSize and Fare. Violin plots show the density of data along the y-axis. The lines on each violin represent the first quartile, median, and third quartile (from lowest to highest).

Figure 4.11. Faceted plot of Survived against Pclass and Sex. Filled bars represent the proportion of passengers at each level of the factors that survived (1 = survival).

Chapter 5. Classifying by maximizing separation with discriminant analysis

Figure 5.1. Discriminant analysis algorithms take the original data and combine continuous predictor variables together into new variables that maximize the separation of the classes.

Figure 5.2. Data becomes more sparse as the number of dimensions increases. Two classes are shown in one-, two-, and three-dimensional feature spaces. The dotted lines in the three-dimensional representation are to clarify the position of the points along the z-axis. Note the increasing empty space with increased dimensions.

Figure 5.3. Learning a discriminant function in two dimensions. LDA learns a new axis such that, when the data is projected onto it (dashed lines), it maximizes the difference between class means while minimizing intra-class variance. and s2 are the mean and variance of each class along the new axis, respectively.

Figure 5.4. Constructing a new axis that only maximizes class centroid separation doesn’t fully resolve the classes (left example). Constructing a new axis that maximizes centroid separation while also minimizing variance within each class results in better separation of the classes (right). and s2 are the mean and variance of each class along the new axis, respectively.

Figure 5.5. When there are more than two classes, LDA maximizes the distance between each class centroid (triangles) and the grand centroid (cross) while minimizing intra-class variance. Once the first discriminant function is found, a second is constructed that is orthogonal to it. The original data can be plotted against these functions.

Figure 5.6. When there are more than two predictors, the cube represents a feature space with three predictor variables (x, y, and z) and three classes (dotted lines help indicate the position of each case along the z-axis). Discriminant function 1 (DF1) is found, and then DF2, which is orthogonal to DF1, is found. Dotted lines indicate “shadows” of DF1 and DF2 to help show their depth along the z-axis. The data can be projected onto DF1 and DF2.

Figure 5.7. Examples of two classes which have equal covariance (the relationship between variable 1 and 2 is the same for both classes) and different covariances. Ovals represent distributions of data within each class. Quadratic and linear DFs (QDF and LDF) are shown. The projection of the classes with different covariances onto each DF is shown.

Figure 5.8. The probability density of each class is assumed to be normally distributed, where the center of each distribution is the centroid of the class. This is shown for one DF (for classes k and j) and for two.

Figure 5.9. Box and whiskers plots of each continuous variable in the data against vineyard number. For the box and whiskers, the thick horizontal line represents the median, the box represents the interquartile range (IQR), the whiskers represent the Tukey range (1.5 times the IQR above and below the quartiles), and the dots represent data outside of the Tukey range.

Figure 5.10. Plotting the DFs against each other. The values for LD1 and LD2 for each case are plotted against each other, shaded by their class.

Chapter 6. Classifying with naive Bayes and support vector machines

Figure 6.1. Using Bayes’ rule to calculate the posterior probability of having unicorn disease, given a positive test result. The priors are the proportion of people with or without the disease. The likelihoods are the probabilities of getting positive or negative test results for each disease status. The evidence is the probability of getting a positive test result (true positives plus the false positives).

Figure 6.2. Filled bar charts showing the proportion of Democrats and Republicans that voted for (y) or against (n) or abstained (NA) on 16 different votes.

Figure 6.3. The SVM algorithm finds a hyperplane (solid line) that passes through the feature space. An optimal hyperplane is one that maximizes the margin around itself (dotted lines). The margin is a region around the hyperplane that touches the fewest cases. Support vectors are shown with double circles.

Figure 6.4. The position of the hyperplane is entirely dependent on the position of support vectors. Moving a support vector moves the hyperplane from its original position (dotted line) to a new position (top two plots). Moving a non-support vector has no impact on the hyperplane (bottom two plots).

Figure 6.5. The SVM algorithm adds an extra dimension to linearly separate the data. The classes in the original data are not linearly separable. The SVM algorithm adds an extra dimension that, in a two-dimensional feature space, can be illustrated as a “stretching” of the data into a third dimension. This additional dimension allows the data to be linearly separated. When this hyperplane is projected back onto the original two dimensions, it appears as a curved decision boundary.

Figure 6.6. Examples of kernel functions. For each example, the solid line indicates the decision boundary (projected back onto the original feature space), and the dashed lines indicate the margin. With the exception of the linear kernel, imagine that the cases of one of the groups are raised off the page in a third dimension.

Figure 6.7. The impact of the cost and gamma hyperparameters. Larger values of the cost hyperparameter give greater penalization for having cases inside the margin. Larger values of the gamma hyperparameter mean individual cases have greater influence on the position of the decision boundary, leading to more complex decision boundaries.

Figure 6.8. One-versus-all and one-versus-one approaches to multiclass SVMs. In the one-versus-all approach, a hyperplane is learned per class, separating it from all the other cases. In the one-versus-one approach, a hyperplane is learned for every pair of classes, separating them while ignoring the data from the other classes.

Figure 6.9. How Platt scaling is used to get probabilities for each hyperplane. This example shows a one-versus-all approach (it also applies to one-versus-one). For each hyperplane, the distance of each case from the hyperplane is recorded (indicated by double-headed arrows). These distances are converted into probabilities using the logistic function.

Chapter 7. Classifying with decision trees

Figure 7.1. The structure of a decision tree. The root node is the node that contains all the data prior to splitting. Nodes are split by a splitting criterion into two branches, each of which leads to another node. Nodes that do not split any further are called leaves.

Figure 7.2. An example decision tree split for 20 cases belonging to classes A and B

Figure 7.3. Calculating the Gini index of the parent node and the left and right leaves

Figure 7.4. How splitting is performed for continuous predictors. Cases belonging to three classes are plotted against two continuous variables. The first node splits the feature space into rectangles based on the value of variable 2. The second node further splits the variable 2 ≥ 20 feature space into rectangles based on the value of variable 1.

Figure 7.5. How the split point is chosen for continuous predictors. Cases (circles) are arranged in order of their value of the continuous predictor. The midpoint between each adjacent pair of cases is considered as a candidate split, and the Gini gain is calculated for each. If one of these splits has the highest Gini gain of any candidate split, it will be used to split the tree at this node.

Figure 7.6. How the split point is chosen for categorical predictors. The Gini index of each factor level is calculated using the proportion of cases from each class with that factor level. The factor levels are arranged in order of their Gini indices, and the Gini gain is evaluated for each split between adjacent levels.

Figure 7.7. Hyperparameters of rpart. Important nodes are highlighted in each example, and the numbers in each node represent the number of cases. The minsplit, maxdepth, cp, and minbucket hyperparameters all simultaneously constrain the splitting of each node.

Figure 7.8. Graphical representation of our decision tree model. The splitting criterion is shown for each node. Each leaf node shows the predicted class, the proportion of each of the classes in that leaf, and the proportion of all cases in that leaf.

Chapter 8. Improving decision trees with random forests and boosting

Figure 8.1. Bootstrap aggregating (bagging) with decision trees. Multiple decision trees are learned in parallel, each one trained on a bootstrap sample of cases from the training set. When predicting new data, each tree makes a prediction, and the modal (most frequent) prediction wins.

Figure 8.2. Adaptive boosting with decision trees. An initial model is trained on a random sample of the training set. Correctly classified cases get lower weights, while incorrectly classified cases get higher weights (indicated by data point size). The probability of subsequent models sampling each case is proportional to the case’s weight. As trees are added, they vote to form an ensemble model, the predictions of which are used to update the weights at every iteration.

Figure 8.3. Plotting the mean out-of-bag error against tree number. For a given forest size during training, the mean out-of-bag error is plotted on the y-axis for each class (different lines) and for the overall out of bag (OOB) error. The out-of-bag error is the mean prediction error for each case, by trees that did not include that case in their bootstrap sample. The y-axis shows the mean out-of-bag error across all cases.

Figure 8.4. Plotting log loss against the number of trees during model building. The curve flattens out after 15 trees, suggesting there is no benefit to adding more trees to the model.

Figure 8.5. Plotting individual trees from our XGBoost model

Chapter 9. Linear regression

Figure 9.1. Imaginary data of how the pH of cider batches changes with apple content

Figure 9.2. Finding the least squares line through the data. Residuals are the vertical distances between the cases and the line. The area of the boxes represents the squared residuals for three of the cases. The intercept (β0) is where the line hits the y-axis when x = 0. The slope is the change in y (Δy) divided by the change in x (Δx).

Figure 9.3. Adding an additional variable: the size of each dot corresponds to the fermentation time of each cider batch.

Figure 9.4. Representing a linear model with two predictors. Combining apple content and fermentation time in our linear model can be represented as a surface. The solid lines show the residual error for each case (its vertical distance from the surface).

Figure 9.5. Finding the slope between two levels of a categorical variable using a dummy variable. The apple types are recoded as 0 and 1 and treated as a continuous variable. The slope now represents the difference in means between the two apple types, and the intercept represents the mean of the reference category (Gala).

Figure 9.6. Recoding a polytomous categorical variable into k – 1 dummy variables. A four-level factor can be represented using three (k – 1) dummy variables. The reference level (Gala) has a value of 0 for each dummy variable. The other levels have a value of 1 for a particular dummy variable. A slope is estimated for each dummy.

Figure 9.7. Plotting each predictor variable in the Ozone dataset against the Ozone variable. The straight lines represent linear regression lines, and the curved lines represent GAM lines.

Figure 9.8. Plotting the correlation of each predictor against the ozone level using plotFilterValues()

Figure 9.9. Combining multiple preprocessing wrappers. Once a learner and preprocessing step (such as imputation) have been combined in a wrapper, this wrapper can be used as the learner in another wrapper.

Figure 9.10. Plotting diagnostic plots for our linear model. The Residuals vs. Fitted and Scale-Location plots help identify patterns that suggest nonlinearity and heteroscedasticity. The Normal Q-Q plot helps identify non-normality of residuals, and the Residuals vs. Leverage plot helps identify influential outliers.

Chapter 10. Nonlinear regression with generalized additive models

Figure 10.1. Shapes of polynomial functions from the first to the fourth degree. When the x variable is raised to the first power, the equation models a straight line. As we increase the power that x is raised to, the equations model lines with varying degrees of flexibility.

Figure 10.2. Comparing linear and quadratic fits to an imaginary nonlinear relationship between apple content and cider acidity

Figure 10.3. Comparing the shape of polynomial functions that do and do not include the first-degree term. Vertical dotted lines indicate the position of each function’s vertex on the x-axis.

Figure 10.4. Fitting a spline to a nonlinear relationship. The solid dots indicate the knots. Individual polynomial functions fit the data between the knots and connect to each other through them.

Figure 10.5. Smoothing functions for continuous variables in GAMs are commonly the sum of a series of basis functions, which are often splines. Three spline basis functions are summed at each value of x to predict the value of y. The dotted line shows the sum of the three basis functions, which models the nonlinear relationship in the data.

Figure 10.6. Plotting the DayOfYear variable against ozone levels

Figure 10.7. Plotting the nonlinear relationships learned by our GAM. The rug at the base of each plot shows the position of each case along the x-axis. The residual vs. fitted plot (middle panel of the second row) shows a pattern suggestive of heteroscedasticity, and the normal Q-Q plot (right panel of the second row) shows the residuals are normally distributed.

Chapter 11. Preventing overfitting with ridge regression, LASSO, and elastic net

Figure 11.1. Examples of underfitting, optimal fitting, and overfitting for a two-class classification problem. The dotted line represents a decision boundary. class classification problem. The dotted line represents a decision boundary.

Figure 11.2. Examples of underfitting, optimal fitting, and overfitting for a singlepredictor regression problem. The dotted line represents the regression line.

Figure 11.3. Sampling error leads to models that don’t generalize well to new data. In the left-side example, a regression line is fit, considering only the more darkly shaded cases. In the right-side example, all the cases are used to construct the regression line. The dotted lines help indicate that the magnitude of the slope is larger on the left side than on the right.

Figure 11.4. Calculating the sum of squares from a model that predicts the number of people in a park based on the temperature

Figure 11.5. Calculating the sum of squares and the L2 norm for the slope between temperature and the number of people at the park.

Figure 11.6. A graphical representation of the ridge regression penalty. The x- and y-axes represent the values of two model parameters. The solid, concentric circles represent the sum of squares value for different combinations of the parameters. The dashed circles represent the L2 norm multiplied by lambda.

Figure 11.7. A graphical representation of the LASSO penalty. The x- and y-axes represent the values of two model parameters. The solid, concentric circles represent the sum of squares value for different combinations of the parameters. The dashed diamonds represent the L2 norm multiplied by lambda.

Figure 11.8. Comparing the ridge regression and LASSO penalties

Figure 11.9. Comparing the shape of the ridge regression, LASSO, and elastic net penalties

Figure 11.10. Plotting each of the predictors against wheat yield for the Iowa dataset. Lines represent linear model fits between each predictor and yield.

Figure 11.11. Plotting the result of generateFilterValuesData(). Bar height represents how much information each predictor contains about wheat yield.

Figure 11.12. Plotting the ridge regression lambda-tuning process. The x-axis represents lambda, and the y-axis represents the mean MSE. Dots represent values of lambda sampled by the random search. The line connects the dots.

Figure 11.13. Comparing the parameter estimates of our ridge regression model to our OLS regression model

Figure 11.14. Plotting the LASSO lambda-tuning process. The x-axis represents lambda, and the y-axis represents the mean MSE. Dots represent values of lambda sampled by the random search. The line connects the dots.

Figure 11.15. Comparing the parameter estimates of our ridge regression model, LASSO model, and OLS regression model

Figure 11.16. Plotting the hyperparameter tuning process for our elastic net model. The x-axis represents lambda, the y-axis represents alpha, and the shading represents mean MSE. The white dot represents the combination of hyperparameters chosen by our tuning process.

Figure 11.17. Comparing the parameter estimates of our ridge regression model, LASSO model, elastic net model, and OLS regression model

Chapter 12. Regression with kNN, random forest, and XGBoost

Figure 12.1. An example relationship for how long your commute to work takes, depending on what time you leave the house

Figure 12.2. The kNN algorithm for classification: identifying the k nearest neighbors and taking the majority vote. Lines connect the unlabeled data with their one, three, and five nearest neighbors. The majority vote in each scenario is indicated by the shape drawn under each cross.

Figure 12.3. How the kNN algorithm predicts continuous variables. The crosses represent new data points for which we wish to predict the journey length. For the one-, three-, and five-nearest neighbor models, the nearest neighbors to each new data point are highlighted in a lighter shade. In each case, the predicted value is the mean journey length of the nearest neighbors.

Figure 12.4. How splitting is performed for classification problems. Cases belonging to three classes are plotted against two continuous variables. The first node splits the feature space into rectangles based on the value of variable 2. The second node further splits the variable 2 ≥ 20 feature space into rectangles based on the value of variable 1.

Figure 12.5. How splitting is performed for regression problems. The feature space is split into shaded regions based on the nodes of the tree next to each plot. The predicted journey length is shown inside each region. The dashed line in the top plot demonstrates how journey length is predicted from departure time based on the tree. The bottom plot shows a two-predictor situation.

Figure 12.6. How candidate splits are chosen for regression problems. The measure of purity is the sum of squares for the split, which is the combined sums of squares for the left and right nodes. Each sum of squares is the vertical distance between each case and the predicted value for the leaf it belongs to.

Figure 12.7. Plotting the relationships in the fuelTib dataset. The topmost plots show absorbance against heatan with separate lines drawn for each wavelength, faceted by near-infrared (NIR) or ultraviolet (UVV) light. The middle plots show wavelength against absorbance shaded by heatan with separate lines drawn for each batch of fuel, faceted by NIR or UVV light. The bottom plot shows h20 against heatan.

Figure 12.8. Plotting our hyperparameter tuning process. The average MSE (mse.test.mean) is shown for each value of k.

Figure 12.9. Plotting the out-of-bag error for our random forest model. The Error y-axis shows the mean square error for all cases, predicted by trees that didn’t include that case in the training set. This is shown for varying numbers of trees in the ensemble. The flattening out of the line suggests we have included enough individual trees in the forest.

Figure 12.10. Plotting the average root mean square error (train_rmse) against the iteration of the boosting process. The curve flattens out just before 30 iterations, suggesting that we have included enough trees in our ensemble.

Chapter 13. Maximizing variance with principal component analysis

Figure 13.1. Data becomes more sparse as the number of dimensions increases. Two classes are shown in one-, two-, and three-dimensional feature spaces. The dotted lines in the three-dimensional representation are to clarify the position of the points along the z-axis. Note the increasing empty space with increased dimensions.

Figure 13.2. The first thing we do before applying the PCA algorithm is (usually) to center the data by subtracting the mean of each variable for each case. This places the origin at the center of the data. The first principal axis is then found: it is the axis that passes through the origin and maximizes the variance of the data when projected onto it.

Figure 13.3. What it means for the first principal axis to “maximize the variance.” The left-side plot shows a sub-optimal candidate principal axis. The right-side plot shows the optimal candidate principal axis. The data is shown projected onto each principal axis below the respective plots. The variance of the data along the axis is greatest on the right side.

Figure 13.4. Calculating the eigenvector for a principal component. The distances along each variable are scaled so that they mark a point that is one unit along the principal axis away from the origin. We can illustrate this graphically by taking a triangle defined by the change in one variable over the change in the other variable, and using Pythagoras’s theorem to find the distance from the origin to divide by.

Figure 13.5. In a two-dimensional feature space, the first principal axis is the one that maximizes the variance (as it always is), and the second principal axis is orthogonal (at a right angle) to the first. In this situation, plotting the principal components simply results in a rotation of the data.

Figure 13.6. In a three-dimensional feature space, the second principal axis is still orthogonal to the first principal axis, but it has freedom to rotate around the first (indicated by the ellipse with arrows in the left-side plot) until it maximizes the remaining variance. The third principal axis is orthogonal to the first and second principal axes and so has no freedom to rotate; it explains the least amount of variance.

Figure 13.7. The result of calling the ggpairs() function on our banknote dataset. Each variable is plotted against every other variable, with different plot types drawn depending on the combination of variable types.

Figure 13.8. Typical exploratory plots for PCA analysis as supplied by the factoextra package. The top-left plot shows a biplot, combining each case’s component scores with arrows to show the variable loadings. The top-right plot shows the variable loading plot with a correlation circle (the boundary within which the variable loadings must lie). The bottom scree plots show the eigenvalue (left) and percentage explained variance (right).

Figure 13.9. The PCA component scores are plotted for each case, shaded by whether they were genuine or counterfeit.

Chapter 14. Maximizing similarity with t-SNE and UMAP

Figure 14.1. t-SNE measures the distance from each case to every other case, converted into a probability by fitting a normal distribution over the current case. These probabilities are scaled by dividing them by their sum, so that they add to 1.

Figure 14.2. The scaled probabilities for each case are stored as a matrix of values. This is visualized here as a heatmap: the closer two cases are, the darker the box is that represents their distance in the heatmap.

Figure 14.3. When converting distances in the lower-dimensional representation into probabilities, t-SNE fits a Student’s t distribution over the current case instead of a normal distribution. The Student’s t distribution has longer tails, meaning dissimilar cases are pushed further away to achieve the same probability as in the high-dimensional representation.

Figure 14.4. Cases are randomly initialized over the new axes (one axis is shown here). The probability matrix is computed for this axis, and the cases are shuffled around to make this matrix resemble the original, high-dimensional matrix by minimizing the Kullback-Leibler (KL) divergence. During shuffling, cases are attracted toward cases that are similar to them (lines with circles) and repulsed away from cases that are dissimilar (lines with triangles).

Figure 14.5. The effect on the final t-SNE embedding of the banknote dataset of changing theta (row facets) and perplexity (column facets) using the default values of eta and max_iter

Figure 14.6. The effect on the final t-SNE embedding of the banknote dataset of changing max_iter (row facets) and eta (column facets) using the default values of theta and perplexity

Figure 14.7. tSNE1 and tSNE2 axes plotted against each other, faceted and shaded by the original variables, and shaped by whether each case was a genuine or counterfeit banknote

Figure 14.8. How UMAP learns a manifold. UMAP expands a search region around each case. A naive form of this is shown in the top row, where the radius of each search region is the same. When cases with overlapping search regions are connected by edges, there are gaps in the manifold. In the bottom row, the search region extends to the nearest neighbor and then extends outward in a fuzzy way, with a radius inversely related to the density of data in that region. This results in a complete manifold.

Figure 14.9. The effect on the final UMAP embedding of the banknote dataset of changing n_neighbors (row facets) and min_dist (column facets) using the default values of metric and n_epochs

Figure 14.10. The effect on the final UMAP embedding of the swissTib dataset of changing metric (row facets) and n_epochs (column facets) using the default values of n_neighbors and min_dist

Figure 14.11. UMAP1 and UMAP2 axes plotted against each other, faceted and shaded by the original variables, and shaped by whether each case was a genuine or counterfeit banknote

Chapter 15. Self-organizing maps and locally linear embedding

Figure 15.1. Placing beads into bowls based on their characteristics. Similar beads are placed in the same or nearby bowls, while dissimilar beads are placed in bowls far away from each other. One bowl didn’t have any beads placed in it, but that’s okay.

Figure 15.2. A sphere is a three-dimensional manifold. We can reconstruct a sphere as a series of linear patches that connect to one another. This three-dimensional manifold of a sphere can be represented in two dimensions by cutting a sheet of paper in a certain way.

Figure 15.3. Common graphical representations of square and hexagonal self-organizing maps. The top two maps show a grid of rectangular nodes that are each connected to four neighbors. The bottom two maps show a grid of hexagonal nodes that are each connected to six neighbors.

Figure 15.4. How distances between each case to each node are calculated. The arrows pointing from each variable to each node represent the weight for that variable on that particular node (for example, the weights of node 1 are 3, 9, and 1). Distance is calculated by finding the difference between a node’s weights and a case’s value for each variable, squaring these differences, and summing them.

Figure 15.5. At each stage of the algorithm, the node whose weights have the smallest distance to a particular case is selected as the best matching unit (BMU) for that case.

Figure 15.6. Between the first and last iteration of the algorithm, both the radius of the neighborhood around a BMU (the darkest node) and the amount by which neighboring node weights are updated get smaller. The radius of a Gaussian neighborhood function is shown as a translucent circle centered over the BMU, and the amount each neighboring node is updated is represented by how dark its shading is. If the bubble neighborhood function was shown, all nodes would be shaded the same (as they’re updated by the same amount).

Figure 15.7. A matrix of plots created using the ggpairs() function, plotting all variables against each other from the flea dataset. Because the individual plots are quite small, I’ve manually zoomed in on one plot with a virtual magnifying glass (much like one you might need to use to see the fleas).

Figure 15.8. Diagnostic plots for our SOM. The Codes fan plot for each node indicates the weight for each variable. The Training Progress plot shows the mean distance between each case and its BMU for each iteration. The Counts plot shows the number of cases per node. The Quality plot shows the mean distance between each case and the weights of its BMU. The Neighbor Distance plot shows the sum of differences between cases in one node and cases in neighboring nodes. The Mapping plot draws the cases inside their assigned nodes.

Figure 15.9. Separate heatmaps showing node weights for each original variable. The scales are in standard deviation units.

Figure 15.10. Showing class membership on the SOM. The left-side mapping plot shows cases drawn inside their assigned nodes, shaded by which flea species they belong to. The right-side plot shows the same information, but nodes are shaded by cluster membership after applying a clustering algorithm to the nodes. The solid black lines separate nodes assigned to different clusters.

Figure 15.11. New data can be mapped onto an existing SOM. This mapping plot shows a graphical representation of the nodes to which the two new cases are assigned.

Figure 15.12. The distance between each case and every other case is calculated, and their k-nearest neighbors are assigned (distance along the z-axis in the top-left plot is indicated by the size of the circle). For each case, the algorithm learns a set of weights, one for each nearest neighbor, that sum to 1. Each neighbor’s variable values are multiplied by its weight (so row 1 becomes x = 3.1 × 0.1, y = 2.0 × 0.1, z = 0.1 × 0.1). The weighted values of each neighbor are summed (the columns are summed) to approximate the original values of the selected case.

Figure 15.13. The S-curve dataset plotted in three dimensions using the scatter3D() function. The shading of the points is mapped to the z variable.

Figure 15.14. Plotting 1 – ρ2 against k to find the optimal value of k. The solid horizontal line indicates the value of k with the lowest 1 – ρ2.

Figure 15.15. Plotting the two-dimensional embedding of the S-curve data. The shading of the points is mapped to the z variable, the same as in figure 15.11.

Figure 15.16. Plotting the output of listing 15.12. The top plot shows 1 – ρ2 for different values of k. The lower plot shows the two-dimensional embedding of the flea data, shaded by species.

Chapter 16. Clustering by finding centers with k-means

Figure 16.1. Five iterations of k-means clustering. In the top-left plot, three initial centers are randomly generated in the feature space (crosses). Cases are assigned to the cluster of their nearest center. At each iteration, each center moves to the mean of the cases in its cluster (indicated by arrows). The feature space can be partitioned up into Voronoi cells (I'll discuss these shortly), indicated by the shaded regions, that show regions of the feature space closest to a particular centroid.

Figure 16.2. ggpairs() plot of each variable against every other variable in our GvHD dataset. Scatterplots are shown below the diagonal, 2D density plots are shown above the diagonal, and 1D density plots are drawn on the diagonal. It appears as if there are multiple clusters in the data.

Figure 16.3. What under- and overfitting looks like for clustering tasks. In the left-side plot, the clusters are underfit (fewer clusters have been identified than actually exist). In the right-side plot, the clusters are overfit (real clusters are broken up into smaller clusters). In the center plot, an optimal clustering model has been found that faithfully represents the structure in the data.

Figure 16.4. The Davies-Bouldin index calculates the intracluster (within-cluster) variance (left-side plot) and the distance between the centroids of each cluster (right-side plot). For each cluster, its nearest neighboring cluster is identified, and the sum of their intracluster variances is divided by the difference between their centroids. This value is calculated for each cluster, and the Davies-Bouldin index is the mean of these values.

Figure 16.5. The Dunn index quantifies the ratio between the smallest distance between cases in different clusters (left-side plot) and the largest distance within a cluster (right-side plot).

Figure 16.6. The pseudo F statistic is the ratio of the between-cluster sum of squares (right-side plot) to the within-cluster sum of squares (left-side plot). The grand centroid is shown as a square in the right-side plot.

Figure 16.7. Plotting our tuning process. Each subplot shows the values of a different internal cluster metric. The different lines indicate the performance of each of the three different algorithms.

Figure 16.8. ggpairs() plot with the k-means cluster membership mapped to the color aesthetic. Box plots and histograms show how the values of the continuous variables vary between clusters.

Chapter 17. Hierarchical clustering

Figure 17.1. Dendrogram showing an imaginary clustering of instruments in an orchestra. Horizontal lines indicate the merging of separate clusters. The height of a merge indicates the similarity between the clusters (lower merge, higher similarity).

Figure 17.2. Agglomerative hierarchical clustering merges clusters that are closest to each other at each iteration. Ellipses indicate the formation of clusters at each iteration, going from top left to bottom right.

Figure 17.3. Different linkage methods to define the distance between clusters. Centroid linkage calculates the distance between cluster centroids. Single linkage calculates the smallest distance between clusters. Complete linkage calculates the largest distance between clusters. Average linkage calculates all pairwise distances between cases in two clusters and finds the mean. Ward’s method calculates the within-cluster sum of squares for each candidate merge and chooses the one with the smallest value.

Figure 17.4. The resulting dendrogram representing our hierarchical clustering model. The y-axis represents the distances between cases. Horizontal lines indicate the positions at which cases/clusters merge with each other. The higher the merge, the less similar the clusters are to each other.

Figure 17.5. Plotting the result of our bootstrap experiment. Each subplot shows the result of a different internal cluster metric. The x-axis shows the cluster number, and the y-axis shows the value of each metric. Faint lines connect the results of each individual bootstrap sample, while the bold line connects the mean. The top and bottom of each crossbar indicate the 95% confidence interval for that particular value, and the horizontal line represents the mean.

Figure 17.6. The same plot as in figure 17.4, but this time with rectangles indicating the clusters resulting from cutting the tree

Figure 17.7. ggpairs() plot showing the result of our hierarchical clustering model. Compare these clusters to the ones obtained by k-means in figure 16.8.

Figure 17.8. The graphical output of the clusterboot() function. The first and last plots show the full, original clusters of data, while the plots in between show the clusters on the bootstrap samples. The cluster membership of each case is indicated by a number. Notice the relatively high stability of the clusters.

Chapter 18. Clustering based on density: DBSCAN and OPTICS

Figure 18.1. A challenging clustering problem. The dataset shown in each facet contains clusters of varying shapes and diameters, with cases that could be considered noise. The three subplots show the data clustered using DBSCAN, hierarchical clustering (complete linkage), and k-means (Hartigan-Wong). Of the three algorithms used, only DBSCAN is able to faithfully represent these shapes as distinct clusters.

Figure 18.2. The DBSCAN algorithm. A case is selected at random, and if its epsilon radius (ϵ) contains minPts cases, it is considered a core point. Reachable cases of this core point are evaluated the same way until there are no more reachable cases. This network of density-connected cases is considered a cluster. Cases that are reachable from core points but are not themselves core points are border points. The algorithm moves on to the next unvisited case. Cases that are neither core nor border points are labeled as noise.

Figure 18.3. Defining the core distance and reachability distance. In OPTICS, epsilon (ϵ) is the maximum search distance. The core distance (ϵ′) is the minimum search distance needed to include minPts cases (including the case in question). The reachability distance for a case is the larger of the core distance and the distance between the case in question, and another case inside its epsilon (maximum search distance).

Figure 18.4. Defining core points in OPTICS. Cases for which the core distance (ϵ′) is less than or equal to the maximum search distance (ϵ) are considered core points. Cases for which the core distance is greater than the maximum search distance are not considered core points.

Figure 18.5. The OPTICS algorithm. A case is selected, and its core distance (ϵ′) is measured. The reachability distance is calculated between this case and all the cases inside this case’s maximum search distance (ϵ). The processing order of the dataset is updated such that the nearest case is visited next. The reachability score and the processing order are recorded for this case, and the algorithm moves on to the next one.

Figure 18.6. Reachability plot of a simulated dataset. The top plot shows the reachability plot learned by the OPTICS algorithm from the data shown in the bottom plot. The plots are shaded to indicate where each cluster in the feature space maps onto the reachability plot.

Figure 18.7. An illustration of different ways clusters can be extracted from a reachability plot. In the top plot, a single reachability score cut-off has been defined, and any troughs bordered by peaks above this cut-off are defined as clusters. In the bottom plot, a hierarchy of clusters is defined, based on the steepness of changes in reachability, allowing for clusters within clusters.

Figure 18.8. Plotting the Swiss banknote dataset with ggpairs(). 2D density plots are shown above the diagonal.

Figure 18.9. K-nearest neighbor distance plot with k = 5. Horizontal lines have been drawn using abline() to highlight the 5-NN distances at the start and end of the knee/elbow in the plot.

Figure 18.10. Visualizing the result of our tuning experiment. Each subplot shows the Right and Diagonal variables plotted against each other for a different permutation of epsilon and minPts. Cases are shaded by their cluster membership.

Figure 18.11. Visualizing the cluster performance experiment. Each subplot shows a heatmap for the number of clusters, Davies-Bouldin index (db), Dunn index (dunn), and pseudo F statistic (G1) returned by the cluster models. Each tile represents the combination of epsilon and minPts, and the depth of shading of the tile indicates its value for each metric. The blank region at the top right in the metric plots indicates a region with no data, and semitransparent tiles indicate fewer than 10 samples.

Figure 18.12. Plotting our final DBSCAN cluster model with ggpairs(). This plot includes the noise cluster.

Figure 18.13. Plotting our final DBSCAN cluster model with ggpairs(). This plot excludes the noise cluster.

Figure 18.14. The reachability plot generated from applying the OPTICS algorithm to our data. The x-axis shows the processing order of the cases, and the y-axis shows the reachability distance for each case. We can see two main troughs in the plot, bordered by peaks of higher reachability distance.

Figure 18.15. Plotting our final OPTICS cluster model with ggpairs()

Chapter 19. Clustering based on distributions with mixture modeling

Figure 19.1. The expectation-maximization algorithm for two, one-dimensional Gaussians. Dots represent cases along a number line. Two Gaussians are randomly initialized along the line. In the expectation step, the posterior probability of each case for each Gaussian is calculated (indicated by shading). In the maximization step, the means, variances, and priors for each Gaussian are updated based on the calculated posteriors. The process continues until the likelihood converges.

Figure 19.2. The expectation-maximization algorithm for two, two-dimensional Gaussians. Two Gaussians are randomly initialized in the feature space. In the expectation step, the posterior probabilities for each case are calculated for each Gaussian. In the maximization step, the centroids, covariance matrices, and priors are updated for each Gaussian, based on the posteriors. The process continues until the likelihood converges.

Figure 19.3. The BIC plot from our mclust model. The x-axis shows the number of clusters, the y-axis shows the Bayesian information criterion (BIC), and each line shows a different model, with the three-letter code indicating which constraints are put on the covariance matrix. In this arrangement of the BIC, higher values indicate better-fitting and/or more parsimonious models.

Figure 19.4. The classification plot from our mclust model. All variables in the original data are plotted against each other in a scatterplot matrix, with cases shaded and shaped according to their cluster. Ellipses indicate the covariances of each Gaussian, and stars indicate their centroids.

Figure 19.5. The uncertainty plot from our mclust model. This plot is similar to the classification plot, except that the size of each case corresponds to its uncertainty under the final model.

Figure 19.6. The density plot from our mclust model. This matrix of plots shows the 2D density of the final model for each combination of variables in the feature space.

Chapter 20. Final notes and further reading

Figure 20.1. Supervised vs. unsupervised machine learning. Supervised algorithms take data already labeled with a ground truth and build a model that can predict the labels of new, unlabeled data. Unsupervised algorithms take unlabeled data and learn patterns within it, such that new data can be mapped onto these patterns.

Figure 20.2. Summary of the algorithms we cover in the book, whether they are supervised or unsupervised learners, and whether they can be used for classification, regression, dimension reduction, or clustering

Figure 20.3. Holdout cross-validation. The data is randomly split into a training set and a test set. The training set is used to train the model, which is then used to make predictions on the test set. The similarity of the predictions to the true values of the test set is used to evaluate model performance.

Figure 20.4. K-fold cross-validation. The data is randomly split into near equally sized folds. Each fold is used as the test set once, and the rest of the data is used as the training set. The similarity of the predictions to the true values of the test set is used to evaluate model performance.

Figure 20.5. Leave-one-out cross-validation is the extreme of k-fold, where we reserve a single case as the test set and train the model on the remaining data. The similarity of the predictions to the true values of the test set is used to evaluate model performance.

Figure 20.6. Nested cross-validation. The dataset is split into folds. For each fold, the training set is used to create sets of inner k-fold cross-validation. Each of these inner sets cross-validates a single hyperparameter value by splitting the data into training and test sets. For each fold in these inner sets, a model is trained using the training set and evaluated on the test set, using that set’s hyperparameter value. The hyperparameter from each inner cross-validation loop that gives the best-performing model is used to train the models on the outer loop.

Appendix Appendix. Refresher on statistical concepts

Figure A.1. The difference between the population and the sample. The population is the set of all units we would like to generalize our results to. The population is often considered nearly infinite in size. The sample is a more manageable subset we measure, which we hope will represent the population.

Figure A.2. An example of a two-dimensional vector at point x = 3, y = 5. The arrow shows how vectors encode magnitude, which we can represent as their distance from the origin (or another vector). The curved line representing the angle between the x-axis and the arrow indicates how vectors encode direction.

Figure A.3. A histogram showing the distribution of an imaginary sample of hippo tusk lengths. The distribution approximates a Gaussian distribution. The curved line represents the probability density function of a Gaussian distribution with the same mean and standard deviation as the sample.

Figure A.4. For a perfect, Gaussian-distributed variable, 68% of the cases lie within one standard deviation’s distance of the mean. 95% and 99.7% of cases lie within two and three standard deviations, respectively.

Figure A.5. Box and whiskers plots for imaginary hippo tusk data. The thick horizontal lines are the medians, the edges of the boxes represent the first and third quartiles, and the vertical whiskers represent the full range.

Figure A.6. The impact of log10 transformation on variables. The left-side plot shows a y variable with both very small and very large values. In the right-side plot, the y variable has been log10-transformed.