On va profiter de la quarantaine pour mettre en ligne un billet sur la courbe ROC, la receiver operating characteristic. Considérons une petite base de données avec n=10 observations, deux variables continues, x_1 et x_2, et la variable d’intérêt binaire y\in\{0,1\}. On peut représenter les points dans le plan (x_1,x_2), et on utilise une couleur différente pour y\in\{0,1\}.
x1 = c(.4,.55,.65,.9,.1,.35,.5,.15,.2,.85) x2 = c(.85,.95,.8,.87,.5,.55,.5,.2,.1,.3) y = c(1,1,1,1,1,0,0,1,0,0) df = data.frame(x1=x1,x2=x2,y=as.factor(y)) plot(x1,x2,col=c("red","blue")[1+y],pch=19,cex=1.5) |
On peut alors faire une régression logistique, de telle sorte que \mathbb{P}(Y=1|x_1,x_2)=\frac{e^{\beta_0+\beta_1x_1+\beta_2x_2}}{1+e^{\beta_0+\beta_1x_1+\beta_2x_2}}On peut visualiser l’ensemble des points (x_1,x_2) pour lesquels \mathbb{P}(Y=0|x_1,x_2)=\mathbb{P}(Y=1|x_1,x_2)(on a alors autant de chance d’être rouge et bleu) soit \beta_0+\beta_1x_1+\beta_2x_2=0 qui correspond à une droite,
reg = glm(y~x1+x2,data=df,family=binomial(link = "logit")) b = coefficients(reg) abline(a=-b[1]/b[3],b=-b[2]/b[3]) |
On peut alors représenter y_i en fonction du score, i.e. l’estimation de \mathbb{P}(Y=1|x_{1,i},x_{2,i}),
Y = df$y S = predict(reg,type="response") plot(S,y,xlab="probabilité prédite",ylab="y") |
On va alors se donner un seuil (par exemple 50\%) : si la probabilité que Y prenne la valeur 1 excède le seuil, on prédit 1 (et sinon 0). Sur la figure ci-dessus, on a alors 4 sortes de points : ceux à gauche du seuil (et qui sont prédits 0), qui sont bien classés s’ils sont en bas, et mal classés en haut; à droite du seuil (et qui sont prédits 1), ils sont bien classés s’ils sont en haut, et mal classés en bas [dans le code ci-dessous, le symbole > désigne l’opérateur “supérieur“, qui malheureusement ne passe pas dans cet éditeur]
seuil = .5 Yhat = (S>seuil)*1 plot(S,y,xlab="probabilité prédite",ylab="y",pch=19, col=c("red","blue")[1+(y==Yhat)]) abline(v=seuil,lty=2) |
Les couleurs reflètent le bon ou mauvais classement : les points rouges correspondent à des erreurs de classement. On peut retrouver tout ça dans le tableau de contingence ci-dessous, qui correspond au tableau standard d’un test d’hypothèse
table(Yhat,Y) Y Yhat 0 1 0 3 1 1 1 5 |
Ce qui va nous intéresser ici à deux grandeurs particulières : le taux de faux positifs et le taux de vrais positifs,
FP=sum((Ps==1)*(Y==0))/sum(Y==0) TP=sum((Ps==1)*(Y==1))/sum(Y==1) |
On a obtenu ce tableau à un seuil donné (ici 50\%) mais on peut regarder ce qui se passe lors que le seuil change, comme sur l’animation ci-dessous, où on trace, à droite, le taux de vrais positifs (sur l’axe y) en fonction du taux de faux positifs (sur l’axe x
L’ensemble des points donne la courbe ROC.
roc.curve=function(s,print=FALSE){ Ps=(S>s)*1 FP=sum((Ps==1)*(Y==0))/sum(Y==0) TP=sum((Ps==1)*(Y==1))/sum(Y==1) if(print==TRUE){ print(table(Observed=Y,Predicted=Ps)) } vect=c(FP,TP) names(vect)=c("FPR","TPR") return(vect)} u = seq(0,1,length=251) V = Vectorize(roc.curve)(u) plot(t(V),type="s",xlab="Faux Positifs",ylab="Vrais Positifs") segments(0,0,1,1,col="light blue") |
On peut vérifier que le point qu’on avait obtenu avec un seuil de 50\% est bien sur la courbe
table(Yhat,Y) Y Yhat 0 1 0 3 1 1 1 5 (FP = sum((Yhat)*(Y==0))/sum(Y==0)) [1] 0.25 (TP = sum((Yhat==1)*(Y==1))/sum(Y==1)) [1] 0.8333333 abline(v=FP,lty=2,col="blue") abline(h=TP,lty=2,col="blue") points(FP,TP,pch=19,cex=1.5) |
Bien entendu, il y a (beaucoup) de packages R qui permettent d’avoir cette courbe,
library(ROCR) pred = prediction(S,Y) plot(performance(pred,"tpr","fpr")) |
Une grandeur intéressante est appelée aire sous la courbe (ou AUC) qu’on peut calculer ici à la main (on a une simple fonction en escalier)
p1 = roc.curve(1/3) p2 = roc.curve(.7) p2[1]*p2[2]+(p1[1]-p2[1])*p1[2]+(1-p1[1]) [1] 0.875 |
mais qu’on peut avoir automatiquement
auc.perf = performance(pred, measure = "auc") auc.perf@y.values[[1]] [1] 0.875 |
Allez, tentons un autre classifieur : toujours une régression logistique, mais sur un facteur obtenu en coupant la seconde variable en deux, \boldsymbol{1}_{[s,\infty)}(x_2)
reg = glm(y~I(x2>.525),data=df,family=binomial(link = "logit")) abline(h=.525) |
La droite horizontale n’est plus la droite qui donne autant de chance d’être rouge que bleu, mais qui coupe la variable x_2. Ici, on prédit juste deux valeurs : 40\% de chance d’être bleu en bas, et 80\% de chance d’être bleu, en haut. Si on représente observations y_i en fonction des probabilités prédites, on obtient
Y = df$y S = predict(reg,type="response") plot(S,y,xlab="probabilité prédite",ylab="y",xlim=0:1) |
Avec un seuil à 50\%, on obtient le tableau de contingence suivant (avec 3 erreurs, contre 2 auparavant)
seuil = .5 Yhat = (S>seuil)*1 table(Yhat,Y) Y Yhat 0 1 0 3 2 1 1 4 |
Si on trace la courbe ROC, on obtient
roc.curve=function(s,print=FALSE){ Ps=(S>s)*1 FP=sum((Ps==1)*(Y==0))/sum(Y==0) TP=sum((Ps==1)*(Y==1))/sum(Y==1) if(print==TRUE){ print(table(Observed=Y,Predicted=Ps)) } vect=c(FP,TP) names(vect)=c("FPR","TPR") return(vect)} u = seq(0,1,length=251) V = Vectorize(roc.curve)(u) plot(t(V),type="l",xlab="Faux Positifs",ylab="Vrais Positifs") segments(0,0,1,1,col="light blue") |
Cette fois, la courbe n’est plus constante par morceaux, mais linéaire par morceaux, et continue… L’interprétation est un peu plus subtile: cette fois, on a deux régions de l’espace, et dans chaque région, on ne sait pas trop comment distinguer (la probabilité est plate partout ici, contraitement à la régression précédante). Autrement dit, dans cette région, on a une proba constante, par exemple 40\% (en bas) : quand on doit prévoir pour un individu dans cette région, on lui attribue les valeurs \{0,1\} respectivement avec les probabilités \{40\%,60\%\}. Quand on a une probabilité constante, on parle de classifieur aléatoire… La diagonale bleue sur la figure ci-dessus est justement un classifieur aléatoire… C’est ce qu’on obtient si on prédit au hasard
pred = prediction(S,Y) plot(performance(pred,"tpr","fpr")) |
Le point est obtenu avec un seuil de 50\% (ou en fait, n’importe quelle valeur entre 40\% et 80\%). On peut là encore calculer l’aire sous la courbe, cette fois à l’aide de trapèzes (ou de triangles)
p1 = roc.curve(.5) p2[1]*p2[2]/2+(1-p1[1])*p1[2]+(1-p1[1])*(1-p1[2])/2 [1] 0.7083333 auc.perf = performance(pred, measure = "auc") auc.perf@y.values[[1]] [1] 0.7083333 |
