1.1 L’estimateur MC

L’estimateur des moindres carrés, \(\widehat{\mathbf{w}}\) du vecteur \(\mathbf{w}\), est défini par \[\begin{eqnarray} \widehat{\mathbf{w}} & := & \underset{(w_0,w_1,\ldots,w_p)\in\mathbb{R}^{1+p}}{\arg\min} \,\, \frac{1}{n} \sum_{i=1}^n\left(Y_i-w_0-w_1\, X_{i,1}-\cdots - w_p\,X_{i,p}\right)^2\nonumber\\ & =: & \underset{\mathbf{w}\in\mathbb{R}^{1+p}}{\arg\min} \,\, \frac{1}{n} \left\| \mathbf{Y} - \mathbb{X}\, \mathbf{w} \right\|^2.\nonumber\\ & = & (\mathbb{X}^\top \, \mathbb{X})^{-1}\, \mathbb{X}^\top \,\, \mathbf{Y}. \end{eqnarray}\] La dernière égalité a lieu si la matrice \(\mathbb{X}^\top \, \mathbb{X}\) est inversible, ce qui est équivalent à supposer que les colonnes de \(\mathbb{X}\) sont linéairement indépendantes (donc le nombre d’observations \(n\) doit être supérieur à \(1+p\), \(p\) étant le nombre de variables explicatives). Cette hypothèse garantit l’existence et l’unicité de l’estimateur MC précédent.

1.2 Lien avec l’estimation par maximum de vraisemblance (MV)

On suppose que les termes d’erreur, \[\varepsilon_1,\, \ldots,\, \varepsilon_n,\] sont i.i.d de même loi normale \(\mathcal{N}(0,\sigma^2)\), \(\mathbb{E}(\varepsilon\,|\, \mathbf{X}=\mathbf{x})=0\), et \(\text{Var}(\varepsilon\,|\, \mathbf{X}=\mathbf{x} )=\sigma^2\) ne dépendant pas de \(\mathbf{x}\) (hypothèse d’homoscédasticité). Par conséquent, la loi de \(Y\), conditionnellement à \(\mathbf{X}=\mathbf{x}\), est une loi normale d’espérance \(w_0+w_1x_1+\ldots +w_px_p\) et de variance \(\sigma^2.\,\) La log-vraisemblance (conditionnelle) s’écrit donc sous la forme \[\begin{eqnarray*} \mathcal{L}_n(\mathbf{w},\sigma^2) = -\frac{n}{2}\log \sigma^2 -\frac{n}{2}\log(2\pi)-\frac{1}{2\sigma^2} \|\mathbf{Y}-\mathbb{X}\,\mathbf{w}\|^2. \end{eqnarray*}\] On peut voir dans ce cas que l’estimateur du maximum de vraisemblance (MV), de \(\mathbf{w}\), coïncide avec l’estimateur des MC : Notons \((\widetilde{\mathbf{w}},\widetilde{\sigma}^2)\) l’EMV de \((\mathbf{w},\sigma^2)\), i.e., \[(\widetilde{\mathbf{w}},\widetilde{\sigma}^2) :=\underset{\mathbf{w}\in\mathbb{R}^{1+p},\, \sigma^2\in\mathbb{R}_+^*}{\arg\max} \,\, \mathcal{L}_n(\mathbf{w},\sigma^2).\] Il est clair que l’EMV \(\widetilde{\mathbf{w}} = \underset{\mathbf{w}}{\arg\min} \,\, \left\|\mathbf{Y}-\mathbb{X}\,\mathbf{w}\right\|^2 =: \widehat{\mathbf{w}} =: EMC.\) L’EMV \(\widetilde{\sigma}^2\), de \(\sigma^2\), vaut \[ \widetilde{\sigma}^2 = \frac{1}{n} \left\|\mathbf{Y}-\mathbb{X}\,\widehat{\mathbf{w}}\right\|^2.\]

La log-vraisemblance du modèle (évaluée à l’EMV) est donnée par \[\begin{eqnarray}\label{la vraisemblance du MRLM} \mathcal{L}_n(\widetilde{\mathbf{w}},\widetilde{\sigma}^2) = - \frac{n}{2}\log \frac{\left\|\widehat{\boldsymbol{\varepsilon}}\right\|^2}{n} - \frac{n}{2} (1+\log (2\pi)), \end{eqnarray}\] où \[\widehat{\boldsymbol{\varepsilon}} := (\widehat{\varepsilon}_1,\ldots, \widehat{\varepsilon}_n)^\top := \mathbf{Y}-\mathbb{X}\,\widehat{\mathbf{w}} =: \mathbf{Y} - \widehat{\mathbf{Y}}\] représentant le vecteur des résidus : les écarts entre les valeurs observées \(\mathbf{Y} :=(Y_1,\ldots,Y_n)^\top\), de la variable réponse \(Y\), et les valeurs ajustées \[\widehat{\mathbf{Y}} := \mathbb{X}\, \widehat{\mathbf{w}} =: (\widehat{Y}_1,\ldots,\widehat{Y}_n)^\top.\]

3.1 Proposition

On a

1. \(\widehat{\mathbf{w}}\) est un vecteur gaussien d’espérance \(\mathbf{w}\) et de matrice de variance-covariance \(\sigma^2 \left(\mathbb{X}^\top\mathbb{X}\right)^{-1}\);

2. La statistique \((n-p-1) \,\, \frac{\widehat{\sigma}^2}{\sigma^2}\) suit la loi du \(\chi^2_{(n-p-1)}\, (\text{la loi du } \chi^2 \text{ à } n-p-1 \text{ degrès de liberté})\);

3. \(\widehat{\mathbf{w}}\) et \(\widehat{\sigma}^2\) sont indépendants.

3.2 Intervalles de confiance et tests

On note \(\widehat{\sigma}^2_j := \widehat{\sigma}^2 \left[\left(\mathbb{X}^\top \mathbb{X}\right)^{-1}\right]_{j,j}\), pour \(j=0,1,\ldots,p.\) On a \[\forall j= 0,\ldots,p,\, \text{ la statistique }\, \frac{\widehat{w}_j-w_j}{\widehat{\sigma}_j} \, \text{ suit la loi }\, \mathcal{T}_{(n-p-1)},\] (la loi de Student à \(n-p-1\) degrès de liberté), ce qui permet de construire des intervalles de confiance, au niveau \(1-\alpha\), pour les paramètres \(w_j\), et de réaliser des tests d’hypothèses du type \(\mathcal{H}_0 : w_j =0\) contre \(\mathcal{H}_1 : w_j \neq 0\).

3.3 Intervalle de confiance, au niveau \(1-\alpha\), pour \(w_j\)

\[\left[ \widehat{w}_j - t_{(n-p-1)}(1-\alpha/2) \, \widehat{\sigma}_j, \, \, \widehat{w}_j + t_{(n-p-1)}(1-\alpha/2) \, \widehat{\sigma}_j \right],\] où \(t_{(n-p-1)}(1-\alpha/2)\) est le quantile d’ordre \(1-\alpha/2\) de la loi de Student à (n-p-1) degrés de liberté.

3.4 Test de Student

Considérons le problème de test de l’hypothèse nulle \(\mathcal{H}_0 \, : \, w_j =0\) contre l’alternative \(\mathcal{H}_1 \, : \, w_j \neq 0\).

Notons \(t := \frac{\widehat{w}_j}{\widehat{\sigma}_j}\) (la statistique de Student). La P-value du test (de Student) est donnée par \[\text{P-value} = \mathbb{P}\left(|T|>|t|\right),\] où \(T\) est une variable aléatoire suivant la loi \(\mathcal{T}_{(n-p-1)}.\)

3.5 Prévision et intervalles de confiance

On dispose d’une nouvelle observation \(\mathbb{X}_{n+1}:=(1,X_{n+1,1},\ldots,X_{n+1,p})\). On prédit la valeur \(Y_{n+1}\) correspondante par \[ \widehat{Y}_{n+1} := \mathbb{X}_{n+1}\, \widehat{\mathbf{w}}.\]

Intervalle d’estimation, au niveau \(1-\alpha\), pour \(Y_{n+1}\) : \[\left[\widehat{Y}_{n+1} - t_{(n-p-1)}(1-\alpha/2) \, \widehat{\sigma} \, \sqrt{\mathbb{X}_{n+1}\, \left(\mathbb{X}^\top \, \mathbb{X}\right)^{-1} \mathbb{X}_{n+1}^\top },\,\,\widehat{Y}_{n+1} + t_{(n-p-1)}(1-\alpha/2) \, \widehat{\sigma} \, \sqrt{\mathbb{X}_{n+1}\, \left(\mathbb{X}^\top \, \mathbb{X}\right)^{-1} \mathbb{X}_{n+1}^\top }\right];\]

Intervalle de prévision, au niveau \(1-\alpha\), pour \(Y_{n+1}\) : \[\left[\widehat{Y}_{n+1} - t_{(n-p-1)}(1-\alpha/2) \, \widehat{\sigma} \, \sqrt{1 + \mathbb{X}_{n+1} \, \left(\mathbb{X}^\top \, \mathbb{X}\right)^{-1} \mathbb{X}_{n+1}^\top}, \,\, \widehat{Y}_{n+1} + t_{(n-p-1)}(1-\alpha/2) \, \widehat{\sigma} \, \sqrt{1 + \mathbb{X}_{n+1} \, \left(\mathbb{X}^\top \, \mathbb{X}\right)^{-1} \mathbb{X}_{n+1}^\top} \right].\]

3.6 Équation de l’analyse de variance

On a d’après le théorème de Pythagore \[ \begin{array}{rcl} \|\mathbf{Y} - \overline{\mathbf{Y}}\mathbf{1}\|^2 & = & \|\widehat{\mathbf{Y}}-\overline{\mathbf{Y}}\mathbf{1}\|^2+ \| \widehat{\boldsymbol{\varepsilon}}\|^2\\ SST & = & SSR + SSE. \end{array} \] (somme totale des carrés \(=\) somme des carrés de la régression \(+\) somme des carrés résiduels)

(total sum of squares \(=\) regression sum of squares \(+\) sum of squared errors)

3.7 Coefficient de détermination \(R^2\)

Le coefficient de détermination \(R^2\) est défini par \[R^2 := \frac{\|\widehat{\mathbf{Y}}-\overline{\mathbf{Y}}\mathbf{1}\|^2}{\|\mathbf{Y} - \overline{\mathbf{Y}}\mathbf{1}\|^2} =:\frac{SSR}{SST}.\]

Il vérifie les propriétés suivantes

1. \(0\leq R^2 \leq 1;\)

2. Si \(R^2=1\), la variabilité de la variable réponse est entièrement expliquée par le modèle;

3. Si \(R^2=0\), toute la variabilité se trouve dans le bruit (le terme d’erreur).

3.8 Test du modèle global

Le modèle de RLM s’écrit \[Y_i = w_0 + w_1\, X_{i,1}+\cdots + w_p\, X_{i,p} + \varepsilon_i, \,\, i=1,\ldots,n,\] où les termes d’erreurs \(\varepsilon_i, i=1,\ldots,n,\) sont supposés ici i.i.d. de même loi \(\mathcal{N}(0,\sigma^2)\).

On veut tester \[\mathcal{H}_0 : w_1=\cdots = w_p=0 \,\,\text{ contre } \,\, \mathcal{H}_1 : \exists j \in \{1,\ldots,p\} \text{ t.q. } w_j\neq 0. \]

Sous \(\mathcal{H}_0\), la statistique de Fisher \[F := \frac{R^2}{1-R^2}\frac{n-p-1}{p} \,\text{ suit la loi }\, \mathcal{F}_{(p,n-p-1)} \,\, (\text{la loi de Fisher à } p \text{ et } n-p-1 \text{ degrés de liberté}).\]

On rejette \(\mathcal{H}_0\) si \(F > F_{(p,n-p-1)}(1-\alpha)\), où \(F_{(p,n-p-1)}(1-\alpha)\) est le quantile d’ordre \((1-\alpha)\) de la loi \(\mathcal{F}_{(p,n-p-1)}\).

La P-value est donc donnée par \(\text{P-value} = \mathbb{P}(Z > F),\) où \(Z\) est une variable aléatoire suivant la loi \(\mathcal{F}_{(p,n-p-1)}\).

3.9 Tests entre modèles emboîtés

On veut tester le modèle réduit \(\mathcal{M}_0\) \[Y_i = w_0 + w_{q+1}\, X_{i,q+1}+\cdots + w_p\, X_{i,p} + \varepsilon_i,\] avec \(1\leq q <p\), à l’intérieur du modèle plus large \(\mathcal{M}\) \[Y_i = w_0 + w_1\, X_{i,1}+\cdots + w_p\, X_{i,p} + \varepsilon_i.\]

Cela revient à tester (à l’intérieur du modèle \(\mathcal{M}\)) l’hypothèse nulle \[\mathcal{H}_0 : w_1=\cdots = w_q=0 \,\,\text{ contre } \,\, \mathcal{H}_1 : \exists j \in \{1,\ldots,q\} \text{ t.q. } w_j\neq 0.\]

Pour réaliser le test précédent, on utilise la statistique de Fisher suivante \[ F := \frac{\|\widehat{\mathbf{Y}}_0 - \widehat{\mathbf{Y}}\|^2/q} {\|\mathbf{Y} - \widehat{\mathbf{Y}}\|^2/(n-p-1)} \] qui suit la loi \(\mathcal{F}_{(q,n-p-1)}\), sous \(\mathcal{H}_0\). (\(\widehat{\mathbf{Y}}_0\) désigne le vecteur des valeurs ajustées du modèle réduit).

On rejette \(\mathcal{H}_0\) si \(F>F_{(q,n-p-1)}(1-\alpha).\)

La P-value est donc définie par \(\text{P-value} := \mathbb{P}\left(Z > F\right),\) où \(Z\) est une variable aléatoire suivant la loi \(\mathcal{F}_{(q,n-p-1)}.\)

Ce test peut se faire sous R à l’aide de la fonction anova() qu’on applique à des modèles linéaires calculés avec la fonction lm().

3.10 Exercice 1

1. Construire le modèle de RLM de la variable Sales en fonction de toutes les variables de la base de données Carseats : utiliser la fonction lm();

2. Commenter tous les résultats;

3. Vérifier graphiquement : (i) la non-corrélation des erreurs (utiliser la fonction acf()); (ii) la relation linéaire entre la variable réponse et les variables explicatives; (iii) l’hypothèse d’homoscédasticité des erreurs (appliquer la fonction plot() au modèle calculé par la fonction lm());

4. Tester l’hypothèse d’homoscédasticité des erreurs: utiliser le test de Breusch-Pagan (fonction bptest() du package lmtest);

5. Vérifier graphiquement l’hypothèse de normalité des erreurs : utiliser la représentation Q-Q plot, et comparer l’histograme des erreurs et la densité gaussienne;

6. Donner les résulats du test de Student de l’hypothèse nulle de non-significativité de chacune des variables explicatives. Ordonner les variables explicatives de la plus significative à la moins significative;

7. Donner les résulats du test de Fisher de l’hypothèse nulle de non-significativité de chacune des variables explicatives. Ordonner les variables explicatives de la plus significative à la moins significative;

8. Comparer les résultats des deux questions précédentes, et commenter.

4.1 Taille d’un modèle et précision

Lorsque la taille du modèle et petite (ou nombre petit de variables explicatives) :

1. variance faible, biais très élevé (erreur quadratique élevée);

2. erreur théorique de prévision élevée;

3. erreur empirique (d’ajustement) élevée.

Lorsque la taille du modèle et grande (ou nombre élevé de variables explicatives) :

1. variance très élevée, biais faible (erreur quadratique élevée);

2. erreur théorique de prévision élevée;

3. erreur empirique (d’ajustement) très faible (problème de sur-ajustement).

D’où la nécessité de développer des procédures de sélection de modèles (de variables).

4.2 Le cadre

On considère un modèle de RLM \[Y = w_0+w_1\, X_1+\ldots +w_p\, X_p +\varepsilon.\] On dispose d’un échantillon \((\mathbf{X}_1,Y_1),\ldots,(\mathbf{X}_n,Y_n)\) du vecteur \((\mathbf{X},Y)=:(X_1,\ldots,X_p,Y)\in \mathbb{R}^p\times \mathbb{R}\).

L’objectif est d’obtenir le sous-ensemble de variables explicatives qui conduit au ``meilleur’’ modèle RLM au sens d’un critère donné.

Avec \(p\) variables explicatives candidates, \(X_1,\ldots,X_p\), on peut construire \(2^p-1\) modèles de régression linéaires différents (les modèles à une variable, à deux variables, …, à \(p\) variables).

4.3 Exemples de critères de sélection de modèles

1. L’Akaike Information Criterion (AIC), d’un modèle de RLM, constitué de \(k\) variables explicatives, est défini par \[ AIC = -2 \, \mathcal{L}_n(\widehat{\mathbf{w}},\widetilde{\sigma}^2) + 2 \, (k+2),\] où \(\mathcal{L}_n(\widehat{\mathbf{w}},\widetilde{\sigma}^2)\) est la log-vraisemblance du modèle, définie ci-dessus, sous les hypothèses d’homoscédasticité et de normalité des erreurs;

2. Le Bayesian Information Criterion (BIC) : \[ BIC = -2\, \mathcal{L}_n(\widehat{\mathbf{w}},\widetilde{\sigma}^2) + \log(n) \, (k+2);\]
3. \(R^2\)-ajusté : \[R^2_a := 1-\frac{n-1}{n-k-1}\left(1-R^2\right);\]

4. Le \(C_p\) de Mallow d’un modèle de régression utilisant \(k\) variables explicatives (\(1\leq k\leq p\)) est donné par : \[C_p:=\frac{1}{n} \left(\|\mathbf{Y}-\widehat{\mathbf{Y}}_0\mathbf{1}\|^2 + 2 (1+k) \, \widehat{\sigma}^2\right),\] où \(\widehat{\mathbf{Y}}_0\) est le vecteur des valeurs ajustées selon le modèle utilisant les \(k\) variables explicatives;

5. Le critère \(F\) de Fisher : Notons \(\widehat{\mathbf{Y}}\) le vecteur des valeurs ajustées selon le modèle de RLM complet à \(p\) variables explicatives, et \(\widehat{\mathbf{Y}}_0\) le vecteur des valeurs ajustées selon le modèle réduit à \(p-q\) variables explicatives (\(1\leq q <p\)). Le critère \(F\) du modèle réduit est défini par \[F := \frac{\|\widehat{\mathbf{Y}}_0 -\widehat{\mathbf{Y}}\|^2/q}{\|\mathbf{Y}-\widehat{\mathbf{Y}}\|^2/(n-p-1)}.\]

4.4 Algorithme de recherche exhaustive

1. Construire les \(2^p-1\) modèles;

2. Choisir celui qui optimise un critère donné.

4.5 Algorithme de recherche pas-à-pas

L’approche exhaustive permet de comparer tous les modèles; l’inconvénient est que le temps de calcul devient très important si le nombre de variables est grand;

Lorsque le nombre de variables est grand, on privilégie souvent les méthodes pas-à-pas qui consistent à construire les modèles de façon récursive, en ajoutant/supprimant une variable explicative à chaque étape.

4.6 Sélection par algorithme génétique

1. On l’utilise quand le nombre de variables devient de plus en plus grand et qu’une recherche exhaustive est impossible et une recherche pas à pas peut mener à une solution qui n’est pas tout à fait optimale;

2. Cet algorithme est implémenté dans la fonction glmulti() du package glmulti en spécifiant l’argument method = "g";

3. Cette méthode est donc supposée trouver le meilleur modèle sans avoir besoin de calculer le critère à considérer sur tous les modèles possibles (recherche exhaustive).

5.1 Méthode 1 : de l’ensemble de validation

Cette approche est simple. Elle consiste à diviser de manière aléatoire l’ensemble d’observations en deux parties : un ensemble d’apprentissage (\(\mathcal{A}\)) et un ensemble de validation (ou de test) (\(\mathcal{V}\)). Le modèle est construit à partir des données d’apprentissage, et il est ensuite utilisé pour prédire la variable réponse des données de l’ensemble de validation. L’erreur théorique du modèle considéré, que l’on note \(\mathcal{E}\), est alors estimée par l’écart moyen entre les valeurs prédites et les valeurs observées de la variable réponse des données de l’ensemble de validation : \[\widehat{\mathbf{w}} :=\arg\min_{\mathbf{w}} \, \, \frac{1}{\text{Card}(\mathcal{A})} \sum_{i\in\mathcal{A}} \left(Y_i- f(\mathbf{X}_i, \mathbf{w})\right)^2\] et \[\widehat{\mathcal{E}} := \frac{1}{\text{Card}(\mathcal{V})} \sum_{i\in\mathcal{V}} \left(Y_i-f(\mathbf{X}_i, \widehat{\mathbf{w}})\right)^2. \] Cette méthode est simple, mais elle a au moins deux inconvénients :

1. L’estimation de l’erreur théorique dépend de la partition, et cette estimation peut varier beaucoup si on considère une partition différente;

2. Dans cette approche seulement une partie des données est utilisée pour ajuster le modèle, ce qui peut mener à une sur-estimation de l’erreur.

5.2 Méthode 2 : leave-one-out cross-validation (LOOCV)

Comme la méthode précédente, la méthode LOOCV consiste à diviser l’ensemble des observations en deux parties. Cependant, au lieu de créer deux parties de tailles ``comparables’’, une seule observation \((\mathbf{X}_1,Y_1)\) est utilisée pour la validation et tout le reste est utilisé pour l’apprentissage (pour estimer les paramètres du modèle). Donc le modèle est ajusté sur la base des \((n-1)\) observations restantes, et une prédiction \(\widehat{Y}_1\) est calculée en fonction de \(\mathbf{X}_1\). Plus précisément, si on note \(\widehat{\mathbf{w}}^{(-1)}\) l’estimateur obtenu en utilisant toutes les observations sauf l’observation \((\mathbf{X}_1,Y_1)\), alors \(\widehat{Y}_1:= f(\mathbf{X}_1,\widehat{\mathbf{w}}^{(-1)})\). Puisque l’observation \((\mathbf{X}_1,Y_1)\) n’a pas était utilisée pour ajuster le modèle, alors
\(e_1:=\left(Y_1-\widehat{Y}_1\right)^2= \left(Y_1-f(\mathbf{X}_1,\widehat{\mathbf{w}}^{(-1)})\right)^2\) est une estimation approximativement sans biais de l’erreur théorique, mais elle n’est pas précise puisque elle est basée sur une seule observation. On répète alors la procédure précédente en choisissant \((\mathbf{X}_2,Y_2)\) pour la validation, et les observations restantes pour ajuster le modèle, et on calcule \(e_2:=(Y_2-\widehat{Y}_2)^2=(Y_2-f(\mathbf{X}_2, \widehat{\mathbf{w}}^{(-2)}))^2\). En répétant la procédure \(n\) fois, on obtient \(n\) estimations, \(e_1,e_2,\ldots, e_n\), de l’erreur théorique \(\mathcal{E}\). Le LOOCV donne l’estimation suivante de l’erreur théorique
\[\begin{equation}\label{estimation de l erreur par LOOCV} \widehat{\mathcal{E}} := \frac{1}{n} \sum_{i=1}^n e_i = \frac{1}{n} \sum_{i=1}^n \left(Y_i-f(\mathbf{X}_i, \widehat{\mathbf{w}}^{(-i)})\right)^2. \end{equation}\]
Cette méthode a au moins deux avantages par rapport à la méthode 1 précédente (l’approche de l’ensemble de validation) :

1. En LOOCV, l’estimation des paramètres du modèle est plus précise, car à chaque fois on utilise toutes les données (sauf une), alors que la méthode 1 utilise seulement une partie des données pour l’estimation des paramètres du modèle;

2. Contrairement à la méthode 1 de l’ensemble de validation - qui donne des résultats différents en changeant la répartition des données entre les deux ensembles apprentissage/validation - le LOOCV donne le même résultat s’il est appliqué plusieurs fois sur la même base de données.

Par ailleurs, le LOOCV est coûteux en temps de calcul en général, en particulier si \(n\) est grand, puisque le calcul des estimateurs des paramètres du modèle est répété \(n\) fois. Cependant, ce problème de recalcul peut être évité dans le cas d’un modèle de régression (linéaire en \(\mathbf{w}\)) et une fonction de perte quadratique. En effet, on peut montrer que l’estimation précédente peut s’écrire sous la forme plus simple \[\begin{eqnarray} \widehat{\mathcal{E}} & := & \frac{1}{n} \sum_{i=1}^n \left(Y_i - f(\mathbf{X}_i, \widehat{\mathbf{w}}^{(-i)})\right)^2\nonumber\\ & = & \frac{1}{n} \sum_{i=1}^n \left(\frac{Y_i-f(\mathbf{X}_i, \widehat{\mathbf{w}})}{1-h_i}\right)^2, \end{eqnarray}\] où \(\widehat{\mathbf{w}}=(\mathbb{X}^\top\, \mathbb{X})^{-1}\, \mathbb{X}^\top \mathbf{Y}\) , est l’estimateur des moindres-carrés utilisant toutes les données, et \(h_i\) est le \(i^{eme}\) élément diagonal de la matrice \(\mathbb{X} \, (\mathbb{X}^\top\, \mathbb{X})^{-1}\, \mathbb{X}^\top,\) pour tout \(i=1,\ldots,n.\)

5.3 Méthode 3 : K-fold CV

Cette méthode est une alternative aux deux méthodes précédentes LOOCV et Apprentissage/Validation. Elle consiste à diviser de manière aléatoire le tableau des données en \(K\) groupes (folds) de même effectif (ou presque). Notons \(G_1,\ldots, G_K\) ces groupes. Le premier groupe est considérer comme un ensemble de validation, les \(K-1\) autres groupes sont utilisés pour ajuster le modèle. Une première approximation de l’erreur théorique est donnée par \[e_1 = \frac{1}{\text{Card}(G_1)}\sum_{i\in G_1} \left(Y_i - f(\mathbf{X}_i,\widehat{\mathbf{w}}^{(-1)})\right)^2.\] Cette procédure est répétée \(K\) fois, et à chaque fois, un groupe différent est traité comme ensemble de validation et les autres comme ensemble d’apprentissage; Ceci permet d’avoir \(K\) ``estimations’’ de l’erreur théorique de prévision \[e_k = \frac{1}{\text{Card}(G_k)}\sum_{i\in G_k} \left(Y_i- f(\mathbf{X}_i,\widehat{\mathbf{w}}^{(-k)})\right)^2, \quad \forall k=1,\ldots,K.\] L’estimation de l’erreur théorique par la K-fold CV est calculée en moyennant les valeurs précédentes \[\begin{eqnarray}\label{estimation de l erreur par K-fold CV} \widehat{\mathcal{E}} & := & \frac{1}{K} \sum_{k=1}^K e_k\nonumber\\ & = & \frac{1}{K} \sum_{k=1}^K \frac{1}{\text{Card}(G_k)} \sum_{i\in G_k} \left(Y_i- f(\mathbf{X}_i,\widehat{\mathbf{w}}^{(-k)})\right)^2. \end{eqnarray}\]

6.1 Idée des méthodes de régression pénalisé

Pour réduire la variance de l’estimateur MC (quitte à augmenter un ``peu’’ le biais), on impose une contrainte sur les paramètres \((w_1,\ldots,w_p)\) du type \(\|(w_1,\ldots,w_p)\|_? \leq b\), où \(b>0\) est un seuil donné, et \(\|\cdot\|_?\) est une norme sur \(\mathbb{R}^p\).

L’estimateur MC ``pénalisé’’ sera alors défini par \[\widehat{\mathbf{w}}_{\text{pen}} := \underset{\{\mathbf{w} \in\mathbb{R}^{1+p} \text{ t.q. } \|(w_1,\ldots,w_p)\|_{?} \leq b\}}{\arg\min}\,\, \frac{1}{n} \left\|\mathbf{Y}-\mathbb{X}\mathbf{w}\right\|^2.\] Plusieurs questions se posent alors : quelle norme utiliser pour la contrainte? existence, unicité et calcul de la solution \(\widehat{\mathbf{w}}_{\text{pen}}\) du problème de minimisation précédent? le choix du seuil \(b\)?

6.2 Régression ridge

La régression ridge consiste à minimiser le critère des moindres carrés pénalisé par la norme \(L_2\) du vecteur \((w_1,\ldots,w_p)\). L’estimateur ridge de \(\mathbf{w}:=(w_0,w_1,\ldots,w_p)^\top\) est défini alors comme suit \[\begin{equation}\label{ridge 1} \widehat{\mathbf{w}}^R := \underset{\mathbf{w} \in\mathbb{R}^{1+p} \text{ t.q. } \sum_{j=1}^p w_j^2\leq b} {\arg\min} \,\, \frac{1}{n}\sum_{i=1}^n \left(Y_i-w_0-\sum_{j=1}^p w_j\, X_{i,j} \right)^2 \end{equation}\] ou de façon équivalente \[\begin{equation}\label{ridge 2} \widehat{\mathbf{w}}^R := \underset{\mathbf{w}\in\mathbb{R}^{1+p}}{\arg\min} \,\, \left\{\frac{1}{n}\sum_{i=1}^n \left(Y_i-w_0-\sum_{j=1}^p w_j\, X_{i,j} \right)^2 + \lambda \, \sum_{j=1}^p w_j^2\right\}. \end{equation}\] Les deux définitions précédentes sont équivalentes dans le sense où pour tout \(b>0\) il existe un unique \(\lambda >0\) tels que les deux solutions coïncident.

6.3 Régression Lasso

La régression lasso consiste à minimiser le critère des moindres carrés pénalisé par la norme \(L_1\) du vecteur \((w_1,\ldots,w_p)\).

L’estimateur lasso \(\widehat{\mathbf{w}}^L\) est défini par \[\begin{equation}\label{lasso 1} \widehat{\mathbf{w}}^L := \underset{\mathbf{w} \in\mathbb{R}^{1+p} \text{ t.q. } \sum_{j=1}^p|w_j|\leq b} {\arg\min} \,\, \frac{1}{n}\sum_{i=1}^n \left(Y_i-w_0-\sum_{j=1}^pw_j\, X_{i,j} \right)^2 \end{equation}\] ou de façon équivalente \[\begin{equation}\label{lasso 2} \widehat{\mathbf{w}}^L := \underset{\mathbf{w}\in\mathbb{R}^{1+p}}{\arg\min} \,\, \left\{\frac{1}{n}\sum_{i=1}^n \left(Y_i-w_0-\sum_{j=1}^pw_j\, X_{i,j} \right)^2 + \lambda \, \sum_{j=1}^p |w_j|\right\}. \end{equation}\]

6.4 Régression group-lasso

Dans certains cas, les variables explicatives appartiennent à des groupes de variables prédéfinis; c’est le cas par exemple des variables indicatrices des modalités d’une même variable qualitative (on veut les sélectionner toutes ou les exclure toutes). En présence de \(d\) variables réparties en \(L\) groupes \(\mathbf{X}_1, \ldots, \mathbf{X}_L\) de cardinal \(d_1,\ldots, d_L\), on note \(\mathbf{w}_\ell :=(w_{\ell,1},\ldots,w_{\ell,d_\ell})^\top\), le vecteur des coefficients associé au groupe \(\mathbf{X}_\ell\), \(\ell = 1,\ldots,L\). Les estimateurs group-lasso s’obtiennent en minimisant, en \((w_0,\mathbf{w}_1,\ldots,\mathbf{w}_L)\), le critère suivant \[(w_0,\mathbf{w}_1,\ldots,\mathbf{w}_L) \mapsto \frac{1}{n}\sum_{i=1}^n\left(Y_i-w_0-\sum_{\ell = 1}^L \mathbf{X}_{i,\ell} \,\, \mathbf{w}_\ell \right)^2 + \lambda \sum_{\ell = 1}^L \sqrt{d_\ell}\, \|\mathbf{w}_\ell\|_2.\]

Puisque \(\|\mathbf{w}_\ell\|_2 = 0 \text{ ssi } w_{\ell,1} = \cdots = w_{\ell,d_\ell}=0,\) la méthode group-lasso encourage la mise à zéro des coefficients d’un même groupe.

Pour calculer les estimateurs group-lasso, on peut utiliser la fonction gglasso() du package gglasso.