Comparing SVM kernels

This vignette is an adaptation to R of the vignette created by Fernando Wittmann in Python avaiable at:

https://gist.github.com/WittmannF/60680723ed8dd0cb993051a7448f7805

WARNING: PHOTOSENSITIVITY/EPILEPSY SEIZURES

A very small percentage of individuals may experience epileptic seizures or blackouts when exposed to certain light patterns or flashing lights. These conditions may trigger previously undetected epileptic symptoms or seizures in persons who have no history of prior seizures or epilepsy. If you, or anyone in your family has an epileptic condition or has had seizures of any kind IMMEDIATELY DISCONTINUE the reading of this vignette.

RESUME THE READING OF THIS VIGNETTE ONLY WITH APPROVAL OF YOUR PHYSICIAN

Packages

Let’s start loading the packages we need:

library(tidyverse) # data wrangling and plotting
library(mlbench) # spirals dataset
library(e1071) # for svm algorithm
library(fdm2id) #two moons and targets datasets
library(kableExtra) #html tables
library(gridExtra) #array of graphical objects (grobs)

Datasets

Let’s define some datasets with different boundaries between classes:

#iris (linear frontier)
iris <- iris %>% rename("X" = "Sepal.Length", "Y" = "Sepal.Width") %>% mutate(Class = as.factor(ifelse(Species == "setosa", "Class 1", "Class 2")))

#two moons
two_moons <- data.twomoons(sigma = 0.1, graph = FALSE)

#target
target <- data.target2(initn = 600, graph = FALSE)

#spirals
spirals <- mlbench.spirals(n=600, cycles=2, sd=0.05)
spirals <- data.frame(X=spirals$x[,1], Y=spirals$x[,2], Class=spirals$classes)
levels(spirals$Class) <- c("Class 1", "Class 2")

SVM accuracy

We’ll classify each dataset with the four kernels available in the svm function of the e1071 package, and obtain the accuracy of each classifyer.

accuracy <- sapply(list(iris, two_moons, target, spirals), function(df){
  
  models <- lapply(list("linear", "poly", "radial", "sigmoid"), function(meth){
    return(svm(Class ~ X + Y, df, kernel=meth, prob=TRUE))
  })
  
  names(models) <- c("linear", "poly", "radial", "sigmoid")
  
  accuracy <- sapply(models, function(x) sum(x$fitted == df$Class)/nrow(df))
  names(accuracy) <- names(models)
  return(accuracy)
  
})

accuracy <- data.frame(accuracy)
colnames(accuracy) <- c("iris", "moons", "target", "spirals")

accuracy %>%
  kable(digits = 3) %>% kable_styling(bootstrap_options = "striped", full_width = F)

	iris	moons	target	spirals
linear	0.993	0.925	0.825	0.500
poly	0.987	0.925	0.825	0.507
radial	0.993	1.000	1.000	0.605
sigmoid	0.980	0.792	0.747	0.403

We see that when boundary between classes is not linear, the best results are obtained with the radial kernel.

SVM plots

Let’s see how each kernel behaves plotting the prediction of each kernel, obtaining the probability of classification for each point in the graph. Bold colors in background correspond to high probabilities of classification (above 0.75). We define the kernel_plot function to obtain the plots for a dataset:

kernel_plot <- function(df){
    grid <- expand.grid(X=seq(min(df$X), max(df$X), 0.05), Y=seq(min(df$Y), max(df$Y), 0.05))
    
    svm <- lapply(list("linear", "poly", "radial", "sigmoid"), function(x){
      return(svm(Class ~ X + Y, df, kernel=x, prob=TRUE))
    })
    names(svm) <- c("linear", "poly", "radial", "sigmoid")
    
    #grid values
    pred <- sapply(svm, function(x){
      pr <- attr(predict(x, grid, probability=TRUE), "probabilities")
      return(pr[ , which(attr(pr, "dimnames")[[2]] == "Class 1")])
    }) 
    pred_cut <- apply(pred, 2, function(x) cut(x, breaks = 4))
    grid <- cbind(grid, pred_cut)
    
    plot_base <- ggplot(df) + 
      geom_point(aes(X, Y, col=Class), show.legend = FALSE) +
      scale_color_manual(values=c("#0000FF", "#FF0000")) +
      theme_minimal()
    
    plots <- lapply(names(svm), function(x){
      plot_base +
        geom_tile(data=grid, aes(X, Y, fill=grid[, x]), show.legend = FALSE, alpha=0.5) + 
        scale_fill_manual(values=c("#FF0000", "#FF9999", "#99CCFF", "#0000FF")) + 
        annotate(geom="label", x=min(grid$X), y=max(grid$Y), label=x, hjust="inward", vjust="inward")
    })
       
    outcome <- list(plot_base)
    outcome <- c(outcome, plots)
    
    names(outcome) <- c("base", "linear", "poly", "radial", "sigmoid")
    return(outcome)
  
}

… and use grid.arrange to show the plots:

grid.arrange(grobs = c(kernel_plot(iris), kernel_plot(target), kernel_plot(two_moons), kernel_plot(spirals)), ncol=5)

No animals were harmed, and no for loops have been used in the making of this page.