I. Hypothèses sur le modèle
1) Cas général
Un modèle de régression linéaire est défini par : \(Y=X\beta + \varepsilon\) où
- \(Y\) est un vecteur aléatoire de dimension \(n\) appelé variable à expliquer
- \(X\) est une matrice (non-aléatoire) de taille \(n\times p\) appelée matrice du plan d’expérience
- \(\beta\) est le vecteur de dimension \(p\) des paramètres inconnus du modèle
- \(\varepsilon\) est le vecteur aléatoire de dimension \(n\) des erreurs (ou bruits)
Remarque : Les colonnes de \(X\) sont aussi appelées variables explicatives du modèle.
De plus, il y a des hypothèses additionnelles au modèle : \[ \begin{cases} \text{rg}(X)=p \\ \mathbb{E}[\varepsilon]=0 \\ \mathbb{V}(\varepsilon)=\sigma ^2I_n\\ \end{cases} \] Si le rang de \(X\) est inférieur à \(p\) alors il y a des vecteurs de la matrice qui sont combinaisons linéaires des autres et donc il y a une redondance de l’information. On peut expliquer un vecteur avec les autres. Les erreurs respectent la propriété d’homoscédasticité : \(\forall i \;\; \mathbb{V}(\varepsilon_i)=\sigma^2\) et elles sont également centrées.
2) Cas où \(p=2\)
Supposons que l’intercepte est dans le modèle. Alors la matrice \(X\) et \(\beta\) s’écrivent :
\[ \left(\begin{matrix} 1 & x_1 \\ 1 & x_2 \\ . & . \\ . & . \\ 1 & x_n \\ \end{matrix}\right) \;\;\; \text{et} \;\;\; \left(\begin{matrix} \beta_1 \\ \beta_2 \\ \end{matrix}\right) \] Et pour tout \(i\) on a \(y_i= \beta_1 +\beta_2 x_i + \varepsilon_i\).
Remarque : Comme rg\((X)=2\), la variable explicative n’est pas constante.
3) Exemple
n <- 100
beta_1 <- 4
beta_2 <- 8
X <- runif(n, -4, 4)
eps <- rnorm(n, 0, 5)
y <- beta_1 + beta_2*X + eps
plot(X, y, type='p', main='Exemple')
Dans cet exemple, \(\beta_1=4\) et
\(\beta_2=8\). La variable \(X\) explicative est de taille \(100\) et est modélisée par une loi uniforme
\(\mathcal{U}[-4,4]\). Les erreurs,
eps, sont modélisées par une loi normale de moyenne \(0\) et de variance \(25\). Ensuite, on trace \(Y=\beta_1 + \beta_2X+\varepsilon\).
II. Estimateurs OLS
1) Définition et expression
On définit l’estimateur des moindres carrés ordinaires (OLS = Ordinary Least Squared) comme : \[\hat{\beta}=\underset{\beta \in \mathbb{R}^p}{\text{argmin}} ||Y-X\beta||^2\] Si on note \(\mathcal{M}(X)\) l’espace engendré par les colonnes de \(X\), \(X\hat{\beta}\) est la projection orthogonale de \(Y\) sur \(\mathcal{M}(X)\) et on note \(\hat{Y}=X\hat{\beta}\). Donc, \(Y-X\hat{\beta}\) est orthogonal à \(\mathcal{M}(X)\) et on a :
\[ \forall i \in \{1,..,p\},~~<Y-X\hat{\beta},X_i>=0 \\ \Leftrightarrow \forall i \in \{1,..,p\},~~ X_i'(Y-X\hat{\beta})=0 \\ \Leftrightarrow X'Y=X'X\hat{\beta} \\ \Leftrightarrow (X'X)^{-1}X'Y=\hat{\beta} \]
On obtient alors l’expression suivante : \[\hat{\beta}= (X'X)^{-1}X'Y\] On remarque bien que cet estimateur est linéaire en \(Y\) car de la forme \(AY\) avec \(A\) une matrice.
Rappel : On note la matrice de projection orthogonale sur \(\mathcal{M}(X)\) : \[P_X=X(X'X)^{-1}X'\] Et on a alors que \(\hat{Y}=P_XY\)
2) Moments
L’espérance et la matrice de variance de \(\hat{\beta}\) sont :
\[ \mathbb{E}[\hat{\beta}]=\beta \;\;\;\;\;\;\;\;\;\;\;\;\;\;\; \mathbb{V}[\hat{\beta}]=\sigma^2(X'X)^{-1} \]
De plus, l’estimateur OLS est de variance minimale parmi les estimateurs sans biais, linéaires en \(Y\) (Théorème de Gauss-Markov).
Note : Il existe un ordre partiel dans l’ensemble des matrices symétriques \(\mathcal{S}_n(\mathbb{R})\) défini par : \[ \forall S_1, S_2 \in \mathcal{S}_n(\mathbb{R}), \;\;S_1\le S_2 \Leftrightarrow S=S_2-S_1 \;\; \text{semi-définie positive} \] Cet ordre est utilisé dans le théorème de Gauss-Markov.
3) Expressions pour le cas \(p=2\)
En utilisant les formules précédentes, on retrouve facilement le cas de la régression linéaire simple :
\[ \hat{\beta}_1=\bar{y}-\hat{\beta}_2\bar{x} \;\;\;\;\; \text{et} \;\;\;\;\; \hat{\beta}_2=\frac{s_{X,Y}^{2}}{s_X^2} \] où \(s_{X,Y}^{2}=\frac{1}{n}\underset{i=1}{\overset{n}{\sum}}(x_i-\bar{x})(y_i-\bar{y})\) et \(s_X^2=\frac{1}{n}\underset{i=1}{\overset{n}{\sum}}(x_i-\bar{x})^2\)
De plus, \(\hat{\beta}_1\) et \(\hat{\beta}_2\) sont des estimateurs sans biais de \(\beta_1\) et \(\beta_2\) respectivement. Donc : \[ \mathbb{E}[\hat{\beta}_1]=\beta_1 \;\;\;\; \text{et} \;\;\;\; \mathbb{E}[\hat{\beta}_2]=\beta_2 \] On retrouve également les résultats des variances : \[ \mathbb{V}[\hat{\beta}_1]=\sigma^2\bigg(\frac{1}{n}+\frac{\bar{x}^2}{\underset{i=1}{\overset{n}{\sum}}(x_i-\bar{x})^2}\bigg) \;\;\;\;\; \text{et} \;\;\;\;\; \mathbb{V}[\hat{\beta}_2]=\frac{\sigma^2}{\underset{i=1}{\overset{n}{\sum}}(x_i-\bar{x})^2} \] Puis la covariance : \[ \mathbb{C}ov(\hat{\beta}_1,\hat{\beta}_2)=\frac{-\sigma^2\bar{x}}{\underset{i=1}{\overset{n}{\sum}}(x_i-\bar{x})^2} \]
4) Exemple
# On utilise pour cet exemple les données de la partie précédente
S_X <- sum((X-mean(X))^2)
S_XY <- sum((X-mean(X))*y)
h_beta2 <- S_XY/S_X
h_beta1 <- mean(y)-h_beta2*mean(X)
plot(X,y, main ="Régression linéaire")
points(X,h_beta1+h_beta2*X,type='l',col='red')
legend("topleft",legend = c("Droite de régression", "y"), col = c("red", "black"),
lty = c(1,0), pch = c(-1,1))
Sous R, il y a la fonction lm() qui donne
plusieurs informations sur la régression, permettant notamment de tracer
la droite :
model <- lm(y~X)
plot(X,y, main='Regression linéaire avec lm')
abline(model, col='blue')
legend("topleft",legend = c("Droite de régression", "y"), col = c("blue", "black"),
lty = c(1,0), pch = c(-1,1))
De plus, on voit que les coefficients correspondent :
## (Intercept)
## 4.07389 4.07389## X
## 8.220286 8.220286III. Résidus et prédictions
1) Résidus et propriétés
Les résidus sont définis par : \[ \hat{\varepsilon}=\left(\begin{matrix} \hat{\varepsilon}_1\\ \hat{\varepsilon}_2\\ ...\\ \hat{\varepsilon}_n\\ \end{matrix}\right) = Y-\hat{Y} = P_{X^\perp}Y=P_{X^\perp}\varepsilon \] Les résidus ne sont donc rien d’autres que la projection du vecteur d’erreur \(\varepsilon\) sur l’orthogonal de \(\mathcal{M}(X)\).
\(\hat{\varepsilon}\) et \(\hat{Y}\) possèdent les propriétés suivantes :
- \(\mathbb{E}[\hat{\varepsilon}]=0\), en moyenne l’erreur est nulle.
- \(\mathbb{V}[\hat{\varepsilon}]=\sigma^2P_{X^\perp}\)
- \(\mathbb{E}[\hat{Y}]=X\beta\)
- \(\mathbb{V}[\hat{Y}]=\sigma^2P_X\)
- \(\mathbb{C}ov(\hat{\varepsilon},\hat{Y})=0\)
2) Estimation de \(\sigma ^2\)
On peut donner un estimateur non biaisé naturel de la variance \(\sigma ^2\) qui s’exprime :
\[ \hat{\sigma}^2 = \frac{1}{n-p}||\hat{\varepsilon}||^2 = \frac{1}{n-p} \sum_{i=1}^n \hat{\varepsilon}_i^2=\frac{RSS}{n-p} \] où \(RSS=\) Residual Sum of Squares.
3) Prévision
Si on se donne une nouvelle observation \(x^{n+1}\), il est possible d’utiliser le modèle de régression dans le but de faire une prédiction sur la valeur de la variable \(y\) correspondante. Si on note \(y_{n+1}\) la vraie valeur de cette variable, et \(\hat{y}_{n+1} = x^{n+1} \hat{\beta}\) la valeur approchée prédite par le modèle, on peut définir alors l’erreur de prévision :
\[ \hat{\varepsilon}_{n+1}= y_{n+1} - \hat{y}_{n+1} \] Si on peut écrire le modèle : \[ y_{n+1} = x^{n+1}\beta+ \varepsilon_{n+1} \] avec \(\mathbb{E}[\varepsilon_{n+1}] = 0\), \(\mathbb{V}[\varepsilon_{n+1}] =\sigma^2\) et pour tout \(j \in \{1,...,n\}, \mathbb{C}ov(\varepsilon_j,\varepsilon_{n+1})=0\), la variable aléatoire \(\hat{\varepsilon}_{n+1}\) possède les propriétés suivantes :
- \(\mathbb{E}[\hat{\varepsilon}_{n+1}] = 0\)
- \(\mathbb{V}[\hat{\varepsilon}_{n+1}] = \sigma^2(1+x^{n+1}(X'X)^{-1}(x^{n+1})')\)
Lorsque \(p=2\), on en déduit la formule suivante :
\[ \mathbb{V}[\hat{\varepsilon}_{n+1}] = \sigma^2\bigg(1+ \frac{1}{n} + \frac{(x_{n+1} - \bar{x})^2}{\underset{i=1}{\overset{n}{\sum}}(x_i-\bar{x})^2}\bigg) \] Illustrons cela en utilisant notre exemple précédent :
## 1
## 28.73475
Il est aussi possible (en faisant une hypothèse gaussienne sur le modèle), de pouvoir construire des intervalles de prédiction !
## fit lwr upr
## 1 28.73475 19.89144 37.57805library(ggplot2)
data <- data.frame(y=y, x=X)
ggplot(data=data,aes(x=x, y=y))+
geom_point()+
geom_point(aes(x=3, y=P[1]), shape=7, col='red')+
geom_smooth(method='lm', se=F)+
geom_segment(aes(x=3, y=P[2], xend=3, yend=P[3]),
color='red', linetype="twodash", show.legend=TRUE)+
labs(title="Régression avec intervalle de prédiction pour x=3")+
scale_linetype_manual("", values=c("Intervalle de prédiction"=2))+
theme(plot.title=element_text(hjust=0.5, face="bold"))
4) Interprétation géométrique : vers le \(R^2\)
Si la constante fait partie du modèle, en utilisant le théorème de Pythagore sur la figure suivante on obtient : \[ ||Y-\bar{Y}\mathbb{1}||^2=||\hat{Y}-\bar{Y}\mathbb{1}||^2+||\hat{\varepsilon}||^2 \] \[ \Leftrightarrow TSS = ESS + RSS \] où \(TSS=\) Total Sum of Squares, \(ESS=\) Explained Sum of Squares et \(RSS=\) Residual Sum of Squares.
De cette remarque nous pouvons en déduire le coefficient de détermination, noté \(R^2\), défini par :
\[ R^2 = cos^2(\theta) = \frac{||\hat{Y}-\bar{Y}\mathbb{1}||^2}{||Y-\bar{Y}\mathbb{1}||^2} = \frac{ESS}{TSS} = 1 - \frac{RSS}{TSS} \]
Si le \(R^2=1\), alors \(Y\) est dans \(\mathcal{M}(X))\) c’est-à-dire \(y_i=\beta_1+\beta_2x_i \;\; \forall i\)
Si le \(R^2=0\), cela veut dire que \(\hat{y}_i=\bar{y} \;\; \forall i\), par conséquent le modèle de régression linéaire est alors inadapté.
Ce coefficient peut être ajusté à l’aide du coefficient de détermination ajusté, noté \(R_{a}^2\), défini par :
\[ R_{a}^2 = 1-\frac{RSS/(n-p)}{TSS/(n-1)} =1-\frac{n-1}{n-p}(1-R^2) \] Le \(R^2\) augmente lorsqu’on ajoute des variables ne décrivant pas le modèle c’est pourquoi on utilise le \(R^2_a\) qui, lui, pénalise les points de données qui ne correspondent pas au modèle de régression et n’augmente donc pas.
On regardera donc principalement le \(R_{a}^2\). Sous R, il suffit
d’utiliser la fonction summary sur notre modèle pour avoir
accès à cette quantité :
##
## Call:
## lm(formula = y ~ X)
##
## Residuals:
## Min 1Q Median 3Q Max
## -10.8545 -2.7869 -0.7722 2.8166 11.5458
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 4.0739 0.4396 9.266 4.73e-15 ***
## X 8.2203 0.1941 42.354 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 4.396 on 98 degrees of freedom
## Multiple R-squared: 0.9482, Adjusted R-squared: 0.9477
## F-statistic: 1794 on 1 and 98 DF, p-value: < 2.2e-16Ici, le \(R_{a}^2\) est représenté par la valeur de “Adjusted R-squared”.