2 La régression linéaire multiple

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

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

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{β}$ est évidente. Calculons sa variance : $\begin{array}{r} V (\hat{β}) = V ((X^{'} X)^{- 1} X^{'} Y) = (X^{'} X)^{- 1} X^{'} V (Y) X (X^{'} X)^{- 1} = σ^{2} (X^{'} X)^{- 1} . \end{array}$ Nous allons montrer que, pour tout autre estimateur $\tilde{β}$ de $β$ linéaire et sans biais, $V (\tilde{β}) \geq V (\hat{β})$ . Décomposons la variance de $\tilde{β}$ $\begin{array}{r} V (\tilde{β}) = V (\tilde{β} - \hat{β} + \hat{β}) = V (\tilde{β} - \hat{β}) + V (\hat{β}) - 2 Cov (\tilde{β} - \hat{β}, \hat{β}) . \end{array}$ Les variances étant définies positives, si nous montrons que $Cov (\tilde{β} - \hat{β}, \hat{β}) = 0$ , nous aurons fini la démonstration.\ Puisque $\tilde{β}$ est linéaire, $\tilde{β} = A Y$ . De plus, nous savons qu’il est sans biais, c’est-à-dire $E (\tilde{β}) = β$ pour tout $β$ , donc $A X = I$ . La covariance devient : $\begin{array}{rcl} Cov (\tilde{β} - \hat{β}, \hat{β}) & = & Cov (A Y, (X^{'} X)^{- 1} X^{'} Y) - V (\hat{β}) \\ = & σ^{2} A X (X^{'} X)^{- 1} - σ^{2} (X^{'} X)^{- 1} = 0. \end{array}$

Exercice 4 (Représentation des variables) Nous représentons les données dans $R^{2}$ pour le premier jeu et dans $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 $\vec{O X}$ . Nous trouvons par le calcul $\hat{β} = 1.47$ , résultat que nous aurions pu trouver graphiquement car $\vec{O \hat{Y}} = \hat{β} . \vec{O X}$ .

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

Exercice 5 (Modèles emboîtés) Nous obtenons $\begin{array}{r} {\hat{Y}}_{p} = X \hat{β} et {\hat{Y}}_{q} = X_{q} \hat{γ} . \end{array}$ Par définition du $R 2$ , 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$ , $M_{X_{q}}$ et $M_{X}$ , nous avons $M_{X_{q}} \subset M_{X}$ . Nous obtenons alors $\begin{array}{rcl} {\hat{Y}}_{p} = P_{X_{p}} Y = (P_{X_{q}} + P_{X_{q}^{⊥}}) P_{X_{p}} Y & = & P_{X_{q}} P_{X_{p}} Y + P_{X_{q}^{⊥}} P_{X_{p}} Y \\ = & P_{X_{q}} Y + P_{X_{q}^{⊥} \cap X_{p}} Y \\ = & {\hat{Y}}_{q} + P_{X_{q}^{⊥} \cap X_{p}} Y . \end{array}$ En utilisant le théorème de Pythagore, nous avons $\begin{array}{rcl} ∥ {\hat{Y}}_{p} ∥^{2} & = & ∥ {\hat{Y}}_{q} ∥^{2} + ∥ P_{X_{q}^{⊥} \cap X_{p}} Y ∥^{2} \geq ∥ {\hat{Y}}_{q} ∥^{2}, \end{array}$ d’où $\begin{array}{r} R 2 (p) = \frac{∥ {\hat{Y}}_{p} ∥^{2}}{∥ Y ∥^{2}} \geq \frac{∥ {\hat{Y}}_{q} ∥^{2}}{∥ Y ∥^{2}} = R 2 (q) . \end{array}$

En conclusion, lorsque les modèles sont emboîtés $M_{X_{q}} \subset M_{X}$ , le $R 2$ du modèle le plus grand (ayant le plus de variables) sera toujours plus grand que le $R 2$ 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 $ρ_{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 .$ Nous avons $\begin{array}{rcl} y_{i} & = & - 2 + x_{i} + z_{i} + {\hat{ε}}_{i} \end{array}$ et la moyenne vaut alors $\begin{array}{rcl} \bar{y} & = & - 2 + \bar{x} + \bar{z} + \frac{1}{n} \sum_{i} {\hat{ε}}_{i} . \end{array}$ La constante étant dans le modèle, la somme des résidus est nulle car le vecteur $\hat{ε}$ est orthogonal au vecteur $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{array}{rcl} ∥ \hat{Y} ∥^{2} & = & \sum_{i = 1}^{30} (- 2 + x_{i} + 2 z_{i})^{2} \\ = & 4 + 10 + 60 + 14 = 88. \end{array}$ Par le théorème de Pythagore, nous concluons que $\begin{array}{r} SCT = SCE + SCR = 88 + 12 = 100. \end{array}$

Exercice 7 (Changement d’échelle des variables explicatives) Nous avons l’estimation sur le modèle avec les variables originales qui minimise $\begin{aligned} MCO (β) & = ∥ Y - \sum_{j = 1}^{p} X_{j} β_{j} ∥^{2} \end{aligned}$ Cette solution est notée $\hat{β}$ .

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

$\begin{aligned} \tilde{MCO} (β) & = ∥ Y - \sum_{j = 1}^{p} {\tilde{X}}_{j} β_{j} ∥^{2} = ∥ Y - \sum_{j = 1}^{p} a_{j} X_{j} β_{j} ∥^{2} \\ = ∥ Y - \sum_{j = 1}^{p} X_{j} γ_{j} ∥^{2} = MCO (γ), \end{aligned}$

en posant en dernière ligne $γ_{j} = a_{j} β_{j}$ . La solution de de $MCO (γ)$ (ou encore $MCO (β)$ ) est $\hat{β}$ . La solution de $\tilde{MCO} (β)$ est alors donnée par ${\hat{β}}_{j} = a_{j} {\tilde{β}}_{j}$ .

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

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{aligned} (X^{'} X) & = (\begin{array}{c} ∥ X ∥^{2} & < X, Z > \\ < X, Z > & ∥ Z ∥^{2} \end{array}) \end{aligned}$ Son déterminant est $Δ = ∥ X ∥^{2} ∥ Z ∥^{2} - 2 < X, Z >$ et son inverse est $\begin{array}{r} \frac{1}{Δ} (\begin{array}{c} ∥ Z ∥^{2} & - < X, Z > \\ - < X, Z > & ∥ X ∥^{2} \end{array}) \end{array}$ Ensuite $X^{'} Y$ est simplement le vecteur colonne de coordonnées $< X, Y >$ et $< Z, Y >$ . En rassemblant le tout nous avons $\begin{aligned} {\hat{β}}_{1} & = \frac{1}{Δ} (∥ Z ∥^{2} < X, Y > - < X, Z >< Z, Y >), \\ {\hat{β}}_{2} & = \frac{1}{Δ} (∥ X ∥^{2} < Z, Y > - < X, Z >< X, Y >) . \end{aligned}$ Si $< X, Z >= 0$ (les deux vecteurs sont orthogonaux) alors cette écriture se simplifie en ${\hat{β}}_{1} = \frac{< X, Y >}{∥ X ∥^{2}}, {\hat{β}}_{2} = \frac{< Z, Y >}{∥ Z ∥^{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{array}{r} {\hat{β}}_{X} = \frac{< X, Y >}{∥ X ∥^{2}} . \end{array}$ Passons maintenant à la matrice qui possède ici une colonne (notée ici $Z$ ) et $n$ lignes. On a donc $\begin{array}{r} {\hat{β}}_{Z} = \frac{< Z, Y >}{∥ Z ∥^{2}} . \end{array}$
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.
Nous avons ici les résidus de la première régression qui sont $\begin{array}{r} \hat{ε} = Y - {\hat{β}}_{X} X . \end{array}$ La deuxième régression (sur les résidus) donne le coefficient $\begin{array}{r} {\hat{β}}_{Z} = \frac{< Z, \hat{ε} >}{∥ Z ∥^{2}} = {\hat{β}}_{Z} - {\hat{β}}_{X} \frac{< Z, X >}{∥ Z ∥^{2}} \end{array}$ 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)

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 $σ = 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

Le graphique est obtenu avec
```
rgl::plot3d(X1,X2,Y)
```
On observe grâce au code ci-dessus que les points sont répartis autour d’un plan (d’équation $z = 2 - 3 x + 4 y$ ).

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)

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.

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

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

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_{1 c}$ et $X_{2 c}$ étant centrées elles sont orthogonales à la constante. Nous avons donc $\begin{array}{r} ℑ (1, X_{1 c}, X_{2 c}) = ℑ (1) \overset{⊥}{\oplus} ℑ (X_{1 c}, X_{2 c}) \end{array}$ Le modèle précédent nous donne grâce à l’orthogonalité ci-dessus ${\hat{Y}}_{1 c} = P_{1, X_{1 c}} Y = P_{1} Y + P_{X_{1 c}} Y$ . Le modèle complet donne de son côté $\hat{Y} = P_{1, X_{1 c}, X_{2 c}} Y = P_{1} Y + P_{X_{1 c}, X_{2 c}} Y$ . Nous retrouvons donc le coefficient constant qui est la coordonnée sur $1$ de $P_{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_{1} Y - P_{X_{1 c}} Y$ et que $X_{2 c} ⊥ 1$ on a la projection des résidus sur $ℑ (1, X_{2 c})$ qui vaut $\begin{aligned} P_{1, X_{2 c}} (Y - {\hat{Y}}_{1 c}) & = P_{1} (Y - {\hat{Y}}_{1 c}) + P_{X_{2 c}} (Y - {\hat{Y}}_{1 c}) \\ = P_{1} Y - P_{1} P_{1, X_{1 c}} Y + P_{X_{2 c}} Y - P_{X_{2 c}} P_{1, X_{1 c}} Y \\ = P_{1} Y - P_{1} Y - 0 + P_{X_{2 c}} Y - 0 - P_{X_{2 c}} P_{X_{1 c}} Y \\ = P_{X_{2 c}} Y - P_{X_{2 c}} P_{X_{1 c}} Y \end{aligned}$ Cette dernière somme de projection est dans $ℑ (X_{1 c}, X_{2 c})$ qui est un sous-espace orthogonal à $ℑ (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)

ozone <- read.table("../donnees/ozone_complet.txt",sep=";",header=TRUE)
nomvar <- names(ozone)

(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"

(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"

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 $M_{X} = M_{U} \overset{⊥}{\oplus} M_{V}$ . Nous pouvons alors écrire $\begin{array}{rcl} {\hat{Y}}_{X} = P_{X} Y & = & (P_{U} + P_{U^{⊥}}) P_{X} Y \\ = & P_{U} P_{X} Y + P_{U^{⊥}} P_{X} Y = P_{U} Y + P_{U^{⊥} \cap X} Y \\ = & {\hat{Y}}_{U} + {\hat{Y}}_{V} . \end{array}$ 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)

Comme la dernière colonne de $X$ , notée $X_{p}$ vaut $1_{n}$ sa moyenne empirique vaut $1$ et la variable centrée issue de $X_{p}$ est donc $X_{p} - 1 \times 1_{n} = 0_{n}$ .
Nous avons le modèle sur variable centrée $\begin{array}{rcl} \tilde{Y} & = & \tilde{X} \tilde{β} + ε \\ Y - \bar{Y} 1_{n} & = & \sum_{j = 1}^{p - 1} (X_{j} - {\bar{X}}_{j} 1_{n}) {\tilde{β}}_{j} + ε \\ Y & = & \sum_{j = 1}^{p - 1} {\tilde{β}}_{j} X_{j} + (\bar{Y} - \sum_{j = 1}^{p - 1} {\bar{X}}_{j} {\tilde{β}}_{j}) 1_{n} + ε . \end{array}$ En identifiant cela donne $\begin{matrix} (1) & \begin{array}{rcl} β_{j} & = & {\tilde{β}}_{j}, \forall j \in {1, \dots, p - 1}, \\ β_{p} & = & \bar{Y} 1_{n} - \sum_{j = 1}^{p - 1} {\bar{X}}_{j} {\tilde{β}}_{j} . \end{array} \end{matrix}$ Si l’on utilise des variables centrées dans le modèle de régression, on ne met pas de colonne $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).
Maintenant les variables explicatives sont centrées et réduites : $\begin{array}{rcl} \tilde{Y} & = & \tilde{X} \tilde{β} + ε \\ Y - \bar{Y} 1_{n} & = & \sum_{j = 1}^{p - 1} \frac{(X_{j} - {\bar{X}}_{j} 1_{n})}{{\hat{σ}}_{X_{j}}} {\tilde{β}}_{j} + ε \\ Y & = & \sum_{j = 1}^{p - 1} \frac{{\tilde{β}}_{j}}{{\hat{σ}}_{X_{j}}} X_{j} + (\bar{Y} - \sum_{j = 1}^{p - 1} {\bar{X}}_{j} \frac{{\tilde{β}}_{j}}{{\hat{σ}}_{X_{j}}}) 1_{n} + ε . \end{array}$ En identifiant cela donne $\begin{array}{rcl} β_{j} & = & \frac{{\tilde{β}}_{j}}{{\hat{σ}}_{X_{j}}}, \forall j \in {1, \dots, p - 1}, \\ β_{p} & = & \bar{Y} 1_{n} - \sum_{j = 1}^{p - 1} {\bar{X}}_{j} \frac{{\tilde{β}}_{j}}{{\hat{σ}}_{X_{j}}} . \end{array}$ 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 $β_{j}$ sera réduit par rapport à ${\tilde{β}}_{j}$ . Le coefficient constant est donné par la formule ci-dessus.
La variable à expliquer $Y$ est elle aussi centrée-réduite~: $\begin{array}{rcl} \tilde{Y} & = & \tilde{X} \tilde{β} + \tilde{ε} \\ \frac{Y - \bar{Y} 1_{n}}{{\hat{σ}}_{Y}} & = & \sum_{j = 1}^{p - 1} \frac{(X_{j} - {\bar{X}}_{j} 1_{n})}{{\hat{σ}}_{X_{j}}} {\tilde{β}}_{j} + \tilde{ε} \\ Y & = & {\hat{σ}}_{Y} \sum_{j = 1}^{p - 1} \frac{{\tilde{β}}_{j}}{{\hat{σ}}_{X_{j}}} X_{j} + (\bar{Y} - {\hat{σ}}_{Y} \sum_{j = 1}^{p - 1} {\bar{X}}_{j} \frac{{\tilde{β}}_{j}}{{\hat{σ}}_{X_{j}}}) 1_{n} + {\hat{σ}}_{Y} \tilde{ε} . \end{array}$ En identifiant cela donne $\begin{array}{rcl} β_{j} & = & {\hat{σ}}_{Y} \frac{{\tilde{β}}_{j}}{{\hat{σ}}_{X_{j}}}, \forall j \in {1, \dots, p - 1}, \\ β_{p} & = & \bar{Y} 1_{n} - {\hat{σ}}_{Y} \sum_{j = 1}^{p - 1} {\bar{X}}_{j} \frac{{\tilde{β}}_{j}}{{\hat{σ}}_{X_{j}}}, \\ ε & = & {\hat{σ}}_{Y} \tilde{ε} . \end{array}$ 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$ .

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 $β = (1, 2, 0.1, - 4)^{'}$ et $σ = 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

Centrons 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]

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

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

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)

L’estimateur des MC vaut $\begin{array}{r} \hat{β} = (X^{'} X)^{- 1} X^{'} Y, \end{array}$
Calculons maintenant l’estimateur contraint. Nous pouvons procéder de deux manières différentes.

La première consiste à écrire le lagrangien $\begin{array}{r} L = S (β) - λ^{'} (R β - r) . \end{array}$ Les conditions de Lagrange permettent d’obtenir un minimum $\begin{array}{r} {\begin{cases} \frac{\partial L}{\partial β} = - 2 X^{'} Y + 2 X^{'} X {\hat{β}}_{c} - R^{'} \hat{λ} = 0, \\ \frac{\partial L}{\partial λ} = R {\hat{β}}_{c} - r = 0, \end{cases} \end{array}$ Multiplions à gauche la première égalité par $R (X^{'} X)^{- 1}$ , nous obtenons $\begin{array}{rcl} - 2 R (X^{'} X)^{- 1} X^{'} Y + 2 R (X^{'} X)^{- 1} X^{'} X {\hat{β}}_{c} - R (X^{'} X)^{- 1} R^{'} \hat{λ} & = & 0 \\ - 2 R (X^{'} X)^{- 1} X^{'} Y + 2 R {\hat{β}}_{c} - R (X^{'} X)^{- 1} R^{'} \hat{λ} & = & 0 \\ - 2 R (X^{'} X)^{- 1} X^{'} Y + 2 r - R (X^{'} X)^{- 1} R^{'} \hat{λ} & = & 0. \end{array}$ Nous obtenons alors pour $\hat{λ}$ $\begin{array}{r} \hat{λ} = 2 {[R (X^{'} X)^{- 1} R^{'}]}^{- 1} [r - R (X^{'} X)^{- 1} X^{'} Y] . \end{array}$ Remplaçons ensuite $\hat{λ}$ $\begin{array}{rcl} - 2 X^{'} Y + 2 X^{'} X {\hat{β}}_{c} - R^{'} \hat{λ} & = & 0 \\ - 2 X^{'} Y + 2 X^{'} X {\hat{β}}_{c} - 2 R^{'} {[R (X^{'} X)^{- 1} R^{'}]}^{- 1} [r - R (X^{'} X)^{- 1} X^{'} Y] & = & 0, \end{array}$ d’où nous calculons ${\hat{β}}_{c}$ $\begin{array}{rcl} {\hat{β}}_{c} & = & (X^{'} X)^{- 1} X^{'} Y + (X^{'} X)^{- 1} R^{'} {[R (X^{'} X)^{- 1} R^{'}]}^{- 1} (r - R \hat{β}) \\ = & \hat{β} + (X^{'} X)^{- 1} R^{'} {[R (X^{'} X)^{- 1} R^{'}]}^{- 1} (r - R \hat{β}) . \end{array}$ La fonction $S (β)$ à 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 β = 0$ . Calculons analytiquement le projecteur orthogonal sur $M_{0}$ . Rappelons que $\dim (M_{0}) = p - q$ , nous avons de plus $\begin{array}{rcl} R β & = & 0 \Leftrightarrow β \in K e r (R) \\ R (X^{'} X)^{- 1} X^{'} X β & = & 0 \\ U^{'} X β & = & 0 où U = X (X^{'} X)^{- 1} R^{'} . \end{array}$ Nous avons donc que $\forall β \in \ker (R)$ , $U^{'} X β = 0$ , c’est-à-dire que $M_{U}$ , l’espace engendré par les colonnes de $U$ , est orthogonal à l’espace engendré par $X β$ , $\forall β \in \ker (R)$ . Nous avons donc que $M_{U} ⊥ M_{0}$ . Comme $U = X [(X^{'} X)^{- 1} R^{'}]$ , $M_{U} \subset M_{X}$ . En résumé, nous avons $\begin{array}{r} M_{U} \subset M_{X} et M_{U} ⊥ M_{0} donc M_{U} \subset (M_{X} \cap M_{0}^{⊥}) . \end{array}$ Afin de montrer que les colonnes de $U$ engendrent $M_{X} \cap M_{0}^{⊥}$ , 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 $M_{U}$ vaut $q$ . De plus, nous avons vu que $\begin{array}{r} M_{X} = M_{0} \overset{⊥}{\oplus} (M_{0}^{⊥} \cap M_{X}) \end{array}$ et donc, en passant aux dimensions des sous-espaces, nous en déduisons que $\dim (M_{0}^{⊥} \cap M_{X}) = q$ . Nous venons de démontrer que $\begin{array}{r} M_{U} = M_{X} \cap M_{0}^{⊥} . \end{array}$ Le projecteur orthogonal sur $M_{U} = M_{X} \cap M_{0}^{⊥}$ s’écrit $\begin{array}{r} P_{U} = U (U^{'} U)^{- 1} U^{'} = X (X^{'} X)^{- 1} R^{'} [R (X^{'} X)^{- 1} R^{'}]^{- 1} R (X^{'} X)^{- 1} X^{'} . \end{array}$ Nous avons alors $\begin{array}{rcl} \hat{Y} - {\hat{Y}}_{0} & = & P_{U} Y \\ X \hat{β} - X {\hat{β}}_{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{β} . \end{array}$ Cela donne $\begin{array}{r} {\hat{β}}_{0} = \hat{β} - (X^{'} X)^{- 1} R [R (X^{'} X)^{- 1} R^{'}]^{- 1} R \hat{β} . \end{array}$ Si maintenant $r \neq 0$ , nous avons alors un sous-espace affine défini par ${β \in R^{p} : R β = 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 $β_{p} \in R^{p}$ tel que $R β_{p} = r$ et le sous-espace vectoriel associé $M_{0}^{v} = {β \in R^{p} : R β = 0}$ . Les points du sous-espace affine sont alors ${β_{0} \in R^{p} : β_{0} = β_{p} + β_{0}^{v}, β_{0}^{v} \in M_{0}^{v} e t β_{p} : R β_{p} = r}$ . La solution qui minimise les moindres carrés, notée ${\hat{β}}_{0}$ , est élément de ce sous-espace affine et est définie par ${\hat{β}}_{0} = β_{p} + {\hat{β}}_{0}^{v}$ où $\begin{array}{r} {\hat{β}}_{0}^{v} = \hat{β} - (X^{'} X)^{- 1} R^{'} [R (X^{'} X)^{- 1} R^{'}]^{- 1} R \hat{β} . \end{array}$ Nous savons que $R β_{p} = r$ donc $\begin{array}{r} R β_{p} = [R (X^{'} X)^{- 1} R^{'}] [R (X^{'} X)^{- 1} R^{'}]^{- 1} r \end{array}$ donc une solution particulière est $β_{p} = (X^{'} X)^{- 1} R^{'} [R (X^{'} X)^{- 1} R^{'}]^{- 1} r$ . La solution ${\hat{β}}_{0}$ qui minimise les moindres carrés sous la contrainte $R β = r$ est alors $\begin{matrix} (2) & \begin{aligned} {\hat{β}}_{0} & = β_{p} + {\hat{β}}_{0}^{v} \\ = (X^{'} X)^{- 1} R^{'} [R (X^{'} X)^{- 1} R^{'}]^{- 1} r + \hat{β} - (X^{'} X)^{- 1} R^{'} [R (X^{'} X)^{- 1} R^{'}]^{- 1} R \hat{β} \\ = \hat{β} + (X^{'} X)^{- 1} R^{'} [R (X^{'} X)^{- 1} R^{'}]^{- 1} (r - R \hat{β}) . \end{aligned} \end{matrix}$
Calculons l’EQM de $\hat{β}$ qui vaut selon la formule classique: $\begin{aligned} EQM & = (E (\hat{β}) - β) (E (\hat{β}) - β)^{'} + V (\hat{β}) \end{aligned}$ avec $E (\hat{β}) = β$ (sous l’hypothèse que $Y$ est généré par le modèle de régression) et $V (\hat{β}) = σ^{2} (X^{'} X)^{- 1}$ .

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

Calculons maintemant la variance: $\begin{aligned} V ({\hat{β}}_{0}) & = V [\hat{β} + (X^{'} X)^{- 1} R^{'} [R (X^{'} X)^{- 1} R^{'}]^{- 1} (r - R \hat{β})] \\ = V (\hat{β}) + V [(X^{'} X)^{- 1} R^{'} [R (X^{'} X)^{- 1} R^{'}]^{- 1} (r - R \hat{β})] + \\ 2 Cov [\hat{β}; (X^{'} X)^{- 1} R^{'} [R (X^{'} X)^{- 1} R^{'}]^{- 1} (r - R \hat{β})] \\ = σ^{2} V_{X}^{- 1} + (II) + (III) \end{aligned}$ Intéressons nous à la seconde partie $(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{aligned} (II) & = V_{X}^{- 1} R^{'} [R V_{X}^{- 1} R^{'}]^{- 1} V (r - R \hat{β}) [R V_{X}^{- 1} R^{'}]^{- 1} R V_{X}^{- 1} \\ = V_{X}^{- 1} R^{'} [R V_{X}^{- 1} R^{'}]^{- 1} V (R \hat{β}) [R V_{X}^{- 1} R^{'}]^{- 1} R V_{X}^{- 1} \\ = V_{X}^{- 1} R^{'} [R V_{X}^{- 1} R^{'}]^{- 1} R σ^{2} V_{X}^{- 1} R^{'} [R V_{X}^{- 1} R^{'}]^{- 1} R V_{X}^{- 1} \\ = σ^{2} V_{X}^{- 1} R^{'} [R V_{X}^{- 1} R^{'}]^{- 1} R V_{X}^{- 1} \end{aligned}$ La troisième partie est $\begin{aligned} (III) & = 2 V_{X}^{- 1} R^{'} [R V_{X}^{- 1} R^{'}]^{- 1} Cov [\hat{β}; (r - R \hat{β})] \\ = - 2 V_{X}^{- 1} R^{'} [R V_{X}^{- 1} R^{'}]^{- 1} R Cov [\hat{β}; \hat{β}] \\ = - 2 V_{X}^{- 1} R^{'} [R V_{X}^{- 1} R^{'}]^{- 1} R σ^{2} V_{X}^{- 1} \end{aligned}$ Ce qui donne au final $\begin{aligned} V ({\hat{β}}_{0}) & = σ^{2} (X^{'} X)^{- 1} - σ^{2} (X^{'} X)^{- 1} R^{'} [R (X^{'} X)^{- 1} R^{'}]^{- 1} R (X^{'} X)^{- 1} \end{aligned}$ L’écart entre les 2 variances est donc de $σ^{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.