Tag Archives: ROC

Ce que la courbe ROC (et l’AUC) ne raconte pas

En préparant une intervention pour mardi prochain, j’épluchais les résultats renvoyés pour un exercice, et j’ai eu un résultat assez étrange avec un modèle de classification. J’avais donné la même base cet automne à l’ensae, et j’avais donc près d’une trentaine d’autres modèles, pour comparer (disons plutôt que sur la même base de test, j’ai une trentaine de prévisions). Les observations noires sont celles obtenues cet automne (le trait correspond aux meilleurs AUC sur la base de test), et les observations rouges sont celles que j’ai obtenu pour l’intervention de mardi (là encore, le trait vertical correspond aux meilleurs modèles, au sens du AUC), sur une observation de la base de test,

Ce sont les probabilités prédites (mais j’ai enlevé l’échelle).

Pour presque toutes mes observations, les poids rouges sont bien au dessus des autres… Mais cela n’enlève en rien le fait que l’AUC obtenu (pour les deux modèles rouges) est très bon. C’est effectivement un résultat important (et connu) : le critère AUC (mais plus généralement la courbe ROC en entier) n’indique en rien si la valeur prédite est bonne, ou pas. Il nous dit juste si l’ordre obtenu est correct. Si les valeurs les plus importantes sont effectivement les valeurs pour lesquelles on a un 1, l’AUC sera très bon.

C’est ce qu’on peut observer sur le petit exemple ci-dessous. Considérons un modèle logistique simulé assez simple,

> n=1e3
> set.seed(1)
> x1=rnorm(n)
> x2=runif(n)
> u=-3+x2+x1
> p=exp(u)/(1+exp(u))
> y=rbinom(n,prob=p,size=1)
> library(ROCR)
> df=data.frame(y,x1,x2)
> mean(df$y)
[1] 0.116
> reg=glm(y~.,data=df,family=binomial)
> p=predict(reg,type="response")
> mean(p)
[1] 0.116
> pred1=prediction(p, df$y)
> L=performance(pred1, "tpr", "fpr")

L’AUC est ici

> auc=performance(pred1, "auc")@y.values[[1]]
> auc
[1] 0.7681191

et la courbe ROC est la suivante,

> plot(unlist(L@x.values),unlist(L@y.values),
+ type="s",col="blue")

Supposons que l’on change la constante du modèle logistique,

> reg$coefficients[1]=0

Dans ce cas, notre prévision est assez mauvaise, car la probabilité moyenne prédite est ici

> u=reg$coefficients[1]+reg$coefficients[2]*
+ df$x1+reg$coefficients[3]*df$x2
> p=exp(u)/(1+exp(u))
> mean(p)
[1] 0.6060676

(on est loin des 11,6% de 1 dans la base). Pourtant le AUC est bon

> pred1=prediction(p, df$y)
> L=performance(pred1, "tpr", "fpr")
> auc=performance(pred1, "auc")@y.values[[1]]
> auc
[1] 0.7681191

(c’est en fait la même valeur qu’auparavant, ce qui a du sens puisque la courbe ROC est identique)

> lines(unlist(L@x.values),unlist(L@y.values),
+ type="s",col="red")

Autrement dit, ces outils, classiquement utilisés pour juger la qualité d’un classifieur, ne permettent en aucun cas de dire que la probabilité prédite a du sens. Ces critères permettent juste de dire qu’on identifie assez bien les personnes qui ont le plus de chance d’avoir la réponse 1. Ce qui n’est pas si mal… mais c’est un autre problème que celui d’avoir une probabilité qui soit pertinente.

Classification on the German Credit Database

In our data science course, this morning, we’ve use random forrest to improve prediction on the German Credit Dataset. The dataset is

> url="http://freakonometrics.free.fr/german_credit.csv"
> credit=read.csv(url, header = TRUE, sep = ",")

Almost all variables are treated a numeric, but actually, most of them are factors,

> str(credit)
'data.frame':	1000 obs. of  21 variables:
 $ Creditability   : int  1 1 1 1 1 1 1 1 1 1 ...
 $ Account.Balance : int  1 1 2 1 1 1 1 1 4 2 ...
 $ Duration        : int  18 9 12 12 12 10 8  ...
 $ Purpose         : int  2 0 9 0 0 0 0 0 3 3 ...

(etc). Let us convert categorical variables as factors,

> F=c(1,2,4,5,7,8,9,10,11,12,13,15,16,17,18,19,20)
> for(i in F) credit[,i]=as.factor(credit[,i])

Let us now create our training/calibration and validation/testing datasets, with proportion 1/3-2/3

> i_test=sample(1:nrow(credit),size=333)
> i_calibration=(1:nrow(credit))[-i_test]

The first model we can fit is a logistic regression, on selected covariates

> LogisticModel <- glm(Creditability ~ Account.Balance + Payment.Status.of.Previous.Credit + Purpose + 
Length.of.current.employment + 
Sex...Marital.Status, family=binomial, 
data = credit[i_calibration,])

Based on that model, it is possible to draw the ROC curve, and to compute the AUC (on ne validation dataset)

> fitLog <- predict(LogisticModel,type="response",
+                   newdata=credit[i_test,])
> library(ROCR)
> pred = prediction( fitLog, credit$Creditability[i_test])
> perf <- performance(pred, "tpr", "fpr")
> plot(perf)
> AUCLog1=performance(pred, measure = "auc")@y.values[[1]]
> cat("AUC: ",AUCLog1,"\n")
AUC:  0.7340997

An alternative is to consider a logistic regression on all explanatory variables

> LogisticModel <- glm(Creditability ~ ., 
+  family=binomial, 
+  data = credit[i_calibration,])

We might overfit, here, and we should observe that on the ROC curve

> fitLog <- predict(LogisticModel,type="response",
+                   newdata=credit[i_test,])
> pred = prediction( fitLog, credit$Creditability[i_test])
> perf <- performance(pred, "tpr", "fpr")
> plot(perf)
> AUCLog2=performance(pred, measure = "auc")@y.values[[1]]
> cat("AUC: ",AUCLog2,"\n")
AUC:  0.7609792

There is a slight improvement here,  compared with the previous model, where only five explanatory variables were considered.

Consider now some regression tree (on all covariates)

> library(rpart)
> ArbreModel <- rpart(Creditability ~ ., 
+  data = credit[i_calibration,])

We can visualize the tree using

> library(rpart.plot)
> prp(ArbreModel,type=2,extra=1)

The ROC curve for that model is

> fitArbre <- predict(ArbreModel,
+                     newdata=credit[i_test,],
+                     type="prob")[,2]
> pred = prediction( fitArbre, credit$Creditability[i_test])
> perf <- performance(pred, "tpr", "fpr")
> plot(perf)
> AUCArbre=performance(pred, measure = "auc")@y.values[[1]]
> cat("AUC: ",AUCArbre,"\n")
AUC:  0.7100323

As expected, a single has a lower performance, compared with a logistic regression. And a natural idea is to grow several trees using some boostrap procedure, and then to agregate those predictions.

> library(randomForest)
> RF <- randomForest(Creditability ~ .,
+ data = credit[i_calibration,])
> fitForet <- predict(RF,
+                     newdata=credit[i_test,],
+                     type="prob")[,2]
> pred = prediction( fitForet, credit$Creditability[i_test])
> perf <- performance(pred, "tpr", "fpr")
> plot(perf)
> AUCRF=performance(pred, measure = "auc")@y.values[[1]]
> cat("AUC: ",AUCRF,"\n")
AUC:  0.7682367

Here this model is (slightly) better than the logistic regression. Actually, if we create many training/validation samples, and compare the AUC, we can observe that – on average – random forests perform better than logistic regressions,

> AUC=function(i){
+   set.seed(i)
+   i_test=sample(1:nrow(credit),size=333)
+   i_calibration=(1:nrow(credit))[-i_test]
+   LogisticModel <- glm(Creditability ~ ., 
+    family=binomial, 
+    data = credit[i_calibration,])
+   summary(LogisticModel)
+   fitLog <- predict(LogisticModel,type="response",
+                     newdata=credit[i_test,])
+   library(ROCR)
+   pred = prediction( fitLog, credit$Creditability[i_test])
+   AUCLog2=performance(pred, measure = "auc")@y.values[[1]] 
+   RF <- randomForest(Creditability ~ .,
+   data = credit[i_calibration,])
+   fitForet <- predict(RF,
+                       newdata=credit[i_test,],
+                       type="prob")[,2]
+   pred = prediction( fitForet, credit$Creditability[i_test])
+   AUCRF=performance(pred, measure = "auc")@y.values[[1]]
+   return(c(AUCLog2,AUCRF))
+ }
> A=Vectorize(AUC)(1:200)
> plot(t(A))

Choosing a Classifier

In order to illustrate the problem of chosing a classification model consider some simulated data,

> n = 500
> set.seed(1)
> X = rnorm(n)
> ma = 10-(X+1.5)^2*2
> mb = -10+(X-1.5)^2*2
> M = cbind(ma,mb)
> set.seed(1)
> Z = sample(1:2,size=n,replace=TRUE)
> Y = ma*(Z==1)+mb*(Z==2)+rnorm(n)*5
> df = data.frame(Z=as.factor(Z),X,Y)

A first strategy is to split the dataset in two parts, a training dataset, and a testing dataset.

> df1 = training = df[1:300,]
> df2 = testing  = df[301:500,]
  • The Holdout Method: Training and Testing Datasets

The two datasets can be visualised below, with the training dataset on top, and the testing dataset below

> plot(df1$X,df1$Y,pch=19,col=c(rgb(1,0,0,.4),
+ rgb(0,0,1,.4))[df1$Z])

Continue reading Choosing a Classifier

Please, never use my codes without checking twice (at least)!

I wanted to get back on some interesting experience, following a discussion I had with Carlos after my class, this morning. Let me simplify the problem, and change also the dataset. Consider the following dataset

> db = read.table("http://freakonometrics.free.fr/db2.txt",header=TRUE,sep=";")

Let me change also one little thing (in the course, we use the age of people as explanatory variables, so let us consider rounded figures to),

> db$X1=round(db$X1*10)
> db$X2=round(db$X2*10)

Assume that you want to work with factors, because you don’t see why there should be some linear model (and you did not look at the awesome posts on smoothing techniques)

> db$X1F=cut(db$X1,c(-12,45,75,120))
> db$X2F=cut(db$X2,c(100,200,300))

Then you run your regression,

> reg = glm(Y~X1F+X2F+X3,family=binomial,data=db)

So far, nothing wrong, you can try, no error, no warning. Then Carlos wanted to use a ROC curve to see how the model was performing… so he did use some code I uploaded on the blog, something like

> reg = glm(Y~X1+X2+X3,family=binomial,data=db)
> S = predict(reg,type="response")
> Y = db$Y
> plot(0:1,0:1,xlab="False Positive Rate",ylab="True Positive Rate",cex=.5)
> for(s in seq(0,1,by=.01)){
+   Ps=(S>s)*1
+   FP=sum((Ps==1)*(nombre==0))/sum(nombre==0)
+   TP=sum((Ps==1)*(nombre==1))/sum(nombre==1)
+   points(FP,TP,cex=.5,col="red")
+ }

To make it nicer, let us use the following code

> ROCcurve=function(s){
+ Ps=(S>s)*1
+ FP=sum((Ps==1)*(Y==0))/sum(Y==0)
+ TP=sum((Ps==1)*(Y==1))/sum(Y==1)
+ return(c(FP,TP))}
> u=seq(0,1,by=.001)
> vectROC=Vectorize(ROCcurve)(u)

Here, I should mention that I got a warning, but to be honest, when you see the graph, you’re so puzzled that you forget about it…

There were 50 or more warnings (use warnings() to see the first 50)

Carlos ran the code, show me the graph, and asked me “can my model be that bad?”,

 My first answer is “no, your model cannot be that bad! you use (almost) all your explanatory variables, you cannot have such a bad model…“. So, where the problem comes from? My first guess is that this is what you get using some random sort of a variable. For instance, if you use the code above on

> reg2 = glm(Y~X1+X2+X3,family=binomial,data=db)
> S = sort(predict(reg2,type="response"))
> Y = db$Y

Here, one variable is sorted, and we get

This confirms the idea that using a random classifier, the ROC curve is on the diagonal. But here, when you look at the code, you do not see any sorting operation.

So, again, what went wrong? The problem is actually very simple. The dataset looks like

> head(db)
  Y X1  X2 X3      X1F       X2F
1 1 33 163  B (-12,45] (100,200]
2 1 64 185  D  (45,75] (100,200]
3 1 53 166  B  (45,75] (100,200]
4 1 55 197  C  (45,75] (100,200]
5 1 41 184  C (-12,45] (100,200]
6 1 78 196  C (75,120] (100,200]

If we look at the factor variable, the lowest part is not included, in the interval. More precisely,

> levels(db$X1F)
[1] "(-12,45]" "(45,75]"  "(75,120]"

and if we look more precisely and the range of the two variables, we get

> range(db$X1)
[1] -12 120

Wait… the minimum is -12, but it does not appear in the factor (and we do not get any warning)? Yes,

> db[which.min(db$X1)+(-2):2,]
    Y  X1  X2 X3      X1F       X2F
429 1  76 227  A (75,120] (200,300]
430 0  35 186  D (-12,45] (100,200]
431 0 -12 109  E     <NA> (100,200]
432 1  76 225  B (75,120] (200,300]
433 1  61 206  A  (45,75] (200,300]

Now we start so see more clearly what went wrong…  there is a missing value in the dataset. And when we get our prediction, there is no missing value: the missing value has be droped.

> length(predict(reg))
[1] 999
> nrow(db)
[1] 1000

So when we compare the prediction and the observed value…. there is a problem, since the vectors don’t match. Actually, it was mentioned in the output of the regression (but we did not look at it)

> summary(reg)

Call:
glm(formula = Y ~ X1F + X2F + X3, family = binomial, data = db)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-2.9432   0.1627   0.2900   0.4823   1.1377  

Coefficients:
             Estimate Std. Error z value Pr(>|z|)    
(Intercept)   0.22696    0.23206   0.978 0.328067    
X1F(45,75]    1.86569    0.23397   7.974 1.53e-15 ***
X1F(75,120]   2.97463    0.65071   4.571 4.85e-06 ***
X2F(200,300]  1.11643    0.32695   3.415 0.000639 ***
X3B          -0.06131    0.31076  -0.197 0.843609    
X3C           0.75013    0.35825   2.094 0.036268 *  
X3D           0.13846    0.31399   0.441 0.659226    
X3E          -0.13277    0.31859  -0.417 0.676853    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 802.34  on 998  degrees of freedom
Residual deviance: 607.61  on 991  degrees of freedom
  (1 observation deleted due to missingness)
AIC: 623.61

Number of Fisher Scoring iterations: 7

Yes, there is a tiny little sentence,

  (1 observation deleted due to missingness)

So that was it? Yes! If you drop that missing value, you get something more realistic,

> S= predict(reg,type="response")
> Y=db$Y[-which.min(db$X1)]

or (probably better), change the left point of the interval

> db$X1F=cut(db$X1,c(-13,45,75,120))
> reg = glm(Y~X1F+X2F+X3,family=binomial,data=db)
> S= predict(reg,type="response")
> Y= db$Y

(the output would have been almost the same).

Observe that if we had used a dedicated package, we would not have encountered this problem. For instance (starting with the initial values of the vectors)

> library(ROCR)
> prediction(S,Y)
Error in prediction(S, Y) : 
  Number of predictions in each run must be equal to the number of labels for each run.

We do have the answer here: both vectors do not have the same length…

So, what is my point? R is great, because of all the packages, but as a teacher, I do not feel comfortable asking my student to use those functions, as black boxes. So I try to write my own codes, to get the same output. So yes, I do write codes to explain what the black box is doing, to simplify the algorithm, and show what’s going on. When working on a (forthcoming) book as the Editor, we had a discussion with Rob Hyndman about that issue. I wanted the contributors to explain the core of the code, with a simplified algorithm, when using a dedicated package. I do truly believe that using simplified codes might help to understand better. Until you start to have problems. Because I write a code to deal with one specific problem, there is no check for possible errors. And once the code is understood, please, please do not use it! use R function that can handle errors…

ROC curves and classification

To get back to a question asked after the last course (still on non-life insurance), I will spend some time to discuss ROC curve construction, and interpretation. Consider the dataset we’ve been using last week,

> db = read.table("http://freakonometrics.free.fr/db.txt",header=TRUE,sep=";")
> attach(db)

The first step is to get a model. For instance, a logistic regression, where some factors were merged together,

> X3bis=rep(NA,length(X3))
> X3bis[X3%in%c("A","C","D")]="ACD"
> X3bis[X3%in%c("B","E")]="BE"
> db$X3bis=as.factor(X3bis)
> reg=glm(Y~X1+X2+X3bis,family=binomial,data=db)

From this model, we can predict a probability, not a  variable,

> S=predict(reg,type="response")

Let https://latex.codecogs.com/gif.latex?\widehat{S} denote this variable (actually, we can use the score, or the predicted probability, it will not change the construction of our ROC curve). What if we really want to predict a  variable. As we usually do in decision theory. The idea is to consider a threshold https://perso.univ-rennes1.fr/arthur.charpentier/latex/ROC-04.png, so that

  • if https://perso.univ-rennes1.fr/arthur.charpentier/latex/ROC-05.png, then  https://perso.univ-rennes1.fr/arthur.charpentier/latex/ROC-02.png will be https://latex.codecogs.com/gif.latex?1, or “positive” (using a standard terminology)
  • si https://perso.univ-rennes1.fr/arthur.charpentier/latex/ROC-06.png, then  https://perso.univ-rennes1.fr/arthur.charpentier/latex/ROC-02.png will be https://latex.codecogs.com/gif.latex?0, or “negative

Then we derive a contingency table, or a confusion matrix

     observed value https://perso.univ-rennes1.fr/arthur.charpentier/latex/ROC-01.png
predicted
value
https://perso.univ-rennes1.fr/arthur.charpentier/latex/ROC-02.png
“positive“ “négative“
“positive“ TP FP
“négative“ FN TN

where TP are the so-called true positive, TN  the true negative, FP are the false positive (or type I error) and FN are the false negative (type II errors). We can get that contingency table for a given threshold https://perso.univ-rennes1.fr/arthur.charpentier/latex/ROC-04.png

> 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)
+ }
> threshold = 0.5
> roc.curve(threshold,print=TRUE)
        Predicted
Observed   0   1
       0   5 231
       1  19 745
      FPR       TPR 
0.9788136 0.9751309

Here, we also compute the false positive rates, and the true positive rates,

  • TPR = TP / P = TP / (TP + FN) also called sentivity, defined as the rate of true positive: probability to be predicted positve, given that someone is positive (true positive rate)
  • FPR = FP / N = FP / (FP + TN) is the rate of false positive: probability to be predicted positve, given that someone is negative (false positive rate)

The ROC curve is then obtained using severall values for the threshold. For convenience, define

> ROC.curve=Vectorize(roc.curve)

First, we can plot https://latex.codecogs.com/gif.latex?(\widehat{S}_i,Y_i) (a standard predicted versus observed graph), and visualize true and false positive and negative, using simple colors

> I=(((S>threshold)&(Y==0))|((S<=threshold)&(Y==1)))
> plot(S,Y,col=c("red","blue")[I+1],pch=19,cex=.7,,xlab="",ylab="")
> abline(v=threshold,col="gray")

And for the ROC curve, simply use

> M.ROC=ROC.curve(seq(0,1,by=.01))
> plot(M.ROC[1,],M.ROC[2,],col="grey",lwd=2,type="l")

This is the ROC curve. Now, to see why it can be interesting, we need a second model. Consider for instance a classification tree

> library(tree)
> ctr <- tree(Y~X1+X2+X3bis,data=db)
> plot(ctr)
> text(ctr)

To plot the ROC curve, we just need to use the prediction obtained using this second model,

> S=predict(ctr)

All the code described above can be used. Again, we can plot https://latex.codecogs.com/gif.latex?(\widehat{S}_i,Y_i) (observe that we have 5 possible values for https://latex.codecogs.com/gif.latex?\widehat{S}_i, which makes sense since we do have 5 leaves on our tree). Then, we can plot the ROC curve,

An interesting idea can be to plot the two ROC curves on the same graph, in order to compare the two models

> plot(M.ROC[1,],M.ROC[2,],type="l")
> lines(M.ROC.tree[1,],M.ROC.tree[2,],type="l",col="grey",lwd=2)

The most difficult part is to get a proper interpretation. The tree is not predicting well in the lower part of the curve. This concerns people with a very high predicted probability. If our interest is more on those with a probability lower than 90%, then, we have to admit that the tree is doing a good job, since the ROC curve is always higher, comparer with the logistic regression.

Construire une courbe ROC

Juste avant les vacances, Jean-Pierre Liégeois, un jeune lecteur du var, me demandais par courriel, “à partir d’une régression logistique (ou d’une matrice de confusion 2×2), comment programmer en R, un programme qui construit la courbe ROC associée“. Avant d’aller plus loin (et de répondre a la question), je vais renvoyer vers un vieux billet sur les matrices de confusion. L’idée est que l’on suppose que l’on dispose d’un prédicteur d’une variable prenant des valeurs 0 et 1 (ou pour reprendre la terminologie classique “positif” et “négatif”), par exemple un modèle logistique. Formellement, pour l’ensemble de nos observations, on a une valeur observée https://perso.univ-rennes1.fr/arthur.charpentier/latex/ROC-01.png et (comme je l’expliquais dans un autre billet) et d’un score https://perso.univ-rennes1.fr/arthur.charpentier/latex/ROC-03.png. Et c’est ce score qu’on va utiliser pour construire la courbe ROC. Ce score sera utilise pour prédire https://perso.univ-rennes1.fr/arthur.charpentier/latex/ROC-02.png . La règle d’affectation est alors simple: on se fixe un seuil https://perso.univ-rennes1.fr/arthur.charpentier/latex/ROC-04.png, et

  • si https://perso.univ-rennes1.fr/arthur.charpentier/latex/ROC-05.png, alors  https://perso.univ-rennes1.fr/arthur.charpentier/latex/ROC-02.png est “positif”
  • si https://perso.univ-rennes1.fr/arthur.charpentier/latex/ROC-06.png, alors  https://perso.univ-rennes1.fr/arthur.charpentier/latex/ROC-02.png est “négatif”

On peut alors construire une matrice dite de confusion, qui est simplement un table de contingence,

                valeur observée https://perso.univ-rennes1.fr/arthur.charpentier/latex/ROC-01.png
valeur prédite
 https://perso.univ-rennes1.fr/arthur.charpentier/latex/ROC-02.png
“positif” “négatif”
“positif” TP FP
“négatif” FN TN

où TP désigne les vrais positifs (true positive), TN  les vrais négatifs (true negative),FP désigne les faux positifs (false positive) ou erreur de type I (dans une terminologie de théorie de la décision, ou de théorie des tests), et FN désigne les faux négatifs (false negative) ou erreur de type II.
Quid de la mise en œuvre sous R ? Commençons par générer des données, et estimons un modèle de régression.

set.seed(1)
n=50
X=rnorm(n)
Y=rbinom(n,size=1,prob=
exp(2*X-1)/(1+exp(2*X-1)))
B=data.frame(Y,X)
reg=glm(Y~X,family=binomial,data=B)
S=predict(reg,type="response")

On a maintenant notre observation (variable prenant les valeurs 0 ou 1) et notre score. On va ensuite pouvoir choisir plusieurs valeurs possible pour le seuil, et visualiser le taux de vrais positifs, en fonction du taux de faux positifs.

plot(0:1,0:1,xlab="False Positive Rate",
ylab="True Positive Rate",cex=.5)
for(s in seq(0,1,by=.01)){
Ps=(S>s)*1
FP=sum((Ps==1)*(Y==0))/sum(Y==0)
TP=sum((Ps==1)*(Y==1))/sum(Y==1)
points(FP,TP,cex=.5,col="red")
}

On a alors le graphique suivant,

Si on relit les points, on a alors la courbe ROC,

FP=TP=rep(NA,101)
plot(0:1,0:1,xlab="False Positive Rate",
ylab="True Positive Rate",cex=.5)
for(s in seq(0,1,by=.01)){
Ps=(S>s)*1
FP[1+s*100]=sum((Ps==1)*(Y==0))/sum(Y==0)
TP[1+s*100]=sum((Ps==1)*(Y==1))/sum(Y==1)
}
lines(c(FP),c(TP),type="s",col="red")

En fait, le code est assez simple, et il traîne dans différents packages de R, e.g.

library(ROCR)
pred=prediction(P,Y)
perf=performance(pred,"tpr", "fpr")
plot(perf,colorize = TRUE)

On peut aussi s’amuser a bootstrapper l’échantillon pour construire des intervalles de confiance, ou ajuster des modèles théoriques,

library(verification)
roc.plot(Y,P, xlab = "False Positive Rate",
ylab = "True Positive Rate", main = "", CI = TRUE,
n.boot = 100, plot = "both", binormal = TRUE)

ou encore (toujours avec des bornes de confiance obtenues par bootstrap)

library(pROC)
PROC=plot.roc(Y,P,main="", percent=TRUE,
ci=TRUE)
SE=ci.se(PROC,specificities=seq(0, 100, 5))
plot(SE, type="shape", col="light blue")

De la qualité d’un score de classification

Un petit mot sur les courbes dites ROC, pour Receiver Operating Characteristic. Pour cela, on suppose que l’on dispose d’un prédicteur d’un variable prenant des valeurs 0 et 1 (pour simplifier), ou mieux encore “positif” et “négatif”. Peu importe le prédicteur, on peut considérer une régression logistique, une analyse discriminante, un classificateur nonparamétrique…. Bref, pour l’ensemble de nos observations, on a une valeur observée https://perso.univ-rennes1.fr/arthur.charpentier/latex/ROC-01.png et une valeur prédite https://perso.univ-rennes1.fr/arthur.charpentier/latex/ROC-02.png. En fait, comme je l’expliquais ici, on dispose plus précisément d’un score https://perso.univ-rennes1.fr/arthur.charpentier/latex/ROC-03.png. La règle d’affectation est alors simple: on se fixe un seuil https://perso.univ-rennes1.fr/arthur.charpentier/latex/ROC-04.png, et

  • si https://perso.univ-rennes1.fr/arthur.charpentier/latex/ROC-05.png, alors  https://perso.univ-rennes1.fr/arthur.charpentier/latex/ROC-02.png est “positif”
  • si https://perso.univ-rennes1.fr/arthur.charpentier/latex/ROC-06.png, alors  https://perso.univ-rennes1.fr/arthur.charpentier/latex/ROC-02.png est “négatif”

On peut alors construire une matrice dite de confusion, qui est simplement un table de contingence,

                valeur observée https://perso.univ-rennes1.fr/arthur.charpentier/latex/ROC-01.png
valeur prédite
 https://perso.univ-rennes1.fr/arthur.charpentier/latex/ROC-02.png
“positif” “négatif”
“positif” TP FP
“négatif” FN TN

où TP désigne les vrais positifs (true positive), TN  les vrais négatifs (true negative),FP désigne les faux positifs (false positive) ou erreur de type I (dans une terminologie de théorie de la décision, ou de théorie des tests), et FN désigne les faux négatifs (false negative) ou erreur de type II.
On peut alors définir toute une batterie d’indicateurs permettant de juger de la qualité de notre prédicteur (ou plutôt de notre score),

  • TPR = TP / P = TP / (TP + FN) appelé sensibilité, correspondant au taux de vrais positifs (true positive rate)
  • FPR = FP / N = FP / (FP + TN) correspondant au taux de faux positifs (false positive rate)
  • ACC = (TP + TN) / (P + N) appelé précision ou accuracy,
  • SPC = TN / N = TN / (FP + TN) = 1 − FPR appelé spécificité ou taux de vrais négatifs (True Negative Rate)
  • PPV = TP / (TP + FP) le taux de positifs prédits (positive predictive value)
  • NPV = TN / (TN + FN) le taux de négatifs prédits (negative predictive value)
  • FDR = FP / (FP + TP) correspondant au false discovery rate

Bref, on convertit cette matrice en probabilités conditionnelles, et beaucoup de notions peuvent être définie en changeant le conditionnement.
On peut aussi essayer de visualiser ces quantités. La représentation graphique la plus connue est probablement la courbe ROC, Receiver Operating Characteristic. L’idée est simple: il s’agit de juger du modèle, indépendmment du seuil s choisi. Ou plutôt de se donner un outils permettant de choisir le seuil s. On définie alors la fonction de sensibilité

https://perso.univ-rennes1.fr/arthur.charpentier/latex/ROC-07.png

et la fonction de spécificité

https://perso.univ-rennes1.fr/arthur.charpentier/latex/ROC-08.png

La courbe ROC est alors la courbe

https://perso.univ-rennes1.fr/arthur.charpentier/latex/ROC-09.pngSi cette courbe coïncide avec la diagonale, c’est que le modèle n’est pas plus performant qu’un modèle aléatoire (où on attribue la classe au hasard). Plus la courbe ROC s’approche du coin supérieur gauche, meilleur est le modèle, car il permet de capturer le plus possible de vrais positifs avec le moins possible de faux positifs. Notons que de part sa construction, la courbe ROC est invariante par toute transformation monotone croissante de la fonction de score (ce qui peut être particulièrement intéressant si on s’amuse à “normaliser” la fonction de score).
De plus, on notera que l’aire sous la courbe ROC doit pouvoir être vu comme une mesure de la qualité de l’ajustement. Le dessin ci-dessous montre la construction de la courbe ROC pour un échantillon discriminé de manière assez simple, par analyse discriminante.

Si on regarde une relecture probabiliste, on obtient la courbe suivante,

Il existe une autre courbe relativement classique, appelée courbe de lift. Cette dernière correspond à la courbe de Lorenz (que j’avais évoquée ici). On pose alors

https://perso.univ-rennes1.fr/arthur.charpentier/latex/ROC-10.png

et la courbe de lift est alors la courbe https://perso.univ-rennes1.fr/arthur.charpentier/latex/ROC-11.png. Là aussi, un calcul d’aire donne un indicateur de la qualité, et on retrouve l’indice de Gini.