{"word_count": 769, "character_count": 6964}
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 [latex]n=10[/latex] observations, deux variables continues, [latex]x_1[/latex] et [latex]x_2[/latex], et la variable d'intérêt binaire [latex]y\in\{0,1\}[/latex]. On peut représenter les points dans le plan [latex](x_1,x_2)[/latex], et on utilise une couleur différente pour [latex]y\in\{0,1\}[/latex].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 [latex display="true"]\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}}[/latex]On peut visualiser l'ensemble des points [latex](x_1,x_2)[/latex] pour lesquels [latex display="true"]\mathbb{P}(Y=0|x_1,x_2)=\mathbb{P}(Y=1|x_1,x_2)[/latex](on a alors autant de chance d'être rouge et bleu) soit [latex]\beta_0+\beta_1x_1+\beta_2x_2=0[/latex] 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 [latex]y_i[/latex] en fonction du score, i.e. l'estimation de [latex]\mathbb{P}(Y=1|x_{1,i},x_{2,i})[/latex],Y = df$yS = predict(reg,type="response")plot(S,y,xlab="probabilité prédite",ylab="y")On va alors se donner un seuil (par exemple [latex]50\%[/latex]) : si la probabilité que [latex]Y[/latex] prenne la valeur [latex]1[/latex] excède le seuil, on prédit [latex]1[/latex] (et sinon [latex]0[/latex]). Sur la figure ci-dessus, on a alors 4 sortes de points : ceux à gauche du seuil (et qui sont prédits [latex]0[/latex]), qui sont bien classés s'ils sont en bas, et mal classés en haut; à droite du seuil (et qui sont prédits [latex]1[/latex]), ils sont bien classés s'ils sont en haut, et mal classés en bas [dans le code ci-dessous, le symbole &#38;gt désigne l'opérateur "supérieur", qui malheureusement ne passe pas dans cet éditeur]seuil = .5Yhat = (S&#62;seuil)*1plot(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èsetable(Yhat,Y)    YYhat 0 1   0 3 1   1 1 5Ce 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 [latex]50\%[/latex]) 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 [latex]y[/latex]) en fonction du taux de faux positifs (sur l'axe [latex]x[/latex]L'ensemble des points donne la courbe ROC.roc.curve=function(s,print=FALSE){  Ps=(S&#62;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 [latex]50\%[/latex] est bien sur la courbetable(Yhat,Y)    YYhat 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.8333333abline(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.875mais qu'on peut avoir automatiquementauc.perf = performance(pred, measure = "auc")auc.perf@y.values[[1]][1] 0.875Allez, tentons un autre classifieur : toujours une régression logistique, mais sur un facteur obtenu en coupant la seconde variable en deux, [latex]\boldsymbol{1}_{[s,\infty)}(x_2)[/latex]reg = glm(y~I(x2&#62;.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 [latex]x_2[/latex]. Ici, on prédit juste deux valeurs : [latex]40\%[/latex] de chance d'être bleu en bas, et [latex]80\%[/latex] de chance d'être bleu, en haut. Si on représente observations [latex]y_i[/latex] en fonction des probabilités prédites, on obtientY = df$yS = predict(reg,type="response")plot(S,y,xlab="probabilité prédite",ylab="y",xlim=0:1)Avec un seuil à [latex]50\%[/latex], on obtient le tableau de contingence suivant (avec 3 erreurs, contre 2 auparavant)seuil = .5Yhat = (S&#62;seuil)*1table(Yhat,Y)    YYhat 0 1   0 3 2   1 1 4Si on trace la courbe ROC, on obtientroc.curve=function(s,print=FALSE){  Ps=(S&#62;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 [latex]40\%[/latex] (en bas) : quand on doit prévoir pour un individu dans cette région, on lui attribue les valeurs [latex]\{0,1\}[/latex] respectivement avec les probabilités [latex]\{40\%,60\%\}[/latex]. 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 hasardpred = prediction(S,Y)plot(performance(pred,"tpr","fpr"))Le point est obtenu avec un seuil de [latex]50\%[/latex] (ou en fait, n'importe quelle valeur entre [latex]40\%[/latex] et [latex]80\%[/latex]). 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.7083333auc.perf = performance(pred, measure = "auc")auc.perf@y.values[[1]][1] 0.7083333