SVM

Object

Les Support Vector Machines (SVM) sont une famille d’algorithmes dédiés aux problèmes de régression et de classification supervisée. Nous présenterons la démarche dans le cas de la classification supervisée avec une variable à expliquer \(Y\) binaire. On désigne par \(X = (X_1, . . . , X_p)\) le vecteur des variables explicatives.

Le succès des SVM dans les années 90 est lié d’une part à la facilité d’implémentation et d’autre part à la possibilité de transporter les données dans des espaces de grandes dimensions où on espère qu’elles seront plus séparables sans perte d’efficacité. Ce second point repose sur le fait que l’algorithme permettant d’obtenir l’hyperplan de marge maximale n’a besoin que des produits scalaires entre données pour fonctionner. Plus précisément, il est possible que les données ne soient pas linéairement séparables mais qu’elles le deviennent après leur avoir appliqué une transformation pertinente. Une telle transformation revient à plonger les données dans un nouvel espace, souvent de plus grande dimension que l’espace original. Il s’avère qu’appliquer la SVM linéaire sur les données transformées revient à appliquer une SVM sur les données initiales en utilisant un noyau à la place du produit scalaire usuel. Cette astuce, dite astuce du noyau, incite à définir directement le noyau (le produit scalaire après transformation) et non la transformation explicite.Plusieurs noyaux existent et ce choix est également primordial. Nous utiliserons les noyaux suivants dans cette fiche : – linéaire : \[k(x, x') = (x;x')\] – polynomial : \[k(x, x')=(scale x, x') + offset)^{degree}\] – gaussien ou radial : \[k(x, x')= exp(−σ||x-x||)^2\]

Ici \(x\) et \(x\)’sont dans \(R^p\),\((x,x')\) est le produit scalaire usuel entre \(x\) et \(x'\) tandis que scale, offset, degree et σ sont des paramètres à calibrer.

1. Importer les données

Le jeu de données se trouve dans le package kernlab, on peut le charger par :

library(kernlab)
data(spam)

On sépare l’échantillon en un échantillon d’apprentissage de taille 3000 et un échantillon de validation de taille 1601. L’échantillon d’apprentissage sera utilisé pour construire les algorithmes de régression sous contraintes et sélectionner leurs paramètres, celui de validation pourestimer leurs performances.

set.seed(5678)
perm <- sample(4601,3000)
app <- spam[perm,]
valid <- spam[-perm,]

2. Construire et analyser un algorithme de SVM

Il existe plusieurs fonctions sur R permettant de faire des SVM. Les deux plus utilisées sont certainement la fonction svm du package e1071 et la fonction ksvm du package kernlab. Nous présentons dans cette fiche la fonction ksvm qui présentel’avantage d’être prise en compte dans le package caret pour calibrer les paramètres d’un SVM.

La syntaxe est similaire aux autres fonctions de régression ou de discrimination : on utilise une formule pour indiquer la variable à expliquer et les variables explicatives. Il faut faire attention aux paramètres : par défaut, la fonction ksvm transforme les données avec le noyau gaussien, et le paramètre de coût \(C\) vaut \(1\). Pour traiter les données initiales, c’est-à-dire non transformées ou non transportées dans un espace plus grand, il suffit de spécifier kernel=“vanilladot” :

mod.svm <- ksvm(type~.,data=app,kernel="vanilladot",C=1)
##  Setting default kernel parameters
mod.svm
## Support Vector Machine object of class "ksvm" 
## 
## SV type: C-svc  (classification) 
##  parameter : cost C = 1 
## 
## Linear (vanilla) kernel function. 
## 
## Number of Support Vectors : 597 
## 
## Objective Function Value : -546.6826 
## Training error : 0.062

La sortie renseigne sur : – le type de SVM ajusté : C-svc (classification), puisque la variable à expliquer est ici qualitative ; – le paramètre cost, qui correspond au paramètre de coût C, vaut 1 ; – la transformation des données effectuées : ici le noyau est linéaire, les variablessont conservées en l’état ; – le nombre de vecteurs supports (Support Vectors), qui est le nombre de données à conserver pour calculer la fonction de classification : ici il vaut 599, cequi signifie qu’il suffit de mémoriser 599 exemples pour calculer la solution et non pas les 3000 exemples de la base d’apprentissage ; – la valeur de la fonction optimisée pour obtenir les vecteurs supports à l’optimum (cette valeur est rarement utile) ; – l’erreur de classification calculée sur les données d’apprentissage, ici 0.067.

3. Sélectionner les paramètres d’un SVM

Le SVM dépend de plusieurs paramètres, notamment le paramètre de coût \(C\) et le noyau, dont le choix se révèle crucial pour la performance de la méthode. La sélection de ces paramètres est la partie difficile de la mise en œuvre des SVM.

L’approche classique consiste à essayer plusieurs noyaux et à se donner une grille de valeurs pour les paramètres quantitatifs. Pour chaque jeu de paramètres et pour chaque noyau, on estime ensuite un critère (par exemple par validation croisée) et on choisit le jeu de paramètres qui optimise ce critère estimé. L’utilisateur doit choisir les valeurs de la grille, il faut en général faire quelques essais sur les données pour se donner les valeurs minimales et maximales de chaque grille de paramètres. Pour notre exemple, nous proposons de confronter 2 noyaux : le noyau polynomial et le noyau radial (le noyau linéaire étant un cas particulier du noyau polynomial). Les grilles pour chaque paramètre sont définies comme suit :

C <- c(0.1,1,10,100)
degree <- c(1,2,3)
scale <- 1
sigma <- c(0.0001,0.001,0.01,0.1,1)

Comme nous l’avons dit, l’utilisateur doit faire des choix. Nous choisissons ici de ne pas optimiser le paramètre scale du noyau polynomial et de le laisser à sa valeur par défaut 1. Ce paramètre n’est généralement pas le plus important. Pour chaque jeu de paramètres, nous allons estimer le taux de bon classement par validation croisée 3 blocs à l’aide du package caret (voir fiche 10.1). Nous commençons par le noyau polynomial :

library(caret)
## Loading required package: ggplot2
## 
## Attaching package: 'ggplot2'
## The following object is masked from 'package:kernlab':
## 
##     alpha
## Loading required package: lattice
gr.poly <- expand.grid(C=C,degree=degree,scale=scale)
ctrl <- trainControl(method="cv",number=3)
set.seed(123)
sel.poly <- train(type~.,data=app,method="svmPoly",trControl=ctrl,
tuneGrid=gr.poly)
sel.poly
## Support Vector Machines with Polynomial Kernel 
## 
## 3000 samples
##   57 predictor
##    2 classes: 'nonspam', 'spam' 
## 
## No pre-processing
## Resampling: Cross-Validated (3 fold) 
## Summary of sample sizes: 2000, 2001, 1999 
## Resampling results across tuning parameters:
## 
##   C      degree  Accuracy   Kappa    
##     0.1  1       0.9279930  0.8480626
##     0.1  2       0.9043359  0.7987107
##     0.1  3       0.9020023  0.7948546
##     1.0  1       0.9316593  0.8561199
##     1.0  2       0.8863373  0.7616311
##     1.0  3       0.8950066  0.7798967
##    10.0  1       0.9333286  0.8595239
##    10.0  2       0.8840009  0.7563813
##    10.0  3       0.8876699  0.7639519
##   100.0  1       0.9306623  0.8539060
##   100.0  2       0.8739983  0.7359073
##   100.0  3       0.8790019  0.7466107
## 
## Tuning parameter 'scale' was held constant at a value of 1
## Accuracy was used to select the optimal model using the largest value.
## The final values used for the model were degree = 1, scale = 1 and C = 10.

Les taux de bien classés estimés se trouvent dans la colonne Accuracy. On observe que ce taux est maximal pour degree=1 et C=1. Pour le noyau radial on a :

gr.radial <- expand.grid(C=C,sigma=sigma)
set.seed(345)
sel.radial <- train(type~.,data=app,method="svmRadial",
trControl=ctrl,tuneGrid=gr.radial)
sel.radial
## Support Vector Machines with Radial Basis Function Kernel 
## 
## 3000 samples
##   57 predictor
##    2 classes: 'nonspam', 'spam' 
## 
## No pre-processing
## Resampling: Cross-Validated (3 fold) 
## Summary of sample sizes: 2000, 1999, 2001 
## Resampling results across tuning parameters:
## 
##   C      sigma  Accuracy   Kappa      
##     0.1  1e-04  0.6086667  0.000000000
##     0.1  1e-03  0.7536652  0.421735478
##     0.1  1e-02  0.8989933  0.781491960
##     0.1  1e-01  0.7696675  0.463890913
##     0.1  1e+00  0.6096667  0.003109888
##     1.0  1e-04  0.7583335  0.434374323
##     1.0  1e-03  0.8996639  0.784189510
##     1.0  1e-02  0.9249983  0.841160325
##     1.0  1e-01  0.9013296  0.787240997
##     1.0  1e+00  0.7616695  0.438967988
##    10.0  1e-04  0.8976643  0.779621855
##    10.0  1e-03  0.9279983  0.847302008
##    10.0  1e-02  0.9309983  0.854409908
##    10.0  1e-01  0.9019946  0.790347214
##    10.0  1e+00  0.7653388  0.450553767
##   100.0  1e-04  0.9226629  0.835987885
##   100.0  1e-03  0.9333296  0.859231698
##   100.0  1e-02  0.9263286  0.844285705
##   100.0  1e-01  0.9046609  0.796227727
##   100.0  1e+00  0.7650058  0.450610779
## 
## Accuracy was used to select the optimal model using the largest value.
## The final values used for the model were sigma = 0.001 and C = 100.

On peut maintenant ajuster les SVM avec les paramètres optimaux pour les deux noyaux concurrents :

mod.poly <- ksvm(type~.,data=app,kernel="polydot",
kpar=list(degree=1,scale=1,offset=1),C=1,prob.model = TRUE)
 mod.radial <- ksvm(type~.,data=app,kernel="rbfdot",
kpar=list(sigma=0.01),C=10,prob.model = TRUE)

4. Faire de la prévision

La fonction predict.ksvm, que l’on peut appeler avec le raccourci predict, permet de prédire le label d’un nouveau courriel. On peut obtenir les labels prédits pour les individus de l’échantillon de validation avec :

prev.class.poly <- predict(mod.poly,newdata=valid)
prev.class.radial <- predict(mod.radial,newdata=valid)
prev.class.poly[1:10]
##  [1] spam    spam    spam    spam    spam    nonspam spam    spam    nonspam
## [10] spam   
## Levels: nonspam spam

Parmi les 10 premiers messages de l’échantillon de validation, 9 sont considérés comme des spams par les deux SVM. La fonction predict permet d’estimer la probabilité qu’un nouveau courriel soit un spam. On peut par exemple obtenir ces estimations pour les courriels de l’échantillon de validation avec les commandes :

prev.prob.poly <- predict(mod.poly,newdata=valid, type="probabilities")
prev.prob.radial <- predict(mod.radial,newdata=valid,
type="probabilities")
 round(head(prev.prob.poly),3)
##      nonspam  spam
## [1,]   0.415 0.585
## [2,]   0.332 0.668
## [3,]   0.397 0.603
## [4,]   0.404 0.596
## [5,]   0.388 0.612
## [6,]   0.729 0.271

La probabilité que le premier courriel de l’échantillon de validation soit un spam est estimée par le SVM polynomial à 0.944, celle pour le second à 0.588, etc. Par défaut, quand les probabilités sont supérieures à 0.5, le courriel est prédit comme spam.

5. Estimer les performances de l’algorithme

À partir des prévisions faites sur les données de validation dans la partie précédente, on peut obtenir une estimation de l’erreur de classification des deux SVM avec :

library(tidyverse)
## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr     1.1.0     ✔ readr     2.1.4
## ✔ forcats   1.0.0     ✔ stringr   1.5.0
## ✔ lubridate 1.9.2     ✔ tibble    3.1.8
## ✔ purrr     1.0.1     ✔ tidyr     1.3.0
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ ggplot2::alpha() masks kernlab::alpha()
## ✖ purrr::cross()   masks kernlab::cross()
## ✖ dplyr::filter()  masks stats::filter()
## ✖ dplyr::lag()     masks stats::lag()
## ✖ purrr::lift()    masks caret::lift()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
prev.class <- data.frame(poly=prev.class.poly,
radial=prev.class.radial,obs=valid$type)
prev.class %>% summarise_all(funs(err=mean(obs!=.))) %>%
select(-obs_err) %>% round(3)
## Warning: `funs()` was deprecated in dplyr 0.8.0.
## ℹ Please use a list of either functions or lambdas:
## 
## # Simple named list: list(mean = mean, median = median)
## 
## # Auto named with `tibble::lst()`: tibble::lst(mean, median)
## 
## # Using lambdas list(~ mean(., trim = .2), ~ median(., na.rm = TRUE))
##   poly_err radial_err
## 1    0.075      0.067

Les taux d’erreur sont relativement proches avec une légère préférence pour le noyau radial. De plus, ces taux d’erreur sont du même ordre que ceux obtenus par validation croisée sur les données d’apprentissage. On peut donc considérer qu’il n’y a pas de sur-ajustement pour ces deux SVM. D’autres indicateurs peuvent bien entendu être utilisés pour mesurer la performance, on peut par exemple comparer les courbes ROC. Cela peut s’effectuer avec le package pROC comme dans la fiche sur les forêts aléatoires . Nous proposons de le faire ici en ggplot à l’aide du package plotROC.

library(plotROC)
## Warning: package 'plotROC' was built under R version 4.2.3
prev.prob <- data.frame(poly=prev.prob.poly[,2],
radial=prev.prob.radial[,2],obs=valid$type)
df.roc <- prev.prob %>% gather(key=Methode,value=score,poly,radial)
ggplot(df.roc)+aes(d=obs,m=score,color=Methode)+geom_roc()+theme_classic()
## Warning in verify_d(data$d): D not labeled 0/1, assuming nonspam = 0 and spam =
## 1!
## Warning: The following aesthetics were dropped during statistical transformation: d, m
## ℹ This can happen when ggplot fails to infer the correct grouping structure in
##   the data.
## ℹ Did you forget to specify a `group` aesthetic or to convert a numerical
##   variable into a factor?