## Archive for the ‘**R**’ Category

## Tackling the John Smith Problem – deduplicating data via fuzzy matching in R

Last week I attended a CRM & data user group meeting for not-for-profits (NFPs), organized by my friend Yael Wasserman from Mission Australia. Following a presentation from a vendor, we broke up into groups and discussed common data quality issues that NFPs (and dare I say most other organisations) face. Number one on the list was the vexing issue of duplicate constituent (donor) records – henceforth referred to as dupes. I like to call this the *John Smith Problem* as it is likely that a typical customer database in a country with a large Anglo population is likely to have a fair number of records for customers with that name. The problem is tricky because one has to identify John Smiths who appear to be distinct in the database but are actually the same person, while also ensuring that one does not inadvertently merge two distinct John Smiths.

The John Smith problem is particularly acute for NFPs as much of their customer data comes in either via manual data entry or bulk loads with less than optimal validation. To be sure, all the NFPs represented at the meeting have some level of validation on both modes of entry, but all participants admitted that dupes tend to sneak in nonetheless…and at volumes that merit serious attention. Yael and his team have had some success in cracking the dupe problem using SQL-based matching of a combination of fields such as first name, last name and address or first name, last name and phone number and so on. However, as he pointed out, this method is limited because:

- It does not allow for typos and misspellings.
- Matching on too few fields runs the risk of false positives – i.e. labelling non-dupes as dupes.

The problems arise because SQL-based matching requires one to pre-specify match patterns. The solution is straightforward: use fuzzy matching instead. The idea behind fuzzy matching is simple: allow for inexact matches, assigning each match a *similarity score* ranging from 0 to 1 with 0 being complete dissimilarity and 1 being a perfect match. My primary objective in this article is to show how one can make headway with the John Smith problem using the fuzzy matching capabilities available in R.

### A bit about fuzzy matching

Before getting down to fuzzy matching, it is worth a brief introduction on how it works. The basic idea is simple: one has to generalise the notion of a match from a binary *“match” / “no match”* to allow for partial matching. To do this, we need to introduce the notion of an edit distance, which is essentially the *minimum number of operations required to transform one string into another*. For example, the edit distance between the strings *boy* and *bay* is 1: there’s only one edit required to transform one string to the other. The Levenshtein distance is the most commonly used edit distance. It is essentially, “*the minimum number of single-character edits (insertions, deletions or substitutions) required to change one word into the other.”*

A variant called the Damerau-Levenshtein distance, which additionally allows for the transposition of two adjacent characters (counted as one operation, not two) is found to be more useful in practice. We’ll use an implementation of this called the *optimal string alignment* (osa) distance. If you’re interested in finding out more about osa, check out the Damerau-Levenshtein article linked to earlier in this paragraph.

Since longer strings will potentially have larger numeric distances between them, it makes sense to normalise the distance to a value lying between 0 and 1. We’ll do this by dividing the calculated osa distance by the length of the larger of the two strings . Yes, this is crude but, as you will see, it works reasonably well. The resulting number is a normalised measure of the dissimilarity between the two strings. To get a similarity measure we simply subtract the dissimilarity from 1. So, a normalised dissimilarity of 1 translates to similarity score of 0 – i.e. the strings are perfectly *dissimilar*. I hope I’m not belabouring the point; I just want to make sure it is perfectly clear before going on.

### Preparation

In what follows, I assume you have R and RStudio installed. If not, you can access the software here and here for Windows and here for Macs; installation for both products is usually quite straightforward.

You may also want to download the Excel file **many_john_smiths** which contains records for ten fictitious John Smiths. At this point I should affirm that as far as the dataset is concerned, *a**ny resemblance to actual John Smiths, living or dead, is purely coincidental! *Once you have downloaded the file you will want to open it in Excel and examine the records and save it as a csv file in your R working directory (or any other convenient place) for processing in R.

As an aside, if you have access to a database, you may also want to load the file into a table called **many_john_smiths** and run the following dupe-detecting SQL statement:

select * from many_john_smiths t1

where exists

(select 'x' from many_john_smiths t2

where

t1.FirstName=t2.FirstName

and

t1.LastName=t2.LastName

and

t1.AddressPostcode=t2.AddressPostcode

and

t1.CustomerID <> t2.CustomerID)

You may also want to try matching on other column combinations such as *First/Last Name and AddressLine1* or F*irst/Last Name and AddressSuburb* for example. The limitations of column-based exact matching will be evident immediately. Indeed, I have deliberately designed the records to highlight some of the issues associated with dirty data: misspellings, typos, misheard names over the phone etc. A quick perusal of the records will show that there are *probably* two distinct John Smiths in the list. The problem is to quantify this observation. We do that next.

### Tackling the John Smith problem using R

We’ll use the following libraries: stringdist and stringi . The first library, stringdist, contains a bunch of string distance functions, we’ll use stringdistmatrix() which returns a matrix of pairwise string distances (osa by default) when passed a vector of strings, and stringi has a number of string utilities from which we’ll use str_length(), which returns the length of string.

OK, so on to the code. The first step is to load the required libraries:

`#load libraries`

library("stringdist")

library("stringr")

We then read in the data, ensuring that we override the annoying default behaviour of R, which is to convert strings to categorical variables – we want strings to remain strings!

`#read data, taking care to ensure that strings remain strings`

df <- read.csv("many_john_smiths.csv",stringsAsFactors = F)

#examine dataframe

str(df)

The output from str(df) (not shown) indicates that all columns barring *CustomerID* are indeed strings (i.e. type=character).

The next step is to find the length of each row:

#find length of string formed by each row (excluding title)

rowlen <- str_length(paste0(df$FirstName,df$LastName,df$AddressLine1,

df$AddressPostcode,df$AddressSuburb,df$Phone))

#examine row lengths

rowlen

> [1] 41 43 39 42 28 41 42 42 42 43

Note that I have excluded the *Title* column as I did not think it was relevant to determining duplicates.

Next we find the distance between every pair of records in the dataset. We’ll use the stringdistmatrix()function mentioned earlier:

#stringdistmatrix - finds pairwise osa distance between every pair of elements in a

#character vector

d <- stringdistmatrix(paste0(df$FirstName,df$LastName,df$AddressLine1,

df$AddressPostcode,df$AddressSuburb,df$Phone))

d

1 2 3 4 5 6 7 8 9

2 7

3 10 13

4 15 21 24

5 19 26 26 15

6 22 21 28 12 18

7 20 23 26 9 21 14

8 10 13 17 20 23 25 22

9 19 22 19 21 24 29 23 22

10 17 22 25 13 22 19 16 22 24

stringdistmatrix() returns an object of type dist (distance), which is essentially a vector of pairwise distances.

For reasons that will become clear later, it is convenient to normalise the distance – i.e. scale it to a number that lies between 0 and 1. We’ll do this by dividing the distance between two strings by the length of the longer string. We’ll use the nifty base R function combn() to compute the maximum length for every pair of strings:

#find the length of the longer of two strings in each pair

pwmax <- combn(rowlen,2,max,simplify = T)

The first argument is the vector from which combinations are to be generated, the second is the group size (2, since we want pairs) and the third argument indicates whether or not the result should be returned as an array (simplify=T) or list (simplify=F). The returned object, pwmax, is a one-dimensional array containing the pairwise maximum lengths. This has the same length and is organised in the same way as the object d returned by stringdistmatrix() (check that!). Therefore, to normalise d we simply divide it by pwmax

#normalised distance

dist_norm <- d/pwmax

The normalised distance lies between 0 and 1 (check this!) so we can define similarity as 1 minus distance:

#similarity = 1 - distance

similarity <- round(1-dist_norm,2)

sim_matrix <- as.matrix(similarity)

sim_matrix

1 2 3 4 5 6 7 8 9 10

1 0.00 0.84 0.76 0.64 0.54 0.46 0.52 0.76 0.55 0.60

2 0.84 0.00 0.70 0.51 0.40 0.51 0.47 0.70 0.49 0.49

3 0.76 0.70 0.00 0.43 0.33 0.32 0.38 0.60 0.55 0.42

4 0.64 0.51 0.43 0.00 0.64 0.71 0.79 0.52 0.50 0.70

5 0.54 0.40 0.33 0.64 0.00 0.56 0.50 0.45 0.43 0.49

6 0.46 0.51 0.32 0.71 0.56 0.00 0.67 0.40 0.31 0.56

7 0.52 0.47 0.38 0.79 0.50 0.67 0.00 0.48 0.45 0.63

8 0.76 0.70 0.60 0.52 0.45 0.40 0.48 0.00 0.48 0.49

9 0.55 0.49 0.55 0.50 0.43 0.31 0.45 0.48 0.00 0.44

10 0.60 0.49 0.42 0.70 0.49 0.56 0.63 0.49 0.44 0.00

The diagonal entries are 0, but that doesn’t matter because we know that every string is perfectly similar to itself! Apart from that, the similarity matrix looks quite reasonable: you can, for example, see that records 1 and 2 (similarity score=0.84) are quite similar while records 1 and 6 are quite dissimilar (similarity score=0.46). Now let’s extract some results more systematically. We’ll do this by printing out the top 5 non-diagonal similarity scores and the associated records for each of them. This needs a bit of work. To start with, we note that the similarity matrix (like the distance matrix) is symmetric so we’ll convert it into an upper triangular matrix to avoid double counting. We’ll also set the diagonal entries to 0 to avoid comparing a record with itself:

#convert to upper triangular to prevent double counting

sim_matrix[lower.tri(sim_matrix)] <- 0

#set diagonals to zero to avoid comparing row to itself

diag(sim_matrix) <- 0

Next we create a function that returns the n largest similarity scores and their associated row and column number – we’ll need the latter to identify the pair of records that are associated with each score:

#adapted from:

#https://stackoverflow.com/questions/32544566/find-the-largest-values-on-a-matrix-in-r

nlargest <- function(m, n) {

res <- order(m, decreasing = T)[seq_len(n)];

pos <- arrayInd(res, dim(m), useNames = TRUE);

list(values = m[res],

position = pos)

}

The function takes two arguments: a matrix m and a number n indicating the top n scores to be returned. Let’s set this number to 5 – i.e. we want the top 5 scores and the associated record indexes. We’ll store the output of nlargest in the variable sim_list:

top_n <- 5

sim_list <- nlargest(sim_matrix,top_n)

Finally, we loop through sim_list printing out the scores and associated records as we go along:

for (i in 1:top_n){

rec <- as.character(df[sim_list$position[i],])

sim_rec <- as.character(df[sim_list$position[i+top_n],])

cat("score: ",sim_list$values[i],"\n")

cat("record 1: ",rec,"\n")

cat ("record 2: ",sim_rec,"\n\n")

}

score: 0.84

record 1: 1 John Smith Mr 12 Acadia Rd Burnton 9671 1234 5678

record 2: 2 Jhon Smith Mr 12 Arcadia Road Bernton 967 1233 5678

score: 0.79

record 1: 4 John Smith Mr 13 Kynaston Rd Burnton 9671 34561234

record 2: 7 Jon Smith Mr. 13 Kinaston Rd Barnston 9761 36451223

score: 0.76

record 1: 1 John Smith Mr 12 Acadia Rd Burnton 9671 1234 5678

record 2: 3 J Smith Mr. 12 Acadia Ave Burnton 867`1 1233 567

score: 0.76

record 1: 1 John Smith Mr 12 Acadia Rd Burnton 9671 1234 5678

record 2: 8 John Smith Dr 12 Aracadia St Brenton 9761 12345666

score: 0.71

record 1: 4 John Smith Mr 13 Kynaston Rd Burnton 9671 34561234

record 2: 6 John S Dr. 12 Kinaston Road Bernton 9677 34561223

As you can see, the method correctly identifies close matches: there appear to be 2 distinct records (1 and 4) – and possibly more, depending on where one sets the similarity threshold. I’ll leave you to explore this further on your own.

### The John Smith problem in real life

As a proof of concept, I ran the following SQL on a real CRM database hosted on SQL Server:

select

FirstName+LastName,

count(*)

from

TableName

group by

FirstName+LastName

having

count(*)>100

order by

count(*) desc

I was gratified to note that John Smith did indeed come up tops – well over 200 records. I suspected there were a few duplicates lurking within, so I extracted the records and ran the above R code (with a few minor changes). I found there indeed were some duplicates! I also observed that the code ran with no noticeable degradation despite the dataset having well over 10 times the number of records used in the toy example above. I have not run it for larger datasets yet, but I suspect one will run into memory issues when the number of records gets into the thousands. Nevertheless, based on my experimentation thus far, this method appears viable for small datasets.

The problem of deduplicating large datasets is left as an exercise for motivated readers 😛

### Wrapping up

Often organisations will turn to specialist consultancies to fix data quality issues only to find that their work, besides being quite pricey, comes with a lot of caveats and cosmetic fixes that do not address the problem fully. Given this, there is a case to be made for doing as much of the exploratory groundwork as one can so that one gets a good idea of what can be done and what cannot. At the very least, one will then be able to keep one’s consultants on their toes. In my experience, the John Smith problem ranks right up there in the list of data quality issues that NFPs and many other organisations face. This article is intended as a starting point to address this issue using an easily available and cost effective technology.

Finally, I should reiterate that the approach discussed here is just one of many possible and is neither optimal nor efficient. Nevertheless, it works quite well on small datasets, and is therefore offered here as a starting point for your own attempts at tackling the problem. If you come up with something better – as I am sure you can – I’d greatly appreciate your letting me know via the contact page on this blog or better yet, a comment.

**Acknowledgements:**

I’m indebted to Homan Zhao and Sree Acharath for helpful conversations on fuzzy matching. I’m also grateful to all those who attended the NFP CRM and Data User Group meetup that was held earlier this month – the discussions at that meeting inspired this piece.

## An intuitive introduction to support vector machines using R – Part 1

About a year ago, I wrote a piece on support vector machines as a part of my *gentle introduction to data science R* series. So it is perhaps appropriate to begin this piece with a few words about my motivations for writing yet another article on the topic.

Late last year, a curriculum lead at DataCamp got in touch to ask whether I’d be interested in developing a course on SVMs for them.

My answer was, obviously, an enthusiastic “Yes!”

Instead of rehashing what I had done in my previous article, I thought it would be interesting to try an approach that focuses on building an intuition for how the algorithm works using examples of increasing complexity, supported by visualisation rather than math. This post is the first part of a two-part series based on this approach.

The article builds up some basic intuitions about support vector machines (abbreviated henceforth as SVM) and then focuses on linearly separable problems. Part 2 (to be released at a future date) will deal with radially separable and more complex data sets. The focus throughout is on developing an understanding *what* the algorithm does rather than the technical details of how it does it.

Prerequisites for this series are a basic knowledge of R and some familiarity with the *ggplot* package. However, even if you don’t have the latter, you should be able to follow much of what I cover so I encourage you to press on regardless.

<advertisement> if you have a DataCamp account, you may want to check out my course on support vector machines using R. Chapters 1 and 2 of the course closely follow the path I take in this article. </advertisement>

### A one dimensional example

A soft drink manufacturer has two brands of their flagship product: *Choke* (sugar content of 11g/100ml) and *Choke-R* (sugar content 8g/100 ml). The actual sugar content can vary quite a bit in practice so it can sometimes be hard to figure out the brand given the sugar content. Given sugar content data for 25 samples taken randomly from both populations (see file sugar_content.xls), our task is to come up with a decision rule for determining the brand.

Since this is one-variable problem, the simplest way to discern if the samples fall into distinct groups is through visualisation. Here’s one way to do this using *ggplot*:

…and here’s the resulting plot:

Note that we’ve simulated a one-dimensional plot by setting all the y values to 0.

From the plot, it is evident that the samples fall into distinct groups: low sugar content, bounded above by the 8.8 g/100ml sample and high sugar content, bounded below by the 10 g/100ml sample.

Clearly, any point that lies between the two points is an acceptable *decision boundary.* We could, for example, pick 9.1g/100ml and 9.7g/100ml. Here’s the R code with those points added in. Note that we’ve made the points a bit bigger and coloured them red to distinguish them from the sample points.

label=d_bounds$sep, size=2.5,

vjust=2, hjust=0.5, colour=”red”)

And here’s the plot:

Now, a bit about the decision rule. Say we pick the first point as the decision boundary, the decision rule would be:

Say we pick 9.1 as the decision boundary, our classifier (in R) would be:

The other one is left for you as an exercise.

Now, it is pretty clear that although either these points define an acceptable decision boundary, neither of them are the best. Let’s try to formalise our intuitive notion as to why this is so.

The *margin* is the distance between the points in both classes that are closest to the decision boundary. In case at hand, the margin is 1.2 g/100ml, which is the difference between the two extreme points at 8.8 g/100ml (Choke-R) and 10 g/100ml (Choke). It should be clear that the best separator is the one that lies halfway between the two extreme points. This is called the *maximum margin separator*. The maximum margin separator in the case at hand is simply the average of the two extreme points:

geom_point(data=mm_sep,aes(x=mm_sep$sep, y=c(0)), colour=”blue”, size=4)

And here’s the plot:

We are dealing with a one dimensional problem here so the decision boundary is a point. In a moment we will generalise this to a two dimensional case in which the boundary is a straight line.

Let’s close this section with some general points.

Remember this is a sample not the entire population, so it is quite possible (indeed likely) that there will be as yet unseen samples of Choke-R and Choke that have a sugar content greater than 8.8 and less than 10 respectively. So, the best classifier is one that lies at the greatest possible distance from both classes. The maximum margin separator is that classifier.

This toy example serves to illustrate the main aim of SVMs, which is to find an optimal separation boundary in the sense described here. However, doing this for real life problems is not so simple because life is not one dimensional. In the remainder of this article and its yet-to-be-written sequel, we will work through examples of increasing complexity so as to develop a good understanding of how SVMs work in addition to practical experience with using the popular SVM implementation in R.

<Advertisement> Again, for those of you who have DataCamp premium accounts, here is a course that covers pretty much the entire territory of this two part series. </Advertisement>

### Linearly separable case

The next level of complexity is a two dimensional case (2 predictors) in which the classes are separated by a straight line. We’ll create such a dataset next.

Let’s begin by generating 200 points with attributes x1 and x2, randomly distributed between 0 and 1. Here’s the R code:

Let’s visualise the generated data using a scatter plot:

And here’s the plot

Now let’s classify the points that lie above the line *x1=x2* as belonging to the class +1 and those that lie below it as belonging to class -1 (the class values are arbitrary choices, I could have chosen them to be anything at all). Here’s the R code:

Let’s modify the plot in Figure 4, colouring the points classified as +1n blue and those classified -1 red. For good measure, let’s also add in the decision boundary. Here’s the R code:

Note that the parameters in *geom_abline()* are derived from the fact that the line* x1=x2* has slope 1 and y intercept 0.

Here’s the resulting plot:

Next let’s introduce a margin in the dataset. To do this, we need to exclude points that lie within a specified distance of the boundary. A simple way to approximate this is to exclude points that have *x1* and *x2* values that differ by less a pre-specified value, *delta*. Here’s the code to do this with *delta* set to 0.05 units.

The check on the number of datapoints tells us that a number of points have been excluded.

Running the *previous* *ggplot* code block yields the following plot which clearly shows the reduced dataset with the depopulated region near the decision boundary:

Let’s add the margin boundaries to the plot. We know that these are parallel to the decision boundary and lie delta units on either side of it. In other words, the margin boundaries have slope=1 and y intercepts *delta* and –*delta*. Here’s the *ggplot* code:

And here’s the plot with the margins:

OK, so we have constructed a dataset that is *linearly separable*, which is just a short code for saying that the classes can be separated by a straight line. Further, the dataset has a margin, i.e. there is a “gap” so to speak, between the classes. Let’s save the dataset so that we can use it in the next section where we’ll take a first look at the *svm()* function in the *e1071* package.

That done, we can now move on to…

### Linear SVMs

Let’s begin by reading in the datafile we created in the previous section:

We then split the data into training and test sets using an 80/20 random split. There are many ways to do this. Here’s one:

The next step is to build the an SVM classifier model. We will do this using the *svm() f*unction which is available in the *e1071* package. The *svm()* function has a range of parameters. I explain some of the key ones below, in particular, the following parameters: *type*, *cost*, *kernel* and *scale*. It is recommended to have a browse of the documentation for more details.

The *type* parameter specifies the algorithm to be invoked by the function. The algorithm is capable of doing both classification and regression. We’ll focus on classification in this article. Note that there are two types of classification algorithms, nu and C classification. They essentially differ in the way that they penalise margin and boundary violations, but can be shown to lead to equivalent results. We will stick with C classification as it is more commonly used. The “C” refers to the cost which we discuss next.

The *cost* parameter specifies the penalty to be applied for boundary violations. This parameter can vary from 0 to infinity (in practice a large number compared to 0, say 10^6 or 10^8). We will explore the effect of varying *cost* later in this piece. To begin with, however, we will leave it at its default value of 1.

The *kernel* parameter specifies the kind of function to be used to construct the decision boundary. The options are linear, polynomial and radial. In this article we’ll focus on linear kernels as we know the decision boundary is a straight line.

The *scale* parameter is a Boolean that tells the algorithm whether or not the datapoints should be scaled to have zero mean and unit variance (i.e. shifted by the mean and scaled by the standard deviation). Scaling is generally good practice to avoid undue influence of attributes that have unduly large numeric values. However, in this case we will avoid scaling as we know the attributes are bounded and (more important) we would like to plot the boundary obtained from the algorithm manually.

Building the model is a simple one-line call, setting appropriate values for the parameters:

We expect a linear model to perform well here since the dataset it is linear by construction. Let’s confirm this by calculating training and test accuracy. Here’s the code:

The perfect accuracies confirm our expectation. However, accuracies by themselves are misleading because the story is somewhat more nuanced. To understand why, let’s plot the predicted decision boundary and margins using *ggplot*. To do this, we have to first extract information regarding these from the svm model object. One can obtain summary information for the model by typing in the model name like so:

kernel = “linear”, scale = FALSE)

Which outputs the following: the function *call*, SVM *type*, *kernel* and *cost* (which is set to its default). In case you are wondering about *gamma, * although it’s set to 0.5 here, it plays no role in linear SVMs. We’ll say more about it in the sequel to this article in which we’ll cover more complex kernels. More interesting are the *support vectors*. In a nutshell, these are training dataset points that *specify the location of the decision boundary*. We can develop a better understanding of their role by visualising them. To do this, we need to know their coordinates and indices (position within the dataset). This information is stored in the SVM model object. Specifically, the *index* element of *svm_model* contains the indices of the training dataset points that are support vectors and the *SV* element lists the coordinates of these points. The following R code lists these explicitly (Note that I’ve not shown the outputs in the code snippet below):

Let’s use the indices to visualise these points in the training dataset. Here’s the ggplot code to do that:

And here is the plot:

We now see that the support vectors are clustered around the boundary and, in a sense, serve to define it. We will see this more clearly by plotting the predicted decision boundary. To do this, we need its slope and intercept. These aren’t available directly available in the *svm_model*, but they can be extracted from the *coefs*, *SV* and *rho* elements of the object.

The first step is to use *coefs* and the support vectors to build the what’s called the *weight vector*. The *weight vector* is given by the product of the *coefs* matrix with the matrix containing the SVs. Note that the fact that only the support vectors play a role in defining the boundary is consistent with our expectation that the boundary should be fully specified by them. Indeed, this is often touted as a feature of SVMs in that it is one of the few classifiers that depends on only a small subset of the training data, i.e. the datapoints closest to the boundary rather than the entire dataset.

Once we have the weight vector, we can calculate the slope and intercept of the predicted decision boundary as follows:

Note that the slope and intercept are quite different from the correct values of 1 and 0 (reminder: the actual decision boundary is the line *x1=x2* by construction). We’ll see how to improve on this shortly, but before we do that, let’s plot the decision boundary using the slope and intercept we have just calculated. Here’s the code:

And here’s the augmented plot:

The plot clearly shows how the support vectors “support” the boundary – indeed, if one draws line segments from each of the points to the boundary in such a way that the intersect the boundary at right angles, the lines can be thought of as “holding the boundary in place”. Hence the term *support vector*.

This is a good time to mention that the *e1071* library provides a built-in plot method for *svm* function. This is invoked as follows:

The svm *plot *function takes a formula specifying the plane on which the boundary is to be plotted. This is not necessary here as we have only two predictors (x1 and x2) which automatically define a plane.

Here is the plot generated by the above code:

Note that the axes are switched (x1 is on the y axis). Aside from that, the plot is reassuringly similar to our *ggplot* version in Figure 9. Also note that that the support vectors are marked by “x”. Unfortunately the built in function does not display the margin boundaries, but this is something we can easily add to our home-brewed plot. Here’s how. We know that the margin boundaries are parallel to the decision boundary, so all we need to find out is their intercept. It turns out that the intercepts are offset by an amount *1/w[2]* units on either side of the decision boundary. With that information in hand we can now write the the code to add in the margins to the plot shown in Figure 9. Here it is:

geom_abline(slope=slope_1,intercept = intercept_1+1/w[2], linetype=”dashed”)

And here is the plot with the margins added in:

Note that the predicted margins are much wider than the actual ones (compare with Figure 7). As a consequence, many of the support vectors lie within the predicted margin – that is, they violate it. The upshot of the wide margin is that the decision boundary is not tightly specified. This is why we get a significant difference between the slope and intercept of predicted decision boundary and the actual one. We can sharpen the boundary by narrowing the margin. How do we do this? We make margin violations more expensive by increasing the *cost*. Let’s see this margin-narrowing effect in action by building a model with *cost = *100 on the same training dataset as before. Here is the code:

I’ll leave you to calculate the training and test accuracies (as one might expect, these will be perfect).

Let’s inspect the *cost=100* model:

kernel = “linear”,cost=100, scale = FALSE)

The number of support vectors is reduced from 55 to 6! We can plot these and the boundary / margin lines using *ggplot* as before. The code is identical to the previous case (see code block preceding Figure 8). If you run it, you will get the plot shown in Figure 12.

Since the boundary is more tightly specified, we would expect the slope and intercept of the predicted boundary to be considerably closer to their actual values of 1 and 0 respectively (as compared to the default cost case). Let’s confirm that this is so by calculating the slope and intercept as we did in the code snippets preceding Figure 9. Here’s the code:

Which nicely confirms our expectation.

The decision boundary and margins for the high cost case can also be plotted with the code shown earlier. Her it is for completeness:

geom_abline(slope=slope_100,intercept = intercept_100+1/w[2], linetype=”dashed”)

And here’s the plot:

SVMs that allow margin violations are called *soft margin classifiers* and those that do not are called *hard*. In this case, the hard margin classifier does a better job because it specifies the boundary more accurately than its soft counterpart. However, this does not mean that hard margin classifier are to be preferred over soft ones in all situations. Indeed, in real life, where we usually do not know the shape of the decision boundary upfront, soft margin classifiers can allow for a greater degree of uncertainty in the decision boundary thus improving generalizability of the classifier.

OK, so now we have a good feel for what the SVM algorithm does in the linearly separable case. We will round out this article by looking at a real world dataset that fortuitously turns out to be almost linearly separable: the famous (notorious?) iris dataset. It is instructive to look at this dataset because it serves to illustrate another feature of the *e1071* SVM algorithm – its capability to handle classification problems that have more than 2 classes.

### A multiclass problem

The iris dataset is well-known in the machine learning community as it features in many introductory courses and tutorials. It consists of 150 observations of 3 *species* of the iris flower – *setosa*, *versicolor* and *virginica*. Each observation consists of numerical values for 4 independent variables (predictors): petal length, petal width, sepal length and sepal width. The dataset is available in a standard installation of R as a built in dataset. Let’s read it in and examine its structure:

Now, as it turns out, petal length and petal width are the key determinants of species. So let’s create a scatterplot of the datapoints as a function of these two variables (i.e. project each data point on the petal length-petal width plane). We will also distinguish between species using different colour. Here’s the ggplot code to do this:

And here’s the plot:

On this plane we see a clear linear boundary between *setosa* and the other two species, *versicolor* and *virginica*. The boundary between the latter two is almost linear. Since there are four predictors, one would have to plot the other combinations to get a better feel for the data. I’ll leave this as an exercise for you and move on with the assumption that the data is nearly linearly separable. If the assumption is grossly incorrect, a linear SVM will not work well.

Up until now, we have discussed binary classification problem, i.e. those in which the predicted variable can take on only two values. In this case, however, the predicted variable, *Species*, can take on 3 values (setosa, versicolor and virginica). This brings up the question as to how the algorithm deals multiclass classification problems – i.e those involving datasets with more than two classes. The SVM algorithm does this using a *one-against-one *classification strategy. Here’s how it works:

- Divide the dataset (assumed to have N classes) into N(N-1)/2 datasets that have two classes each.
- Solve the binary classification problem for each of these subsets
- Use a simple voting mechanism to assign a class to each data point.

Basically, each data point is assigned the most frequent classification it receives from all the binary classification problems it figures in.

With that said, let’s get on with building the classifier. As before, we begin by splitting the data into training and test sets using an 80/20 random split. Here is the code to do this:

Then we build the model (default cost) and examine it:

The main thing to note is that the function call is identical to the binary classification case. We get some basic information about the model by typing in the model name as before:

kernel = “linear”)

And the train and test accuracies are computed in the usual way:

This looks good, but is potentially misleading because it is for a particular train/test split. Remember, in this case, unlike the earlier example, we do not know the shape of the actual decision boundary. So, to get a robust measure of accuracy, we should calculate the average test accuracy over a number of train/test partitions. Here’s some code to do that:

Which is not too bad at all, indicating that the dataset is indeed nearly linearly separable. If you try different values of *cost* you will see that it does not make much difference to the average accuracy.

This is a good note to close this piece on. Those who have access to DataCamp premium courses will find that the content above is covered in chapters 1 and 2 of the course on support vector machines in R. The next article in this two-part series will cover chapters 3 and 4.

## Summarising

My main objective in this article was to help develop an intuition for how SVMs work in simple cases. We illustrated the basic principles and terminology with a simple 1 dimensional example and then worked our way to linearly separable binary classification problems with multiple predictors. We saw how the latter can be solved using a popular svm implementation available in R. We also saw that the algorithm can handle multiclass problems. All through, we used visualisations to see what the algorithm does and how the key parameters affect the decision boundary and margins.

In the next part (yet to be written) we will see how SVMs can be generalised to deal with complex, nonlinear decision boundaries. In essence, the use a mathematical trick to “linearise” these boundaries. We’ll delve into details of this trick in an intuitive, visual way as we have done here.

Many thanks for reading!

## A gentle introduction to data visualisation using R

Data science students often focus on machine learning algorithms, overlooking some of the more routine but important skills of the profession. I’ve lost count of the number of times I have advised students working on projects for industry clients to curb their keenness to code and work on understanding the data first. This is important because, as people (ought to) know, data doesn’t speak for itself, it has to be given a voice; and as data-scarred professionals know from hard-earned experience, one of the best ways to do this is through visualisation.

Data visualisation is sometimes (often?) approached as a bag of tricks to be learnt individually, with no little or no reference to any underlying principles. Reading Hadley Wickham’s paper on the grammar of graphics was an epiphany; it showed me how different types of graphics can be constructed in a consistent way using common elements. Among other things, the grammar makes visualisation a logical affair rather than a set of tricks. This post is a brief – and hopefully logical – introduction to visualisation using ggplot2, Wickham’s implementation of a grammar of graphics.

In keeping with the practical bent of this series we’ll focus on worked examples, illustrating elements of the grammar as we go along. We’ll first briefly describe the elements of the grammar and then show how these are used to build different types of visualisations.

### A grammar of graphics

Most visualisations are constructed from common elements that are pieced together in prescribed ways. The elements can be grouped into the following categories:

**Data**– this is obvious, without data there is no story to tell and definitely no plot!**Mappings**– these are correspondences between data and display elements such as spatial location, shape or colour. Mappings are referred to as*aesthetics*in Wickham’s grammar.**Scales**– these are transformations (conversions) of data values to numbers that can be displayed on-screen. There should be one scale per mapping. ggplot typically does the scaling transparently, without users having to worry about it. One situation in which you might need to mess with default scales is when you want to zoom in on a particular range of values. We’ll see an example or two of this later in this article.**Geometric objects**– these specify the geometry of the visualisation. For example, in ggplot2 a scatter plot is specified via a*point*geometry whereas a fitting curve is represented by a*smooth*geometry. ggplot2 has a range of geometries available of which we will illustrate just a few.**Coordinate system –**this specifies the system used to position data points on the graphic. Examples of coordinate systems are Cartesian and polar. We’ll deal with Cartesian systems in this tutorial. See this post for a nice illustration of how one can use polar plots creatively.**Facets**– a facet specifies how data can be split over multiple plots to improve clarity. We’ll look at this briefly towards the end of this article.

The basic idea of a *layered* grammar of graphics is that each of these elements can be combined – literally added layer by layer – to achieve a desired visual result. Exactly how this is done will become clear as we work through some examples. So without further ado, let’s get to it.

### Hatching (gg)plots

In what follows we’ll use the NSW Government Schools dataset, made available via the state government’s open data initiative. The data is in csv format. If you cannot access the original dataset from the aforementioned link, you can download an Excel file with the data **here** (remember to save it as a csv before running the code!).

The first task – assuming that you have a working R/RStudio environment – is to load the data into R. To keep things simple we’ll delete a number of columns (as shown in the code) and keep only rows that are complete, i.e. those that have no missing values. Here’s the code:

A note regarding the last line of code above, a couple of schools have “np” entered for the* student_number* variable. These are coerced to NA in the numeric conversion. The last line removes these two schools from the dataset.

Apart from *student numbers* and location data, we have retained *level of schooling* (primary, secondary etc.) and* ICSEA ranking*. The location information includes attributes such as *suburb*, *postcode*, *region*, *remoteness* as well as *latitude* and *longitude*. We’ll use only *remoteness* in this post.

The first thing that caught my eye in the data was was the ICSEA ranking. Before going any further, I should mention that the Australian Curriculum Assessment and Reporting Authority (the organisation responsible for developing the ICSEA system) emphasises that the score is *not* a school ranking, but a measure of *socio-educational advantage* of the student population in a school. Among other things, this is related to family background and geographic location. The average ICSEA score is set at an average of 1000, which can be used as a reference level.

I thought a natural first step would be to see how ICSEA varies as a function of the other variables in the dataset such as *student number*s and *location* (*remoteness*, for example). To begin with, let’s plot ICSEA rank as a function of student number. As it is our first plot, let’s take it step by step to understand how the layered grammar works. Here we go:

This displays a blank plot because we have not specified a *mapping* and *geometry *to go with the data. To get a plot we need to specify both. Let’s start with a scatterplot, which is specified via a point geometry. Within the geometry function, variables are mapped to visual properties of the using aesthetic mappings. Here’s the code:

The resulting plot is shown in Figure 1.

At first sight there are two points that stand out: 1) there are fewer number of large schools, which we’ll look into in more detail later and 2) larger schools seem to have a higher ICSEA score on average. To dig a little deeper into the latter, let’s add a linear trend line. We do that by adding another layer (geometry) to the scatterplot like so:

The result is shown in Figure 2.

The *lm* method does a linear regression on the data. The shaded area around the line is the 95% confidence level of the regression line (i.e that it is 95% certain that the true regression line lies in the shaded region). Note that *geom_smooth * provides a range of smoothing functions including generalised linear and local regression (loess) models.

You may have noted that we’ve specified the aesthetic mappings in both *geom_point* and *geom_smooth*. To avoid this duplication, we can simply specify the mapping, once in the top level ggplot call (the first layer) like so:

geom_point()+

geom_smooth(method=”lm”)

From Figure 2, one can see a clear positive correlation between student numbers and ICSEA scores, let’s zoom in around the average value (1000) to see this more clearly…

The *coord_cartesian* function is used to zoom the plot to without changing any other settings. The result is shown in Figure 3.

To make things clearer, let’s add a reference line at the average:

The result, shown in Figure 4, indicates quite clearly that larger schools tend to have higher ICSEA scores. That said, there is a twist in the tale which we’ll come to a bit later.

As a side note, you would use *geom_vline* to zoom in on a specific range of x values and *geom_abline* to add a reference line with a specified slope and intercept. See this article on ggplot reference lines for more.

OK, now that we have seen how *ICSEA scores* vary with *student number*s let’s switch tack and incorporate another variable in the mix. An obvious one is *remoteness*. Let’s do a scatterplot as in Figure 1, but now colouring each point according to its remoteness value. This is done using the colour aesthetic as shown below:

geom_point()

The resulting plot is shown in Figure 5.

Aha, a couple of things become apparent. First up, large schools tend to be in metro areas, which makes good sense. Secondly, it appears that metro area schools have a distinct socio-educational advantage over regional and remote area schools. Let’s add trendlines by remoteness category as well to confirm that this is indeed so:

The plot, which is shown in Figure 6, indicates clearly that ICSEA scores decrease on the average as we move away from metro areas.

Moreover, larger schools metropolitan areas tend to have higher than average scores (above 1000), regional areas tend to have lower than average scores overall, with remote areas being markedly more disadvantaged than both metro and regional areas. This is no surprise, but the visualisations show just how stark the differences are.

It is also interesting that, in contrast to metro and (to some extent) regional areas, there *negative correlation* between student numbers and scores for remote schools. One can also use local regression to get a better picture of how ICSEA varies with student numbers and remoteness. To do this, we simply use the *loess* method instead of *lm:*

geom_point() + geom_hline(yintercept=1000) + geom_smooth(method=”loess”)

The result, shown in Figure 7, has a number of interesting features that would have been worth pursuing further were we analysing this dataset in a real life project. For example, why do small schools tend to have lower than benchmark scores?

From even a casual look at figures 6 and 7, it is clear that the confidence intervals for remote areas are huge. This suggests that the number of datapoints for these regions are a) small and b) very scattered. Let’s quantify the number by getting counts using the *table* function (I know, we could plot this too…and we will do so a little later). We’ll also transpose the results using *data.frame* to make them more readable:

The number of datapoints for remote regions is *much* less than those for metro and regional areas. Let’s repeat the loess plot with only the two remote regions. Here’s the code:

geom_point() + geom_hline(yintercept=1000) + geom_smooth(method=”loess”)

The plot, shown in Figure 8, shows that there is indeed a huge variation in the (small number) of datapoints, and the confidence intervals reflect that. An interesting feature is that some small remote schools have above average scores. If we were doing a project on this data, this would be a feature worth pursuing further as it would likely be of interest to education policymakers.

Note that there is a difference in the x axis scale between Figures 7 and 8 – the former goes from 0 to 2000 whereas the latter goes up to 400 only. So for a fair comparison, between remote and other areas, you may want to re-plot Figure 7, zooming in on student numbers between 0 and 400 (or even less). This will also enable you to see the complicated dependence of scores on student numbers more clearly across all regions.

We’ll leave the scores vs student numbers story there and move on to another geometry – the well-loved bar chart. The first one is a visualisation of the remoteness category count that we did earlier. The relevant geometry function is *geom_bar*, and the code is as easy as:

The plot is shown in Figure 9.

The category labels on the x axis are too long and look messy. This can be fixed by tilting them to a 45 degree angle so that they don’t run into each other as they most likely did when you ran the code on your computer. This is done by modifying the *axis.text* element of the plot theme. Additionally, it would be nice to get counts on top of each category bar. The way to do that is using another geometry function, geom_text. Here’s the code incorporating the two modifications:

theme(axis.text.x=element_text(angle=45, hjust=1))

The result is shown in Figure 10.

Some things to note: : *stat=count* tells ggplot to compute counts by category and the aesthetic *label = ..count.. *tells ggplot to access the internal variable that stores those counts. The the vertical justification setting, *vjust=-1*, tells ggplot to display the counts on top of the bars. Play around with different values of *vjust* to see how it works. The code to adjust label angles is self explanatory.

It would be nice to reorder the bars by frequency. This is easily done via *fct_infreq* function in the *forcats* package like so:

geom_bar(mapping = aes(x=fct_infreq(ASGS_remoteness)))+

theme(axis.text.x=element_text(angle=45, hjust=1))

The result is shown in Figure 11.

To reverse the order, invoke *fct_rev,* which reverses the sort order:

geom_bar(mapping = aes(x=fct_rev(fct_infreq(ASGS_remoteness))))+

theme(axis.text.x=element_text(angle=45, hjust=1))

The resulting plot is shown in Figure 12.

If this is all too grey for us, we can always add some colour. This is done using the *fill* aesthetic as follows:

geom_bar(mapping = aes(x=ASGS_remoteness, fill=ASGS_remoteness))+

theme(axis.text.x=element_text(angle=45, hjust=1))

The resulting plot is shown in Figure 13.

Note that, in the above, that we have mapped fill and x to the *same* variable, *remoteness *which makes the legend superfluous. I will leave it to you to figure out how to suppress the legend – Google is your friend.

We could also map fill to another variable, which effectively adds another dimension to the plot. Here’s how:

geom_bar(mapping = aes(x=ASGS_remoteness, fill=level_of_schooling))+

theme(axis.text.x=element_text(angle=45, hjust=1))

The plot is shown in Figure 14. The new variable, level of schooling, is displayed via proportionate coloured segments stacked up in each bar. The default stacking is one on top of the other.

Alternately, one can stack them up side by side by setting the *position *argument to *dodge* as follows:

geom_bar(mapping = aes(x=ASGS_remoteness,fill=level_of_schooling),position =”dodge”)+

theme(axis.text.x=element_text(angle=45, hjust=1))

The plot is shown in Figure 15.

Finally, setting the *position* argument to *fill* normalises the bar heights and gives us the proportions of* level of schooling* for each *remoteness* category. That sentence will make more sense when you see Figure 16 below. Here’s the code, followed by the figure:

geom_bar(mapping = aes(x=ASGS_remoteness,fill=level_of_schooling),position = “fill”)+

theme(axis.text.x=element_text(angle=45, hjust=1))

Obviously, we lose frequency information since the bar heights are normalised.

An interesting feature here is that the proportion of central and community schools increases with remoteness. Unlike primary and secondary schools, central / community schools provide education from Kindergarten through Year 12. As remote areas have smaller numbers of students, it makes sense to consolidate educational resources in institutions that provide schooling at all levels .

Finally, to close the loop so to speak, let’s revisit our very first plot in Figure 1 and try to simplify it in another way. We’ll use faceting to split it out into separate plots, one per remoteness category. First, we’ll organise the subplots horizontally using *facet_grid*:

facet_grid(~ASGS_remoteness)

The plot is shown in Figure 17 in which the different remoteness categories are presented in separate plots (facets) against a common y axis. It shows, the sharp differences between student numbers between remote and other regions.

To get a vertically laid out plot, switch the faceted variable to other side of the formula (left as an exercise for you).

If one has too many categories to fit into a single row, one can wrap the facets using *facet_wrap* like so:

geom_point(mapping = aes(x=student_number,y=ICSEA_Value))+

facet_wrap(~ASGS_remoteness, ncol= 2)

The resulting plot is shown in Figure 18.

One can specify the number of rows instead of columns. I won’t illustrate that as the change in syntax is quite obvious.

…and I think that’s a good place to stop.

### Wrapping up

Data visualisation has a reputation of being a dark art, masterable only by the visually gifted. This may have been partially true some years ago, but in this day and age it definitely isn’t. Versatile packages such as ggplot, that use a consistent syntax have made the art much more accessible to visually ungifted folks like myself. In this post I have attempted to provide a brief and (hopefully) logical introduction to ggplot. In closing I note that although some of the illustrative examples violate the principles of good data visualisation, I hope this article will serve its primary purpose which is pedagogic rather than artistic.

**Further reading**:

Where to go for more? Two of the best known references are Hadley Wickham’s books:

I highly recommend his R for Data Science , available online here. Apart from providing a good overview of ggplot, it is an excellent introduction to R for data scientists. If you haven’t read it, do yourself a favour and buy it now.

People tell me his ggplot book is an excellent book for those wanting to learn the ins and outs of ggplot . I have not read it myself, but if his other book is anything to go by, it should be pretty damn good.

## A gentle introduction to logistic regression and lasso regularisation using R

In this day and age of artificial intelligence and deep learning, it is easy to forget that simple algorithms can work well for a surprisingly large range of practical business problems. And the simplest place to start is with the granddaddy of data science algorithms: linear regression and its close cousin, logistic regression. Indeed, in his acclaimed MOOC and accompanying textbook, Yaser Abu-Mostafa spends a good portion of his time talking about linear methods, and with good reason too: linear methods are not only a good way to learn the key principles of machine learning, they can also be remarkably helpful in zeroing in on the most important predictors.

My main aim in this post is to provide a beginner level introduction to logistic regression using R and also introduce LASSO (Least Absolute Shrinkage and Selection Operator), a powerful feature selection technique that is very useful for regression problems. Lasso is essentially a regularization method. If you’re unfamiliar with the term, think of it as a way to reduce overfitting using less complicated functions (and if that means nothing to you, check out my prelude to machine learning). One way to do this is to toss out less important variables, after checking that they aren’t important. As we’ll discuss later, this can be done manually by examining p-values of coefficients and discarding those variables whose coefficients are not significant. However, this can become tedious for classification problems with many independent variables. In such situations, lasso offers a neat way to model the dependent variable while *automagically* *selecting significant variables by* *shrinking the coefficients of unimportant predictors to zero. A*ll this without having to mess around with p-values or obscure information criteria. How good is that?

### Why not linear regression?

In linear regression one attempts to model a dependent variable (i.e. the one being predicted) using the best straight line fit to a set of predictor variables. The best fit is usually taken to be one that minimises the root mean square error, which is the sum of square of the differences between the actual and predicted values of the dependent variable. One can think of logistic regression as the equivalent of linear regression for a classification problem. In what follows we’ll look at binary classification – i.e. a situation where the dependent variable takes on one of two possible values (Yes/No, True/False, 0/1 etc.).

First up, you might be wondering why one can’t use linear regression for such problems. The main reason is that classification problems are about determining *class membership* rather than predicting *variable values*, and linear regression is more naturally suited to the latter than the former. One could, in principle, use linear regression for situations where there is a natural ordering of categories like *High*, *Medium* and *Low* for example. However, one then has to map sub-ranges of the predicted values to categories. Moreover, since predicted values are potentially unbounded (in data as yet unseen) there remains a degree of arbitrariness associated with such a mapping.

Logistic regression sidesteps the aforementioned issues by modelling class probabilities instead. Any input to the model yields a number lying between 0 and 1, representing the probability of class membership. One is still left with the problem of determining the threshold probability, i.e. the probability at which the category flips from one to the other. By default this is set to p=0.5, but in reality it should be settled based on how the model will be used. For example, for a marketing model that identifies potentially responsive customers, the threshold for a positive event might be set low (much less than 0.5) because the client does not really care about mailouts going to a non-responsive customer (the negative event). Indeed they may be more than OK with it as there’s always a chance – however small – that a non-responsive customer will actually respond. As an opposing example, the cost of a false positive would be high in a machine learning application that grants access to sensitive information. In this case, one might want to set the threshold probability to a value closer to 1, say 0.9 or even higher. The point is, the setting an appropriate threshold probability is a business issue, not a technical one.

### Logistic regression in brief

So how does logistic regression work?

For the discussion let’s assume that the outcome (predicted variable) and predictors are denoted by Y and X respectively and the two classes of interest are denoted by + and – respectively. We wish to model the conditional probability that the outcome Y is +, given that the input variables (predictors) are X. The conditional probability is denoted by p(Y=+|X) which we’ll abbreviate as p(X) since we know we are referring to the positive outcome Y=+.

As mentioned earlier, we are after the *probability of class membership* so we must ensure that the hypothesis function (a fancy word for the model) always lies between 0 and 1. The function assumed in logistic regression is:

You can verify that does indeed lie between 0 and 1 as varies from to . Typically, however, the values of that make sense are bounded as shown in the example (stolen from Wikipedia) shown in Figure 1. The figure also illustrates the typical S-shaped curve characteristic of logistic regression.

As an aside, you might be wondering where the name *logistic* comes from. An equivalent way of expressing the above equation is:

The quantity on the left is the logarithm of the odds. So, the model is a linear regression of the *log-odds*, sometimes called* logit*, and hence the name *logistic*.

The problem is to find the values of and that results in a that most accurately classifies all the observed data points – that is, those that belong to the positive class have a probability as close as possible to 1 and those that belong to the negative class have a probability as close as possible to 0. One way to frame this problem is to say that we wish to maximise the product of these probabilities, often referred to as the likelihood:

Where represents the products over i and j, which run over the +ve and –ve classed points respectively. This approach, called maximum likelihood estimation, is quite common in many machine learning settings, especially those involving probabilities.

It should be noted that in practice one works with the *log* likelihood because it is easier to work with mathematically. Moreover, one *minimises* the *negative * log likelihood which, of course, is the same as maximising the log likelihood. The quantity one minimises is thus:

However, these are technical details that I mention only for completeness. As you will see next, they have little bearing on the practical use of logistic regression.

### Logistic regression in R – an example

In this example, we’ll use the logistic regression option implemented within the *glm* function that comes with the base R installation. This function fits a class of models collectively known as *generalized linear models*. We’ll apply the function to the *Pima Indian Diabetes* dataset that comes with the mlbench package. The code is quite straightforward – particularly if you’ve read earlier articles in my “gentle introduction” series – so I’ll just list the code below noting that the logistic regression option is invoked by setting *family=”binomial”* in the *glm* function call.

Here we go:

Although this seems pretty good, we aren’t quite done because there is an issue that is lurking under the hood. To see this, let’s examine the information output from the model summary, in particular the coefficient estimates (i.e. estimates for ) and their significance. Here’s a summary of the information contained in the table:

- Column 2 in the table lists coefficient estimates.
- Column 3 list s the standard error of the estimates (the larger the standard error, the less confident we are about the estimate)
- Column 4 the z statistic (which is the coefficient estimate (column 2) divided by the standard error of the estimate (column 3)) and
- The last column (Pr(>|z|) lists the p-value, which is the probability of getting the listed estimate assuming the predictor has no effect. In essence, the
*smaller*the p-value, the*more significant*the estimate is likely to be.

From the table we can conclude that only 4 predictors are significant – *pregnant*, *glucose*, *mass *and *pedigree *(and possibly a fifth – *pressure*). The other variables have little predictive power and worse, may contribute to overfitting. They should, therefore, be eliminated and we’ll do that in a minute. However, there’s an important point to note before we do so…

In this case we have only 9 variables, so are able to identify the significant ones by a manual inspection of p-values. As you can well imagine, such a process will quickly become tedious as the number of predictors increases. Wouldn’t it be be nice if there were an algorithm that could somehow automatically shrink the coefficients of these variables or (better!) set them to zero altogether? It turns out that this is precisely what lasso and its close cousin, ridge regression, do.

### Ridge and Lasso

Recall that the values of the logistic regression coefficients and are found by minimising the negative log likelihood described in equation (3). Ridge and lasso regularization work by adding a *penalty term* to the log likelihood function. In the case of ridge regression, the penalty term is and in the case of lasso, it is (Remember, is a vector, with as many components as there are predictors). The quantity to be minimised in the two cases is thus:

– for ridge regression,

and

– for lasso regression.

Where is a free parameter which is usually selected in such a way that the resulting model minimises the out of sample error. Typically, the optimal value of is found using grid search with cross-validation, a process akin to the one described in my discussion on cost-complexity parameter estimation in decision trees. Most canned algorithms provide methods to do this; the one we’ll use in the next section is no exception.

In the case of ridge regression, the effect of the penalty term is to shrink the coefficients that contribute most to the error. Put another way, it reduces the magnitude of the coefficients that contribute to *increasing* . In contrast, in the case of lasso regression, the effect of the penalty term is to set the these coefficients *exactly to zero! *This is cool because what it mean that lasso regression works like a feature selector that picks out the most important coefficients, i.e. those that are most predictive (and have the lowest p-values).

Let’s illustrate this through an example. We’ll use the glmnet package which implements a combined version of ridge and lasso (called elastic net). Instead of minimising (4) or (5) above, glmnet minimises:

where controls the “mix” of ridge and lasso regularisation, with being “pure” ridge and being “pure” lasso.

### Lasso regularisation using glmnet

Let’s reanalyse the Pima Indian Diabetes dataset using glmnet with (pure lasso). Before diving into code, it is worth noting that glmnet:

- does not have a formula interface, so one has to input the predictors as a matrix and the class labels as a vector.
- does not accept categorical predictors, so one has to convert these to numeric values before passing them to glmnet.

The glmnet function model.matrix creates the matrix and also converts categorical predictors to appropriate dummy variables.

Another important point to note is that we’ll use the function cv.glmnet, which automatically performs a grid search to find the optimal value of .

OK, enough said, here we go:

The plot is shown in Figure 2 below:

The plot shows that the *log *of the optimal value of lambda (i.e. the one that minimises the root mean square error) is approximately -5. The exact value can be viewed by examining the variable *lambda_min* in the code below. In general though, the objective of regularisation is to balance accuracy *and *simplicity. In the present context, this means a model with the* smallest number of coefficients that also gives a good accuracy*. To this end, the cv.glmnet function finds the value of lambda that gives the simplest model but also lies within one standard error of the optimal value of lambda. This value of lambda (*lambda.1se*) is what we’ll use in the rest of the computation. Interested readers should have a look at this article for more on* lambda.1se* vs* lambda.min*.

The output shows that *only* those variables that we had determined to be significant on the basis of p-values have non-zero coefficients. The coefficients of all other variables have been set to zero by the algorithm! Lasso has reduced the complexity of the fitting function massively…and you are no doubt wondering what effect this has on accuracy. Let’s see by running the model against our test data:

Which is a bit less than what we got with the more complex model. So, we get a similar out-of-sample accuracy as we did before, and we do so using a way simpler function (4 non-zero coefficients) than the original one (9 nonzero coefficients). What this means is that the simpler function does at least as good a job fitting the signal in the data as the more complicated one. The bias-variance tradeoff tells us that the simpler function should be preferred because it is less likely to overfit the training data.

Paraphrasing William of Ockham: *all other things being equal, a simple hypothesis should be preferred over a complex one*.

### Wrapping up

In this post I have tried to provide a detailed introduction to logistic regression, one of the simplest (and oldest) classification techniques in the machine learning practitioners arsenal. Despite it’s simplicity (or I should say, because of it!) logistic regression works well for many business applications which often have a simple decision boundary. Moreover, because of its simplicity it is less prone to overfitting than flexible methods such as decision trees. Further, as we have shown, variables that contribute to overfitting can be eliminated using lasso (or ridge) regularisation, without compromising out-of-sample accuracy. Given these advantages and its inherent simplicity, it isn’t surprising that logistic regression remains a workhorse for data scientists.

## A gentle introduction to support vector machines using R

### Introduction

Most machine learning algorithms involve minimising an error measure of some kind (this measure is often called an objective function or loss function). For example, the error measure in linear regression problems is the famous mean squared error – i.e. the averaged sum of the squared differences between the predicted and actual values. Like the mean squared error, most objective functions depend on *all* points in the training dataset. In this post, I describe the support vector machine (SVM) approach which focuses instead on finding the *optimal separation boundary between* *datapoints that have different classifications. * I’ll elaborate on what this means in the next section.

Here’s the plan in brief. I’ll begin with the rationale behind SVMs using a simple case of a binary (two class) dataset with a simple separation boundary (I’ll clarify what “simple” means in a minute). Following that, I’ll describe how this can be generalised to datasets with more complex boundaries. Finally, I’ll work through a couple of examples in R, illustrating the principles behind SVMs. In line with the general philosophy of my “*Gentle Introduction to Data Science Using R*” series, the focus is on developing an intuitive understanding of the algorithm along with a practical demonstration of its use through a toy example.

### The rationale

The basic idea behind SVMs is best illustrated by considering a simple case: a set of data points that belong to one of two classes, red and blue, as illustrated in figure 1 below. To make things simpler still, I have assumed that the boundary separating the two classes is a straight line, represented by the solid green line in the diagram. In the technical literature, such datasets are called *linearly separable*.

In the linearly separable case, there is usually a fair amount of freedom in the way a separating line can be drawn. Figure 2 illustrates this point: the two broken green lines are also valid separation boundaries. Indeed, because there is a non-zero distance between the two closest points between categories, there are an infinite number of possible separation lines. This, quite naturally, raises the question as to whether it is possible to choose a separation boundary that is optimal.

The short answer is, yes there is. One way to do this is to select a boundary line that *maximises* the *margin*, i.e. the *distance between the separation boundary and the points that are closest to it*. Such an optimal boundary is illustrated by the black brace in Figure 3. The really cool thing about this criterion is that the *location of the separation boundary depends only on the points that are closest to it*. This means, unlike other classification methods, the classifier does not depend on any other points in dataset. The directed lines between the boundary and the closest points on either side are called *support vectors* (these are the solid black lines in figure 3). A direct implication of this is that the fewer the support vectors, the better the generalizability of the boundary.

Although the above sounds great, it is of limited practical value because *real data sets are seldom (if ever) linearly separable*.

So, what can we do when dealing with real (i.e. non linearly separable) data sets?

A simple approach to tackle small deviations from linear separability is to allow a small number of points (those that are close to the boundary) to be misclassified. The number of possible misclassifications is governed by a free parameter C, which is called the cost. The cost is essentially the penalty associated with making an error: the higher the value of C, the less likely it is that the algorithm will misclassify a point.

This approach – which is called *soft margin classification* – is illustrated in Figure 4. Note the points on the wrong side of the separation boundary. We will demonstrate soft margin SVMs in the next section. (*Note*: At the risk of belabouring the obvious, the purely linearly separable case discussed in the previous para is simply is a special case of the soft margin classifier.)

Real life situations are much more complex and cannot be dealt with using soft margin classifiers. For example, as shown in Figure 5, one could have widely separated clusters of points that belong to the same classes. Such situations, which require the use of multiple (and nonlinear) boundaries, can sometimes be dealt with using a clever approach called the kernel trick.

### The kernel trick

Recall that in the linearly separable (or soft margin) case, the SVM algorithm works by finding a separation boundary that maximises the margin, which is the *distance* between the boundary and the points closest to it. The distance here is the usual straight line distance between the boundary and the closest point(s). This is called the Euclidean distance in honour of the great geometer of antiquity. The point to note is that this process results in a separation boundary that is a straight line, which as Figure 5 illustrates, does not always work. In fact in most cases it won’t.

So what can we do? To answer this question, we have to take a bit of a detour…

What if we were able to generalize the notion of distance in a way that generates nonlinear separation boundaries? It turns out that this is possible. To see how, one has to first understand how the notion of distance can be generalized.

The key properties that any measure of distance must satisfy are:

**Non-negativity**– a distance cannot be negative, a point that needs no further explanation I reckon 🙂**Symmetry**– that is, the distance between point A and point B is the same as the distance between point B and point A.**Identity**– the distance between a point and itself is zero.**Triangle inequality**– that is the sum of distances between point A and B and points B and C must be less than or equal to the distance between A and C (equality holds only if all three points lie along the same line).

Any mathematical object that displays the above properties is akin to a distance. Such *generalized distances* are called metrics and the mathematical space in which they live is called a metric space. Metrics are defined using special mathematical functions designed to satisfy the above conditions. These functions are known as *kernels*.

The essence of the kernel trick lies in mapping the classification problem to a metric space in which the problem is rendered separable via a separation boundary that is simple in the new space, but complex – as it has to be – in the original one. Generally, the transformed space has a higher dimensionality, with each of the dimensions being (possibly complex) combinations of the original problem variables. However, this is not necessarily a problem because in practice one doesn’t actually mess around with transformations, one just tries different kernels (the transformation being implicit in the kernel) and sees which one does the job. The check is simple: we simply test the predictions resulting from using different kernels against a held out subset of the data (as one would for any machine learning algorithm).

It turns out that a particular function – called the radial basis function kernel (RBF kernel) – is very effective in many cases. The RBF kernel is essentially a Gaussian (or Normal) function with the Euclidean distance between pairs of points as the variable (see equation 1 below). The basic rationale behind the RBF kernel is that it creates separation boundaries that it tends to classify points close together (in the Euclidean sense) in the original space in the same way. This is reflected in the fact that the kernel decays (i.e. drops off to zero) as the Euclidean distance between points increases.

The rate at which a kernel decays is governed by the parameter – the higher the value of , the more rapid the decay. This serves to illustrate that the RBF kernel is extremely flexible….but the flexibility comes at a price – the danger of overfitting for large values of . One should choose appropriate values of C and so as to ensure that the resulting kernel represents the best possible balance between flexibility and accuracy. We’ll discuss how this is done in practice later in this article.

Finally, though it is probably obvious, it is worth mentioning that the separation boundaries for arbitrary kernels are also defined through support vectors as in Figure 3. To reiterate a point made earlier, this means that a solution that has fewer support vectors is likely to be more robust than one with many. Why? Because the data points defining support vectors are ones that are most sensitive to noise- therefore the fewer, the better.

There are many other types of kernels, each with their own pros and cons. However, I’ll leave these for adventurous readers to explore by themselves. Finally, for a much more detailed….and dare I say, better… explanation of the kernel trick, I highly recommend this article by Eric Kim.

### Support vector machines in R

In this demo we’ll use the *svm* interface that is implemented in the *e1071* R package. This interface provides R programmers access to the comprehensive *libsvm* library written by Chang and Lin. I’ll use two toy datasets: the famous iris dataset available with the base R package and the sonar dataset from the *mlbench* package. I won’t describe details of the datasets as they are discussed at length in the documentation that I have linked to. However, it is worth mentioning the reasons why I chose these datasets:

- As mentioned earlier, no real life dataset is linearly separable, but the iris dataset is almost so. Consequently, it is a good illustration of using linear SVMs. Although one almost never uses these in practice, I have illustrated their use primarily for pedagogical reasons.
- The sonar dataset is a good illustration of the benefits of using RBF kernels in cases where the dataset is hard to visualise (60 variables in this case!). In general, one would almost always use RBF (or other nonlinear) kernels in practice.

With that said, let’s get right to it. I assume you have R and RStudio installed. For instructions on how to do this, have a look at the first article in this series. The processing preliminaries – loading libraries, data and creating training and test datasets are much the same as in my previous articles so I won’t dwell on these here. For completeness, however, I’ll list all the code so you can run it directly in R or R studio (a complete listing of the code can be found **here**):

The output from the SVM model show that there are 24 support vectors. If desired, these can be examined using the *SV *variable in the model – i.e via* svm_model$SV*.

The test prediction accuracy indicates that the linear performs quite well on this dataset, confirming that it is indeed near linearly separable. To check performance by class, one can create a confusion matrix as described in my post on random forests. I’ll leave this as an exercise for you. Another point is that we have used a soft-margin classification scheme with a cost C=1. You can experiment with this by explicitly changing the value of C. Again, I’ll leave this for you an exercise.

Before proceeding to the RBF kernel, I should mention a point that an alert reader may have noticed. The predicted variable, *Species*, can take on 3 values (setosa, versicolor and virginica). However, our discussion above dealt with a binary (2 valued) classification problem. This brings up the question as to how the algorithm deals multiclass classification problems – i.e those involving datasets with more than two classes. The *libsvm* algorithm (which *svm* uses) does this using a one-against-one classification strategy. Here’s how it works:

- Divide the dataset (assumed to have N classes) into N(N-1)/2 datasets that have two classes each.
- Solve the binary classification problem for each of these subsets
- Use a simple voting mechanism to assign a class to each data point.

Basically, each data point is assigned the most frequent classification it receives from all the binary classification problems it figures in.

With that said for the unrealistic linear classifier, let’s move to the real world. In the code below, I build SVM models using three different kernels

- Linear kernel (this is for comparison with the following 2 kernels).
- RBF kernel with default values for the parameters and .
- RBF kernel with optimal values for and . The optimal values are obtained using the
*tune.svm*function (also available in*e1071*), which essentially builds models for multiple combinations of parameter values and selects the best.

OK, lets go:

I’ll leave you to examine the contents of the model. The important point to note here is that the performance of the model with the test set is quite dismal compared to the previous case. This simply indicates that the linear kernel is not appropriate here. Let’s take a look at what happens if we use the RBF kernel with default values for the parameters:

That’s a pretty decent improvement from the linear kernel. Let’s see if we can do better by doing some parameter tuning. To do this we first invoke *tune.svm* and use the parameters it gives us in the call to *svm*:

**tune_out <- tune.svm(x=trainset[,-typeColNum],y=trainset[,typeColNum],gamma=10^(-3:3),cost=c(0.01,0.1,1,10,100,1000),kernel="radial")**

**tune_out$best.parameters$cost**

**tune_out$best.parameters$gamma**

Which is fairly decent improvement on the un-optimised case.

### Wrapping up

This bring us to the end of this introductory exploration of SVMs in R. To recap, the distinguishing feature of SVMs in contrast to most other techniques is that they attempt to construct optimal separation boundaries between different categories.

SVMs are quite versatile and have been applied to a wide variety of domains ranging from chemistry to pattern recognition. They are best used in binary classification scenarios. This brings up a question as to where SVMs are to be preferred to other binary classification techniques such as logistic regression. The honest response is, “it depends” – but here are some points to keep in mind when choosing between the two. A general point to keep in mind is that SVM algorithms tend to be expensive both in terms of memory and computation, issues that can start to hurt as the size of the dataset increases.

Given all the above caveats and considerations, the best way to figure out whether an SVM approach will work for your problem may be to do what most machine learning practitioners do: try it out!

## A gentle introduction to random forests using R

### Introduction

In a previous post, I described how decision tree algorithms work and demonstrated their use via the rpart library in R. Decision trees work by splitting a dataset *recursively*. That is, subsets arising from a split are further split until a predetermined termination criterion is reached. At each step, a split is made based on the independent variable that results in the *largest possible reduction in heterogeneity* *of the dependent variable*.

(*Note*: readers unfamiliar with decision trees may want to read that post before proceeding)

The main drawback of decision trees is that they are prone to overfitting. The reason for this is that trees, if grown deep, are able to fit all kinds of variations in the data, including noise. Although it is possible to address this partially by pruning, the result often remains less than satisfactory. This is because the algorithm makes a locally optimal choice at each split without any regard to whether the choice made is the best one overall. A poor split made in the initial stages can thus doom the model, a problem that cannot be fixed by post-hoc pruning.

In this post I describe random forests, a tree-based algorithm that addresses the above shortcoming of decision trees. I’ll first describe the intuition behind the algorithm via an analogy and then do a demo using the R randomForest library.

### Motivating random forests

One of the reasons for the popularity of decision trees is that they reflect the way humans make decisions: by weighing up options at each stage and choosing the best one available. The analogy is particularly useful because it also suggests how decision trees can be improved.

One of the lifelines in the game show, Who Wants to be A Millionaire, is “Ask The Audience” wherein a contestant can ask the audience to vote on the answer to a question. The rationale here is that the majority response from a large number of independent decision makers is more likely to yield a correct answer than one from a randomly chosen person. There are two factors at play here:

- People have different experiences and will therefore draw upon different “data” to answer the question.
- People have different knowledge bases and preferences and will therefore draw upon different “variables” to make their choices at each stage in their decision process.

Taking a cue from the above, it seems reasonable to build many decision trees using:

- Different sets of training data.
- Randomly selected subsets of variables at each split of every decision tree.

Predictions can then made by taking the majority vote over all trees (for classification problems) or averaging results over all trees (for regression problems). This is essentially how the random forest algorithm works.

The net effect of the two strategies is to reduce overfitting by a) averaging over trees created from different samples of the dataset and b) decreasing the likelihood of a small set of strong predictors dominating the splits. The price paid is reduced interpretability as well as increased computational complexity. But then, there is no such thing as a free lunch.

### The mechanics of the algorithm

Although we will not delve into the mathematical details of the algorithm, it is important to understand how two points made above are implemented in the algorithm.

#### Bootstrap aggregating… and a (rather cool) error estimate

A key feature of the algorithm is the use of multiple datasets for training individual decision trees. This is done via a neat statistical trick called **b**ootstrap **agg**regat**ing** (also called *bagging*).

Here’s how bagging works:

Assume you have a dataset of size *N*. From this you create a sample (i.e. a subset) of size *n* (*n less than or equal to N*) by choosing n data points randomly with replacement. “Randomly” means every point in the dataset is equally likely to be chosen and “with replacement” means that a specific data point can appear more than once in the subset. Do this *M* times to create *M* equally-sized samples of size *n* each. It can be shown that this procedure, which statisticians call bootstrapping, is legit when samples are created from large datasets – that is, when *N* is large.

Because a bagged sample is created by selection with replacement, there will generally be some points that are not selected. In fact, it can be shown that, on the average, each sample will use about two-thirds of the available data points. This gives us a clever way to estimate the error as part of the process of model building.

Here’s how:

For every data point, obtain predictions for trees in which the point was *out of bag*. From the result mentioned above, this will yield approximately *M/3* predictions per data point (because a third of the data points are out of bag). Take the majority vote of these *M/3 *predictions as the predicted value for the data point. One can do this for the entire dataset. From these out of bag predictions for the whole dataset, we can estimate the overall error by computing a classification error (Count of correct predictions divided by N) for classification problems or the root mean squared error for regression problems. This means there is no need to have a separate test data set, which is kind of cool. However, if you have enough data, it is worth holding out some data for use as an independent test set. This is what we’ll do in the demo later.

#### Using subsets of predictor variables

Although bagging reduces overfitting somewhat, it does not address the issue completely. The reason is that in most datasets a small number of predictors tend to dominate the others. These predictors tend to be selected in early splits and thus influence the shapes and sizes of a significant fraction of trees in the forest. That is, strong predictors enhance correlations between trees which tends to come in the way of variance reduction.

A simple way to get around this problem is to use a *random subset of variables at each split*. This avoids over-representation of dominant variables and thus creates a more diverse forest. This is precisely what the random forest algorithm does.

### Random forests in R

In what follows, I use the famous Glass dataset from the mlbench library. The dataset has 214 data points of six types of glass with varying metal oxide content and refractive indexes. I’ll first build a decision tree model based on the data using the rpart library (recursive partitioning) that I covered in an earlier article and then use then show how one can build a random forest model using the randomForest library. The rationale behind this is to compare the two models – single decision tree vs random forest. In the interests of space, I won’t explain details of the rpart here as I’ve covered it at length in the previous article. However, for completeness, I’ll list the demo code for it before getting into random forests.

#### Decision trees using rpart

Here’s the code listing for building a decision tree using rpart on the Glass dataset (please see my previous article for a full explanation of each step). Note that I have not used pruning as there is little benefit to be gained from it (*Exercise for the reader*: try this for yourself!).

Now, we know that decision tree algorithms tend to display high variance so the hit rate from any one tree is likely to be misleading. To address this we’ll generate a bunch of trees using different training sets (via random sampling) and calculate an average hit rate and spread (or standard deviation).

The decision tree algorithm gets it right about 69% of the time with a variation of about 5%. The variation isn’t too bad here, but the accuracy has hardly improved at all (*Exercise for the reader*: why?). Let’s see if we can do better using random forests.

#### Random forests

As discussed earlier, a random forest algorithm works by averaging over multiple trees using bootstrapped samples. Also, it reduces the correlation between trees by splitting on a random subset of predictors at each node in tree construction. The key parameters for randomForest algorithm are the number of trees (*ntree*) and the number of variables to be considered for splitting (*mtry*). The algorithm sets a *default of 500 for ntree* and sets *mtry* to the *square root of the the number of predictors for classification problems *or *one-third the total number of predictors for regression.* * *These defaults can be overridden by explicitly providing values for these variables.

The preliminary stuff – the creation of training and test datasets etc. – is much the same as for decision trees but I’ll list the code for completeness.

randomForest(formula = Type ~ ., data = trainGlass, importance = TRUE, xtest = testGlass[, -typeColNum], ntree = 1001)

1 | 2 | 3 | 5 | 6 | 7 | class.error | |

1 | 40 | 7 | 2 | 0 | 0 | 0 | 0.1836735 |

2 | 8 | 49 | 1 | 2 | 2 | 1 | 0.2222222 |

3 | 6 | 3 | 6 | 0 | 0 | 0 | 0.6000000 |

5 | 0 | 1 | 0 | 11 | 0 | 1 | 0.1538462 |

6 | 1 | 2 | 0 | 1 | 6 | 0 | 0.5000000 |

7 | 1 | 2 | 0 | 1 | 0 | 21 | 0.1600000 |

The first thing to note is the out of bag error estimate is ~ 24%. Equivalently the hit rate is 76%, which is better than the 69% for decision trees. Secondly, you’ll note that the algorithm does a terrible job identifying type 3 and 6 glasses correctly. This could possibly be improved by a technique called boosting, which works by iteratively improving poor predictions made in earlier stages. I plan to look at boosting in a future post, but if you’re curious, check out the gbm package in R.

Finally, for completeness, let’s see how the test set does:

1 | 2 | 3 | 5 | 6 | 7 | |

1 | 19 | 2 | 0 | 0 | 0 | 0 |

2 | 1 | 9 | 1 | 0 | 0 | 0 |

3 | 1 | 1 | 1 | 0 | 0 | 0 |

5 | 0 | 1 | 0 | 0 | 0 | 0 |

6 | 0 | 0 | 0 | 0 | 3 | 0 |

7 | 0 | 0 | 0 | 0 | 0 | 4 |

The test accuracy is better than the out of bag accuracy and there are some differences in the class errors as well. However, overall the two compare quite well and are significantly better than the results of the decision tree algorithm.

### Variable importance

Random forest algorithms also give measures of *variable importance*. Computation of these is enabled by setting *importance, *a boolean parameter, to TRUE. The algorithm computes two measures of variable importance: *mean decrease in Gini* and *mean decrease in accuracy*. Brief explanations of these follow.

#### Mean decrease in Gini

When determining splits in individual trees, the algorithm looks for the largest class (in terms of population) and attempts to isolate it first. If this is not possible, it tries to do the best it can, always focusing on isolating the largest remaining class in every split.This is called the Gini splitting rule (see this article for a good explanation of the rule).

The “goodness of split” is measured by the Gini Impurity, . For a set containing K categories this is given by:

where is the fraction of the set that belongs to the ith category. Clearly, is 0 when the set is homogeneous or pure (1 class only) and is maximum when classes are equiprobable (for example, in a two class set the maximum occurs when and are 0.5). At each stage the algorithm chooses to split on the predictor that leads to the largest decrease in . The algorithm tracks this decrease for each predictor for all splits and all trees in the forest. The average is reported as the *mean decrease in Gini.*

#### Mean decrease in accuracy

The mean decrease in accuracy is calculated using the out of bag data points for each tree. The procedure goes as follows: when a particular tree is grown, the out of bag points are passed down the tree and the prediction accuracy (based on all out of bag points) recorded . The predictors are then randomly permuted and the out of bag prediction accuracy recalculated. The decrease in accuracy for a given predictor is the difference between the accuracy of the original (unpermuted) tree and the those obtained from the permuted trees in which the predictor was excluded. As in the previous case, the decrease in accuracy for each predictor can be computed and tracked as the algorithm progresses. These can then be averaged by predictor to yield a mean decrease in accuracy.

#### Variable importance plot

From the above, it would seem that the mean decrease in accuracy is a more global measure as it uses fully constructed trees in contrast to the Gini measure which is based on individual splits. In practice, however, there could be other reasons for choosing one over the other…but that is neither here nor there, if you set importance to TRUE, you’ll get both. The numerical measures of importance are returned in the randomForest object (Glass.rf in our case), but I won’t list them here. Instead, I’ll just print out the variable importance plots for the two measures as these give a good visual overview of the relative importance of variables. The code is a simple one-liner:

The plot is shown in Figure 1 below.

In this case the two measures are pretty consistent so it doesn’t really matter which one you choose.

### Wrapping up

Random forests are an example of a general class of techniques called ensemble methods. These techniques are based on the principle that averaging over a large number of not-so-good models yields a more reliable prediction than a single model. This is true only if models in the group are independent of each other, which is precisely what bootstrap aggregation and predictor subsetting are intended to achieve.

Although considerably more complex than decision trees, the logic behind random forests is not hard to understand. Indeed, the intuitiveness of the algorithm together with its ease of use and accuracy have made it very popular in the machine learning community.

## A gentle introduction to decision trees using R

### Introduction

Most techniques of predictive analytics have their origins in probability or statistical theory (see my post on Naïve Bayes, for example). In this post I’ll look at one that has more a commonplace origin: the way in which humans make decisions. When making decisions, we typically identify the options available and then evaluate them based on criteria that are important to us. The intuitive appeal of such a procedure is in no small measure due to the fact that it can be easily explained through a visual. Consider the following graphic, for example:

(Original image: https://www.flickr.com/photos/dullhunk/7214525854, Credit: Duncan Hull)

The tree structure depicted here provides a neat, easy-to-follow description of the issue under consideration and its resolution. The decision procedure is based on asking a series of questions, each of which serve to further reduce the domain of possibilities. The predictive technique I discuss in this post,classification and regression trees (CART), works in much the same fashion. It was invented by Leo Breiman and his colleagues in the 1970s.

In what follows, I will use the open source software, R. If you are new to R, you may want to follow this link for more on the basics of setting up and installing it. Note that the R implementation of the CART algorithm is called RPART (Recursive Partitioning And Regression Trees). This is essentially because Breiman and Co. trademarked the term CART. As some others have pointed out, it is somewhat ironical that the algorithm is now commonly referred to as RPART rather than by the term coined by its inventors.

### A bit about the algorithm

The *rpart* algorithm works by splitting the dataset recursively, which means that the subsets that arise from a split are further split until a predetermined termination criterion is reached. At each step, the split is made based on the independent variable that results in the *largest possible reduction in heterogeneity* *of the dependent (predicted) variable*.

Splitting rules can be constructed in many different ways, all of which are based on the notion of *impurity- *a measure of the degree of heterogeneity of the leaf nodes. Put another way, a leaf node that contains a single class is homogeneous and has impurity=0. There are three popular impurity quantification methods: Entropy (aka information gain), Gini Index and Classification Error. Check out this article for a simple explanation of the three methods.

The *rpart* algorithm offers the entropy and Gini index methods as choices. There is a fair amount of fact and opinion on the Web about which method is better. Here are some of the better articles I’ve come across:

https://www.garysieling.com/blog/sklearn-gini-vs-entropy-criteria

http://www.salford-systems.com/resources/whitepapers/114-do-splitting-rules-really-matter

The answer as to which method is the best is: *it depends*. Given this, it may be prudent to try out a couple of methods and pick the one that works best for your problem.

Regardless of the method chosen, the splitting rules partition the decision space (a fancy word for the entire dataset) into *rectangular regions* each of which correspond to a split. Consider the following simple example with two predictors x1 and x2. The first split is at x1=1 (which splits the decision space into two regions x11), the second at x2=2, which splits the (x1>1) region into 2 sub-regions, and finally x1=1.5 which splits the (x1>1,x2>2) sub-region further.

It is important to note that the algorithm works by making the* best possible choice at each particular stage, without any consideration of whether those choices remain optimal in future stages*. That is, the algorithm makes *a locally optimal decision *at each stage. It is thus quite possible that such a choice at one stage turns out to be sub-optimal in the overall scheme of things. In other words, the algorithm does not find a globally optimal tree.

Another important point relates to well-known bias-variance tradeoff in machine learning, which in simple terms is a tradeoff between the degree to which a model fits the training data and its predictive accuracy. This refers to the general rule that beyond a point, it is counterproductive to improve the fit of a model to the training data as this increases the likelihood of overfitting. It is easy to see that deep trees are more likely to overfit the data than shallow ones. One obvious way to control such overfitting is to construct shallower trees by stopping the algorithm at an appropriate point based on whether a split significantly improves the fit. Another is to grow a tree unrestricted and then prune it back using an appropriate criterion. The *rpart *algorithm takes the latter approach.

Here is how it works in brief:

Essentially one minimises the *cost*, , a quantity that is a linear combination of the *error* (essentially, the fraction of misclassified instances, or variance in the case of a continuous variable), and the *number of leaf nodes in the tree, *:

First, we note that when , this simply returns the original fully grown tree. As increases, we incur a penalty that is proportional to the number of leaf nodes. This tends to cause the minimum cost to occur for a tree that is a subtree of the original one (since a subtree will have a smaller number of leaf nodes). In practice we vary and pick the value that gives the subtree that results in the smallest cross-validated prediction error. One does not have to worry about programming this because the *rpart* algorithm actually computes the errors for different values of for us. All we need to do is pick the value of the coefficient that gives the lowest cross-validated error. I will illustrate this in detail in the next section.

An implication of their tendency to overfit data is that decision trees tend to be sensitive to relatively minor changes in the training datasets. Indeed, small differences can lead to radically different looking trees. Pruning addresses this to an extent, but does not resolve it completely. A better resolution is offered by the so-called ensemble methods that *average over many differently constructed trees*. I’ll discuss one such method at length in a future post.

Finally, I should also mention that decision trees can be used for both classification and regression problems (i.e. those in which the predicted variable is discrete and continuous respectively). I’ll demonstrate both types of problems in the next two sections.

### Classification trees using rpart

To demonstrate classification trees, we’ll use the *Ionosphere* dataset available in the *mlbench* package in R. I have chosen this dataset because it nicely illustrates the points I wish to make in this post. In general, you will almost always find that algorithms that work fine on classroom datasets do not work so well in the real world…but of course, you know that already!

We begin by setting the working directory, loading the required packages (*rpart* and *mlbench*) and then loading the Ionosphere dataset.

Next we separate the data into training and test sets. We’ll use the former to build the model and the latter to test it. To do this, I use a simple scheme wherein I randomly select 80% of the data for the training set and assign the remainder to the test data set. This is easily done in a single R statement that invokes the uniform distribution (*runif*) and the vectorised function, *ifelse*. Before *invoking* runif, I set a seed integer to my favourite integer in order to ensure reproducibility of results.

In the above, I have also removed the training flag from the training and test datasets.

Next we invoke *rpart*. I strongly recommend you take some time to go through the documentation and understand the parameters and their defaults values. Note that we need to remove the predicted variable from the dataset before passing the latter on to the algorithm, which is why we need to find the column index of the predicted variable (first line below). Also note that we set the *method* parameter to “class*“,* which simply tells the algorithm that the predicted variable is discrete.* *Finally, *rpart* uses Gini rule for splitting by default, and we’ll stick with this option.

The resulting plot is shown in Figure 3 below. It is quite self-explanatory so I won’t dwell on it here.

Next we see how good the model is by seeing how it fares against the test data.

pred true | bad | good |

bad | 17 | 2 |

good | 9 | 43 |

Note that we need to verify the above results by doing multiple runs, each using different training and test sets. I will do this later, after discussing pruning.

Next, we prune the tree using the cost complexity criterion. Basically, the intent is to see if a shallower subtree can give us comparable results. If so, we’d be better of choosing the shallower tree because it reduces the likelihood of overfitting.

As described earlier, we choose the appropriate pruning parameter (aka cost-complexity parameter) by picking the value that results in the lowest prediction error. Note that all relevant computations have already been carried out by R when we built the original tree (the call to *rpart* in the code above). All that remains now is to pick the value of :

CP | nsplit | rel error | xerror | xstd | |

1 | 0.57 | 0 | 1.00 | 1.00 | 0.080178 |

2 | 0.20 | 1 | 0.43 | 0.46 | 0.062002 |

3 | 0.02 | 2 | 0.23 | 0.26 | 0.048565 |

4 | 0.01 | 4 | 0.19 | 0.35 |

It is clear from the above, that the lowest cross-validation error (*xerror* in the table) occurs for (this is CP in the table above). One can find CP programatically like so:

Next, we prune the tree based on this value of CP:

Note that *rpart* will use a default CP value of 0.01 if you don’t specify one in *prune*.

The pruned tree is shown in Figure 4 below.

Let’s see how this tree stacks up against the fully grown one shown in Fig 3.

This seems like an improvement over the unpruned tree, but one swallow does not a summer make. We need to check that this holds up for different training and test sets. This is easily done by creating multiple random partitions of the dataset and checking the efficacy of pruning for each. To do this efficiently, I’ll create a function that takes the training fraction, number of runs (partitions) and the name of the dataset as inputs and outputs the proportion of correct predictions for each run. It also optionally prunes the tree. Here’s the code:

Note that in the above, I have set the default value of the prune_tree to FALSE, so the function will execute the first branch of the if statement unless the default is overridden.

OK, so let’s do 50 runs with and without pruning, and check the mean and variance of the results for both sets of runs.

So we see that there is an improvement of about 3% with pruning. Also, if you were to plot the trees as we did earlier, you would see that this improvement is achieved with shallower trees. Again, I point out that this is not always the case. In fact, it often happens that pruning results in worse predictions, albeit with better reliability – a classic illustration of the bias-variance tradeoff.

### Regression trees using rpart

In the previous section we saw how one can build decision trees for situations in which the predicted variable is discrete. Let’s now look at the case in which the predicted variable is continuous. We’ll use the Boston Housing dataset from the *mlbench* package. Much of the discussion of the earlier section applies here, so I’ll just display the code, explaining only the differences.

Next we invoke *rpart*, noting that the predicted variable is *medv* (median value of owner-occupied homes in $1000 units) and that we need to set the *method * parameter to “*anova*“. The latter tells *rpart* that the predicted variable is continuous (i.e that this is a regression problem).

The plot of the tree is shown in Figure 5 below.

Next, we need to see how good the predictions are. Since the dependent variable is continuous, we cannot compare the predictions directly against the test set. Instead, we calculate the root mean square (RMS) error. To do this, we request *rpart* to output the predictions as a vector – one prediction per record in the test dataset. The RMS error can then easily be calculated by comparing this vector with the *medv* column in the test dataset.

Here is the relevant code:

Again, we need to do multiple runs to check on the reliability of the predictions. However, you already know how to do that so I will leave it to you.

Moving on, we prune the tree using the cost complexity criterion as before. The code is exactly the same as in the classification problem.

The tree is unchanged so I won’t show it here. This means, as far as the cost complexity pruning is concerned, the optimal subtree is the same as the original tree. To confirm this, we’d need to do multiple runs as before – something that I’ve already left as as an exercise for you :). Basically, you’ll need to write a function analogous to the one above, that computes the root mean square error instead of the proportion of correct predictions.

### Wrapping up

This brings us to the end of my introduction to classification and regression trees using R. Unlike some articles on the topic I have attempted to describe each of the steps in detail and provide at least some kind of a rationale for them. I hope you’ve found the description and code snippets useful.

I’ll end by reiterating a couple points I made early in this piece. The nice thing about decision trees is that they are easy to explain to the users of our predictions. This is primarily because they reflect the way we think about how decisions are made in real life – via a set of binary choices based on appropriate criteria. That said, in many practical situations decision trees turn out to be unstable: small changes in the dataset can lead to wildly different trees. It turns out that this limitation can be addressed by building a variety of trees using different starting points and then averaging over them. This is the domain of the so-called random forest algorithm.We’ll make the journey from decision trees to random forests in a future post.

**Postscript, 20th September 2016**: I finally got around to finishing my article on random forests.