{"word_count": 705, "character_count": 8425}
Cette fin de semaine, je mettais en ligne un court billet du deuxième effet kiss-cool où je rappelais que quand on fait une régression sur plusieurs variables explicatives corrélées, ce n'est pas équivalent à faire plusieurs régressions simple. C'est ce que disait le théorème de Frisch-Waugh, publié il y a près de 90 ans. Mais je glissais qu'il était possible de faire des régressions en cascade, sans aller beaucoup plus loin. Il faut relire l'article de Michael Lovell qui donner une interprétation géométrique de ce dessin. Mais programmons le, pour le voir... Considérons une base avec 3 variables explicatives, correspondant à des nombres d’incendies à Chicagochicago=read.table("http://freakonometrics.free.fr/chicago.txt", header=TRUE,sep=";")La régression multiple s'écrit ici[latex display="true"]y_i=\beta_0+\beta_1 x_{1,i}+\beta_2 x_{2,i}+\beta_3 x_{3,i}+\varepsilon_i[/latex]On peut estimer ces paramètres par moindre carrés,reg=lm(Fire~.,data=chicago)summary(reg)Coefficients:            Estimate Std. Error t value  Pr(|t|)    (Intercept) 22.07525    6.19447   3.564 0.000910 ***X_1         -0.62764    5.28130  -0.119 0.905953    X_2          0.22378    0.06161   3.632 0.000744 ***X_3         -1.55059    0.38195  -4.060 0.000204 ***et on a alors la prévision suivante[latex display="true"]\widehat{y}_i=\widehat{\beta}_0+\widehat{\beta}_1 x_{1,i}+\widehat{\beta}_2 x_{2,i}+\widehat{\beta}_3 x_{3,i}[/latex]En fait, je peux même visualiser cette prévision.BC = cbind(x0=coefficients(reg)[1],           x1=coefficients(reg)[2]*chicago[,2],           x2=coefficients(reg)[3]*chicago[,3],           x3=coefficients(reg)[4]*chicago[,4],           yp=predict(reg),y=chicago$Fire)colr = 1:4dessin_fire = function(i,BC=BC){plot(0:3,0:3,xlim=c(0,35),col="white",ylim=c(-1,3),xlab="",ylab="",axes=FALSE)abline(v=BC[i,1],lty=2)rect(BC[i,1],3,BC[i,1]+BC[i,2],2,col=colr[1],border=NA)arrows(BC[i,1],2.5,BC[i,1]+BC[i,2],2.5,lwd=2,length=.1,col="white")rect(BC[i,1]+BC[i,2],2,BC[i,1]+BC[i,2]+BC[i,3],1,col=colr[2],border=NA)arrows(BC[i,1]+BC[i,2],1.5,BC[i,1]+BC[i,2]+BC[i,3],1.5,lwd=2,length=.1,col="white") rect(BC[i,1]+BC[i,2]+BC[i,3],1,BC[i,1]+BC[i,2]+BC[i,3]+BC[i,4],0,col=colr[3],border=NA)segments(BC[i,5],0,BC[i,5],1,lwd=3)arrows(BC[i,1]+BC[i,2]+BC[i,3],.5,BC[i,1]+BC[i,2]+BC[i,3]+BC[i,4],.5,lwd=2,length=.1,col="white") abline(v=BC[i,1],lty=2)rect(BC[i,5],-1,BC[i,6],0,col=colr[4],density=20,border=NA)segments(BC[i,6],0,BC[i,6],-1,lwd=3,col=colr[4])text(26,2.5,expression(X[1]),pos=4,col=colr[1])text(26,1.5,expression(X[2]),pos=4,col=colr[2])text(26,.5,expression(X[3]),pos=4,col=colr[3])text(26,-.5,expression(epsilon),pos=4,col=colr[4])axis(1)}Par exemple, pour la douzième observation de ma basedessin_fire(12)on obtientOn lit de haut en bas : on commence par constante, vers 22. Puis on rajoute trois fois rien, parce que [latex]x_1[/latex] était non-significative dans notre régression. Puis on rajoute un petit quelque chose à cause de [latex]x_2[/latex], pour passer à 25, et [latex]x_3[/latex] nous fait plonger de plus de 20 points, pour finir vers 3.67.Le trait noir est la prévision que l'on obtient, à partir des trois variables. En bleu, on peut même voir l'erreur.On peut noter que [latex]x_1[/latex] est non-significative, tout en étant relativement corrélée avec [latex]y[/latex]cor(chicago[,1:2])          Fire       X_1Fire 1.0000000 0.3773486X_1  0.3773486 1.0000000Car si [latex]x_1[/latex] est non-significative, c'est parce que cette variable est très corrélée avec une autre variable explicative. En fait, on peut faire une série de régressions en cascade, en commençant par [latex]x_1[/latex], et plus spécifiquement, on peut écrire[latex display="true"]\widehat{y}_i=\underbrace{\widehat{b}_0+\widehat{b}_1 x_{1,i}}_{(1)}+\underbrace{\widehat{b}_{0,2}+\widehat{b}_2 \tilde{x}_{2,i}}_{(2)}+\underbrace{\widehat{b}_{0,3}+\widehat{b}_3 \tilde{x}_{3,i}}_{(3)}[/latex]où le premier terme [latex](1)[/latex] est obtenu en régression simplement [latex]y[/latex] sur [latex]x_1[/latex] (oui, on fait juste une régression simple),  [latex display="true"]y_i=b_0+b_1 x_{1,i}+\eta_i[/latex]puis on apporte une petite correction, pour tenir compte de ce que nous apprend [latex]x_2[/latex] une fois la première régression effectuée. C'est exactement ce que raconte le théorème de Frish-Waugh : on va projeter [latex]y[/latex] sur [latex]x_1^{\perp}[/latex] (car c'est justement ce que [latex]x_1[/latex] n'explique pas), et [latex]x_2[/latex] sur [latex]x_1^{\perp}[/latex]  et faire la régression de ces deux projections, soit, pour [latex](2)[/latex] [latex display="true"]\Pi_{x_1^{\perp}}y_i=b_{0,2}+b_2 \Pi_{x_1^{\perp}}x_{2,i}+\eta_{2,i}[/latex]puis pour [latex](3)[/latex] [latex display="true"]\Pi_{(x_1,x_2)^{\perp}}y_i=b_{0,3}+b_3 \Pi_{(x_1,x_2)^{\perp}}x_{3,i}+\eta_{2,i}[/latex]En terme informatique, cela donnereg1=lm(Fire~X_1,data=chicago)reg2=lm(residuals(lm(Fire~X_1,data=chicago))~residuals(lm(X_2~X_1,data=chicago)))reg3=lm(residuals(lm(Fire~X_1+X_2,data=chicago))~residuals(lm(X_3~X_1+X_2,data=chicago)))soit iciBC123=cbind(x0=coefficients(reg1)[1],            x1=coefficients(reg1)[2]*chicago[,2],            x2=coefficients(reg2)[1]+coefficients(reg2)[2]*residuals(lm(X_2~X_1,data=chicago)),            x3=coefficients(reg3)[1]+coefficients(reg3)[2]*residuals(lm(X_3~X_1+X_2,data=chicago)),            yp=predict(reg),y=chicago$Fire)En adaptant la fonction de dessin précédante on obtientdessin_fire_123(12)Autrement dit, on commence avec une régression simple, en rouge: on part d'une constante vers 3.5, puis on augemente notre prévision, en tenant compte de [latex]x_1[/latex]. Puis on corrige avec [latex]x_2[/latex] (ou plutôt la partie de [latex]x_2[/latex] non expliquée par [latex]x_1[/latex]). sur ce qui n'avait pas été expliqué auparavant, ce qui pousse à revoir un peu à la baisse notre prévision, puis on tient compte de [latex]x_3[/latex] (ou là encore, sa projection sur l'orthogonal de [latex]x_1[/latex] et [latex]x_2[/latex]). On note que la prévision est très exactement la même qu'auparavant.Mais on peut aller plus loin... pourquoi commencer par [latex]x_1[/latex]? On peut aussi commencer par [latex]x_2[/latex]. On considère alors[latex display="true"]\widehat{y}_i=\underbrace{\widehat{b}_0+\widehat{b}_2 x_{2,i}}_{(2)}+\underbrace{\widehat{b}_{0,1}+\widehat{b}_1 \tilde{x}_{1,i}}_{(1)}+\underbrace{\widehat{b}_{0,3}+\widehat{b}_3 \tilde{x}_{3,i}}_{(3)}[/latex]où le premier terme [latex](2)[/latex] est obtenu en régression simplement [latex]y[/latex] sur [latex]x_2[/latex],  [latex display="true"]y_i=b_0+b_2 x_{2,i}+\eta_i[/latex]puis on apporte une petite correction, pour tenir compte de ce que nous apprend [latex]x_1[/latex] une fois la première régression effectuée, [latex](1)[/latex] [latex display="true"]\Pi_{x_2^{\perp}}y_i=b_{0,1}+b_1\Pi_{x_2^{\perp}}x_{1,i}+\eta_{1,i}[/latex]puis pour [latex](3)[/latex] [latex display="true"]\Pi_{(x_1,x_2)^{\perp}}y_i=b_{0,3}+b_3 \Pi_{(x_1,x_2)^{\perp}}x_{3,i}+\eta_{2,i}[/latex] Soitreg1=lm(Fire~X_2,data=chicago)reg2=lm(residuals(lm(Fire~X_2,data=chicago))~residuals(lm(X_1~X_2,data=chicago)))reg3=lm(residuals(lm(Fire~X_1+X_2,data=chicago))~residuals(lm(X_3~X_1+X_2,data=chicago)))BC213=cbind(x0=coefficients(reg1)[1],            x1=coefficients(reg1)[2]*chicago[,3],            x2=coefficients(reg2)[1]+coefficients(reg2)[2]*residuals(lm(X_1~X_2,data=chicago)),            x3=coefficients(reg3)[1]+coefficients(reg3)[2]*residuals(lm(X_3~X_1+X_2,data=chicago)),            yp=predict(reg),y=chicago$Fire)et visuellement, on obtientdessin_fire_213(12)Cette fois, la constante est un peu plus élevée, et on commence par faire une régression seulement avec [latex]x_2[/latex], pour arriver un peu en dessous de 10. Puis on corrige. On note que cette seconde correction nous ramène... très exactement comme dans le cas précédant. Mais si on pense en projection successives, on ne doit pas être surpris.On peut récapituler avec une autre prévision, soit avec le modèle de régresion multiple (mais là encore, [latex]x_1[/latex] n'explique pas grand chose)dessin_fire(15)soit en commencant par une régression simple sur [latex]x_1[/latex] puis en faisant des régression en cascade,dessin_fire_123(15)soit, alternativement, en commencant par une régression simple sur [latex]x_2[/latex]dessin_fire_213(15)Comme annoncé, cette trois approches sont équivalentes, et donnent très exactement la même prévision.