Tag Archives: ROC

What it the interpretation of the diagonal for a ROC curve

Last Friday, we discussed the use of ROC curves to describe the goodness of a classifier. I did say that I will post a brief paragraph on the interpretation of the diagonal. If you look around some say that it describes the “strategy of randomly guessing a class“, that it is obtained with “a diagnostic test that is no better than chance level“, even obtained by “making a prediction by tossing of an unbiased coin“.

Let us get back to ROC curves to illustrate those points. Consider a very simple dataset with 10 observations (that is not linearly separable)

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))

here we can check that, indeed, it is not separable


Consider a logistic regression (the course is on linear models)

reg = glm(y~x1+x2,data=df,family=binomial(link = "logit"))

but any model here can be used… We can use our own function


or any R package actually



We can plot the two simultaneously here

V=Vectorize(roc.curve)(seq(-5,5,length=251))points(V[1,],V[2,])segments(0,0,1,1,col="light blue")

So our code works just fine, here. Let us consider various strategies that should lead us to the diagonal.

The first one is : everyone has the same probability (say 50%)



Indeed, we have the diagonal. But to be honest, we have only two points here : (0,0) and (1,1). Claiming that we have a straight line is not very satisfying… Actually, note that we have this situation whatever the probability we choose



We can try another strategy, like “making a prediction by tossing of an unbiased coin“. This is what we obtain



V=Vectorize(roc.curve)(seq(0,1,length=251))points(V[1,],V[2,])segments(0,0,1,1,col="light blue")

We can also try some sort of “random classifier”, where we choose the score randomly, say uniform on the unit interval



V=Vectorize(roc.curve)(seq(0,1,length=251))points(V[1,],V[2,])segments(0,0,1,1,col="light blue")

Let us try to go further on that one. For convenience, let us consider another function to plot the ROC curve


roc_curve=Vectorize(function(x) max(V[2,which(V[1,]<=x)]))

We have the same line as previously



But now, consider many scoring strategies, all randomly chosen

MY=matrix(NA,500,length(y))for(i in 1:500){
MY[i,]=roc_curve(x)}plot(performance(prediction(S,df$y),"tpr","fpr"),col="white")for(i in 1:500){lines(x,MY[i,],col=rgb(0,0,1,.3),type="s")}lines(c(0,x),c(0,apply(MY,2,mean)),col="red",type="s",lwd=3)segments(0,0,1,1,col="light blue")

The red line is the average of all random classifiers. It is not a straight line, be we observe oscillations around the diagonal.

Consider a dataset with more observations

myocarde = read.table("http://freakonometrics.free.fr/myocarde.csv",head=TRUE, sep=";")

myocarde$PRONO = (myocarde$PRONO=="SURVIE")*1

reg = glm(PRONO~.,data=myocarde,family=binomial(link = "logit"))



V=Vectorize(roc.curve)(seq(-5,5,length=251))points(V[1,],V[2,])segments(0,0,1,1,col="light blue")

Here is a “random classifier” where we draw scores randomly on the unit interval


V=Vectorize(roc.curve)(seq(-5,5,length=251))points(V[1,],V[2,])segments(0,0,1,1,col="light blue")

And if we do that 500 times, we obtain, on average

MY=matrix(NA,500,length(y))for(i in 1:500){
MY[i,]=roc_curve(x)}plot(performance(prediction(S,Y),"tpr","fpr"),col="white")for(i in 1:500){lines(x,MY[i,],col=rgb(0,0,1,.3),type="s")}lines(c(0,x),c(0,apply(MY,2,mean)),col="red",type="s",lwd=3)segments(0,0,1,1,col="light blue")

So, it looks like me might say that the diagonal is what we have, on average, when drawing randomly scores on the unit interval…

I did mention that an interesting visual tool could be related to the use of the Kolmogorov Smirnov statistic on classifiers. We can plot the two empirical cumulative distribution functions of the scores, given the response Y




we can also look at the distribution of the score, with the histogram (or density estimates)



The underlying idea is the following : we do have a “perfect classifier” (top left corner)

is the supports of the scores do not overlap

otherwise, we should have errors. That the case below

we in 10% of the cases, we might have misclassification

or even more missclassification, with overlapping supports

Now, we have the diagonal

when the two conditional distributions of the scores are identical

Of course, that only valid when n is very large, otherwise, it is only what we observe on average….

On the poor performance of classifiers in insurance models

Each time we have a case study in my actuarial courses (with real data), students are surprised to have hard time getting a “good” model, and they are always surprised to have a low AUC, when trying to model the probability to claim a loss, to die, to fraud, etc. And each time, I keep saying, “yes, I know, and that’s what we expect because there a lot of ‘randomness’ in insurance”. To be more specific, I decided to run some simulations, and to compute AUCs to see what’s going on. And because I don’t want to waste time fitting models, we will assume that we have each time a perfect model. So I want to show that the upper bound of the AUC is actually quite low ! So it’s not a modeling issue, it is a fondamental issue in insurance !

By ‘perfect model’ I mean the following : \Omega denotes the heterogeneity factor, because people are different. We would love to get \mathbb{P}[Y=1|\Omega]. Unfortunately, \Omega  is unobservable ! So we use covariates (like the age of the driver of the car in motor insurance, or of the policyholder in life insurance, etc). Thus, we have data (y_i,\boldsymbol{x}_i)‘s and we use them to train a model, in order to approximate \mathbb{P}[Y=1|\boldsymbol{X}]. And then, we check if our model is good (or not) using the ROC curve, obtained from confusion matrices, comparing y_i‘s and \widehat{y}_i‘s where \widehat{y}_i=1 when \mathbb{P}[Y_i=1|\boldsymbol{x}_i] exceeds a given threshold. Here, I will not try to construct models. I will predict \widehat{y}_i=1 each time the true underlying probability \mathbb{P}[Y_i=1|\omega_i] exceeds a threshold ! The point is that it’s possible to claim a loss (y=1) even if the probability is 3% (and most of the time \widehat{y}=0), and to not claim one (y=0) even if the probability is 97% (and most of the time \widehat{y}=1). That’s the idea with randomness, right ?

So, here p(\omega_1),\cdots,p(\omega_n) denote the probabilities to claim a loss, to die, to fraud, etc. There is heterogeneity here, and this heterogenity can be small, or large. Consider the graph below, to illustrate,

In both cases, there is, on average, 25% chance to claim a loss. But on the left, there is more heterogeneity, more dispersion. To illustrate, I used the arrow, which is a classical 90% interval : 90% of the individuals have a probability to claim a loss in that interval. (here 10%-40%), 5% are below 10% (low risk), and 5% are above 40% (high risk). Later on, we will say that we have 25% on average, with a dispersion of 30% (40% minus 10%). On the right, it’s more 25% on average, with a dispersion of of 15%. What I call dispersion is the difference between the 95% and the 5% quantiles.

Consider now some dataset, with Bernoulli variables y, drawn with those probabilities p(\omega). Then, let us assume that we are able to get a perfect model : I do not estimate a model based on some covariates, here, I assume that I know perfectly the probability (which is true, because I did generate those data). More specifically, to generate a vector of probabilities, here I use a Beta distribution with a given mean, and a given variance (to capture the heterogeneity I mentioned above)


from those probabilities, I generate occurences of claims, or deaths,

Y=rbinom(n,size = 1,prob = p)

Then, I compute the AUC of my “perfect” model,


And then, I will generate many samples, to compute the average value of the AUC. And actually, we can do that for many values of the mean and the variance of the Beta distribution. Here is the code

ab_beta = function(m,inter){
  a=uniroot(function(a) qbeta(.95,a,a/m-a)-qbeta(.05,a,a/m-a)-inter,
  essai = try(ab<-ab_beta(m,i),TRUE) if(inherits(essai,what="try-error")) a=-1 if(!inherits(essai,what="try-error")){ a=ab[1] b=ab[2] } if((a>=0)&(b>=0)){
    for(s in 1:ns){
      Y=rbinom(n,size = 1,prob = p)
V=outer(X = Vm,Y = Vi, Vectorize(function(x,y) 
      xlab="Probability (Average)",
      ylab="Dispersion (Q95-Q5)",
        colorRampPalette(brewer.pal(n = 9, name = "YlGn"))(101))

On the x-axis, we have the average probability to claim a loss. Of course, there is a symmetry here. And on the y-axis, we have the dispersion : the lower, the less heterogeneity in the portfolio. For instance, with a 30% chance to claim a loss on average, and 20% dispersion (meaning that in the portfolio, 90% of the insured have between 20% and 40% chance to claim a loss, or 15% and 35% chance), we have on average a 60% AUC. With a perfect model ! So with only a few covariates, having 55% should be great !

My point here is that with a low dispersion, we cannot expect to have a great AUC (again, even with a perfect model). In motor insurance, from my experience, 90% of the insured are between 3% chance and 20% chance to claim a loss ! That’s less than 20% dispersion ! and in that case, even if the (average) probability is rather small, it is very difficult to expect an AUC above 60% or 65% !

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)

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  

             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
“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)
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 http://freakonometrics.hypotheses.org/files/2018/02/ROC-01.png et (comme je l’expliquais dans un autre billet) et d’un score \widehat{S}. Et c’est ce score qu’on va utiliser pour construire la courbe ROC. Ce score sera utilise pour prédire http://freakonometrics.hypotheses.org/files/2018/02/ROC-02.png . La règle d’affectation est alors simple: on se fixe un seuil http://freakonometrics.hypotheses.org/files/2018/02/ROC-04.png, et

  • si http://freakonometrics.hypotheses.org/files/2018/02/ROC-05.png, alors http://freakonometrics.hypotheses.org/files/2018/02/ROC-02.png est “positif”
  • si http://freakonometrics.hypotheses.org/files/2018/02/ROC-06.png, alors http://freakonometrics.hypotheses.org/files/2018/02/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 http://freakonometrics.hypotheses.org/files/2018/02/ROC-01.png
valeur prédite
“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.


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)){

On a alors le graphique suivant,

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

plot(0:1,0:1,xlab="False Positive Rate",
ylab="True Positive Rate",cex=.5)
for(s in seq(0,1,by=.01)){

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

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,

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)

PROC=plot.roc(Y,P,main="", percent=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
“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é


et la fonction de spécificité


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


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.