Freakonometrics

Juste avant les vacances, Jean-Pierre Liégeois, un jeune lecteur du var, me demandais par courriel, "A partir d'une régression logistique(ou d'une matrice de confusion 2x2), 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 . 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"

	valeur observée
valeur prédite		"positif"	"négatif"
	"positif"	TP	FP
	"négatif"	FN	TN

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

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

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

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

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

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