Dans un modèle de régression, on veut écrire
Quand on se limite à un modèle linéaire, on écrit
ou encore
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
qui peut s’étendre d’avantage, à l’ordre 3,
voire davantage.
Supposons que nos variables 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….