TP2 : classification
Nous testerons dans ce TP 3 méthodes de classification supervisée : la méthode des \(k\) plus proches voisins, la régression logistique et les SVM. Les parties et sous-parties sont indépendantes et vous pouvez les traiter dans l’ordre que vous souhaitez.
Ce TP est largement inspiré des exercices d’applications du livre :
James, G., Witten, D., Hastie, T. and Tibshirani, R. An introduction to statistical learning.
Présentation des données
Nous nous intéressons à des données de tumeurs correspondant à des échantillons de tissus cancéreux de différents types. L’objectif est d’identifier le type de tumeurs à partir de données d’expressions géniques.
- La variable à prédire est donc le type de cancer (quatre types possibles, codés 1, 2, 3 ou 4).
- Les variables d’entrée sont les niveaux d’expression de gènes (2308 gènes ont été retenus pour l’étude).
Nous commençons par charger le jeu de données sous R.
library(ISLR)
data(Khan)
summary(Khan)## Length Class Mode
## xtrain 145404 -none- numeric
## xtest 46160 -none- numeric
## ytrain 63 -none- numeric
## ytest 20 -none- numericattach(Khan) # pour pouvoir utiliser directement les variables xtrain,...
ytrain = factor(ytrain) # pour coder les variables qualitatives comme facteurs
ytest = factor(ytest)
ntrain = length(ytrain)
ntest = length(ytest)Nous observons que le jeu de données contient déjà un échantillon d’apprentissage et un échantillon de test. Quelle est la taille de l’échantillon d’apprentissage ? Celle de l’échantillon de test ? Dans l’échantillon d’apprentissage, combien y’a t-il de tissus correspondant à chacun des 4 types de cancer ?
Dans un premier temps, nous nous intéresserons à la classification binaire, c’est-à-dire que nous nous intéresserons au type 2 versus tout les autres types.
Classification binaire
Nous commençons par générer une nouvelles variable z en regroupant les classes 1, 3 et 4 dans une seule classe (ayant pour label 0).
ztrain = ifelse(ytrain==2,T,F)
ztest = ifelse(ytest==2,T,F)
summary(ztrain)## Mode FALSE TRUE
## logical 40 23k-plus proches voisins
Tester la méthode des \(k\)-plus proches voisins pour différentes valeurs de \(k\). Nous pouvons utiliser par exemple la fonction knn() du package ‘class’. Calculer son taux d’erreur. On pourra également sélectionner \(k\) par validation croisée.
knnCV <- function(xtrain, ytrain, kvect, plot=TRUE){
ntrain = nrow(xtrain)
nbk = length(kvect)
risque = rep(NA,nbk)
for (j in 1:nbk){
preds <- sapply(1:ntrain, function(i)
knn(xtrain[-i,],xtrain[i,],ytrain[-i],k=kvect[j]))
risque[j] <- mean(preds != ytrain)
}
khat <- kvect[which.min(risque)]
print(cbind(k=kvect, risk=risque))
if(plot)
{
plot(kvect, risque, type="b", pch=16)
points(khat, min(risque), col="red")
}
khat
}Régression logistique
Nous rappelons que le modèle logistique consiste à modéliser la fonction de régression \(\eta(X)=\mathbb P(Y=1|X=x)\) comme une fonction dépendant linéairement des observations, c’est-à-dire qu’il existe \(\beta_0,...,\beta_d\) tels que : \[ \eta(x) = g(\beta_0+\beta_1x_1+...+\beta_dx_d),\quad \forall x=(x_1,...,x_d)'\in\mathbb R^d, \] avec \[ g(t) = \frac{e^t}{1+e^t}. \]
La fonction glm, qui fonctionne comme la fonction lm permet d’estimer le paramètre \(\beta=(\beta_0,\beta_1,...,\beta_d)\) par maximisation de la vraisemblance.
log1 = glm(ztrain~xtrain,family='binomial')
summary(log1)Essayez d’interpréter la sortie du logiciel. Il semble qu’il y ait un problème lié à la dimension importante des données. Pour remédier à cela, nous pouvons implémenter le GLM Lasso ou le GLM Ridge.
# GLM Lasso
library('glmnet')## Loading required package: Matrix## Loaded glmnet 4.1-1cv.out = cv.glmnet(xtrain,ztrain,alpha=1)
plot(cv.out)
lambda_opt.lasso = cv.out$lambda.min
lasso.mod.opt = glmnet(xtrain,ztrain,alpha = 1, lambda = lambda_opt.lasso,family="binomial")
lasso.pred.opt = predict(lasso.mod.opt,newx = xtest,type='response') # prédit les probabilités P(Y=1|X)=eta(X)
ztest.pred = (lasso.pred.opt>0.5)
table(ztest.pred,ztest)## ztest
## ztest.pred FALSE TRUE
## FALSE 14 1
## TRUE 0 5Faire de même avec le GLM Ridge.
# GLM Ridge
library('glmnet')
cv.out.Ridge = cv.glmnet(xtrain,ztrain,alpha=0)
plot(cv.out.Ridge)
lambda_opt.Ridge = cv.out.Ridge$lambda.min
Ridge.mod.opt = glmnet(xtrain,ztrain,alpha = 0, lambda = lambda_opt.Ridge,family="binomial")
Ridge.pred.opt = predict(Ridge.mod.opt,newx = xtest,type='response')
ztest.pred.Ridge = (Ridge.pred.opt>0.5)
table(ztest.pred.Ridge,ztest)## ztest
## ztest.pred.Ridge FALSE TRUE
## FALSE 14 0
## TRUE 0 6SVM
Nous testons maintenant les SVM à l’aide de la fonction svm du package e1071.
library("e1071")
# Estimation des parametres
dat1=data.frame(x=xtrain , z=as.factor(ztrain))
svmlin = svm(z~., data=dat1, kernel="linear")
table(svmlin$fitted,ztrain)## ztrain
## FALSE TRUE
## FALSE 40 0
## TRUE 0 23# Prediction
dat1.test=data.frame(x=xtest , y=as.factor(ztest))
ztest.pred.svmlin=predict (svmlin , newdata =dat1.test)
table(ztest.pred.svmlin,ztest)## ztest
## ztest.pred.svmlin FALSE TRUE
## FALSE 14 0
## TRUE 0 6Nous remarquons que les SVM linéaires réalisent une classification parfaite, que ce soit sur l’échantillon d’apprentissage ou l’échantillon de test. Le fait que la classification soit parfaite sur l’échantillon d’apprentissage ne caractérise pas les bonnes performances du classifieurs car nous ne pouvons pas être certains de ne pas faire du sur-apprentissage. Cependant, dans le cas des SVM, cela indique que l’échanillon est linéairement séparable, ce qui est souvent le cas en grande dimension lorsque l’on a peu d’observations. Il n’est donc pas utile d’aller plus loin dans l’examen des SVM dans ce cas-là. Pour d’autres jeu de données, il pourra être utile de tester d’autres noyaux, par exemple le noyau polynomial (kernel="polynomial") ou le noyau gaussien (kernel="radial").
Le paramètre \(C\) apparaissant dans la formulation vue en cours du problème d’optimisation est codée dans le paramètre cost. Nous ne l’avons pas précisé ici. Regarder dans l’aide de la fonction pour savoir comment il est fixé par défaut. Vous pouvez tester différentes valeurs. Cela change-t’il le résultat ?
Classification multi-label
Nous revenons maintenant au jeu de données initial, c’est-à-dire où nous devons prédire le type de cancer parmi 4 possibilités. Etant donné qu’il existe des extensions de la régression logistique au cas multi-label mais qu’elles ne sont pas vraiment utilisées en pratique, nous nous limiterons dans cette partie à la méthode des \(k\)-plus proches voisins et aux SVM.
\(k\) plus proches voisins
La fonction knn prend en charge le cas multiclasse.
library(class)
# k=20
res.multi20 = knn(xtrain,xtest,ytrain,k=20)
table(res.multi20,ytest)## ytest
## res.multi20 1 2 3 4
## 1 1 0 0 0
## 2 2 6 4 0
## 3 0 0 2 0
## 4 0 0 0 5sum(res.multi20!=ytest)/length(ytest)## [1] 0.3# Pour k=20 les k-PPV predisent trop souvent le type 2.
# k=1
res.multi1 = knn(xtrain,xtest,ytrain,k=1)
table(res.multi1,ytest)## ytest
## res.multi1 1 2 3 4
## 1 3 0 0 0
## 2 0 4 0 0
## 3 0 0 2 0
## 4 0 2 4 5sum(res.multi1!=ytest)/length(ytest)## [1] 0.3# Pour k=1, les k-PPV predisent trop souvent le type 4.
# k selectionne par validation croisee
k.CV = knnCV(xtrain,ytrain,kvect=1:20) ## k risk
## [1,] 1 0.07936508
## [2,] 2 0.14285714
## [3,] 3 0.06349206
## [4,] 4 0.07936508
## [5,] 5 0.06349206
## [6,] 6 0.12698413
## [7,] 7 0.09523810
## [8,] 8 0.12698413
## [9,] 9 0.14285714
## [10,] 10 0.14285714
## [11,] 11 0.22222222
## [12,] 12 0.17460317
## [13,] 13 0.20634921
## [14,] 14 0.17460317
## [15,] 15 0.19047619
## [16,] 16 0.19047619
## [17,] 17 0.23809524
## [18,] 18 0.20634921
## [19,] 19 0.25396825
## [20,] 20 0.33333333
res.multi.CV = knn(xtrain,xtest,ytrain,k=k.CV)
table(res.multi.CV,ytest)## ytest
## res.multi.CV 1 2 3 4
## 1 3 0 0 0
## 2 0 6 0 0
## 3 0 0 4 0
## 4 0 0 2 5sum(res.multi.CV!=ytest)/length(ytest)## [1] 0.1SVM
La version des SVM vue en cours ne permet pas de faire de la classification multi-classe. Cependant, plusieurs méthodes existent pour étendre les SVM au cas multiclasse. Celle implémentée dans la fonction svm()que nous allons utiliser est la classification un contre un. Supposons que notre variable à prédire \(Y\) puisse prendre \(K\) valeurs différentes. Cette méthode consiste à entrainer un classifieur linéaire pour chaque couple de labels possible.
Par exemple pour notre cas où \(K=4\), nous entrainons 6 classifieurs :
- un classifieur pour différencier le type 1 du type 2;
- un pour différencier le type 2 du type 3;
- un pour type 1 versus type 4;
- un pour type 2 versus type 3;
- un pour type 2 versus type 4;
- et un pour type 3 versus type 4.
Chaque classifieur est entrainé sur le sous-échantillon de l’échantillon d’apprentissage contenant les individus \(i\) ayant les valeurs de \(Y_i\) correspondantes (par exemple les paramètres du classifieur permettant de différencier le type 1 du type 2 sont déterminés uniquement grâce aux individus de l’échantillon d’apprentissage tels que \(Y_i\in\{1,2\}\)).
La classification finale d’un point \(x\) est attribuée à la majorité (la classe 1 est attribuée par exemple si, parmi les 6 classifieurs, la classe 1 est plus souvent attribuée au point \(x\) que toutes les autres). La fonction ‘svm()’ que nous utilisons est basée sur cette méthode.
library("e1071")
# Estimation des parametres
dat2=data.frame(x=xtrain , y=as.factor(ytrain))
svmlin.multi = svm(y~., data=dat2, kernel="linear")
table(svmlin.multi$fitted,ytrain)## ytrain
## 1 2 3 4
## 1 8 0 0 0
## 2 0 23 0 0
## 3 0 0 12 0
## 4 0 0 0 20# Prediction
dat2.test=data.frame(x=xtest , y=as.factor(ytest))
ytest.pred.svmlin.multi=predict (svmlin.multi , newdata =dat2.test)
table(ytest.pred.svmlin.multi,ytest)## ytest
## ytest.pred.svmlin.multi 1 2 3 4
## 1 3 0 0 0
## 2 0 6 2 0
## 3 0 0 4 0
## 4 0 0 0 5sum(ytest.pred.svmlin.multi!=ytest)/length(ytest)## [1] 0.1La classifieur fait quelques erreurs sur l’échantillon de test. Qu’indique le fait qu’il ne fasse aucune erreur sur l’échantillon d’apprentissage ? Dans ce cas-là, est-il utile d’essayer de mieux calibrer le paramètre \(C\) (option cost de la fonction) ?