Régression multilinéaire


Le “Cours du prof” ci-dessous explique en détails les bases de la régression multilinéaire par une approche algébrique. Il reste cependant beaucoup à approfondir, notamment sur les sujets suivants : la variance résiduelle et les résidus, la mesure de l’adéquation du modèle aux données et le coefficient de détermination \(R^2\), les intervalles de confiance, interprétation de l’output d’un programme statistique : les tests de Student et de Fisher et la sélection de variables, et le problème plus général de la sélection d’un modèle (dont la sélection de variable est un cas particulier) parmi plusieurs modèles possibles en utilisant le critère d’information bayesienne (BIC).
Chacun de ces aspects est traité en profondeur dans la partie “Code et explications pratiques”.

Modèle général

Donnons-nous deux variables explicatives \(X_1\) et \(X_2\) pour décrire une variable \(Y\), toutes trois variables à valeurs réelles. Les variables explicatives sont déterministes tandis que la variable expliquée \(Y\) est considérée comme aléatoire.

Le modèle général de la régression multilinéaire est alors donné par Les coefficients \(\beta_i\) sont des nombres réels et \(\varepsilon\) est une variable aléatoire réelle d’espérance nulle.

Il s’agit d’une manière, parmi d’autres possibles, de décrire la grandeur aléatoire \(Y\). C’est ce qu’on appelle un modèle probabiliste. Ce modèle particulier qui est celui de la régression multilinéaire fait l’hypothèse que l’espérance de la variable aléatoire \(Y\) est déterministe et se laisse décrire par une fonction linéaire des variables explicatives. Par conséquent ce modèle décrit \(Y\) comme l’espérance de \(Y\) additionnée d’un terme aléatoire d’espérance nulle : \(Y=E[Y]+\varepsilon\) avec \(E[\varepsilon]=0\).

Tenir compte de possibles interactions

Ce modèle suppose qu’il n’y a pas d’interactions entre les variables explicatives, i.e. que \(X_1\) n’influe pas sur \(X_2\). Parfois cela ne reflète pas la réalité. Si l’on souhaite que notre modèle ne fasse pas cette hypothèse mais décrive une dépendance entre les variables explicatives, il y a plusieurs possibilités. La plus simple est de considérer le modèle suivant. La première ligne ajoute, par rapport au modèle sans interaction, un terme mixte. Les lignes suivantes reforment l’expression du modèle avec interaction de façon à pouvoir le comparer avec le modèle sans interaction en le traçant dans un même système d’axes \((X_1, X_2)\).

Le modèle sans interaction revient à supposer que la fonction liant \(Y\) aux variables explicatives est un hyperplan, représenté ci-dessous :

Figure 1 : Comparaison de modèle de régression linéaire multiple sans interaction des variables explicatives. À gauche se trouve un graphique 3D représentant le modèle et à droite un contour plot.
A gauche : graphique 3D représentant le modèle. A droite : sa représentation en lignes de niveau.

Le modèle avec interaction suppose lui que la fonction liant \(Y\) aux variables explicatives est une surface du type de celle représentée ci-dessous :

On n’est pas obligé de modéliser ainsi une interaction entre les variables explicatives. On peut l’exprimer par n’importe quel modèle polynomial, de degré au plus égal au nombre de variables explicatives, et dont la modélisation que nous avons vue n’est qu’un cas particulier.

Ce type de modélisation tenant compte des interactions rentre parfaitement dans le cadre de la régression multiple “sans interaction”. Les variables d’interaction sont en effet des produits de variables connues et sont donc connues. Dans notre modélisation tenant compte des interactions, on définit simplement une troisième variable explicative \(X_3\) égale au produit \(X_1X_2\) et nous retrouvons ainsi la modélisation sans interaction mais avec la variable explicative \(X_3\) en plus.

En conclusion nous pouvons considérer que n’importe quelle transformation connue et fixée des variables explicatives (logarithme, exponentielle, produit, etc.) rentre dans le modèle de régression multiple. Ainsi la transformée d’une variable explicative \(X_1\) par la fonction log par exemple devient \(\tilde{X}_1 = log(X_1)\) et le modèle reste donc un modèle de régression multiple. Un modèle linéaire ne suppose par conséquent pas forcément que c’est le lien entre les variables explicatives et la variable à expliquer qui est linéaire mais que le modèle est linéaire en les paramètres \(\beta_i\).

Calcul de l’estimateur du vecteur des coefficients \(\underline{\hat{\beta}}\)

La première relation importante concernant l’estimateur \(\hat{Y}\) de \(Y\) est

Ici, \(X\) est la matrice des données (the design matrix), dont chaque ligne est associée à une observation et chaque colonne à une variable explicative \(X_1\) ou \(X_2\), sauf la première qui est un vecteur dont chaque coefficient vaut 1. De cette façon, on peut noter une application affine comme une application linéaire, i.e. comme le produit d’une matrice avec un vecteur \(X\cdot\underline{\beta}\), ce qui simplifie toute la théorie.

Quant à \(\underline{\hat{\beta}}\), c’est l’estimateur du vecteur des coefficients du modèle, estimateur pour lequel nous recherchons une formule.

Intéressons-nous à présent à la géométrie du modèle de régression linéaire :

La seconde relation importante concernant le vecteur \(\hat{Y}\) est la suivante.

L’estimateur \(\hat{Y}\) de \(Y\) est en effet l’élément de l’espace engendré par les vecteurs-colonnes de la matrice \(X\) qui est le plus proche de \(Y\). Par définition cet élément, unique, est la projection de \(Y\) sur cet espace. Or toute projection est une application linéaire et en tant que telle se laisse écrire par une multiplication matricielle. La matrice \(H\) est cette matrice de projection qui projette \(Y\) sur \(\hat{Y}\). Quant à l’espace engendré par les vecteurs-colonnes de \(X\), nous le notons \(Spa(X)\), et son orthogonal, \(Spa(X)^{\perp}\).

Nous remarquons que toute matrice de projection est caractérisée par la relation suivante :

Nous calculons à présent l’estimateur du vecteur des coefficients du modèle de régression linéaire. Nous commençons par remarquer que

Cette relation importante montre que la matrice de projection \(H_{\perp}\) sur l’espace orthogonal à celui engendré par les vecteurs-colonnes de \(X\) est donnée par \(H_{\perp}=I-H\), où \(I\) est la matrice unitaire de dimension égale à celle de \(H\), matrice elle-même de dimension égale au nombre de vecteurs-colonnes de \(X\).

Le calcul suit.

Puisque le vecteur \(v\) est choisi tel que \(v\in Spa(X)\), il se laisse en effet exprimer comme une combinaison linéaire des vecteurs-colonnes de la matrice \(X\). C’est ceci que l’on a exprimé ci-dessus par \(v=X\cdot\underline{\alpha}\), pour \(\underline{\alpha}\) un vecteur de coefficients quelconque.

Nous obtenons finalement que :

Calcul de la variance de \(\underline{\hat{\beta}}\)

On fait usage du fait que \(Var(M\cdot v)=M\cdot Var(v)\cdot M^T\) pour \(v\) un vecteur et \(M\) une matrice, tous deux de dimensions telles que leur produit est défini :

On a fait l’hypothèse d’homoscédasticité, i.e. que l’on a supposé que la variance de \(Y\) -et donc celle des résidus- est constante et égale à \(\sigma^2\cdot I\), où \(\sigma^2\) est un réel positif qui ne dépend d’aucune variable explicative. En outre, affirmer que \(Var(Y)=\sigma^2\cdot I= \begin{pmatrix} \sigma^2 & 0 & 0 \\ 0 & \sigma^2 & 0 \\ 0 & 0 & \sigma^2 \end{pmatrix}\) signifie que les projections de \(Y\) sur chacun des vecteurs-colonnes de \(X\) ne sont pas corrélées, puisque cette matrice de variance-covariance est diagonale.

Calcul de l’estimateur des résidus

Nous faisons usage de la relation entre les matrices de projection que nous avons vue ci-dessus.

Les résidus appartiennent donc à \(Spa(X)^{\perp}\), l’espace dit des résidus. Les résidus sont par conséquent toujours orthogonaux à \(\hat{Y}\).

Nous avons les propriétés suivantes.

Espérance et variance des résidus

La relation suivante motive le fait que dans un modèle de régression simple, la somme des résidus est nulle.

Quant à la variance :

Nous avons fait usage de la propriété des matrices de projection que nous avons indiquée plus haut.

Espérance et variance de \(\hat{Y}\)

Le développement algébrique tracé d’une croix rouge qui consiste à développer \((X^T\cdot X)^{-1}=X^{-1}\cdot (X^T)^{-1}\) est tantant mais inapplicable ici puisque \(X\) n’est pas une matrice carrée.

Covariance

En effet, \[ \begin{eqnarray*} Cov(\hat{\varepsilon},\hat{Y}) &=& Cov(\hat{\varepsilon},Y-\hat{\varepsilon}) \\ &=& Cov(\hat{\varepsilon},Y)-Var(\hat{\varepsilon}) \\ &=& \underbrace{Cov(H_{\perp}Y,Y)}_{H_{\perp}\cdot Var(Y)}-\sigma^2\cdot H_{\perp}Y \\ &=& 0 \end{eqnarray*} \]

Calcul de l’estimateur de \(\sigma^2\)

On sait que \(\sigma^2\cdot I=Var(\varepsilon):=E\left[(\varepsilon-E[\varepsilon])\cdot(\varepsilon-E[\varepsilon])^T\right]\) est une matrice de variance-covariance, et qu’elle est égale à \(E\left[\varepsilon\cdot\varepsilon^T\right]\) en raison du fait que \(E[\varepsilon]=0\).

On a donc \[ \begin{pmatrix} \sigma^2 & 0 & 0 \\ 0 & \sigma^2 & 0 \\ 0 & 0 & \sigma^2 \end{pmatrix} =E\left[\varepsilon\cdot\varepsilon^T\right]. \]

Or en appliquant de part et d’autre l’opérateur trace, nous transformons cette équation matricielle en une équation scalaire :

\[ \begin{eqnarray*} tr\begin{pmatrix} \sigma^2 & 0 & 0 \\ 0 & \sigma^2 & 0 \\ 0 & 0 & \sigma^2 \end{pmatrix} &=& tr\left(E\left[\varepsilon\cdot\varepsilon^T\right]\right) \\ &=& tr\left(E\left[\varepsilon^T\cdot\varepsilon\right]\right) =tr\left(E[||\varepsilon||^2]\right), \end{eqnarray*} \]

où l’on a usé du fait que la trace d’une matrice carrée est égale à celle de sa transposée. Ce raisonnement est la raison pour laquelle nous nous intéressons particulièrement à la grandeur suivante :

Par la suite, nous cherchons à exprimer cette grandeur scalaire de façon à en faire ressortir le \(\sigma^2\).

Nous employons à nouveau le fait que la trace d’une matrice carrée est égale à celle de sa transposée.

Explications :

  • Le passage de la 3e à la 4e ligne : \(E[\hat{\varepsilon}^T\cdot\hat{\varepsilon}]\) est égal à \(Var(\hat{\varepsilon})\), puisque \(E[\hat{\varepsilon}]=0\). Et on sait que \(Var(\hat{\varepsilon})=\sigma^2\cdot H_{\perp}\)

  • A la dernière ligne, on a utilisé le fait que la trace d’une matrice de projection est égale à la dimension de l’espace de projection. Or la matrice de projection orthogonale \(H_{\perp}\) projette dans un espace de dimension \(n-p\), où \(p\) est le nombre de vecteurs-colonnes de la matrice \(X\) et \(n\) est le nombre de lignes de \(X\).

Nous pouvons à présent définir un estimateur de \(\sigma^2\).

Cet estimateur est non biaisé, i.e. \(E[\hat{\sigma}^2]=\sigma^2\).

Prévision

Nous pouvons à présent exprimer pleinement l’estimateur du vecteur des coefficients du modèle de régression multilinéaire :

Position du problème

Considérons une nouvelle observation, non déjà contenue parmi les \(n\) lignes de \(X\) la matrice des données. C’est un vecteur de dimension \(p=\)nombre de vecteurs-colonnes de \(X\), vecteur que nous notons \(x_{n+1}=(x_{n+1}^1,\ldots,x_{n+1}^p)^T\). Nous voulons prédire la valeur associée \(y_{n+1}\).

Nous savons que \(y_{n+1}=x_{n+1}\cdot\underline{\beta}+\varepsilon_{n+1}\),

\(E[\varepsilon_{n+1}]=0\), \(Var(\varepsilon_{n+1})=\sigma^2\) et \(Cov(\varepsilon_{n+1},\varepsilon_i)=0\quad\forall i=1,\ldots,n\).

La prédiction est donc \(\hat{y}_{n+1}=x_{n+1}\cdot\underline{\hat{\beta}}\).

L’erreur de prévision \(y_{n+1}-\hat{y}_{n+1}\)

Nous calculons la variance de l’erreur de prévision :

\[ \begin{eqnarray*} Var(y_{n+1}-\hat{y}_{n+1}) &=& Var(x_{n+1}\cdot\underline{\beta}+\varepsilon_{n+1} -x_{n+1}\cdot\underline{\hat{\beta}}) \\ &=& \sigma^2+x_{n+1}^T\cdot Var(\underline{\hat{\beta}})\cdot x_{n+1} \\ &=& \sigma^2\cdot\left(1+x_{n+1}^T\cdot\left(X^T\cdot X\right)^{-1}\cdot x_{n+1}\right) \end{eqnarray*} \]

C’est l’incertitude sur \(\sigma^2\) plus l’incertitude sur \(\hat{y}_{n+1}\). C’est bien sûr également \(Var(\hat{\varepsilon}_{n+1})\). Cette expression permet d’associer graphiquement à la droite de régression une bande d’erreur de prévision.