# Talk at CIMAT, Guanajuato, Mexico

I will be back in Guanajuato, Mexico, this week, to visit Victor Rivero. And I will give a talk at the Centro de Investigacion en Matematicas (CIMAT) this Wednesday on “Multivariate Archimax Copulas“. The slides are already online.

(there is a lot of material on copulas, as requested, to provide an introduction for students not familiar with this concept).

# Bivariate Densities with N(0,1) Margins

This Monday, in the ACT8595 course, we came back on elliptical distributions and conditional independence (here is an old post on de Finetti’s theorem, and the extension to Hewitt-Savage’s). I have shown simulations, to illustrate those two concepts of dependent variables, but I wanted to spend some time to visualize densities. More specifically what could be the joint density is we assume that margins are $\mathcal{N}(0,1)$ distributions.

• The Bivariate Gaussian distribution

Here, we consider a Gaussian random vector, with margins $\mathcal{N}(0,1)$, and with correlation $r\in[-1,+1]$. This is the standard graph, with elliptical isodensity curves

r=.5
library(mnormt)
S=matrix(c(1,r,r,1),2,2)
f=function(x,y) dmnorm(cbind(x,y),varcov=S)
vx=seq(-3,3,length=201)
vy=seq(-3,3,length=201)
z=outer(vx,vy,f)
set.seed(1)
X=rmnorm(1500,varcov=S)
xhist <- hist(X[,1], plot=FALSE)
yhist <- hist(X[,2], plot=FALSE)
top <- max(c(xhist$density, yhist$density,dnorm(0)))
nf <- layout(matrix(c(2,0,1,3),2,2,byrow=TRUE), c(3,1), c(1,3), TRUE)
par(mar=c(3,3,1,1))
image(vx,vy,z,col=rev(heat.colors(101)))
points(X,cex=.2)
par(mar=c(0,3,1,1))
barplot(xhist$density, axes=FALSE, ylim=c(0, top), space=0,col="light green") lines((density(X[,1])$x-xhist$breaks[1])/diff(xhist$breaks)[1],
dnorm(density(X[,1])$x),col="red") par(mar=c(3,0,1,1)) barplot(yhist$density, axes=FALSE, xlim=c(0, top), space=0,
horiz=TRUE,col="light green")
lines(dnorm(density(X[,2])$x),(density(X[,2])$x-yhist$breaks[1])/ diff(yhist$breaks)[1],col="red")

That was the simple part.

• The Bivariate Student-t distribution

Consider now another elliptical distribution. But we want here to normalize the margins. Thus, instead of a pair $(X,Y)$, we would like to consider the pair $(\Phi^{-1}(T_\nu(X)),\Phi^{-1}(T_\nu(Y)))$, so that the marginal distributions are $\mathcal{N}(0,1)$. The new density is obtained simply since the transformation is a one-to-one increasing transformation. Here, we have

k=3
r=.5
G=function(x) qnorm(pt(x,df=k))
dg=function(x) dt(x,df=k)/dnorm(qnorm(pt(x,df=k)))
Ginv=function(x) qt(pnorm(x),df=k)
S=matrix(c(1,r,r,1),2,2)
f=function(x,y) dmt(cbind(Ginv(x),Ginv(y)),S=S,df=k)/(dg(x)*dg(y))
vx=seq(-3,3,length=201)
vy=seq(-3,3,length=201)
z=outer(vx,vy,f)
set.seed(1)
Z=rmt(1500,S=S,df=k)
X=G(Z)

Because we considered a nonlinear transformation of the margins, the level curves are no longer elliptical. But there is still some kind of symmetry.

• The Exchangeable Case with Conditionally Independent Random Variables

We did consider the case where $X$ and $Y$ with independent random variables, given $\Theta$, and that both variables are exponentially distributed, with parameter $\Theta$. As we’ve seen in class, it might be difficult to visualize that sample, unless we have log scales on both axis. But instead of a log transformation, why not consider a transformation so that margins will be $\mathcal{N}(0,1)$. The only technical problem is that we do not have the (nonconditional) distributions of the margins. Well, we have them, but they are integral based. From a computational point of view, that’s not a bit deal… Computations might take a while, but we can visualize the density using the following code (here, we assume that  is Gamma distributed)

a=.6
b=1
h=.0001
G=function(x) qnorm(ifelse(x<0,0,integrate(function(z) pexp(x,z)*
dgamma(z,a,b),lower=0,upper=Inf)$value)) Ginv=function(x) uniroot(function(z) G(z)-x,lower=-40,upper=1e5)$root
dg=function(x) (Ginv(x+h)-Ginv(x-h))/2/h
H=function(xy) integrate(function(z) dexp(xy[2],z)*dexp(xy[1],z)*
dgamma(z,a,b),lower=0,upper=Inf)$value f=function(x,y) H(c(Ginv(x),Ginv(y)))*(dg(x)*dg(y)) vx=seq(-3,3,length=151) vy=seq(-3,3,length=151) z=matrix(NA,length(vx),length(vy)) for(i in 1:length(vx)){ for(j in 1:length(vy)){ z[i,j]=f(vx[i],vy[j])}} set.seed(1) Theta=rgamma(1500,a,b) Z=cbind(rexp(1500,Theta),rexp(1500,Theta)) X=cbind(Vectorize(G)(Z[,1]),Vectorize(G)(Z[,2])) There is a small technical problem, but no big deal. Here, the joint distribution is quite different. Margins are – one more time – standard Gaussian, but the shape of the joint distribution is quite different, with an asymmetry from the lower (left) tail to the upper (right) tail. More details when we’ll introduce copulas. The only difference will be that the margins will be uniform on the unit interval, and not standard Gaussian. # Modèle de régression et interaction(s) entre facteurs Dans un modèle de régression, on veut écrire $\mathbb{E}(Y\vert \boldsymbol{X}=(X_1,\ldots,X_d))=\varphi(\boldsymbol{X})=\varphi(X_1,\ldots,X_d)$ Quand on se limite à un modèle linéaire, on écrit $\mathbb{E}(Y\vert \boldsymbol{X}=(X_1,\ldots,X_d))=\boldsymbol{X}^{\text{\sffamily T}}\boldsymbol{\beta}$ ou encore $\mathbb{E}(Y\vert \boldsymbol{X}=(X_1,\ldots,X_d))=\beta_0+\sum_{i=1}^d \beta_i X_i$ Mais on de doute que l’on rate quelque chose… en particulier, on va rater toutes les interactions possibles. On peut croiser les variables, et supposer que $\mathbb{E}(Y\vert \boldsymbol{X}=(X_1,\ldots,X_d))=\beta_0+\sum_{i=1}^d \beta_i X_i+\sum_{i=1}^{d-1}\sum_{j=i+1}^d \gamma_{i,j} X_i X_j$ qui peut s’étendre d’avantage, à l’ordre 3, $\beta_0+\sum_{i=1}^d \beta_i X_i+\sum_{i=1}^{d-1}\sum_{j=i+1}^d \gamma_{i,j} X_i X_j+\sum_{i=1}^{d-2}\sum_{j=i+1}^{d-1}\sum_{k=j+1}^{d} \delta_{i,j,k} X_i X_j X_k$ voire davantage. Supposons que nos variables $X_i$ soient ici qualitatives, et plus précisément binaires. Prenons un exemple simple, avec des données (classiques) en risque de crédit1. On peut trouver la base via library(evtree) db=GermanCredit ou encore directement myVariableNames = c("checking_status","duration","credit_history", "purpose","credit_amount","savings","employment","installment_rate", "personal_status","other_parties","residence_since","property_magnitude", "age","other_payment_plans","housing","existing_credits","job", "num_dependents","telephone","foreign_worker","class") GermanCredit = read.table( "http://archive.ics.uci.edu/ml/machine-learning-databases/statlog/german/german.data", header=FALSE,col.names=myVariableNames) Retenons pour commencer trois variables explicatives, db=data.frame(Y=GermanCredit$class-1,
X1=GermanCredit$checking_status%in%c("A12","A13"), X2=GermanCredit$credit_history%in%c("A30","A31"),
X3=GermanCredit$savings%in%c("A61","A62")) reg=glm(Y~X1+X2+X3,data=db,family=binomial) summary(reg) La régression sans interaction donne ici Call: glm(formula = Y ~ X1 + X2 + X3, family = binomial, data = db) Deviance Residuals: Min 1Q Median 3Q Max -1.5431 -0.8421 -0.6295 1.3994 1.9999 Coefficients: Estimate Std. Error z value Pr(>|z|) (Intercept) -1.8544 0.1699 -10.915 < 2e-16 *** X1TRUE 0.3363 0.1496 2.249 0.0245 * X2TRUE 1.3462 0.2347 5.735 9.76e-09 *** X3TRUE 1.0001 0.1787 5.596 2.19e-08 *** --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 (Dispersion parameter for binomial family taken to be 1) Null deviance: 1221.7 on 999 degrees of freedom Residual deviance: 1143.6 on 996 degrees of freedom AIC: 1151.6 Number of Fisher Scoring iterations: 4 Il existe plusieurs interactions possibles ici (limitons nous aux paires). C’est ce que l’on observe quand on fait la régression reg=glm(Y~X1+X2+X3+X1:X2+X1:X3+X2:X3,data=db,family=binomial) summary(reg) Call: glm(formula = Y ~ X1 + X2 + X3 + X1:X2 + X1:X3 + X2:X3, family = binomial, data = db) Deviance Residuals: Min 1Q Median 3Q Max -1.5369 -0.8281 -0.6439 1.3954 1.9638 Coefficients: Estimate Std. Error z value Pr(>|z|) (Intercept) -1.77109 0.20070 -8.825 < 2e-16 *** X1TRUE 0.30296 0.33737 0.898 0.369186 X2TRUE 0.88353 0.54255 1.628 0.103421 X3TRUE 0.87709 0.22583 3.884 0.000103 *** X1TRUE:X2TRUE -0.37917 0.49343 -0.768 0.442225 X1TRUE:X3TRUE 0.09178 0.37278 0.246 0.805522 X2TRUE:X3TRUE 0.80923 0.58185 1.391 0.164293 --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 (Dispersion parameter for binomial family taken to be 1) Null deviance: 1221.7 on 999 degrees of freedom Residual deviance: 1141.0 on 993 degrees of freedom AIC: 1155 Number of Fisher Scoring iterations: 4 On peut faire un dessin pour visualiser les interactions : on a trois sommets (nos trois variables), et on visualiser les interactions indices=cbind(c(1,2,3),c(1,1,2),c(2,3,3)) k=3 theta=pi/2+2*pi*(0:(k-1))/k sommetX=cos(theta) sommetY=sin(theta) plot(sommetX,sommetY,cex=1,axes=FALSE,xlab="",ylab="", xlim=c(-1.5,1.5),ylim=c(-1.5,1.5)) for(i in 1:nrow(indices)){ segments(sommetX[indices[i,2]],sommetY[indices[i,2]], sommetX[indices[i,3]],sommetY[indices[i,3]],col="grey") text(mean(sommetX[indices[i,2:3]]),mean(sommetY[indices[i,2:3]]), trunc(10000*coefficients(reg)[1+k+i])/10000) } points(sommetX,sommetY,cex=6,pch=19,col="yellow") points(sommetX,sommetY,cex=6,pch=1) text(sommetX,sommetY,1:k) ce qui donne ici, pour nos trois variables Ce modèle pourrait sembler incomplet, car on ne regarde que les interactions entre les modalités, par paires. En fait, c’est parce qu’il manque (visuellement) les variables non-croisées. On peut les rajouter si on veut (au risque de surcharger le dessin) cercle=function(c,r,cl) lines(c[1]+r*cos(seq(0,2*pi,length=501)), c[2]+r*sin(seq(0,2*pi,length=501)),col=cl) reg=glm(Y~X1+X2+X3+X1:X2+X1:X3+X2:X3,data=db,family=binomial) indices=cbind(c(1,2,3),c(1,1,2),c(2,3,3)) k=3 theta=pi/2+2*pi*(0:(k-1))/k sommetX=cos(theta) sommetY=sin(theta) plot(sommetX,sommetY,cex=1,axes=FALSE,xlab="",ylab="",xlim=c(-1.5,1.5),ylim=c(-1.5,1.5)) for(i in 1:nrow(indices)){ segments(sommetX[indices[i,2]],sommetY[indices[i,2]], sommetX[indices[i,3]],sommetY[indices[i,3]],col="grey") text(mean(sommetX[indices[i,2:3]]),mean(sommetY[indices[i,2:3]]), trunc(10000*coefficients(reg)[1+k+i])/10000) } for(i in 1:k){ cercle(c(cos(theta)[i]*1.18,sin(theta)[i]*1.18),.18,"grey") text(cos(theta)[i]*1.35,sin(theta)[i]*1.35, trunc(10000*coefficients(reg)[1+i])/10000) } points(sommetX,sommetY,cex=6,pch=19,col="yellow") points(sommetX,sommetY,cex=6,pch=1) text(sommetX,sommetY,1:k) soit ici Si on change le ‘sens‘ de nos variables (en recodant a l’envers, en permutant les vrais et les faux), on obtient le graphique suivant dbinv=db dbinv[,2:k]=1-dbinv[,2:k] reg=glm(Y~X1+X2+X3+X1:X2+X1:X3+X2:X3,data=dbinv,family=binomial) indices=cbind(c(1,2,3),c(1,1,2),c(2,3,3)) k=3 theta=pi/2+2*pi*(0:(k-1))/k sommetX=cos(theta) sommetY=sin(theta) plot(sommetX,sommetY,cex=1,axes=FALSE,xlab="",ylab="",xlim=c(-1.5,1.5),ylim=c(-1.5,1.5)) for(i in 1:nrow(indices)){ segments(sommetX[indices[i,2]],sommetY[indices[i,2]], sommetX[indices[i,3]],sommetY[indices[i,3]],col="grey") text(mean(sommetX[indices[i,2:3]]),mean(sommetY[indices[i,2:3]]), trunc(10000*coefficients(reg)[1+k+i])/10000) } for(i in 1:k){ cercle(c(cos(theta)[i]*1.18,sin(theta)[i]*1.18),.18,"grey") text(cos(theta)[i]*1.35,sin(theta)[i]*1.35, trunc(10000*coefficients(reg)[1+i])/10000) } points(sommetX,sommetY,cex=6,pch=19,col="yellow") points(sommetX,sommetY,cex=6,pch=1) text(sommetX,sommetY,1:k) qui peut alors être comparé au graphique précédant Avec 5 variables, on augmente les interactions possibles… même si beaucoup risquent d’être non-significatifs. On peut déjà se focaliser sur les paires possibles d’interactions croisées. Pour simplifier le code, on va utiliser deux fonctions locales, vrepeach=function(x,e){ v=NULL for(i in 1:length(e)){v=c(v,rep(x[i],each=e[i]))} return(v)} vreplength=function(x,l){ v=NULL for(i in 1:length(l)){v=c(v,x[l[i]:length(x)])} return(v)} et ensuite, on adapte le code précédant indices=cbind(1:(k*(k-1)/2),vrepeach(1:(k-1),(k-1):1),vreplength(2:k,1:(k-1))) formule="Y~1" for(i in 1:k) formule=paste(formule,"+X",i,sep="") for(i in 1:nrow(indices)) formule=paste(formule,"+X",indices[i,2],":X",indices[i,3],sep="") reg=glm(formule,data=db,family=binomial) theta=pi/2+2*pi*(0:(k-1))/k sommetX=cos(theta) sommetY=sin(theta) plot(sommetX,sommetY,cex=1,axes=FALSE,xlab="",ylab="",xlim=c(-1.5,1.5),ylim=c(-1.5,1.5)) for(i in 1:nrow(indices)){ segments(sommetX[indices[i,2]],sommetY[indices[i,2]], sommetX[indices[i,3]],sommetY[indices[i,3]],col="grey") text(mean(sommetX[indices[i,2:3]]),mean(sommetY[indices[i,2:3]]), trunc(10000*coefficients(reg)[1+k+i])/10000) } for(i in 1:k){ cercle(c(cos(theta)[i]*1.18,sin(theta)[i]*1.18),.18,"grey") text(cos(theta)[i]*1.35,sin(theta)[i]*1.35, trunc(10000*coefficients(reg)[1+i])/10000) } points(sommetX,sommetY,cex=6,pch=19,col="yellow") points(sommetX,sommetY,cex=6,pch=1) text(sommetX,sommetY,1:k) ce qui donne un schéma plus complexe, On peut aussi prendre juste 2 variables, prenant 3 et 4 modalités respectivement. On va extraire deux variables indicatrices pour la première (la modalité restante sera la modalité de référence) et trois pour la seconde, db=data.frame(Y=GermanCredit$class-1,
X1=GermanCredit$checking_status=="A12", X2=GermanCredit$checking_status=="A13",
X3=GermanCredit$checking_status=="A14", X4=GermanCredit$employment%in%c("A72","A73"),
X5=GermanCredit$employment%in%c("A74","A75")) k=5 indices=cbind(1:(k*(k-1)/2),vrepeach(1:(k-1),(k-1):1),vreplength(2:k,1:(k-1))) formule="Y~1" for(i in 1:k) formule=paste(formule,"+X",i,sep="") for(i in 1:nrow(indices)) formule=paste(formule,"+X",indices[i,2],":X",indices[i,3],sep="") reg=glm(formule,data=db,family=binomial) theta=pi/2+2*pi*(0:(k-1))/k sommetX=cos(theta) sommetY=sin(theta) plot(sommetX,sommetY,cex=1,axes=FALSE,xlab="",ylab="",xlim=c(-1.5,1.5),ylim=c(-1.5,1.5)) for(i in 1:nrow(indices)){ if(!is.na(coefficients(reg)[1+k+i])){ segments(sommetX[indices[i,2]],sommetY[indices[i,2]], sommetX[indices[i,3]],sommetY[indices[i,3]],col="grey") text(mean(sommetX[indices[i,2:3]]),mean(sommetY[indices[i,2:3]]), trunc(10000*coefficients(reg)[1+k+i])/10000) }} for(i in 1:k){ cercle(c(cos(theta)[i]*1.18,sin(theta)[i]*1.18),.18,"grey") text(cos(theta)[i]*1.35,sin(theta)[i]*1.35, trunc(10000*coefficients(reg)[1+i])/10000) } points(sommetX,sommetY,cex=6,pch=19,col="yellow") points(sommetX,sommetY,cex=6,pch=1) text(sommetX,sommetY,1:k) On voit que plusieurs interactions ne sont alors plus possibles, sur la partie gauche (les trois modalités de la même variable) et sur la partie droite On peut d’ailleurs simplifier les graphs, en ne visualisant que les interactions significatives. indices=cbind(1:(k*(k-1)/2),vrepeach(1:(k-1),(k-1):1),vreplength(2:k,1:(k-1))) formule="Y~1" for(i in 1:k) formule=paste(formule,"+X",i,sep="") for(i in 1:nrow(indices)) formule=paste(formule,"+X",indices[i,2],":X",indices[i,3],sep="") reg=glm(formule,data=db,family=binomial) theta=pi/2+2*pi*(0:(k-1))/k sommetX=cos(theta) sommetY=sin(theta) plot(sommetX,sommetY,cex=1,axes=FALSE,xlab="",ylab="",xlim=c(-1.5,1.5),ylim=c(-1.5,1.5)) for(i in 1:nrow(indices)){ if(!is.na(coefficients(reg)[1+k+i])){ if(summary(reg)$coefficients[1+k+i,4]<.1){
segments(sommetX[indices[i,2]],sommetY[indices[i,2]],
sommetX[indices[i,3]],sommetY[indices[i,3]],col="grey")
text(mean(sommetX[indices[i,2:3]]),mean(sommetY[indices[i,2:3]]),
trunc(10000*coefficients(reg)[1+k+i])/10000)
}}}
for(i in 1:k){
if(summary(reg)$coefficients[1+i]<.1){ cercle(c(cos(theta)[i]*1.18,sin(theta)[i]*1.18),.18,"grey") text(cos(theta)[i]*1.35,sin(theta)[i]*1.35, trunc(10000*coefficients(reg)[1+i])/10000) }} points(sommetX,sommetY,cex=6,pch=19,col="yellow") points(sommetX,sommetY,cex=6,pch=1) text(sommetX,sommetY,1:k) soit ici Ici, une seule interactions croisée est significative, et presque toutes les variables le sont. Et si on reprend le modèle avec 5 facteurs, db=data.frame(Y=GermanCredit$class-1,X1=GermanCredit$checking_status%in%c("A12","A13"), X2=GermanCredit$credit_history%in%c("A30","A31"),
X3=GermanCredit$savings%in%c("A61","A62"), X4=GermanCredit$employment%in%c("A71","A72"),
X5=GermanCredit$other_payment_plans=="A143") indices=cbind(1:(k*(k-1)/2),vrepeach(1:(k-1),(k-1):1),vreplength(2:k,1:(k-1))) formule="Y~1" for(i in 1:k) formule=paste(formule,"+X",i,sep="") for(i in 1:nrow(indices)) formule=paste(formule,"+X",indices[i,2],":X",indices[i,3],sep="") reg=glm(formule,data=db,family=binomial) theta=pi/2+2*pi*(0:(k-1))/k sommetX=cos(theta) sommetY=sin(theta) plot(sommetX,sommetY,cex=1,axes=FALSE,xlab="",ylab="",xlim=c(-1.5,1.5),ylim=c(-1.5,1.5)) for(i in 1:nrow(indices)){ if(!is.na(coefficients(reg)[1+k+i])){ if(summary(reg)$coefficients[1+k+i,4]<.1){
segments(sommetX[indices[i,2]],sommetY[indices[i,2]],
sommetX[indices[i,3]],sommetY[indices[i,3]],col="grey")
text(mean(sommetX[indices[i,2:3]]),mean(sommetY[indices[i,2:3]]),
trunc(10000*coefficients(reg)[1+k+i])/10000)
}}}
for(i in 1:k){
if(summary(reg)\$coefficients[1+i]<.1){
cercle(c(cos(theta)[i]*1.18,sin(theta)[i]*1.18),.18,"grey")
text(cos(theta)[i]*1.35,sin(theta)[i]*1.35,
trunc(10000*coefficients(reg)[1+i])/10000)
}}
points(sommetX,sommetY,cex=6,pch=19,col="yellow")
points(sommetX,sommetY,cex=6,pch=1)
text(sommetX,sommetY,1:k)

on obtient

Je ne sais pas si mes graphiques sont pertinents, ou pas. Mais je trouve ça joli. En fait, je suis tombé un peu par hasard2 sur les Tables de Taguchi, développées par Gen’ichi Taguchi (田口 玄一). Le soucis est que je n’ai rien compris… Enfin, disons que je croyais comprendre, puis j’ai continué à faire des dessins… Si quelqu’un pourrait m’expliquer sur mon exemple les graphiques de Taguchi, je suis preneur ! car je doute que ce soit ce que je fais depuis tout à l’heure…

1. Cette base est largement utilisée dans le quatrième chapitre de Computational Actuarial Science with R, à paraître dans les mois à venir.

2.En l’occurence, le hasard est @Benavent qui a suscité ma curiosité ce matin en me parlant de ces tables, dont je n’avais alors jamais entendu parlé ! J’avais même lu rapidement Taniguchi (谷口 ジロー) et je ne voyais pas le rapport avec les statistiques….

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

# Ecole d’été en actuariat à Brest

Exposé sur mesures de risques et dépendance, à Brest, dans le cadre de l’école d’été, la première partie est en ligne

ainsi que la seconde