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
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.
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
> 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
(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).
> 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 dans le modèle multiple est légèrement différent de
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
et à droite celle de
)
 |
 |
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 et à droite les
, 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 -value: plus la
-value est faible (rejet de l’hypothèse d’indépendance) plus la couleur sera vive (et rouge), et plus la
-value est forte (acceptation de l’hypothèse d’indépendance) plus la couleur sera terne (et jaune, voire blanche),
 |
 |
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 -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
> 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
> 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