9  Classical supervised machine-learning

9.1 Regression trees

This is our first ML algorithm where we’re going to be predicting (continuous) data, and so it’s our first supervised machine learning algorithm. Because of this, I’m going to walk you through how we validate machine learning algorithms by setting aside some of our data for testing. Machine learning algorithms have no coherent statistical philosophy or definition—they are literally just computational instructions that happen to work well—and, as such, verifying your models using independent data is a vital thing to do. If you remember one thing from this course, let it be this: woe betide those who fit machine learning algorithms but do not test their performance with independent data.

9.1.1 An informal introduction to regression trees

When we try to make decisions, we often use rules of thumb or heuristics. An example is “I’m going to boil that egg for five minutes”: we know that boiling it for four minutes and fifty-nine seconds would give much the same result, as would boiling it for five minutes and one second, but we pick a rule (“stop at five minutes”) and stick with it. Similarly, while you’re used to regression methods that work with continuous data (\(y = 2 \times x\)), regression trees split continuous data into discrete decisions. For small, continuous datasets, that may seem like a bad idea, but for datasets with lots of explanatory variables it’s often quicker to split things like this, because it makes it clear what variables are important and makes it very easy to visualise interactions.

Regression trees work exactly like this, and they look remarkably similar to the hierarchical clustering outputs I believe you’re already familiar with from another course (look in Appendix E if you’re not and you’re curious). Starting at the top of the tree, you can follow the course of it all the way down, asking yourself at each step whether the particular row of data you’re trying to predict you’re in has a lesser (go to the left branch) or greater (go to the right branch) value than the decision in front of you. Once you get to the bottom of the tree—you’re finished! You’ve made your prediction. The process for categorical data is much the same: the only difference is you split things according to the categories they all fall into.

The details of regression trees we’re going to learn through experience. But there are four main kinds of regression tree we’re going to cover (briefly) today, that I list in both the order you’ll encounter them and relative complexity:

  1. Regression trees. Exactly what I’ve just described above.
  2. Bagged regression trees. Split your training data into lots of little subsets, fit across them all, and then average across all those trees.
  3. Random forests. Exactly the same as bagged trees, but each random subset also has a subset of the total set of explanatory variables.
  4. Boosted regression trees. Fit a regression tree to data, then fit a regression tree to the residuals of that model. Average across the two, then fit another model, take those residuals, and fit another model. Keep repeating until you get bored.

Regression trees can also be run with a categorical response variable, in which case they’re called decision trees. I’m sure you’ll agree, however, that four new kinds of machine learning technique are more than enough for your first lesson, so don’t worry about them!

9.1.2 Hands-on with regression trees

First of all, let’s simulate some data that we’re going to work with. We’re going to continue our example1 of plant diversity across a series of (a)biotic gradients—temperature, humidity, soil Carbon, and herbivore diversity—and simulate some data under a reasonable biological model (greater diversity in tropics, and a sort-of trophic cascade effect of herbivory). Notice that we’re going to model diversity of plants using a Poisson distribution to give some noise to the data; if you’re not overly familiar with the Poisson, just treat it as a way of getting some (very noisy) variation around the relationships we’re simulating2. Also don’t worry if you don’t track what’s going on with the model of plant diversity; just use boxplot to examine what’s going on with the relationships of all the variables, and notice that you would have had a great deal of difficulty picking out most of these relationships by eye.

# Build model of species diversity
data <- expand.grid(temperature=seq(0,40,4), humidity=seq(0,100,10),
  carbon=seq(1,10,1), herbivores=seq(0,20,2))
data$plants <- runif(nrow(data), 3, 5)
data$plants <- with(data, plants + temperature * .1)
data$plants[data$humidity > 50] <- with(data[data$humidity > 50,],
  plants + humidity * .05)
data$plants[data$carbon < 2] <- with(data[data$carbon < 2,], plants - carbon)
data$plants <- with(data, plants + herbivores * .1)
data$plants[data$herbivores > 5 & data$herbivores < 15] <-
  with(data[data$herbivores > 5 & data$herbivores < 15,], plants - herbivores * .2)
# Draw random data from Poisson based on this
for(i in seq_len(nrow(data)))
  data$plants[i] <- rpois(1, data$plants[i])

Now let’s fit a regression tree to that data, and then plot it out. Once you see the plots, suddenly the structure of a regression tree is going to make much more sense to you. It’s a series of ‘decisions’ you make, working from the top of the tree to the tips, where each decision maps onto the value of the explanatory variables. At the very bottom, we have the estimate of the response variable (plant richness) that we would expect. Notice that you can get the residual variation out of the model, just as you would with a standard linear regression.

library(tree)
# Pick some training data and then fit a model to it
training <- sample(nrow(data), nrow(data)/2)
model <- tree(plants~., data=data[training,])
# Examine the model
plot(model)
text(model)

# Look at the statistics of the model
model
#> node), split, n, deviance, yval
#>       * denotes terminal node
#> 
#>  1) root 6655 104200  7.867  
#>    2) humidity < 55 3621  35220  6.016  
#>      4) herbivores < 15 2618  21010  5.317  
#>        8) temperature < 22 1441   8239  4.323 *
#>        9) temperature > 22 1177   9603  6.535 *
#>      5) herbivores > 15 1003   9590  7.840  
#>       10) temperature < 22 536   4034  6.836 *
#>       11) temperature > 22 467   4395  8.994 *
#>    3) humidity > 55 3034  41790 10.070  
#>      6) herbivores < 15 2193  27250  9.460  
#>       12) temperature < 26 1394  14530  8.669 *
#>       13) temperature > 26 799  10320 10.840 *
#>      7) herbivores > 15 841  11550 11.680 *
summary(model)
#> 
#> Regression tree:
#> tree(formula = plants ~ ., data = data[training, ])
#> Variables actually used in tree construction:
#> [1] "humidity"    "herbivores"  "temperature"
#> Number of terminal nodes:  7 
#> Residual mean deviance:  9.428 = 62680 / 6648 
#> Distribution of residuals:
#>     Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
#> -10.6800  -2.3230  -0.3227   0.0000   2.0060  16.3200

However, if we really want to be sure of whether our model is doing well, we should test its performance on the data we didn’t fit to it. Remember: this is important to do in machine learning methods because they’re not necessarily based on some fundamental, deep aspects of statistics: they just happen to work very well in the right circumstances. So we need to be careful. We can also cross-validate our model to see how it performs under different tree depths (numbers of nodes), seeing how the Mean Squared Error (calculated just as for a normal regression) changes. It’s your job to decide where to draw the line and what constitutes a good fit, but in this case it’s clear that we’re doing a pretty good job (the line goes down), albeit we could fit a simpler model without much change.

# Check performance outside training set
cor.test(predict(model, data[-training,]), data$plants[-training])
#> 
#>  Pearson's product-moment correlation
#> 
#> data:  predict(model, data[-training, ]) and data$plants[-training]
#> t = 68.1, df = 6653, p-value < 2.2e-16
#> alternative hypothesis: true correlation is not equal to 0
#> 95 percent confidence interval:
#>  0.6265179 0.6548398
#> sample estimates:
#>       cor 
#> 0.6408969
# Check cross-validation of model
plot(cv.tree(model))

Maybe you’re not so impressed with the roughly 0.60 sort-of \(r\) of this dataset. Do you think our use of rpois could have had anything to do with it? How could you find out?…

9.1.3 Bagged regression trees and random forests

The problem with fitting regression trees is that they fit your data so well: that’s why we’ve been splitting our data in half and only training on one half of it to check we’ve not over-fit and so built a bad model. One way, potentially, of dealing with this problem is so randomly take subsets of our training data, fit regression trees to those bagged bootstrap replicates, and use the average of all those regression trees. It might not seem obvious, but by averaging across all of our trees like that, we reduce the variance among them and so (hopefully) reduce the variance associated with only working with a sample of data (i.e., the degree of over-fitting in our data). Maybe that makes sense, maybe it doesn’t, but the take-home is that bagged regression trees, where you work with the average of many regression trees fit to many bootstrapped subsets of your data, are a good solution to the problem of over-fitting in data. Fitting them is simple, and they can be tested just as you tested regression trees.

library(randomForest)
#> randomForest 4.7-1.2
#> Type rfNews() to see new features/changes/bug fixes.
model <- randomForest(plants~., data=data[training,], mtry=ncol(data)-1)
cor.test(predict(model, data[-training,]), data$plants[-training])
#> 
#>  Pearson's product-moment correlation
#> 
#> data:  predict(model, data[-training, ]) and data$plants[-training]
#> t = 67.054, df = 6653, p-value < 2.2e-16
#> alternative hypothesis: true correlation is not equal to 0
#> 95 percent confidence interval:
#>  0.6204808 0.6491617
#> sample estimates:
#>       cor 
#> 0.6350401

You might be wondering what that mtry argument to randomForest is all about. This tells the randomForest package to make sure, each time it’s trying to make a new split in the tree, to consider all the variables available to it. If we go with the default (or set any number less than the total number of variables), we restrict the number of options available randomly each time. When we randomly change the available explanatory variables each time a split is being considered in each of the bootstrapped trees, we are building a random forest model. It might sound a bit strange, but this actually improves things because it means the bootstrapped trees resemble each other less. When you take the average of things that are correlated with one-another (resemble each other) more than you expect, you don’t tend to reduce variance as much as you would hope, and so by ‘decorrelating’ your trees like this you improve estimates.

Below I show you how to fit random forest models, and also show off variable importance. This is simply the average decrease in Mean Squared Error each time a split in a regression tree is fit to a particular explanatory variable: it’s essentially the \(r^2\) of each variable. Since we can’t look at a single regression tree anymore (we’ve fit thousands of them!), this is perhaps the simplest way to understand what variable is doing what in your model.

model <- randomForest(plants~., data=data[training,], importance=TRUE)
importance(model)
#>                %IncMSE IncNodePurity
#> temperature 112.747420     11480.872
#> humidity    145.753770     27871.296
#> carbon        5.369292      4209.613
#> herbivores  101.404893     10794.867
cor.test(predict(model, data[-training,]), data$plants[-training])
#> 
#>  Pearson's product-moment correlation
#> 
#> data:  predict(model, data[-training, ]) and data$plants[-training]
#> t = 75.252, df = 6653, p-value < 2.2e-16
#> alternative hypothesis: true correlation is not equal to 0
#> 95 percent confidence interval:
#>  0.6648914 0.6908564
#> sample estimates:
#>       cor 
#> 0.6780854

9.1.4 Boosted regression trees

Boosted regression trees are perhaps the most ‘meta’ of the regression tree family. They also fit a series of models, but each time to the residuals of the previous model, which is then added in to the set of predictions from the previous model, and the process repeated again. So you end up with a model that’s a hybrid of a regression tree fitted to the original data, and a series of models that are trying to fit the variation that the first model didn’t fit very well. The intention is to avoid over-fitting the data through this approach, but the disadvantage is it makes the meaning of the model a little more obscure. Of course, as with all machine learning algorithms, our intention here isn’t necessarily to be easy to interpret, but to perform well!

There are a lot more parameters to play with in this model-type, but there are two I want to draw your attention to. The first is that gbm allows you to fit on a link function, much as you can in a Generalised Linear Model, so here I’ve told it that we’re dealing with count data (distribution="poisson"). If you’re familiar with this, then that’s all great, but if you’re not then you’ll be fine not setting this option (as I do in the second example). Worry about it once you’ve covered Generalised Linear Models. Second, there is a shrinkage parameter that controls the relative importance of earlier vs. later fitted models in the boosting process. A smaller value means trees later in the process (residuals of the residuals of the residuals of the…) are given relatively more weight than if the parameter is greater. Generally, smaller values give better results, because the whole purpose of the exercise is to allow those later trees to matter, but as with everything experimentation pays dividends.

library(gbm)
#> Loaded gbm 2.3.1
#> This version of gbm is no longer under development. Consider transitioning to gbm3, https://github.com/gbm-developers/gbm3
model <- gbm(plants~., data=data[training,], distribution="poisson")
summary(model)

#>                     var   rel.inf
#> humidity       humidity 58.644554
#> temperature temperature 20.813490
#> herbivores   herbivores 19.274843
#> carbon           carbon  1.267113
# A plot of variable importance should also have appeared now
faster.model <- gbm(plants~., data=data[training,], distribution="poisson", shrinkage=.1)

9.2 Lasso regression and Least Angle Regression (LAR)

Least angle regression (LAR) is essentially something called the lasso regression on steroids, so we’re going to start off by learning how to use the lasso and then move onto its somewhat more unwieldy younger brother LAR. What’s interesting about these two is that, unlike regression trees, they were designed to help us pick between very large numbers of candidate explanatory variables. Whereas regression trees are quite happy to make use of all the factors, and relative variable importance is something of a side-note, that’s not the case in LAR. Here the whole purpose is to figure out which explanatory variables you can throw away and ignore, and so it’s a favourite of bioinformaticians, busy MSc students, and lazy lecturers. I think you’ll like it: I do!

9.2.1 A wordy introduction to lasso regression and LAR

All the regression techniques you learnt before this course involved least squares: you were trying to find an equation that minimised the squared error of the difference between a response variable and a regression line. But the method by which you found that line wasn’t really something that was discussed too much… In lasso regression, we are still trying to minimise the modulus of the error, but that is subject to the constraint that we don’t want the sum of absolute value of our coefficients to be too great. The sum of the absolute value of our coefficients is called the ‘\(L_1\) arc length’, or the lasso penalty, and we want to minimise it because we want to keep the complexity of our model as low as possible: a greater \(L_1\) means we’ve got more ‘stuff’ in our model, and more stuff means more complexity and we hate complexity. This is exactly the same as a standard linear regression with multiple variables: the only difference is we can’t use any of the normal calculus3 to estimate what the ‘significant’ variables are, and so instead we must try and minimise some sort of penalty. The great news is, though, that this means we can use a new method to find our best-fits and that method happens to be very efficient…

To figure out what the best model is, we can make this penalty greater or lesser, and plot out how much better or worse our model predictions get. The one big advantage we have here is that, when the arc length is infinitely small, only one variable will have a non-zero coefficient, and as we increase it variables will slowly add in one-by-one. So if we can figure out the arc-length at which the model fit isn’t really much better, we’ve picked our best model! Breathe for a minute (because that was quite a lot!), and if you get lost focus on the basic principle which will become much clearer once you start fitting these models: we are minimising the absolute error of a model, penalising ourselves to have as simple a model as possible and so finding the best set of predictors in our model. Figure 9.1 shows, on the left-hand-side, an example of how this process works. Lasso starts at the left-hand side of the figure, where all the coefficients are equal to zero. It then slowly increases the arc length (the same thing as decreasing the penalty), slowly increasing the absolute value of a coefficient. Each dashed line represents the point where, now the arc length is allowed to be a certain size, another variable is able to ‘jump in’ and making its coefficient non-zero increases model fit. Somewhere between the far left-hand and right-hand sides of the plot lies a model that optimally trades complexity for explanatory power.

Least angle regression (LAR) is so similar to the lasso that the statistical machinery underlying it is often used to more efficiently estimate lasso regression coefficients. As with lasso, there is a penalty term that we can increase from zero; as we do so, we find the explanatory variable most correlated with the response variable and increase its coefficient as much as the penalty/arc length will allow. Eventually, another variable will correlate more with the left-over variation, and so increasing its coefficient would give us a better fir than increasing the coefficient of the explanatory variable we’ve been working with: so we then increase both variables’ coefficients (still constrained by the penalty). This process continues as we increase the arc length (decrease the penalty), adding in more and more explanatory variables as we do so until everything has been added in (the far right-hand side of Figure 9.1). This is called a least angle regression because, geometrically, when we’re adding all these variables in we’re moving the ‘angle’ of the coefficients in a direction that matches the residual variation in the response variable at each step. Don’t worry about the name, but do notice that both panels in Figure 9.1 look incredibly similar—in practice, lasso and LAR are extremely similar.

Figure 9.1: Model coefficients and arc length in lasso and LAR. Taken from Friedman et al. (2001; The elements of statistical learning.).

9.2.2 Going to the rodeo

Enough talk, let’s do it! Let’s start out by simulating a dataset with 1000 explanatory variables, only two of which (columns 123 and 678) are significantly related to our data. Consider how you would approach the problem of determining what ‘significantly’ explains our response variable with this much data using something like a linear regression.

explanatory <- replicate(1000, rnorm(1000))
response <- explanatory[,123]*1.5 -explanatory[,678]*.5

Now let’s fit a lasso regression. Notice how we’re using the package lars: because the LAR approach can be used to generate lasso estimates more efficiently, this package does everything we need for today.

library(lars)
#> Loaded lars 1.3
model <- lars(explanatory, response, type="lasso")
plot(model)

That plot should look remarkably familiar, as it’s basically the same as what you see in Figure 9.1, only this time with some real (simulated) data. Which is all well and good, but (as you’ve probably noticed if you’ve played around a bit) it’s not that easy to get estimates of coefficients out of this model. That’s because, while lasso (and LAR) has its own stopping criteria for when a model is ‘good enough’, it gives you the coefficients for each step in that search. Luckily, however, it’s reasonably trivial to write your own ‘wrapper’ function that will grab the coefficient estimates in a format you can work with. If you disagree, then you’re in luck, because I’ve written one for you:

signif.coefs <- function(model, threshold=0.001){
  coefs <- coef(model)
  signif <- which(abs(coefs[nrow(coefs),]) > threshold)
  return(setNames(coefs[nrow(coefs),signif], signif))
}
signif.coefs(model)
#>  123  678 
#>  1.5 -0.5

Hopefully the above code isn’t complete gibberish to you, but you are quite welcome to treat it as a magic function that will tell you the final coefficients whose coefficients are interesting in the final model. You can, and should, play with the threshold for what defines an ‘interesting’ coefficient (see the next section for the importance of variable scaling, which will affect this too).

In case you’ve missed it, lasso regression works amazingly well, particularly given how quick it is. From one thousand input explanatory variables, it correctly estimates the two variables that mattered. If that’s not impressive, then I don’t know what is.

9.2.3 This isn’t your first rodeo

Least angle regression (LAR) operates in much the same way as the lasso, because it’s the same package.

model <- lars(explanatory, response, type="lar")
plot(model)

signif.coefs(model)
#>  123  678 
#>  1.5 -0.5

So now is an excellent opportunity to show you the importance of scaling your variables. Everyone remembers that standard linear regression assumes independent, Normally-distributed variables, but what people often forget is that the optimisation routines inherent in essentially every single statistical solver (and often the mathematics itself) assume that your variables are on the same scale and centred. This means that each variable should have roughly the same standard deviation, and a mean of 0—your explanatory variables should be z-transformed. There are many reasons for this, but perhaps the most intuitive is that the maths/computing is going to be driven by a variable that varies over a much greater range, because all its coefficients are going to be much larger. Numerically, a computer might miss a coefficient whose value is \(10^{-23}\), even if when the variable is z-transformed it would be much “more significant” than a variable with a greater (but less important) coefficient. If you understood the importance of scaling variables in a PCA, then it might be helpful to know that the underlying logic and reasoning is the same in both cases. Anyway, you can see the effect for yourself below when I turn off lars’s default to automatically normalise all the variables. Notice that, because I simulated the data from a Standard Normal distribution, I have to set the threshold to 0 to see the effect.

bad.model <- lars(explanatory, response, type="lar", normalize=FALSE)
signif.coefs(bad.model, thresh=0) # Wow what a lot of coefficients!
#>  123  678 
#>  1.5 -0.5
signif.coefs(model, thresh=0)     # Nothing wrong here :D
#>  123  678 
#>  1.5 -0.5

…So maybe you should standardise those variables before you fit that regression your advisor is desperate for you to run, eh?

9.3 Support Vector Machines (SVM)

Support Vector Machines (SVMs) are useful for splitting data into groups. Thus, while they can be re-purposed for continuous response variables, you can think of them as the supervised complement to the clustering algorithms you already know. They are tremendously flexible and efficient, which makes them very useful in odd circumstances that tend to crop up a lot. For example, one-class variants of them are very popular in ecology because they can be used to predict where species should be found without any data on where they aren’t found. Equally, because they are so efficient, until quite recently they were used in smartphones to detect when someone wanted their phone to pay attention to them 4

9.3.1 An informal introduction to SVMs

Imagine you have data that are split into two categories (classes): an example is drawn in Figure 9.2 with light-blue (on the left) and dark (on the right) categories. If you can draw a single line—which we’ll call a separating hyperplane5—separating the two classes, then you could probably draw an infinite number of them by just tweaking the line a little up/down or altering its slope somewhat. You could, however, find a single line that has the furthest minimum distance to the points in the data; let’s agree, for the sake of argument, that such a line is the “best” line and is called the maximal margin hyperplane, and that minimal distance is the margin. The crazy thing about this line/hyperplane/whatever, is that it doesn’t really depend on any of the data other than the points lying on the margin6, which is wonderful because it means that no matter how big our dataset is we can just focus on these and we’ll be alright. We call these important points the support vectors of our classifying hyperplane, and identifying these points—and so the hyperplane that we can use to classify our data—is the sole goal of a support vector machine.

What’s really exciting about support vector machines (SVMs) is that they don’t just have to be fit with straight lines. Such linear SVMs are, of course, quite common, but we can fit essentially any kind of fancy complex classifier we want. Instead of drawing a straight line, we could draw a polynomial—or even a circle—by specifying a different equation for that line, which we call a different kernel7. These non-linear (i.e., not straight!) kernels are really useful and let us fit all classify all kinds of whacky datasets, as I give an example of in Figure 9.3.

I’m not going to give you anything more profound a statistical insight about this than what I’ve already given you. SVMs are extremely complex, and it would be a disservice to claim that I can teach your their theoretical underpinings in a single class. You do now, however, know enough about them to use them and, critically, to understand when they perform well.

Figure 9.2: A simple overview of Support Vector Machines. If we want to find a way to classify these two classes of data, we need to draw a line—a separating hyperplane—somewhere through this graph. To find the best line, we pick the one that has the smallest minimum distance from all the points—that distance defines the margin—and the points that lie on that margin are the support vectors. So a support vector machine is something that classifies data using some kind of hyperplane, and that hyperplane is defined using the support vectors themselves. This figure is taken from the vignette for the R package e1071, which you will be using in class.
Figure 9.3: SVMs with non-linear kernels. Not everything needs to be drawn on a straight line, and it’s possible to have polynomial (left) and radial (right) kernels. These figures are taken from An introduction to statistical learning by James et al. (2013).

9.3.2 Practical examples

To get a feeling for how SVMs work, let’s simulate some data. We’re going to make a dataset with two random variables, and then set make a set of points “in the middle” as somehow different. If that doesn’t make a great deal of sense to you right now, don’t worry: just look at the plot you’ll get out at the end and things will make sense.

data <- data.frame(temp=rnorm(1000), humid=rnorm(1000))
data$y = (rowSums(data) > (median(rowSums(data)) - 1)) &
  (rowSums(data) < (median(rowSums(data)) + 1))
with(data, plot(humid ~ temp, pch=20, col=ifelse(y, "red", "black")))

You could pretend, if you wanted, that the y variable here represented whether or not a particular species tends to be found in a particular climatic envelope. Pretend the x-axis in your plot represents the temperature of a series of sites, and the y-axis represents the humidity.

Now let’s see if we can fit an SVM, in this case using the default radial kernel but it could be with anything else you chose (e.g., a simple linear kernel), and see if we can classify between these two different classes.

library(e1071)
training <- sample(nrow(data), nrow(data)/2)
model <- svm(y~., data=data[training,], type="C")
plot(model, data[training,])

Success! In our training dataset of 50% of the data, it seems our model does a good job. Note that I’m stating that on the basis of the plot, which you should generate for yourself in R. Note how the SVM has identified the region in the center as different, and, according to the colors and the symbols, seems to have predicted the training data well. See if you can figure out, looking at the plot for yourself, what is going on. There are many kinds of SVM, and I’m only teaching you one of them in this class. Notice that I can make double-sure that svm fits a classification SVM to our data (type="c"): this function can be a bit picky, and you may find (for example, in the exercises) that you need to force svm to fit the kind of model you want—this is how to do that. Of course, to fit our model properly, we would want to test our model on some independent data that were not used to train it. So let’s see what happens when we do that:

table(predict(model, data[-training,]), data$y[-training])
#>        
#>         FALSE TRUE
#>   FALSE   238    1
#>   TRUE     10  251

Success! Even with data that wasn’t used to fit the SVM, we correctly identified the overwhelming majority of cases. I can tell this because, in this contingency table, most of the data are in the diagonals: the model predicted TRUE most often when the data were TRUE, and FALSE when the data were FALSE.

If you’re interested, you might want to read about the parameters that actually underlie the way the SVMs are fit. All SVMs include a cost parameter, which represents how much ‘wiggle room’ the SVM is allowed when trying to find the margin in comparison with error in the data. Another parameter is the \(\gamma\) (gamma) parameter, which you can think of as a non-linearity ‘fudge’ parameter. Larger values mean the kernel becomes more non-linear, which is a good thing if the data really are non-linear but a bad thing if they’re not as it opens you to over-fitting problems. Perhaps the best way to see if you \(\gamma\) parameter is leading you astray is to check the performance of your parameter on non-training data—but is that cheating?…

tune.svm(factor(y)~., data=data[-training,], gamma=c(.5,1,10), cost=c(1,10))
#> 
#> Parameter tuning of 'svm':
#> 
#> - sampling method: 10-fold cross validation 
#> 
#> - best parameters:
#>  gamma cost
#>    0.5   10
#> 
#> - best performance: 0.008

9.3.3 Details and other applications

There are two extensions to SVMs as presented here that are worth mentioning, but I don’t want you to worry about the details right now. The first is SVMs that can deal with multiple classes of data to predict: a military application of such data might be predicting whether a person is an enemy combatant, a friendly combatant, or a civilian (three classes). One approach for such data is to fit models that classify data into all the pairwise combinations of classes (e.g., enemy–friendly, enemy–civilian, and friendly–civilian)—this is one-versus-one classification. Whichever class is most frequently predicted is treated as the prediction. The second approach is to fit models predicting each class vs. some other class (e.g., enemy–other, friendly–other, civilian–other), measure which of the models most confidently predicted each particular piece of data under observation, and treat the most confident estimate as the correct one. This is one-versus-all prediction. These two options work better than you might imagine, although we’d probably all agree we’d rather not be anywhere near something like this being applied in the real world!

The second extension, which is quite common in ecology, is the one-class SVM. This has become an extremely popular method in ecology, where it is used to predict where a species might be found on the basis of environmental data. This is often called niche modeling, and while there are many statistical techniques that can be used when we have confirmed presences and absences of species on a landscape, dealing with the more common case where we only know where a species is, but do not know with much certainty where it isn’t, is much more difficult. This approach was pioneered by Guo et al. (2005; Ecological Modelling 182; 75–90) and is well-worth a read if you’re interested.

9.4 Exercises

Remember, you have two tasks. (1) To read through the material above, make sure you understand it, and that you can interpret and understand everything that is going on in the code examples. (2) Carry out the exercises below. Note also that they involve answering questions and thinking; you could write the code in seconds (as is the case for much of the work in machine learning), and so your task is to make sure you understand what is going on.

  1. The following questions require you to use a (slightly simplified) version of data from: Quinlan (1993) “Combining Instance-Based and Model-Based Learning” in Proceedings on the Tenth International Conference of Machine Learning, 236–243. This dataset describes the fuel efficiency of several kinds of car (mpg—miles per gallon), and can be described as a function of various other properties of the cars. The dataset is available on your course website as auto-mpg.txt.
  1. Fit a regression tree to a training subset of data.
  2. Plot out the resulting model, and explain in a few short sentences what the output shows.
  3. Validate your regression tree using independent data (i.e., not your training data).
  4. Fit a lasso regression to these data. Plot out your model and explain, in a few sentences, what the fitting process is showing you as it progresses.
  5. Now fit a LAR to these data. Compare your results with those of parts (a), (c), and (d) above. Which model do you find the easiest to interpret? Why?
  1. The following questions concern a dataset of forest fires (see http://www3.dsi.uminho.pt/pcortez/fires.pdf). Your task is to model the area of forest that was on fire as a function of various weather factors and some forest-fire indicators that the Canadian government uses. As your response variable is an area, you might want to employ a good-ole’-fashioned log10(area+1) transformation to it8.
  1. Fit a regression tree to this data. Explain what it means intuitively.
  2. Fit a bagged tree, a random forests model, and a boosted tree to a training subset of these data.
  3. Using your independent data, determine which of these models has performed the best.
  4. Fit both a lasso and a LAR model to these data. Explain their results, and contrast their differences.
  5. Fit a PCA to these data, extract the most significant terms, and fit a LAR to these data. How different is this model from the one you fit in (d)? Why do you think that is, and which model do you prefer?
  1. Many of the classical statistical methods we use today come from the brewing and distillation industries9. This dataset, which comes from a tasting tour of Scottish distillaries, was published by Lapointe & Legendre (1994; Applied Statistics, 237–257). You can find it on your course site as whisky.csv.
  1. Load the whisky dataset into R. Make sure the data is correctly formatted for an analysis in R (not column names that are names, etc.)
  2. The variable spey describes whether a whisky is, or is not, a “speyside” whisky. Run an SVM and, using training data, validate whether you can predict a whisky’s type on the basis of its characteristics.
  3. Train your model’s \(\gamma\) parameter, and see if you do a better job.
  4. How useful is this analysis? Is it something you could easily present to a policy-maker? Why (not)?
  5. If you were to do another analysis on this data (don’t!), what would it be and why?
  1. These questions work with the iris dataset that you might well already be very familiar with, and you can load using data(iris). If you’ve got ‘iris-fatigue’ then I would encourage you to run the exercises below using the Palmer Panguins dataset instead (see https://allisonhorst.github.io/palmerpenguins/); you will obviously have to make trivial changes to the example code but the dataset’s website should help with that.
  1. Perform an SVM on the iris dataset. Validate your model using training data.
  2. Train your model’s \(\gamma\) parameter. Does your model fit any better?
  3. It’s difficult to plot data with more than two dimensions, such as this iris dataset. See if you can figure out what the code below is doing, and do remember that you can always look at the help file for how to plot an svm. Make an informative plot of how your SVM is performing.
plot(model, iris, Sepal.Length ~ Sepal.Width, slice=
list(Petal.Width=median(iris$Petal.Width),Petal.Length=median(iris$Petal.Length))
)

  1. …well, erm, we’re continuing if you did the extension exercises. If you didn’t, don’t worry about it, you don’t need to have done those to follow this example.↩︎

  2. …and perhaps take a quick look over Chapter 1 again.↩︎

  3. Actually, we can’t use any standard calculus, and that’s the problem. See? I told you using the absolute value made things harder. Those of you who know about Lagrange multipliers will recognise that the lasso penalty is often called \(\lambda\) for a reason…↩︎

  4. “Hey Siri, record everything I’m saying and relay it to the government please”. For more details in the context of the evolution of bird song, see Pearse et al. (2018; Evolution 72(4):944–960).↩︎

  5. Separating because we’ve separated the data, and hyperplane is just a fancy word for a line↩︎

  6. If you think about it, you’ll have to agree that there must be at least two such points↩︎

  7. This is something of a simplification, but if you are familiar with how Generalized Linear Models use link functions to fit equations to transform data in order to fit it in a different parameter space, then you understand what’s going on. If you don’t: don’t worry about it, we’re just drawing lines with different shapes.↩︎

  8. Although, of course, as this is a statistics class I am required to remind you that you shouldn’t really do so ‘in the wild’: O’hara & Kotze (2010) Methods in Ecology and Evolution, 1(2), 118–122.↩︎

  9. If you’re interested, look up the history of “‘Student’s’ t-test”: https://en.wikipedia.org/wiki/William_Sealy_Gosset.↩︎