2 La régression linéaire multiple

Exercice 1 (Question de cours) A, A, B, B, B, C.

Exercice 2 (Covariance de \(\hat\varepsilon\) et \(\hat Y\)) Nous allons montrer que, pour tout autre estimateur \(\tilde{\beta}\) de \(\beta\) linéaire et sans biais, \(\mathop{\mathrm{V}}(\tilde{\beta}) \geq \mathop{\mathrm{V}}(\hat \beta)\). Décomposons la variance de \(\tilde{\beta}\) \[\begin{eqnarray*} \mathop{\mathrm{V}}(\tilde{\beta}) = \mathop{\mathrm{V}}(\tilde{\beta} - \hat \beta+\hat \beta) =\mathop{\mathrm{V}}(\tilde{\beta} - \hat \beta)+\mathop{\mathrm{V}}(\hat \beta) - 2 \mathop{\mathrm{Cov}}(\tilde{\beta} - \hat \beta,\hat \beta). \end{eqnarray*}\] Les variances étant définies positives, si nous montrons que \(\mathop{\mathrm{Cov}}(\tilde{\beta}- \hat \beta,\hat \beta)=0\), nous aurons fini la démonstration.\ Puisque \(\tilde{\beta}\) est linéaire, \(\tilde{\beta} = A Y\). De plus, nous savons qu’il est sans biais, c’est-à-dire \(\mathbf E(\tilde{\beta}) = \beta\) pour tout \(\beta\), donc \(A X = I\). La covariance devient : \[\begin{eqnarray*} \mathop{\mathrm{Cov}}(\tilde{\beta} - \hat \beta,\hat \beta) &=& \mathop{\mathrm{Cov}}(A Y,(X'X)^{-1}X'Y) - \mathop{\mathrm{V}}(\hat \beta)\\ &=& \sigma^2 A X (X'X)^{-1} - \sigma^2 (X'X)^{-1}=0. \end{eqnarray*}\]

Exercice 3 (Théorème de Gauss Markov) Nous devons montrer que, parmi tous les estimateurs linéaires sans biais, l’estimateur de MC est celui qui a la plus petite variance. La linéarité de \(\hat \beta\) est évidente. Calculons sa variance : \[\begin{eqnarray*} \mathop{\mathrm{V}}(\hat \beta) = \mathop{\mathrm{V}}((X'X)^{-1}X'Y) = (X'X)^{-1}X'\mathop{\mathrm{V}}(Y)X(X'X)^{-1}=\sigma^2 (X'X)^{-1}. \end{eqnarray*}\] Nous allons montrer que, pour tout autre estimateur \(\tilde{\beta}\) de \(\beta\) linéaire et sans biais, \(\mathop{\mathrm{V}}(\tilde{\beta}) \geq \mathop{\mathrm{V}}(\hat \beta)\). Décomposons la variance de \(\tilde{\beta}\) \[\begin{eqnarray*} \mathop{\mathrm{V}}(\tilde{\beta}) = \mathop{\mathrm{V}}(\tilde{\beta} - \hat \beta+\hat \beta) =\mathop{\mathrm{V}}(\tilde{\beta} - \hat \beta)+\mathop{\mathrm{V}}(\hat \beta) - 2 \mathop{\mathrm{Cov}}(\tilde{\beta} - \hat \beta,\hat \beta). \end{eqnarray*}\] Les variances étant définies positives, si nous montrons que \(\mathop{\mathrm{Cov}}(\tilde{\beta}- \hat \beta,\hat \beta)=0\), nous aurons fini la démonstration.\ Puisque \(\tilde{\beta}\) est linéaire, \(\tilde{\beta} = A Y\). De plus, nous savons qu’il est sans biais, c’est-à-dire \(\mathbf E(\tilde{\beta}) = \beta\) pour tout \(\beta\), donc \(A X = I\). La covariance devient : \[\begin{eqnarray*} \mathop{\mathrm{Cov}}(\tilde{\beta} - \hat \beta,\hat \beta) &=& \mathop{\mathrm{Cov}}(A Y,(X'X)^{-1}X'Y) - \mathop{\mathrm{V}}(\hat \beta)\\ &=& \sigma^2 A X (X'X)^{-1} - \sigma^2 (X'X)^{-1}=0. \end{eqnarray*}\]

Exercice 4 (Représentation des variables) Nous représentons les données dans \(\mathbb R^2\) pour le premier jeu et dans \(\mathbb R^3\) pour le second.

Dans le premier modèle, nous projetons \(Y\) sur l’espace engendré par \(X\), soit la droite de vecteur directeur \(\overrightarrow{OX}\). Nous trouvons par le calcul \(\hat \beta = 1.47\), résultat que nous aurions pu trouver graphiquement car \(\overrightarrow{O \hat Y}= \hat \beta . \overrightarrow{OX}\).

Considérons \(\mathbb R^3\) muni de la base orthonormée \((\vec{i},\vec{j},\vec{k})\). Les vecteurs \(\overrightarrow{OX}\) et \(\overrightarrow{OZ}\) engendrent le même plan que celui engendré par \((\vec{i},\vec{j})\). La projection de \(Y\) sur ce plan donne \(\overrightarrow{O \hat Y}\). Il est quasiment impossible de trouver \(\hat \beta\) et \(\hat \gamma\) graphiquement mais nous trouvons par le calcul \(\hat \beta = -3.33\) et \(\hat \gamma =5\).

Exercice 5 (Modèles emboîtés) Nous obtenons \[\begin{eqnarray*} \hat Y_p = X \hat \beta \quad \hbox{et} \quad \hat Y_q= X_q \hat \gamma. \end{eqnarray*}\] Par définition du \(\mathbb R2\), il faut comparer la norme au carré des vecteurs \(\hat Y_p\) et \(\hat Y_q\). Notons les espaces engendrés par les colonnes de \(X_q\) et \(X\), \(\mathcal M_{X_q}\) et \(\mathcal M_{X}\), nous avons \(\mathcal M_{X_q} \subset \mathcal M_{X}\). Nous obtenons alors \[\begin{eqnarray*} \hat Y_p = P_{X_p}Y = (P_{X_q} + P_{X^{\perp}_q})P_{X_p}Y &=& P_{X_q}P_{X_p}Y + P_{X^{\perp}_q}P_{X_p}Y\\ &=& P_{X_q}Y + P_{X^{\perp}_q \cap X_p} Y\\ &=& \hat Y_q + P_{X^{\perp}_q \cap X_p} Y. \end{eqnarray*}\] En utilisant le théorème de Pythagore, nous avons \[\begin{eqnarray*} \| \hat Y_p \|^2 &=& \|\hat Y_q \|^2 + \| P_{X^{\perp}_q \cap X_p} Y \|^2 \geq \|\hat Y_q \|^2, \end{eqnarray*}\] d’où \[\begin{eqnarray*} \mathbb R2(p)=\frac{\| \hat Y_p \|^2}{\| Y \|^2} \geq \frac{\| \hat Y_q \|^2}{\| Y \|^2} =\mathbb R2(q). \end{eqnarray*}\]

En conclusion, lorsque les modèles sont emboîtés \(\mathcal M_{X_q} \subset \mathcal M_{X}\), le \(\mathbb R2\) du modèle le plus grand (ayant le plus de variables) sera toujours plus grand que le \(\mathbb R2\) du modèle le plus petit.

Exercice 6 La matrice \(X'X\) est symétrique, \(n\) vaut 30 et \(\bar x= \bar z=0\). Le coefficient de corrélation \[\begin{equation*} \rho_{x,z} = \frac{\sum_{i=1}^{30} (x_i -\bar x)(z_i - \bar z)} {\sqrt{\sum_{i=1}^{30} (x_i -\bar x)^2\sum_{i=1}^{30} (z_i - \bar z)^2}} =\frac{\sum_{i=1}^{30} x_i z_i} {\sqrt{\sum_{i=1}^{30} x_i^2 \sum_{i=1}^{30} z_i^2}} =\frac{7}{\sqrt{150}}=0.57. \end{equation*}\] Nous avons \[\begin{eqnarray*} y_i &=& -2 +x_i+z_i+\hat \varepsilon_i \end{eqnarray*}\] et la moyenne vaut alors \[\begin{eqnarray*} \bar y &=& -2 + \bar x +\bar z + \frac{1}{n}\sum_i \hat \varepsilon_i. \end{eqnarray*}\] La constante étant dans le modèle, la somme des résidus est nulle car le vecteur \(\hat \varepsilon\) est orthogonal au vecteur \(\mathbf{1}\). Nous obtenons donc que la moyenne de \(Y\) vaut 2 car \(\bar x=0\) et \(\bar z=0\). Nous obtenons en développant \[\begin{eqnarray*} \|\hat Y \|^2 &=& \sum_{i=1}^{30}(-2+x_i+2z_i)^2\\ &=& 4+10+60+14=88. \end{eqnarray*}\] Par le théorème de Pythagore, nous concluons que \[\begin{eqnarray*} \mathop{\mathrm{SCT}}=\mathop{\mathrm{SCE}}+\mathop{\mathrm{SCR}}=88+12=100. \end{eqnarray*}\]

Exercice 7 (Changement d’échelle des variables explicatives) Nous avons l’estimation sur le modèle avec les variables originales qui minimise \[\begin{align*} \mathop{\mathrm{MCO}}(\beta)&=\|Y - \sum_{j=1}^{p} X_j \beta_j \|^{2} \end{align*}\] Cette solution est notée \(\hat \beta\).

Nous avons l’estimation sur le modèle avec les variables prémultipliées par \(a_{j}\) (changement d’échelle) qui minimise

\[ \begin{align*} \tilde{\mathop{\mathrm{MCO}}}(\beta)&=\|Y - \sum_{j=1}^{p} \tilde X_j \beta_j \|^{2} = \|Y - \sum_{j=1}^{p} a_{j} X_j \beta_j \|^{2}\\ &= \|Y - \sum_{j=1}^{p} X_j \gamma_j \|^{2}=\mathop{\mathrm{MCO}}(\gamma), \end{align*} \]

en posant en dernière ligne \(\gamma_{j}=a_{j} \beta_j\). La solution de de \(\mathop{\mathrm{MCO}}(\gamma)\) (ou encore \(\mathop{\mathrm{MCO}}(\beta)\)) est \(\hat \beta\). La solution de \(\tilde{\mathop{\mathrm{MCO}}}(\beta)\) est alors donnée par \(\hat \beta_{j}=a_{j} \tilde \beta_j\).

Exercice 8 (Différence entre régression multiple et régressions simples)  

  1. Calculons l’estimateur des MCO noté traditionnellement \((X'X)^{-1}X'Y\) avec la matrice \(X\) qui possède ici deux colonnes (notées ici \(X\) et \(Z\)) et \(n\) lignes. On a donc \[\begin{align*} (X'X)&= \begin{pmatrix} \|X\|^{2} & <X,Z>\\ <X,Z> & \|Z\|^{2} \\ \end{pmatrix} \end{align*}\] Son déterminant est \(\Delta= \|X\|^{2}\|Z\|^{2} - 2 <X,Z>\) et son inverse est \[\begin{align*} \frac{1}{\Delta} \begin{pmatrix} \|Z\|^{2} & -<X,Z>\\ -<X,Z> & \|X\|^{2} \\ \end{pmatrix} \end{align*}\] Ensuite \(X'Y\) est simplement le vecteur colonne de coordonnées \(<X,Y>\) et \(<Z,Y>\). En rassemblant le tout nous avons \[\begin{align*} \hat \beta_{1}&=\frac{1}{\Delta}(\|Z\|^{2} <X,Y> - <X,Z><Z,Y>),\\ \hat \beta_{2}&=\frac{1}{\Delta}(\|X\|^{2} <Z,Y> - <X,Z><X,Y>). \end{align*}\] Si \(<X,Z>=0\) (les deux vecteurs sont orthogonaux) alors cette écriture se simplifie en \[ \hat \beta_{1}=\frac{<X,Y>}{\|X\|^{2}},\quad \hat \beta_{2}=\frac{<Z,Y>}{\|Z\|^{2}}. \]

  2. Calculons l’estimateur des MCO noté traditionnellement \((X'X)^{-1}X'Y\) avec la matrice \(X\) qui possède ici une colonne (notée ici \(X\)) et \(n\) lignes. On a donc \[\begin{align*} \hat \beta_{X} = \frac{<X,Y>}{\|X\|^{2}}. \end{align*}\] Passons maintenant à la matrice qui possède ici une colonne (notée ici \(Z\)) et \(n\) lignes. On a donc \[\begin{align*} \hat \beta_{Z} = \frac{<Z,Y>}{\|Z\|^{2}}. \end{align*}\]

  3. En général les coefficients des régressions simples ne sont pas ceux obtenus par régression multiple sauf si les variables sont orthogonales.

  4. Nous avons ici les résidus de la première régression qui sont \[\begin{align} \hat \varepsilon = Y - \hat \beta_{X} X. \end{align}\] La deuxième régression (sur les résidus) donne le coefficient \[\begin{align*} \hat \beta_{Z} = \frac{<Z,\hat \varepsilon>}{\|Z\|^{2}} = \hat \beta_{Z} - \hat \beta_{X}\frac{<Z,X>}{\|Z\|^{2}} \end{align*}\] La régression séquentielle donne des coefficients différents des régressions univariées ou bivariées sauf si les variables sont orthogonales. \end{enumerate}

Exercice 9 (TP : différence entre régression multiple et régressions simples)  

  1. Simulons 2 variables explicatives avec GNU-R pour \(n=100\) individus selon selon deux loi uniforme \([0,1]\) :

    n <- 100
set.seed(4321) # pour fixer les simulations
X1 <- runif(n)
X2 <- runif(n)

    Ensuite simulons \(Y\) selon le modèle avec \(\sigma=0.5\)

    sigma <- 0.5
set.seed(321) # pour fixer les simulations
Y <- 2 -3*X1 + 4*X2  + rnorm(n, sd=sigma)
don <- cbind.data.frame(Y,X1,X2)
head(don)
             Y         X1        X2
1 4.213477 0.33477802 0.5913398
2 1.602748 0.90913948 0.6715464
3 2.958646 0.41152969 0.5830570
4 3.215113 0.04384097 0.3516151
5 3.367004 0.76350011 0.9298713
6 3.555247 0.75043889 0.9181180

  2. Le graphique est obtenu avec

    rgl::plot3d(X1,X2,Y)

  3. On observe grâce au code ci-dessus que les points sont répartis autour d’un plan (d’équation \(z=2-3x + 4y\)).

  4. Effectuons la régression multiple et stockons \(\hat Y\) dans Yhat :

    (regmult <- lm(Y~1+X1+X2, data=don))
    
Call:
lm(formula = Y ~ 1 + X1 + X2, data = don)

Coefficients:
(Intercept)           X1           X2  
      2.080       -2.972        3.827

    Nous sommes assez proches des coefficients recherchés (2, -3 et 4)

  5. Effectuons les régression simples

    (regX1 <- lm(Y~1+X1, data=don))
    
Call:
lm(formula = Y ~ 1 + X1, data = don)

Coefficients:
(Intercept)           X1  
      3.803       -2.441

    Nous voyons que le paramètre ne correspond pas à celui de \(X_{1}\) dans la régression multiple.

    (regX2 <- lm(Y~1+X2, data=don))
    
Call:
lm(formula = Y ~ 1 + X2, data = don)

Coefficients:
(Intercept)           X2  
     0.7221       3.4424

    Nous voyons que le paramètre ne correspond pas à celui de \(X_{2}\) dans la régression multiple.

  6. Enchainons après la régression simple sur \(X_{1}\) la régression des résidus (ce que \(X_{1}\) n’a pas réussi à modéliser) sur \(X_{2}\):

    residX1 <- residuals(regX1)
don <- cbind(don, residX1)
regX2residX1 <- lm(residX1~1+X2, data=don)
regX2residX1
    
Call:
lm(formula = residX1 ~ 1 + X2, data = don)

Coefficients:
(Intercept)           X2  
     -1.965        3.758

    Là encore le coefficient sur \(X_{2}\) ne correspond pas à celui de \(X_{2}\) dans la régression multiple. :La prévision est différente

    (predict(regX2residX1) + predict(regX1))[1:5]
           1        2        3        4        5 
3.242781 2.142450 3.024335 3.051894 3.468733

  7. Changeons l’ordre

    residX2 <- residuals(regX2)
don <- cbind(don, residX2)
regX1residX2 <- lm(residX2~1+X1, data=don)
regX1residX2
    
Call:
lm(formula = residX2 ~ 1 + X1, data = don)

Coefficients:
(Intercept)           X1  
      1.531       -2.918

    On a là encore des différences

    (predict(regX1residX2) + predict(regX2))[1:5]
           1        2        3        4        5 
3.311895 1.911850 3.059400 3.335703 3.226125

  8. Envisageons la covariance empirique

    cov(X1,X2)    
    [1] 0.01050429

    et leur produit scalaire

    sum(X1*X2)
    [1] 28.47696

    Les deux variables sont loin d’être orthogonales. Centrons les

    X1c <- X1 - mean(X1)
X2c <- X2 - mean(X2)
sum(X1c*X2c) 
    [1] 1.039925

    Refaisons les régressions: la multiple

    donc <- cbind.data.frame(Y,X1c,X2c)
(regmultc <- lm(Y~1+X1c+X2c, data=donc))
    
Call:
lm(formula = Y ~ 1 + X1c + X2c, data = donc)

Coefficients:
(Intercept)          X1c          X2c  
      2.522       -2.972        3.827

    Nous retrouvons que seul la constante change puis l’enchainement

    (regX1c <- lm(Y~1+X1c, data=donc))
    
Call:
lm(formula = Y ~ 1 + X1c, data = donc)

Coefficients:
(Intercept)          X1c  
      2.522       -2.441

    Ici les variables \(X_{1c}\) et \(X_{2c}\) étant centrées elles sont orthogonales à la constante. Nous avons donc \[\begin{align} \Im(\mathbf{1}, X_{1c}, X_{2c}) = \Im(\mathbf{1}) \stackrel{\perp}{\oplus} \Im(X_{1c}, X_{2c}) \end{align}\] Le modèle précédent nous donne grâce à l’orthogonalité ci-dessus \(\hat Y_{1c} = P_{\mathbf{1}, X_{1c}}Y=P_{\mathbf{1}}Y + P_{X_{1c}}Y\). Le modèle complet donne de son côté \(\hat Y = P_{\mathbf{1}, X_{1c},X_{2c}}Y=P_{\mathbf{1}}Y + P_{X_{1c},X_{2c}}Y\). Nous retrouvons donc le coefficient constant qui est la coordonnée sur \(\mathbf{1}\) de \(P_{\mathbf{1}}Y\) dans les deux cas.

    Enchainons sur la régression sur résidus

    residX1c <- residuals(regX1c)
donc <- cbind(donc, residX1c)
(regX2residX1c <- lm(residX1c~1+X2c, data=don))
    
Call:
lm(formula = residX1c ~ 1 + X2c, data = don)

Coefficients:
(Intercept)          X2c  
  1.669e-16    3.758e+00

    Puisque le vecteur est \(Y- \hat Y_{1}=Y-P_{\mathbf{1}}Y-P_{X_{1c}}Y\) et que \(X_{2c}\perp\mathbf{1}\) on a la projection des résidus sur \(\Im(\mathbf{1}, X_{2c})\) qui vaut \[\begin{align*} P_{\mathbf{1}, X_{2c}}(Y- \hat Y_{1c})&= P_{\mathbf{1}}(Y- \hat Y_{1c}) + P_{X_{2c}}(Y- \hat Y_{1c})\\ &=P_{\mathbf{1}}Y - P_{\mathbf{1}}P_{\mathbf{1}, X_{1c}}Y + P_{X_{2c}}Y - P_{X_{2c}}P_{\mathbf{1}, X_{1c}}Y\\ &= P_{\mathbf{1}}Y - P_{\mathbf{1}}Y - 0 + P_{X_{2c}}Y - 0 - P_{X_{2c}}P_{X_{1c}}Y\\ &= P_{X_{2c}} Y - P_{X_{2c}}P_{X_{1c}}Y \end{align*}\] Cette dernière somme de projection est dans \(\Im(X_{1c}, X_{2c})\) qui est un sous-espace orthogonal à \(\Im(\mathbf{1})\), d’où le coefficient 0 pour la constante.

    Pour les prévisions rien ne change :

    predict(regmultc)[1:5]
           1        2        3        4        5 
3.348336 1.948483 3.088559 3.295485 3.369866

Exercice 10 (TP : régression multiple et code R)  

  1. ozone <- read.table("../donnees/ozone_complet.txt",sep=";",header=TRUE)
nomvar <- names(ozone)
  2. (ch <- paste(nomvar[-1],collapse = "+"))
    [1] "T6+T9+T12+T15+T18+Ne6+Ne9+Ne12+Ne15+Ne18+Vdir6+Vvit6+Vdir9+Vvit9+Vdir12+Vvit12+Vdir15+Vvit15+Vdir18+Vvit18+Vx+maxO3v"
  3. (ch2 <- paste("maxO3~",ch,sep=""))
    [1] "maxO3~T6+T9+T12+T15+T18+Ne6+Ne9+Ne12+Ne15+Ne18+Vdir6+Vvit6+Vdir9+Vvit9+Vdir12+Vvit12+Vdir15+Vvit15+Vdir18+Vvit18+Vx+maxO3v"
  4. form <- formula(ch2)
lm(form,data=ozone)
    
Call:
lm(formula = form, data = ozone)

Coefficients:
(Intercept)           T6           T9          T12          T15          T18  
  31.628418    -1.937194     0.121834     1.522143     0.609047     0.091871  
        Ne6          Ne9         Ne12         Ne15         Ne18        Vdir6  
   0.091969    -0.692891    -0.941925    -0.033847    -0.107068    -0.001344  
      Vvit6        Vdir9        Vvit9       Vdir12       Vvit12       Vdir15  
   1.586427    -0.009728    -1.163728    -0.007144     0.182232    -0.003468  
     Vvit15       Vdir18       Vvit18           Vx       maxO3v  
   0.246885     0.006251     0.560660     0.255474     0.465541

Exercice 11 (Régression orthogonale) Les vecteurs étant orthogonaux, nous avons \(\mathcal M_X = \mathcal M_U \stackrel{\perp}{\oplus} \mathcal M_V\). Nous pouvons alors écrire \[\begin{eqnarray*} \hat Y_X = P_X Y &=& (P_U + P_{U^{\perp}})P_X Y \\ &=& P_U P_X Y + P_{U^{\perp}}P_X Y = P_U Y + P_{U^{\perp}\cap X} Y \\ &=& \hat Y_U + \hat Y_V. \end{eqnarray*}\] La suite de l’exercice est identique. En conclusion, effectuer une régression multiple sur des variables orthogonales revient à effectuer \(p\) régressions simples.

Exercice 12 (Centrage, centrage-réduction et coefficient constant)  

  1. Comme la dernière colonne de \(X\), notée \(X_p\) vaut \(\mathbf{1}_n\) sa moyenne empirique vaut \(1\) et la variable centrée issue de \(X_p\) est donc \(X_p -1\times\mathbf{1}_n=\boldsymbol{0}_n\).

  2. Nous avons le modèle sur variable centrée \[ \begin{eqnarray*} \tilde Y&=&\tilde X\tilde \beta+ \varepsilon\\ Y-\bar Y\mathbf{1}_n&=&\sum_{j=1}^{p-1}{(X_j -\bar X_j\mathbf{1}_n)\tilde\beta_j}+\varepsilon\\ Y&=&\sum_{j=1}^{p-1}{\tilde\beta_j X_j}+ \Bigl(\bar Y -\sum_{j=1}^{p-1}{\bar X_j\tilde\beta_j}\Bigr)\mathbf{1}_n+\varepsilon. \end{eqnarray*} \] En identifiant cela donne \[ \begin{eqnarray} \beta_j&=&\tilde\beta_j,\ \forall j\in\{1,\dotsc,p-1\},\nonumber\\ \beta_p&=&\bar Y\mathbf{1}_n-\sum_{j=1}^{p-1}{\bar X_j\tilde\beta_j}. \end{eqnarray} \tag{1}\] Si l’on utilise des variables centrées dans le modèle de régression, on ne met pas de colonne \(\mathbf{1}\) (pas de coefficient constant - intercept). Les coefficients du modèle sur les variables originales sont égaux à ceux sur les variables centrées et le coefficient constant est donné par la formule (équation 1).

  3. Maintenant les variables explicatives sont centrées et réduites : \[\begin{eqnarray*} \tilde Y&=&\tilde X\tilde \beta+ \varepsilon\\ Y-\bar Y\mathbf{1}_n&=&\sum_{j=1}^{p-1}{\frac{(X_j -\bar X_j\mathbf{1}_n)}{\hat\sigma_{X_j}}\tilde\beta_j}+\varepsilon\\ Y&=&\sum_{j=1}^{p-1}{\frac{\tilde\beta_j}{\hat\sigma_{X_j}}X_j}+ \Bigl(\bar Y-\sum_{j=1}^{p-1}{\bar X_j\frac{\tilde\beta_j}{\hat\sigma_{X_j}}}\Bigr)\mathbf{1}_n+\varepsilon. \end{eqnarray*}\] En identifiant cela donne \[\begin{eqnarray*} \beta_j&=&\frac{\tilde\beta_j}{\hat\sigma_{X_j}},\ \forall j\in\{1,\dotsc,p-1\},\\ \beta_p&=&\bar Y\mathbf{1}_n-\sum_{j=1}^{p-1}{\bar X_j\frac{\tilde\beta_j}{\hat\sigma_{X_j}}}. \end{eqnarray*}\] Nous obtenons ici que les coefficients du modèle sur les variables originales sont égaux à ceux sur les variables centrées-réduites divisés par l’écart-type empirique des variables explicatives. Plus la variable explicative \(X_j\) est dispersée, plus son coefficient \(\beta_j\) sera réduit par rapport à \(\tilde\beta_j\). Le coefficient constant est donné par la formule ci-dessus.

  4. La variable à expliquer \(Y\) est elle aussi centrée-réduite~: \[\begin{eqnarray*} \tilde Y&=&\tilde X\tilde \beta+ \tilde\varepsilon\\ \frac{Y-\bar Y\mathbf{1}_n}{\hat\sigma_Y}&=&\sum_{j=1}^{p-1}{\frac{(X_j -\bar X_j\mathbf{1}_n)}{\hat\sigma_{X_j}}\tilde\beta_j}+\tilde\varepsilon\\ Y&=&\hat\sigma_Y\sum_{j=1}^{p-1}{\frac{\tilde\beta_j}{\hat\sigma_{X_j}}X_j}+ \Bigl(\bar Y-\hat\sigma_Y\sum_{j=1}^{p-1}{\bar X_j\frac{\tilde\beta_j}{\hat\sigma_{X_j}}}\Bigr)\mathbf{1}_n+\hat\sigma_Y\tilde\varepsilon. \end{eqnarray*}\] En identifiant cela donne \[\begin{eqnarray*} \beta_j&=&\hat\sigma_Y\frac{\tilde\beta_j}{\hat\sigma_{X_j}},\ \forall j\in\{1,\dotsc,p-1\},\\ \beta_p&=&\bar Y\mathbf{1}_n-\hat\sigma_Y\sum_{j=1}^{p-1}{\bar X_j\frac{\tilde\beta_j}{\hat\sigma_{X_j}}},\\ \varepsilon&=&\hat\sigma_Y\tilde\varepsilon. \end{eqnarray*}\] L’écart-type empirique de \(Y\) entre en jeu et nous constatons que les résidus du modèle “centré-réduit” sont égaux à ceux initiaux divisés par l’écart-type empirique de \(Y\).

  5. Simulons 3 variables explicatives avec GNU-R pour \(n=100\) individus selon une loi uniforme \([0,1]\) et équirépartie sur \([1;10]\):

    n <- 100
set.seed(1234) # pour fixer les simulations
X1 <- runif(n)
X2 <- seq(1, 10, length=100)
X3 <- rep(1,n)
X <- cbind(X1, X2, X3)
head(X)
                X1       X2 X3
[1,] 0.1137034 1.000000  1
[2,] 0.6222994 1.090909  1
[3,] 0.6092747 1.181818  1
[4,] 0.6233794 1.272727  1
[5,] 0.8609154 1.363636  1
[6,] 0.6403106 1.454545  1

    Simulons maintenant \(Y\) selon le modèle en fixant par exemple \(\beta=(1, 2, 0.1, -4)'\) et \(\sigma=0.1\)

    beta <- c(1, 0.1, -4)
sigma <- 0.1
set.seed(12345) # pour fixer les simulations
Y <- X%*%beta + rnorm(n, sd=sigma)
Y[1:5]
    [1] -3.727744 -3.197663 -3.283474 -3.294698 -2.942132

  6. Centrons maintenant les variables avec 

    Ytilde <- scale(Y, scale = FALSE)
Ytilde[1:5]
    [1] -0.73976069 -0.20968007 -0.29549076 -0.30671453  0.04585078
    Xtilde <- scale(X, scale = FALSE)
head(Xtilde)
                 X1        X2 X3
[1,] -0.3237939 -4.500000  0
[2,]  0.1848021 -4.409091  0
[3,]  0.1717775 -4.318182  0
[4,]  0.1858822 -4.227273  0
[5,]  0.4234181 -4.136364  0
[6,]  0.2028133 -4.045455  0

    Et éliminons la dernière colonne de \(\tilde X\)

    Xtilde <- Xtilde[, -3]

  7. Effectuons les deux régressions en constituant deux data-frames. Nous n’ajoutons pas la constante dans les modèles (-1) car elle est déjà dans la dernière colonne de \(X\) pour le premier modèle et on ne la souhaite pas pour le second.

    don <- cbind.data.frame(Y,X)
dontilde <- cbind.data.frame(Ytilde,Xtilde)
names(dontilde) <- c("Ytilde", "X1tilde", "X2tilde")
mod <- lm(Y~-1+X1+X2+X3, data=don)
modtilde <- lm(Ytilde~-1+X1tilde+X2tilde, data=dontilde)

    On retrouve bien que les coefficients des variables centrées ou non sont identiques

    coef(mod)
            X1         X2         X3 
 1.0225207  0.1030809 -4.0022782

    et que l’intercept est simplement est bien celui donné dans la formule (équation 1)

    mean(don$Y)- sum( coef(modtilde)*colMeans(X[,1:2]))
    [1] -4.002278

  8. Centrons et réduisons les variables explicatives (les 2 premières colonnes de \(X\))

    Xtilde <- scale(X[, 1:2], scale = TRUE)
head(Xtilde)
                 X1        X2
[1,] -1.1617874 -1.706220
[2,]  0.6630787 -1.671751
[3,]  0.6163456 -1.637282
[4,]  0.6669539 -1.602813
[5,]  1.5192439 -1.568344
[6,]  0.7277037 -1.533875

    Reconstuisons le data-frame qui a changé et le modèle correspondant

    dontilde <- cbind.data.frame(Ytilde,Xtilde)
names(dontilde) <- c("Ytilde", "X1tilde", "X2tilde")
modtilde <- lm(Ytilde~-1+X1tilde+X2tilde, data=dontilde)

    Retrouvons les coefficients de départ pour les variables explicatives et pour la constante

    mean(don$Y)- sum( coef(modtilde)/apply(X[,-3], 2, sd) *
                  colMeans(X[,1:2]) )
    [1] -4.002278

  9. Centrons et réduisons les variables explicatives et la variable à expliquer

    Ytilde <- scale(Y, scale = TRUE)
Ytilde[1:5]
    [1] -1.8943177 -0.5369313 -0.7566682 -0.7854091  0.1174109

    Reconstuisons le data-frame qui a changé et le modèle correspondant

    dontilde <- cbind.data.frame(Ytilde, Xtilde)
names(dontilde) <- c("Ytilde", "X1tilde", "X2tilde")
modtilde <- lm(Ytilde~-1+X1tilde+X2tilde, data=dontilde)

    Retrouvons les coefficients de départ pour les variables explicatives

    coef(modtilde)/apply(X[,-3],2,sd)*sd(Y)
      X1tilde   X2tilde 
1.0225207 0.1030809

    et pour la constante

    mean(don$Y)- sum( coef(modtilde)/apply(X[,-3], 2, sd) *
                  sd(Y)*colMeans(X[, 1:2]) )
    [1] -4.002278

Exercice 13 (Moindres carrés contraints)  

  1. L’estimateur des MC vaut \[\begin{eqnarray*} \hat \beta = (X'X)^{-1}X'Y, \end{eqnarray*}\]

  2. Calculons maintenant l’estimateur contraint. Nous pouvons procéder de deux manières différentes.

    La première consiste à écrire le lagrangien \[\begin{eqnarray*} \mathcal{L} = S(\beta) - \lambda'(R\beta-r). \end{eqnarray*}\] Les conditions de Lagrange permettent d’obtenir un minimum \[\begin{eqnarray*} \left\{ \begin{array}{l} \displaystyle\frac{\partial \mathcal{L}}{\partial \beta} = -2X'Y+2X'X\hat{\beta}_c- R'\hat{\lambda}=0,\\ \displaystyle\frac{\partial \mathcal{L}}{\partial \lambda} = R\hat{\beta}_c-r=0, \end{array} \right. \end{eqnarray*}\] Multiplions à gauche la première égalité par \(R(X'X)^{-1}\), nous obtenons \[\begin{eqnarray*} -2 R(X'X)^{-1}X'Y+2R(X'X)^{-1}X'X \hat{\beta}_c-R(X'X)^{-1}R'\hat{\lambda}&=&0\\ -2 R(X'X)^{-1}X'Y+2R\hat{\beta}_c-R(X'X)^{-1}R'\hat{\lambda}&=&0\\ -2 R(X'X)^{-1}X'Y+2r-R(X'X)^{-1}R'\hat{\lambda}&=&0. \end{eqnarray*}\] Nous obtenons alors pour \(\hat \lambda\) \[\begin{eqnarray*} \hat \lambda = 2 \left[R(X'X)^{-1}R'\right]^{-1}\left[r-R(X'X)^{-1}X'Y\right]. \end{eqnarray*}\] Remplaçons ensuite \(\hat \lambda\) \[\begin{eqnarray*} -2X'Y+2X'X\hat{\beta}_c-R'\hat{\lambda}&=&0\\ -2X'Y+2X'X\hat{\beta}_c-2R'\left[R(X'X)^{-1}R'\right]^{-1}\left[r-R(X'X)^{-1}X'Y\right]&=& 0, \end{eqnarray*}\] d’où nous calculons \(\hat \beta_c\) \[\begin{eqnarray*} \hat \beta_c &=& (X'X)^{-1}X'Y+(X'X)^{-1}R'\left[R(X'X)^{-1}R'\right]^{-1} (r-R\hat \beta)\\ &=& \hat \beta + (X'X)^{-1}R'\left[R(X'X)^{-1}R'\right]^{-1}(r-R\hat \beta). \end{eqnarray*}\] La fonction \(S(\beta)\) à minimiser est une fonction convexe sur un ensemble convexe (contraintes linéaires), le minimum est donc unique.

    Une autre façon de procéder consiste à utiliser les projecteurs. Supposons pour commencer que \(r=0\), la contrainte vaut donc \(R\beta=0\). Calculons analytiquement le projecteur orthogonal sur \(\mathcal M_0\). Rappelons que \(\dim(\mathcal M_0)=p-q\), nous avons de plus \[\begin{eqnarray*} R \beta &=& 0 \quad \quad \Leftrightarrow \quad \beta \in Ker(R)\\ R (X'X)^{-1}X'X \beta &=& 0\\ U' X \beta &=& 0\quad \quad \hbox{où} \quad \quad U = X (X'X)^{-1}R'. \end{eqnarray*}\] Nous avons donc que \(\forall \beta \in \ker(R)\), \(U' X \beta = 0\), c’est-à-dire que \(\mathcal M_U\), l’espace engendré par les colonnes de \(U\), est orthogonal à l’espace engendré par \(X\beta\), \(\forall \beta \in \ker(R)\). Nous avons donc que \(\mathcal M_U \perp \mathcal M_0\). Comme \(U=X[(X'X)^{-1}R']\), \(\mathcal M_U \subset \mathcal M_X\). En résumé, nous avons \[\begin{eqnarray*} \mathcal M_U \subset \mathcal M_X \quad \hbox{et} \quad \mathcal M_U \perp \mathcal M_0 \quad \hbox{donc} \quad \mathcal M_U \subset (\mathcal M_X \cap \mathcal M_0^{\perp}). \end{eqnarray*}\] Afin de montrer que les colonnes de \(U\) engendrent \(\mathcal M_X \cap \mathcal M_0^{\perp}\), il faut démontrer que la dimension des deux sous-espaces est égale. Or le rang de \(U\) vaut \(q\) (\(R'\) est de rang \(q\), \((X'X)^{-1}\) est de rang \(p\) et \(X\) est de rang \(p\)) donc la dimension de \(\mathcal M_U\) vaut \(q\). De plus, nous avons vu que \[\begin{eqnarray*} \mathcal M_X = \mathcal M_0 \stackrel{\perp}{\oplus}\left(\mathcal M_0^{\perp} \cap \mathcal M_X \right) \end{eqnarray*}\] et donc, en passant aux dimensions des sous-espaces, nous en déduisons que \(\dim(\mathcal M_0^{\perp} \cap \mathcal M_X )=q\). Nous venons de démontrer que \[\begin{eqnarray*} \mathcal M_U = \mathcal M_X \cap \mathcal M_0^{\perp}. \end{eqnarray*}\] Le projecteur orthogonal sur \(\mathcal M_U=\mathcal M_X \cap \mathcal M_0^{\perp}\) s’écrit \[\begin{eqnarray*} P_{U} = U (U'U)^{-1} U'= X (X'X)^{-1} R' [R(X'X)^{-1}R']^{-1}R(X'X)^{-1}X'. \end{eqnarray*}\] Nous avons alors \[\begin{eqnarray*} \hat Y - \hat Y_0 &=& P_U Y\\ X \hat \beta - X \hat \beta_0 &=& X (X'X)^{-1} R' [R(X'X)^{-1}R']^{-1}R(X'X)^{-1}X'Y\\ &=& X (X'X)^{-1} R [R(X'X)^{-1}R']^{-1}R \hat \beta. \end{eqnarray*}\] Cela donne \[\begin{eqnarray*} \hat \beta_0 = \hat \beta - (X'X)^{-1} R [R(X'X)^{-1}R']^{-1}R \hat \beta. \end{eqnarray*}\] Si maintenant \(r\neq 0\), nous avons alors un sous-espace affine défini par \(\{\beta\in \mathbb R^p : R\beta=r\}\) dans lequel nous cherchons une solution qui minimise les moindres carrés. Un sous-espace affine peut être défini de manière équivalente par un point particulier \(\beta_p \in \mathbb R^p\) tel que \(R\beta_p=r\) et le sous-espace vectoriel associé \(\mathcal M_0^v=\{\beta\in \mathbb R^p : R\beta=0\}\). Les points du sous-espace affine sont alors \(\{\beta_0 \in \mathbb R^p : \beta_0=\beta_p+\beta_0^v, \beta_0^v \in \mathcal M_0^v \quad et \quad \beta_p : R\beta_p=r\}\). La solution qui minimise les moindres carrés, notée \(\hat \beta_0\), est élément de ce sous-espace affine et est définie par \(\hat \beta_0=\beta_p+\hat \beta_0^v\) où \[\begin{eqnarray*} \hat \beta_0^v = \hat \beta - (X'X)^{-1}R'[R(X'X)^{-1}R']^{-1}R\hat \beta. \end{eqnarray*}\] Nous savons que \(R\beta_p=r\) donc \[\begin{eqnarray*} R\beta_p = [R(X'X)^{-1}R'][R(X'X)^{-1}R']^{-1}r \end{eqnarray*}\] donc une solution particulière est \(\beta_p = (X'X)^{-1}R'[R(X'X)^{-1}R']^{-1}r\). La solution \(\hat \beta_0\) qui minimise les moindres carrés sous la contrainte \(R\beta=r\) est alors \[ \begin{align} \hat \beta_0 &= \beta_p+\hat \beta_0^v\\\nonumber &=(X'X)^{-1}R'[R(X'X)^{-1}R']^{-1}r + \hat \beta - (X'X)^{-1}R'[R(X'X)^{-1}R']^{-1}R\hat \beta\\\nonumber &=\hat \beta + (X'X)^{-1}R'[R(X'X)^{-1}R']^{-1}(r-R\hat \beta). \end{align} \tag{2}\]

  3. Calculons l’EQM de \(\hat \beta\) qui vaut selon la formule classique: \[\begin{align*} \mathop{\mathrm{EQM}}&= (\mathbf E(\hat \beta) - \beta) (\mathbf E(\hat \beta) - \beta)' + \mathop{\mathrm{V}}(\hat\beta) \end{align*}\] avec \(\mathbf E(\hat \beta)=\beta\) (sous l’hypothèse que \(Y\) est généré par le modèle de régression) et\(\mathop{\mathrm{V}}(\hat\beta) =\sigma^{2} (X'X)^{-1}\).

    Pour l’EQM de \(\hat \beta_0\) qui vaut selon la formule classique: \[\begin{align*} \mathop{\mathrm{EQM}}&= (\mathbf E(\hat \beta_0) - \beta_0) (\mathbf E(\hat \beta_0) - \beta_0)' + \mathop{\mathrm{V}}(\hat\beta_0) \end{align*}\] calculons d’abord \(\mathbf E(\hat \beta_0)\) en utilisant l’équation (équation 2): \[\begin{align*} \mathbf E(\hat \beta_0)&=\mathbf E(\hat \beta) + \mathbf E[(X'X)^{-1}R'[R(X'X)^{-1}R']^{-1}(r-R\hat \beta)]\\ &=\beta + (X'X)^{-1}R'[R(X'X)^{-1}R']^{-1}(r-R\beta) \end{align*}\] puisque \(\hat \beta\) est sans biais. Si nous supposons que le modèle satisfait la contrainte (\(R\beta=r\)) alors là encore le biais est nul.

    Calculons maintemant la variance: \[\begin{align*} \mathop{\mathrm{V}}(\hat\beta_0) &= \mathop{\mathrm{V}}[\hat \beta + (X'X)^{-1}R'[R(X'X)^{-1}R']^{-1}(r-R\hat \beta)]\\ &=\mathop{\mathrm{V}}(\hat\beta) + \mathop{\mathrm{V}}[(X'X)^{-1}R'[R(X'X)^{-1}R']^{-1}(r-R\hat \beta)] +\\ &\quad \quad 2\mathop{\mathrm{Cov}}[\hat\beta ; (X'X)^{-1}R'[R(X'X)^{-1}R']^{-1}(r-R\hat \beta)]\\ &= \sigma^{2}V_{X}^{-1} + (\mathrm{II}) + (\mathrm{III}) \end{align*}\] Intéressons nous à la seconde partie \((\mathrm{II})\). Comme \(X\) est déterministe ainsi que \(R\) et \(r\) on a en posant pour alléger \(V_{X}=(X'X)\) (matrice symétrique) \[\begin{align*} (\mathrm{II}) &= V_{X}^{-1}R'[RV_{X}^{-1}R']^{-1} \mathop{\mathrm{V}}(r-R\hat \beta) [RV_{X}^{-1}R']^{-1} R V_{X}^{-1}\\ &= V_{X}^{-1}R'[RV_{X}^{-1}R']^{-1} \mathop{\mathrm{V}}(R\hat \beta)[RV_{X}^{-1}R']^{-1} R V_{X}^{-1}\\ &= V_{X}^{-1}R'[RV_{X}^{-1}R']^{-1} R\sigma^{2}V_{X}^{-1}R'[RV_{X}^{-1}R']^{-1} R V_{X}^{-1}\\ &= \sigma^{2}V_{X}^{-1}R'[RV_{X}^{-1}R']^{-1} R V_{X}^{-1} \end{align*}\] La troisième partie est \[\begin{align*} (\mathrm{III}) &= 2V_{X}^{-1}R'[RV_{X}^{-1}R']^{-1}\mathop{\mathrm{Cov}}[\hat\beta ; (r-R\hat \beta)]\\ &=-2V_{X}^{-1}R'[RV_{X}^{-1}R']^{-1}R\mathop{\mathrm{Cov}}[\hat\beta ;\hat\beta]\\ &= -2V_{X}^{-1}R'[RV_{X}^{-1}R']^{-1}R\sigma^{2}V_{X}^{-1} \end{align*}\] Ce qui donne au final \[\begin{align*} \mathop{\mathrm{V}}(\hat\beta_0) &=\sigma^{2}(X'X)^{-1} - \sigma^{2}(X'X)^{-1}R'[R(X'X)^{-1}R']^{-1}R(X'X)^{-1} \end{align*}\] L’écart entre les 2 variances est donc de \(\sigma^{2}(X'X)^{-1}R'[R(X'X)^{-1}R']^{-1}R(X'X)^{-1}\) qui est une matrice de la forme \(A'A\) donc semi-définie positive.