Etant donné un échantillon séparable \((x_i,y_i),i=1,\dots,n\) où \(x_i\in\mathbb R^2\) et \(y_i\in\{0,1\}\), on rappelle que l’algorithme SVM consiste à trouver un hyperplan \[\langle w,x\rangle+b=0\] qui sépare les \(x_i\) en fonction des \(y_i\). On génère des données selon
library(tidyverse)
n <- 20
set.seed(123)
X1 <- scale(runif(n))
set.seed(567)
X2 <- scale(runif(n))
Y <- rep(0,n)
Y[X1>X2] <- 1
Y <- as.factor(Y)
donnees <- data.frame(X1=X1,X2=X2,Y=Y)
Et on considère la svm suivante :
library(e1071)
mod.svm <- svm(Y~.,data=donnees,kernel="linear",cost=10000000000)
p <- ggplot(donnees)+aes(x=X2,y=X1,color=Y)+geom_point()+theme_classic()
p
Les vecteurs supports se trouvent dans l’élément index de la fonction svm :
ind.svm <- mod.svm$index
sv <- donnees %>% slice(ind.svm)
sv
p1 <- p+geom_point(data=sv,aes(x=X2,y=X1),color="blue",size=2)
On peut ainsi représenter la marge en traçant les droites qui passent par ces points.
sv1 <- sv[,2:1]
b <- (sv1[1,2]-sv1[2,2])/(sv1[1,1]-sv1[2,1])
a <- sv1[1,2]-b*sv1[1,1]
a1 <- sv1[3,2]-b*sv1[3,1]
p1+geom_abline(intercept = c(a,a1),slope=b,col="blue",size=1)
plot(mod.svm,data=donnees,grid=250)
L’hyperplan séparateur est d’équation
\[\langle w^\star,x\rangle+b^\star=0\] avec
\[w^\star=\sum_{i=1}^n\alpha_i^\star y_ix_i\] et \(b^\star\) solution de \(y_i<w^\star,x_i>+b=0\) (pour \(\alpha_i^\star\neq 0\)). La règle s’écrit donc \[g(x)=1_{<w^\star,x>+b^\star\leq 0}-1_{<w^\star,x>+b^\star> 0}.\] L’élément mod.svm$coefs contient les coefficients \(\alpha_i^\star y_i\) pour chaque vecteur support. On peut ainsi récupérer l’équation de l’hyperplan et faire la prévision avec
w <- apply(mod.svm$coefs*donnees[mod.svm$index,1:2],2,sum)
w
## X1 X2
## -1.745100 2.136029
b <- -mod.svm$rho
L’hyperplan séparateur a donc pour équation : \[-1.74x_1+2.12x_2-0.40=0.\]
Il suffit de calculer \(<w^\star,x>+b\) et de prédire en fonction du signe de cette valeur :
newX <- data.frame(X1=-0.5,X2=0.5)
sum(w*newX)+b
## [1] 1.537053
On prédira le groupe 0 pour ce nouvel individu.
decision.values = TRUE
.predict(mod.svm,newX,decision.values = TRUE)
## 1
## 0
## attr(,"decision.values")
## 0/1
## 1 1.537053
## Levels: 0 1
Plus cette valeur est élevée, plus on est loin de l’hyperplan. On peut donc l’interpréter comme un score. Comme souvent, il est possible d’obtenir une estimation des probabilités d’être dans les groupes 0 et 1 à partir de ce score, il “suffit” de ramener ce score sur l’échelle \([0,1]\) avec des transformations de type logit par exemple. Pour la svm, ces probabilités sont obtenues en ajustant un modèle logistique sur les scores \(S(x)\) : \[P(Y=1|X=x)=\frac{1}{1+\exp(aS(x)+b)}.\]
probability=TRUE
dans la fonction svm.mod.svm1 <- svm(Y~.,data=donnees,kernel="linear",cost=10000000000,probability=TRUE)
predict(mod.svm1,newX,decision.values=TRUE,probability=TRUE)
## 1
## 0
## attr(,"decision.values")
## 0/1
## 1 1.537053
## attr(,"probabilities")
## 0 1
## 1 0.8294474 0.1705526
## Levels: 0 1
On peut retrouver ces probabilités avec :
score.newX <- sum(w*newX)+b
1/(1+exp(-(mod.svm1$probB+mod.svm1$probA*score.newX)))
## [1] 0.1705526
On considère le jeu de données suivant où le problème est d’expliquer \(Y\) par \(V1\) et \(V2\).
n <- 750
set.seed(1)
X <- matrix(runif(n*2,-2,2),ncol=2) %>% as.data.frame()
Y <- rep(0,n)
cond <- (X$V1^2+X$V2^2)<=2.8
Y[cond] <- rbinom(sum(cond),1,0.9)
Y[!cond] <- rbinom(sum(!cond),1,0.1)
df <- X %>% mutate(Y=as.factor(Y))
ggplot(df)+aes(x=V1,y=V2,color=Y)+geom_point()+theme_classic()
set.seed(1234)
ind.train <- sample(nrow(df),500)
train <- df %>% slice(ind.train)
test <- df %>% slice(-ind.train)
mod.svm0 <- svm(Y~.,data=train,kernel="linear",cost=1)
plot(mod.svm0,train,grid=250)
La svm linéaire ne permet pas de séparer les groupes (on pouvait s’y attendre).
L’astuce du noyau consiste à envoyer les données dans un espace de représentation (feature space) dans lequel on espère que les données soient linéairement séparables.
mod.svm1 <- svm(Y~.,data=train,kernel="radial",gamma=1,cost=1)
plot(mod.svm1,train,grid=250)
Le noyau radial permet de mettre en évidence une séparation non linéaire.
mod.svm2 <- svm(Y~.,data=train,kernel="radial",gamma=1,cost=0.0001)
mod.svm3 <- svm(Y~.,data=train,kernel="radial",gamma=1,cost=1)
mod.svm4 <- svm(Y~.,data=train,kernel="radial",gamma=1,cost=10000)
plot(mod.svm2,train,grid=250)
plot(mod.svm3,train,grid=250)
plot(mod.svm4,train,grid=250)
mod.svm2$nSV
## [1] 241 241
mod.svm3$nSV
## [1] 124 122
mod.svm4$nSV
## [1] 86 84
Le nombre de vecteurs supports diminue lorsque \(C\) augmente. Une forte valeur de \(C\) autorise moins d’observations à être dans la marge, elle a donc tendance à diminuer (risque de surapprentissage).
set.seed(1234)
tune.out <- tune(svm,Y~.,data=train,kernel="radial",ranges=list(cost=c(0.1,1,10,100,1000),gamma=c(0.5,1,2,3,4)))
summary(tune.out)
##
## Parameter tuning of 'svm':
##
## - sampling method: 10-fold cross validation
##
## - best parameters:
## cost gamma
## 10 1
##
## - best performance: 0.122
##
## - Detailed performance results:
## cost gamma error dispersion
## 1 1e-01 0.5 0.204 0.04299871
## 2 1e+00 0.5 0.158 0.04848826
## 3 1e+01 0.5 0.130 0.06128259
## 4 1e+02 0.5 0.130 0.05517648
## 5 1e+03 0.5 0.132 0.05432413
## 6 1e-01 1.0 0.206 0.04623611
## 7 1e+00 1.0 0.140 0.04988877
## 8 1e+01 1.0 0.122 0.04756282
## 9 1e+02 1.0 0.140 0.06110101
## 10 1e+03 1.0 0.144 0.05947922
## 11 1e-01 2.0 0.196 0.05059644
## 12 1e+00 2.0 0.134 0.04993329
## 13 1e+01 2.0 0.142 0.05452828
## 14 1e+02 2.0 0.150 0.06411795
## 15 1e+03 2.0 0.146 0.05660781
## 16 1e-01 3.0 0.190 0.05354126
## 17 1e+00 3.0 0.134 0.05420127
## 18 1e+01 3.0 0.152 0.06477311
## 19 1e+02 3.0 0.156 0.06310485
## 20 1e+03 3.0 0.162 0.07083627
## 21 1e-01 4.0 0.184 0.05316641
## 22 1e+00 4.0 0.136 0.05641119
## 23 1e+01 4.0 0.156 0.06719788
## 24 1e+02 4.0 0.154 0.07121173
## 25 1e+03 4.0 0.168 0.05181162
La sélection est faite en minimisant l’erreur de classification par validation croisée 10 blocs.
Pour caret il faut utiliser la méthode svmRadial du package kernlab.
library(caret)
C <- c(0.001,0.01,1,10,100,1000)
sigma <- c(0.5,1,2,3,4)
gr <- expand.grid(C=C,sigma=sigma)
ctrl <- trainControl(method="cv")
res.caret1 <- train(Y~.,data=train,method="svmRadial",trControl=ctrl,tuneGrid=gr,prob.model=TRUE)
res.caret1
## Support Vector Machines with Radial Basis Function Kernel
##
## 500 samples
## 2 predictor
## 2 classes: '0', '1'
##
## No pre-processing
## Resampling: Cross-Validated (10 fold)
## Summary of sample sizes: 450, 450, 450, 451, 450, 450, ...
## Resampling results across tuning parameters:
##
## C sigma Accuracy Kappa
## 1e-03 0.5 0.7199040 0.4478997
## 1e-03 1.0 0.7279040 0.4634977
## 1e-03 2.0 0.7599872 0.5256372
## 1e-03 3.0 0.7719888 0.5490030
## 1e-03 4.0 0.7839496 0.5722969
## 1e-02 0.5 0.7179040 0.4441092
## 1e-02 1.0 0.7219040 0.4520669
## 1e-02 2.0 0.7579872 0.5217628
## 1e-02 3.0 0.7719480 0.5489356
## 1e-02 4.0 0.7839496 0.5722977
## 1e+00 0.5 0.8479560 0.6967004
## 1e+00 1.0 0.8698760 0.7395469
## 1e+00 2.0 0.8698351 0.7396646
## 1e+00 3.0 0.8698351 0.7394599
## 1e+00 4.0 0.8678743 0.7355593
## 1e+01 0.5 0.8699152 0.7393388
## 1e+01 1.0 0.8758760 0.7513048
## 1e+01 2.0 0.8758743 0.7512531
## 1e+01 3.0 0.8658727 0.7313577
## 1e+01 4.0 0.8578319 0.7153009
## 1e+02 0.5 0.8659544 0.7315761
## 1e+02 1.0 0.8679136 0.7350568
## 1e+02 2.0 0.8599136 0.7195116
## 1e+02 3.0 0.8360312 0.6713023
## 1e+02 4.0 0.8419120 0.6832832
## 1e+03 0.5 0.8818743 0.7631046
## 1e+03 1.0 0.8618319 0.7229477
## 1e+03 2.0 0.8259904 0.6504446
## 1e+03 3.0 0.8158255 0.6303879
## 1e+03 4.0 0.8197471 0.6385549
##
## Accuracy was used to select the optimal model using the largest value.
## The final values used for the model were sigma = 0.5 and C = 1000.
On peut également répéter plusieurs fois la validation croisée pour stabiliser les résultats (on parallélise avec doParallel) :
library(doParallel) ## pour paralléliser
cl <- makePSOCKcluster(4)
registerDoParallel(cl)
set.seed(12345)
ctrl <- trainControl(method="repeatedcv",number=10,repeats=5)
res.caret2 <- train(Y~.,data=train,method="svmRadial",trControl=ctrl,tuneGrid=gr,prob.model=TRUE)
on.exit(stopCluster(cl))
res.caret2
## Support Vector Machines with Radial Basis Function Kernel
##
## 500 samples
## 2 predictor
## 2 classes: '0', '1'
##
## No pre-processing
## Resampling: Cross-Validated (10 fold, repeated 5 times)
## Summary of sample sizes: 450, 450, 450, 450, 451, 450, ...
## Resampling results across tuning parameters:
##
## C sigma Accuracy Kappa
## 1e-03 0.5 0.7215622 0.4506883
## 1e-03 1.0 0.7203459 0.4484636
## 1e-03 2.0 0.7559566 0.5174585
## 1e-03 3.0 0.7711654 0.5469831
## 1e-03 4.0 0.7799579 0.5641196
## 1e-02 0.5 0.7163704 0.4406961
## 1e-02 1.0 0.7195541 0.4469284
## 1e-02 2.0 0.7563485 0.5181859
## 1e-02 3.0 0.7723654 0.5493151
## 1e-02 4.0 0.7795657 0.5633499
## 1e+00 0.5 0.8528283 0.7061116
## 1e+00 1.0 0.8680293 0.7359375
## 1e+00 2.0 0.8696221 0.7391926
## 1e+00 3.0 0.8680381 0.7360596
## 1e+00 4.0 0.8672224 0.7343738
## 1e+01 0.5 0.8688133 0.7371799
## 1e+01 1.0 0.8712142 0.7421786
## 1e+01 2.0 0.8643976 0.7283885
## 1e+01 3.0 0.8508434 0.7013784
## 1e+01 4.0 0.8448264 0.6893241
## 1e+02 0.5 0.8703982 0.7403436
## 1e+02 1.0 0.8651816 0.7297946
## 1e+02 2.0 0.8431950 0.6857540
## 1e+02 3.0 0.8344192 0.6683523
## 1e+02 4.0 0.8351301 0.6696499
## 1e+03 0.5 0.8736384 0.7466156
## 1e+03 1.0 0.8504126 0.6999760
## 1e+03 2.0 0.8299539 0.6587048
## 1e+03 3.0 0.8166636 0.6322624
## 1e+03 4.0 0.8147216 0.6287175
##
## Accuracy was used to select the optimal model using the largest value.
## The final values used for the model were sigma = 0.5 and C = 1000.
scale <- 1
C <- c(1,10,100)
degree <- c(1:3)
gr <- expand.grid(C=C,scale=scale,degree=degree)
ctrl <- trainControl(method="repeatedcv",number=10,repeats=5)
cl <- makePSOCKcluster(4)
registerDoParallel(cl)
set.seed(12345)
res.caret3 <- train(Y~.,data=train,method="svmPoly",trControl=ctrl,tuneGrid=gr,prob.model=TRUE)
on.exit(stopCluster(cl))
mod.svm0 <- svm(Y~.,data=train,kernel="linear",cost=1,probability=TRUE)
prev.lin <- predict(mod.svm0,newdata=test,probability=TRUE)
prev.lin <- attr(prev.lin, "probabilities")[,2]
prev.rad <- predict(res.caret2,newdata=test,type="prob")[,2]
prev.poly <- predict(res.caret3,newdata=test,type="prob")[,2]
score <- data.frame(lin=prev.lin,rad=prev.rad,poly=prev.poly,obs=as.numeric(as.character(test$Y))) %>%
gather(key="Method",value="score",-obs) %>%
mutate(prev=recode(as.character(score>0.5),"TRUE"=1,"FALSE"=0))
Courbe ROC
library(plotROC)
ggplot(score)+aes(d=obs,m=score,color=Method)+geom_roc()+theme_classic()
AUC
score %>% group_by(Method) %>% summarize(AUC=pROC::auc(obs,score)) %>% arrange(desc(AUC))
Erreur de classification
score %>% group_by(Method) %>% summarize(Erreur=mean(prev!=obs)) %>% arrange(Erreur)
On peut utliser l’option metric de la fonction train :
ctrl <- trainControl(method="cv",classProbs = TRUE,summary = twoClassSummary)
cl <- makePSOCKcluster(3)
registerDoParallel(cl)
train1 <- train %>% mutate(Y=fct_recode(Y,G0="0",G1="1"))
res.caret4 <- train(Y~.,data=train1,method="svmPoly",trControl=ctrl,tuneGrid=gr,prob.model=TRUE,metric="ROC")
on.exit(stopCluster(cl))
res.caret4
## Support Vector Machines with Polynomial Kernel
##
## 500 samples
## 2 predictor
## 2 classes: 'G0', 'G1'
##
## No pre-processing
## Resampling: Cross-Validated (10 fold)
## Summary of sample sizes: 450, 450, 450, 450, 450, 450, ...
## Resampling results across tuning parameters:
##
## C degree ROC Sens Spec
## 1 1 0.5156154 0.1916667 0.9340000
## 1 2 0.8896410 0.8795000 0.8573846
## 1 3 0.8870897 0.8836667 0.8535385
## 10 1 0.4603333 0.1400000 0.9269231
## 10 2 0.8898013 0.8795000 0.8573846
## 10 3 0.8877308 0.8795000 0.8573846
## 100 1 0.4943077 0.2155000 0.8487692
## 100 2 0.8898013 0.8836667 0.8612308
## 100 3 0.8880641 0.8795000 0.8612308
##
## Tuning parameter 'scale' was held constant at a value of 1
## ROC was used to select the optimal model using the largest value.
## The final values used for the model were degree = 2, scale = 1 and C = 10.