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

> 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)  "(-12,45]" "(45,75]" "(75,120]" and if we look more precisely and the range of the two variables, we get > range(db$X1)
 -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))  999 > nrow(db)  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…

# 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 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 . La règle d’affectation est alors simple: on se fixe un seuil , et

• si , alors est “positif”
• si , alors est “négatif”

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

 valeur observée 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.

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