La régression (comme tout calcul d’espérance conditionnelle) est un problème de projections
Un théorème intéressant est le théorème dit de Frisch-Waugh, permettant de comprendre la différence fondamentale entre un modèle de régression multiple, et les modèles de régression simple. On reviendra sur ce point en cours en évoquant rapidement le paradoxe de Simpson. La formulation est la suivante: on veut régresser 
sur 
et 
, deux ensembles (a priori disjoints) de variables explicatives,
Les équations normales (associées au problème de minimisation de la somme des carrés des erreurs) sont
de telle sorte qu’à l’optimum on peut relier les estimateurs des deux jeux de paramètres par
La première partie correspond à la régression de sur 
, mais il reste un second terme dès lors que et ne sont pas orthogonales. Notons 
la matrice de projection (orthogonale) sur
Alors on peut écrire
En posant
et 
,
on retrouve un modèle linéaire classique (à condition de travailler sur la projection des variables sur le sous-espace engendré par , i.e. en transformant les variables),
Pour aller plus loin sur la géométrie des moindres carrées, et sur le théorème de Frisch-Waugh, je peux renvoyer à des notes de cours, et à quelques transparents
Le graphique ci-dessous correspond au cas où les variables explicatives sont orthogonales, et dans ce cas, la régression multiple est équivalente à deux régressions simples
Le graphique ci-dessous au cas non orthogonal,
Le code qui permet de vérifier ces histoires de projections successives est relativement simple. Tout d’abord on importe les données, et on regarde le modèle global,
> chicago=read.table( + "http://freakonometrics.free.fr/chicago.txt", + header=TRUE,sep=";") > Y=chicago$Fire > X1=chicago$X_1 > X2=chicago$X_2 > X3=chicago$X_3 > base=data.frame(Y,X1,X2,X3) > tail(base) Y X1 X2 X3 42 4.8 0.152 19 13.323 43 10.4 0.408 25 12.960 44 15.6 0.578 28 11.260 45 7.0 0.114 3 10.080 46 7.1 0.492 23 11.428 47 4.9 0.466 27 13.731 > regression=lm(Y~X1+X2+X3) > summary(regression) Call: lm(formula = Y ~ X1 + X2 + X3) Residuals: Min 1Q Median 3Q Max -9.737 -4.565 -1.479 3.751 16.079 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 22.07525 6.19447 3.564 0.000910 *** X1 -0.62764 5.28130 -0.119 0.905953 X2 0.22378 0.06161 3.632 0.000744 *** X3 -1.55059 0.38195 -4.060 0.000204 *** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Residual standard error: 6.527 on 43 degrees of freedom Multiple R-squared: 0.4417, Adjusted R-squared: 0.4027 F-statistic: 11.34 on 3 and 43 DF, p-value: 1.314e-05 > n=length(Y) > X=matrix(c(rep(1,n),X1,X2,X3),n,4) > X[1:5,] [,1] [,2] [,3] [,4] [1,] 1 0.604 29 11.744 [2,] 1 0.765 44 9.323 [3,] 1 0.735 36 9.948 [4,] 1 0.669 37 10.656 [5,] 1 0.814 53 9.730
Ensuite, on va projeter seulement sur les deux premières variables (et la constante)
> FWX1=X[,1:3] > regression12=lm(Y~X1+X2) > solve(t(FWX1)%*%FWX1)%*%t(FWX1)%*%Y [,1] [1,] 0.08069764 [2,] 11.56913900 [3,] 0.15108490 > summary(regression12)$coefficients Estimate Std. Error t value Pr(>|t|) (Intercept) 0.08069764 3.49182154 0.02311047 0.98166664 X1 11.56913900 5.05014782 2.29085156 0.02681947 X2 0.15108490 0.06854408 2.20420043 0.03278652 > > > FWX2=X[,4] > H1=FWX1%*%solve(t(FWX1)%*%FWX1)%*%t(FWX1) > M1=diag(rep(1,n))-H1 > FWX2s=M1%*%FWX2 > FWYs =M1%*%Y > (beta2=solve(t(FWX2s)%*%FWX2s)%*%t(FWX2s)%*%FWYs) [,1] [1,] -1.550594 > summary(regression)$coefficients[4] [1] -1.550594 > > (beta1=solve(t(FWX1)%*%FWX1)%*%t(FWX1)%*%Y- + solve(t(FWX1)%*%FWX1)%*%t(FWX1)%*%FWX2%*%beta2) [,1] [1,] 22.0752495 [2,] -0.6276442 [3,] 0.2237765 > summary(regression)$coefficients[1:3] [1] 22.0752495 -0.6276442 0.2237765
On retrouve bien l’estimateur du dernier paramètre à l’aide du théorème de Frisch-Waugh, puis en utilisant les équations normales, on en déduit les trois premiers.