{"word_count": 771, "character_count": 12619}
Dans un modèle de régression, on veut écrireQuand on se limite à un modèle linéaire, on écritou encoreMais 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 quequi 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 vialibrary(evtree)db=GermanCreditou encore directementmyVariableNames = 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 iciCall: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(&#62;|z|)    (Intercept)  -1.8544     0.1699 -10.915  &#60; 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 freedomResidual deviance: 1143.6  on 996  degrees of freedomAIC: 1151.6Number of Fisher Scoring iterations: 4Il existe plusieurs interactions possibles ici (limitons nous aux paires). C'est ce que l'on observe quand on fait la régressionreg=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(&#62;|z|)    (Intercept)   -1.77109    0.20070  -8.825  &#60; 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 freedomResidual deviance: 1141.0  on 993  degrees of freedomAIC: 1155Number of Fisher Scoring iterations: 4On peut faire un dessin pour visualiser les interactions : on a trois sommets (nos trois variables), et on visualiser les interactionsindices=cbind(c(1,2,3),c(1,1,2),c(2,3,3))k=3theta=pi/2+2*pi*(0:(k-1))/ksommetX=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 variablesCe 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=3theta=pi/2+2*pi*(0:(k-1))/ksommetX=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 iciSi on change le 'sens' de nos variables (en recodant a l'envers, en permutant les vrais et les faux), on obtient le graphique suivantdbinv=dbdbinv[,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=3theta=pi/2+2*pi*(0:(k-1))/ksommetX=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édantAvec 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=NULLfor(i in 1:length(e)){v=c(v,rep(x[i],each=e[i]))}return(v)}vreplength=function(x,l){v=NULLfor(i in 1:length(l)){v=c(v,x[l[i]:length(x)])}return(v)}et ensuite, on adapte le code précédantindices=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))/ksommetX=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=5indices=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))/ksommetX=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 droiteOn 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))/ksommetX=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]&#60;.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]&#60;.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 iciIci, 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))/ksommetX=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]&#60;.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]&#60;.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 obtientJe 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....