Tag Archives: Frisch

Régression sur des variables orthogonales

Hier, on a évoqué en cours le théorème de Frisch-Waugh (évoqué sur ce blog il y a quelques mois). Derrière, il y a l’estimation de modèles sur des sous-ensembles de valeurs. Formellement, on a un modèle de la forme

https://latex.codecogs.com/gif.latex?Y=\boldsymbol{X}\boldsymbol{\beta}+\varepsilon=\boldsymbol{X}_1\boldsymbol{\beta}_1+\boldsymbol{X}_2\boldsymbol{\beta}_2+\varepsilon
Les équations normales sont ici (cf transparents du cours, ou équations écrites au tableau) https://latex.codecogs.com/gif.latex?\boldsymbol{X}%27\boldsymbol{X}\widehat{\boldsymbol{\beta}}=\boldsymbol{X}%27Y, que l’on peut écrire sous la forme d’un système,

https://latex.codecogs.com/gif.latex?\left\{\begin{array}{l}\boldsymbol{X}_1%27\boldsymbol{X}_1%20\widehat{\boldsymbol{\beta}}_1+\boldsymbol{X}_1%27\boldsymbol{X}_2%20\widehat{\boldsymbol{\beta}}_2=\boldsymbol{X}_1%27Y%20\\%20\boldsymbol{X}_2%27\boldsymbol{X}_1%20\widehat{\boldsymbol{\beta}}_1+\boldsymbol{X}_2%27\boldsymbol{X}_2%20\widehat{\boldsymbol{\beta}}_2=\boldsymbol{X}_2%27Y%20%20\end{array}\right.
Le point intéressant est que si les variables dans le premier ensemble et les variables dans le second sont orthogonales, i.e. https://latex.codecogs.com/gif.latex?\boldsymbol{X}_1%27\boldsymbol{X}_2%20=\boldsymbol{0} et https://latex.codecogs.com/gif.latex?\boldsymbol{X}_2%27\boldsymbol{X}_1%20=\boldsymbol{0}, alors les termes diagonaux du système disparaissent, et on retrouve les équations normales des régressions simples, i.e.

https://latex.codecogs.com/gif.latex?\left\{\begin{array}{l}\boldsymbol{X}_1%27\boldsymbol{X}_1%20\widehat{\boldsymbol{\beta}}_1=\boldsymbol{X}_1%27Y%20\\%20\boldsymbol{X}_2%27\boldsymbol{X}_2%20\widehat{\boldsymbol{\beta}}_2=\boldsymbol{X}_2%27Y%20%20\end{array}\right.
On a ici une conséquence du théorème de Frisch-Waugh: si les variables sont orthogonales, les estimateurs dans la régression multiple coïncident avec ceux obtenus avec des régressions simples.

Afin de mieux clarifier ce point, écrivons ce qui se passerait avec seulement deux variables. On va associer la constante à la première, i.e.

https://latex.codecogs.com/gif.latex?Y=\boldsymbol{X}\boldsymbol{\beta}+\varepsilon=(1%20\%20\boldsymbol{X}_1)\begin{pmatrix}\boldsymbol{\beta}_0%20\\%20\boldsymbol{\beta}_1\end{pmatrix}+\boldsymbol{X}_2\boldsymbol{\beta}_2+\varepsilon
Dans ce cas, les équations normales sont

https://latex.codecogs.com/gif.latex?\left\{\begin{array}{l}\begin{pmatrix}1%20\\%20\boldsymbol{X}%27_1\end{pmatrix}(1%20\%20\boldsymbol{X}_1)\begin{pmatrix}\widehat{\boldsymbol{\beta}}_0%20\\%20\widehat{\boldsymbol{\beta}}_1\end{pmatrix}+\begin{pmatrix}1%20\\%20\boldsymbol{X}%27_1\end{pmatrix}\boldsymbol{X}_2%20\widehat{\boldsymbol{\beta}}_2=\begin{pmatrix}1%20\\%20\boldsymbol{X}%27_1\end{pmatrix}Y%20\\%20\boldsymbol{X}_2%27(1\%20\boldsymbol{X}_1)%20\begin{pmatrix}\widehat{\boldsymbol{\beta}}_0%20\\%20\widehat{\boldsymbol{\beta}}_1\end{pmatrix}+\boldsymbol{X}_2%27\boldsymbol{X}_2%20\widehat{\boldsymbol{\beta}}_2=\boldsymbol{X}_2%27Y%20%20\end{array}\right.
Bref, il faut que https://latex.codecogs.com/gif.latex?\boldsymbol%20X_1 et https://latex.codecogs.com/gif.latex?\boldsymbol%20X_2 soient orthogonales pour que le terme sur la seconde diagonale disparaisse. C’est ce qu’on retrouve dans la plupart des livres d’économétrie, il suffit de chercher le théorème de Frisch-Waugh dans l’index (même si le théorème de Frisch-Waugh va beaucoup plus loin puisqu’il explique quels ajustement faire pour lier les estimateurs dans une régression simple, et dans le cas multiple).

Mais est-ce un résultat purement théorique ou peut-on l’invoquer dans un cas pratique ? Autrement dit, si deux variables semblent orthogonales, que se passent-il numériquement ? Si je simule des variables explicatives indépendantes, les estimateurs coïncident-ils (comme le suggère la théorie) ?

Considérons un modèle linéaire avec deux variables explicative que l’on supposera indépendantes, mais ayant toutes deux un impact sur la variable explicative. Le plus simple est de considérer un triple de variables telles que

https://latex.codecogs.com/gif.latex?\boldsymbol{Z}=\begin{pmatrix}%20Y%20\\%20X_1\\X_2\end{pmatrix}\sim\mathcal{N}\left(\begin{pmatrix}%20\mu_Y%20\\%20\mu_{X_1}\\\mu_{X_2}\end{pmatrix},\begin{pmatrix}%201%20&%20\star%20&%20\star%20\\%20\star%20&%201&0\\\star&0&1\end{pmatrix}\right)
où https://latex.codecogs.com/gif.latex?\star désigne des valeurs non nulles. Par exemple,

> set.seed(1)
> n=100
> S=matrix(c(1,.4,.3,.4,1,0,.3,0,1),3,3)
> S
[,1] [,2] [,3]
[1,]  1.0  0.4  0.3
[2,]  0.4  1.0  0.0
[3,]  0.3  0.0  1.0
> M=c(1,3,6)
> Z=rmnorm(n,M,varcov=S)
> Y=Z[,1]; X1=Z[,2]; X2=Z[,3]
> base=data.frame(Y,X1,X2)
> cor(base)
Y          X1          X2
Y  1.0000000  0.37796580  0.27859381
X1 0.3779658  1.00000000 -0.05060539
X2 0.2785938 -0.05060539  1.00000000

On a simulé ici deux variables indépendantes, mais toutes deux corrélées avec la variable d’intérêt. On notera que numériquement

https://latex.codecogs.com/gif.latex?\frac{1}{n}\sum_{i=1}^n%20X_{1,i}X_{2,i}%20\neq%20\frac{1}{n}\sum_{i=1}^n%20X_{1,i}\cdot\frac{1}{n}\sum_{i=1}^n%20X_{2,i}

(autrement dit, en moyenne, on simulera des données indépendantes, mais dans ce cas précis, les données ne sont pas rigoureusement orthogonales).

> sum(X1*X2)/n
[1] 18.20194
> sum(X1)*sum(X2)/n^2
[1] 18.25097
Mais un test rapide nous pousse à accepter l’hypothèse d’indépendance entre les deux variables (sans trop d’hésitations).

> cor.test(X1,X2)

Pearsons product-moment correlation

data:  X1 and X2
t = -0.5016, df = 98, p-value = 0.6171
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
-0.2445923  0.1472766
sample estimates:
cor
-0.05060539

Si on regarde maintenant la régression,

> R=lm(Y~X1+X2,data=base)
> summary(R)

Call:
lm(formula = Y ~ X1 + X2, data = base)

Residuals:
Min       1Q   Median       3Q      Max
-1.87885 -0.57656  0.07957  0.52719  1.77407

Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept)  4.71180    0.28246  16.682  < 2e-16 ***
X1           0.37232    0.08441   4.411 2.66e-05 ***
X2           0.25975    0.07755   3.350  0.00115 **
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.7954 on 97 degrees of freedom
Multiple R-squared: 0.2317,	Adjusted R-squared: 0.2159
F-statistic: 14.63 on 2 and 97 DF,  p-value: 2.803e-06

Quel est le lien avec les régressions simples, sur chacune des variables ? Si on régresse sur la première, on obtient

> R1=lm(Y~X1,data=base)
> summary(R1)

Call:
lm(formula = Y ~ X1, data = base)

Residuals:
Min       1Q   Median       3Q      Max
-1.78973 -0.54060  0.04036  0.44138  1.94389

Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept)  5.03165    0.27934  18.013  < 2e-16 ***
X1           0.35802    0.08859   4.041 0.000106 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.8358 on 98 degrees of freedom
Multiple R-squared: 0.1429,	Adjusted R-squared: 0.1341
F-statistic: 16.33 on 1 and 98 DF,  p-value: 0.0001058

et si on régresse sur la seconde

> R2=lm(Y~X2,data=base)
> summary(R2)

Call:
lm(formula = Y ~ X2, data = base)

Residuals:
Min      1Q  Median      3Q     Max
-1.8652 -0.5507  0.1104  0.6109  2.1435

Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept)  5.85053    0.12495  46.824  < 2e-16 ***
X2           0.24244    0.08443   2.872  0.00501 **
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.867 on 98 degrees of freedom
Multiple R-squared: 0.07761,	Adjusted R-squared: 0.0682
F-statistic: 8.246 on 1 and 98 DF,  p-value: 0.005005

Si les données ne sont pas rigoureusement orthogonales, alors https://latex.codecogs.com/gif.latex?\widehat{\beta}_1 dans le modèle multiple est légèrement différent de https://latex.codecogs.com/gif.latex?\widehat{\beta}_1 obtenu dans le modèle simple. Si on regarde la distribution de nos estimateurs (n’oublions pas qu’il y a une marge d’erreur), ils sont proches, ou au moins statistiquement proches (avec à gauche la distribution de https://latex.codecogs.com/gif.latex?\widehat{\beta}_1 et à droite celle de https://latex.codecogs.com/gif.latex?\widehat{\beta}_2)

Si on génère 100 bases de données avec ce modèle (avec 500 observations à chaque fois) on obtient les estimateurs suivants (à gauche les https://latex.codecogs.com/gif.latex?\widehat{\beta}_1 et à droite les https://latex.codecogs.com/gif.latex?\widehat{\beta}_2, avec en abscisse le modèle multiple et en ordonnées, le modèle simple).

(la valeur centrale étant la vraie valeur, utilisée pour les simulations). On pourrait se demander si les valeurs éloignées du barycentre (ou de cette première bissectrice) correspondent à des cas où l’hypothèse d’indépendance entre les deux variables est rejetée ? On peut mettre des couleurs différentes en fonction de la https://latex.codecogs.com/gif.latex?p-value: plus la https://latex.codecogs.com/gif.latex?p-value est faible (rejet de l’hypothèse d’indépendance) plus la couleur sera vive (et rouge), et plus la https://latex.codecogs.com/gif.latex?p-value est forte (acceptation de l’hypothèse d’indépendance) plus la couleur sera terne (et jaune, voire blanche),

sur 1,000 simulations, on obtient

Le code est ici

> set.seed(1)
> library(mnormt)
> S=matrix(c(1,.4,.3,.4,1,0,.3,0,1),3,3)
> BX1=BX2=matrix(NA,1000,2)
> P=rep(NA,1000)
> for(s in 1:1000){
+   n=500
+   Z=rmnorm(n,c(6,3,1),varcov=S)
+   Y=Z[,1]; X1=Z[,2]; X2=Z[,3]
+   base=data.frame(Y,X1,X2)
+   R=lm(Y~X1+X2,data=base)
+   BX1[s,1]=coefficients(R)[2]
+   BX2[s,1]=coefficients(R)[3]
+   BX1[s,2]=coefficients(lm(Y~X1,data=base))[2]
+   BX2[s,2]=coefficients(lm(Y~X2,data=base))[2]
+   P[s]=cor.test(X1,X2)$p.value
+ }
> COL=(heat.colors(100))
> plot(BX1[,1],BX1[,2],pch=19,cex=.6,col=COL[trunc(100*P)])
> abline(a=0,b=1,lty=2)

Autrement dit le fait que l’estimation dans le modèle multiple et les modèles simples diffèrent beaucoup n’a rien a voir avec la https://latex.codecogs.com/gif.latex?p-value d’un test d’indépendance entre les deux variables… étonnant non ? En revanche, si on rajoute une observation pour rendre le modèle parfaitement non-corrélé, par exemple,

> difference=sum(X1*X2)*(n+1)/sum(X1)-sum(X2)
> X1=c(X1,0)
> X2=c(X2,difference)
> n=n+1
> sum(X1*X2)/n
[1] 18.02172
> sum(X1)*sum(X2)/n^2
[1] 18.02172
On retrouve ici une corrélation parfaitement nulle,
> Y=c(Y,0)
>  base=data.frame(Y,X1,X2)
>  cor(base)
Y            X1            X2
Y  1.0000000  3.946091e-01  2.923412e-01
X1 0.3946091  1.000000e+00 -4.817519e-16
X2 0.2923412 -4.817519e-16  1.000000e+00
Si on regarde maintenant la régression multiple
>  R=lm(Y~X1+X2,data=base)
>  coefficients(R)
(Intercept)          X1          X2
-1.5028358   0.3589776   0.2531358
et les deux régressions simples,
>  R1=lm(Y~X1,data=base)
>  coefficients(R1)
(Intercept)          X1
0.02847362  0.35897764
>  R2=lm(Y~X2,data=base)
>  coefficients(R2)
(Intercept)          X2
-0.4334011   0.2531358
on obtient exactement les mêmes coefficient.
Moralité ? l’hypothèse d’orthogonalité entre les variables explicatives n’est pas à prendre à la légère: il faut une orthogonalité stricte pour que numériquement, la régression multiple soit équivalente aux régressions simples (comme le prédit la théorie). Mais de manière un peu surprenante, si on est très proche de l’indépendance, les valeurs numériques des estimateurs peut différer. Et ce, parfois de manière non négligeable, même si statistiquement, on valide l’hypothèse d’orthogonalité entre les variables explicatives.

Régression, variables explicatives et géométrie

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 http://freakonometrics.hypotheses.org/files/2016/05/FW15.gif sur http://freakonometrics.hypotheses.org/files/2016/05/FW16.gif et http://freakonometrics.hypotheses.org/files/2016/05/FW17.gif, deux ensembles (a priori disjoints) de variables explicatives,

http://freakonometrics.hypotheses.org/files/2016/05/FW02.gif

Les équations normales (associées au problème de minimisation de la somme des carrés des erreurs) sont

http://freakonometrics.hypotheses.org/files/2016/05/FW11.gif

de telle sorte qu’à l’optimum on peut relier les estimateurs des deux jeux de paramètres par

http://freakonometrics.hypotheses.org/files/2016/05/FW12.gif

La première partie correspond à la régression de http://freakonometrics.hypotheses.org/files/2016/05/FW15.gif sur http://freakonometrics.hypotheses.org/files/2016/05/FW13.gif, mais il reste un second terme dès lors que http://freakonometrics.hypotheses.org/files/2016/05/FW16.gif et http://freakonometrics.hypotheses.org/files/2016/05/FW17.gif ne sont pas orthogonales. Notons http://freakonometrics.hypotheses.org/files/2016/05/FW04.gifla matrice de projection (orthogonale) sur http://freakonometrics.hypotheses.org/files/2016/05/FW16.gif

http://freakonometrics.hypotheses.org/files/2016/05/FW06.gif

Alors on peut écrire

http://freakonometrics.hypotheses.org/files/2016/05/FW03.gif

En posant

http://freakonometrics.hypotheses.org/files/2016/05/FW07.gif et http://freakonometrics.hypotheses.org/files/2016/05/FW08.gif,

on retrouve un modèle linéaire classique (à condition de travailler sur la projection des variables sur le sous-espace engendré par http://freakonometrics.hypotheses.org/files/2016/05/FW16.gif , i.e. en transformant les variables),

http://freakonometrics.hypotheses.org/files/2016/05/FW10.gif

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

http://freakonometrics.hypotheses.org/files/2016/05/FW1.gif

Le graphique ci-dessous au cas non orthogonal,

http://freakonometrics.hypotheses.org/files/2016/05/FW2.gif

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.