Chapitre 5 Modèles singuliers, orthogonalité et importance des hypothèses sur les erreurs
5.1 Quand H1-H4 ne sont pas respectées…
L’hypothèse de gaussianité des erreurs est la plus difficile à vérifier en pratique. Les tests classiques de normalité (test de Kolmogorov-Smirnov, Cramer-Von Mises, Anderson-Darling ou de Shapiro-Wilks) demanderaient l’observation des erreurs \(\varepsilon_i\) elles-mêmes ; ils perdent beaucoup de leur puissance quand ils sont appliqués sur les résidus \(\widehat{\varepsilon_i}=Y_i-\widehat{Y_i}\), notamment en raison du fait que ces résidus ne sont pas indépendants. Nous pouvons cependant toujours faire des droites de Henry ou des QQ-plots pour mettre en évidence des écarts évidents. Il n’en reste pas moins que l’hypothèse de gaussianité sera le plus souvent un credo que nous ne pourrons pas vraiment vérifier expérimentalement. Fort heureusement, il existe une théorie asymptotique (donc de grands échantillons) du modèle linéaire qui n’a pas besoin de cette hypothèse. Comme il est dit dans la section 2.1, c’est dans cette optique là qu’il faut réellement penser le modèle linéaire.
5.1.1 Propriétés de l’estimateur des moindres carrés \(\widehat{\theta}\)
Proposition 5.1 Soit \(\widehat{\theta} =(X'X)^{-1}X'Y.\)
- \(\widehat{\theta}\) reste sans biais, \(\mathbb{E}[\widehat{\theta}]=\theta\), sous l’hypothèse H1.
- la matrice de variance-covariance de \(\widehat{\theta}\) reste égale à \(\sigma^2 (X'X)^{-1}\) sous les hypothèses H2 et H3, mais si H1 n’est pas vraie cette propriété a peu d’intérêt.
- \(\widehat{\theta}\) n’est plus un estimateur optimal parmi les estimateurs sans biais, mais il le reste parmi les estimateurs linéaires sans biais sous H1-H3.
- \(\widehat{\theta}\) est gaussien sous H3 et H4. Si H4 n’est pas vraie, alors il tend à être gaussien pour de grands échantillons. On dit qu’il est asymptotiquement gaussien.
5.1.2 Propriétés de l’estimateur des moindres carrés \(\widehat{\sigma}^2\)
Cette étude n’a bien sûr d’intérêt que si \(\sigma^2\) est bien définie ce qui nécessite l’hypothèse H2. Nous considérons \[\widehat{\sigma^2} = \frac{1}{n-k}\|Y-X\widehat{\theta}\|^2 \textrm{ avec } \widehat{\theta} =(X'X)^{-1}X'Y.\]
- Sous les hypothèses H1-H3, \(\widehat{\sigma}^2\) reste un estimateur sans biais de \(\sigma^2\) même si l’hypothèse H4 n’est pas vérifiée : \(\mathbb{E}[\widehat{\sigma}^2]=\sigma^2.\)
- Il est clair que \((n-k)\widehat{\sigma}^2\) ne suit plus une loi \(\sigma^2 \chi^2(n-k)\) dès que l’hypothèse H4 n’est pas vérifiée.
- Nous montrons facilement que sous les hypothèses H1-H3, \(\widehat{\sigma}^2\) converge en probabilité vers \(\sigma^2\) quand le nombre d’observations devient grand, même si l’hypothèse H4 n’est pas vérifiée.
- Enfin, sous les seules hypothèses H1-H3, dès que la loi de \(\varepsilon_i\) admet un moment d’ordre 4, alors \(\widehat{\sigma}^2\) converge à la vitesse \(\sqrt{n}\) vers \(\sigma^2\) mais sa vitesse exacte de convergence dépend du type de loi, plus précisément du coefficient de Kurtosis (Azaïs and Bardet 2005).
5.1.3 Modèles avec corrélations
Il est possible de modéliser des corrélations entre erreurs, par exemple en supposant que ces erreurs sont issues d’un processus ARMA, ce qui permet de ne plus avoir besoin de l’hypothèse H3, voir Guyon (2001).
Il est également possible de modéliser les liaisons par des modèles à effets aléatoires et poser un modèle mixte, voir Pinheiro et Bates (Pinheiro and Bates 2006).
5.2 Modèles singuliers
Nous nous sommes jusqu’à présent cantonnés à l’étude des modèles linéaires réguliers. Or certains modèles ne peuvent être paramétrés de façon régulière : ils sont naturellement sur-paramétrés. Un exemple simple est celui du modèle additif en analyse de la variance à 2 facteurs. Considérons le cas où les 2 facteurs ont chacun 2 niveaux et que les 4 combinaisons sont observées une fois et une seule. On a donc, avec les notations vues précédemment : \[ Y_{i,j} = \mu + a-i + b_j + \varepsilon_{i,j},\ i\in\{1,2\},\ j\in\{1,2\}. \]
Le vecteur \(\theta=(\mu,a_1,a_2,b_1,b_2)'\) et la matrice \(X\) du modèle vaut : \[ X=\left( \begin{array}{ccccc} 1 & 1 & 0 & 1 & 0 \\ 1 & 1 & 0 & 0 & 1 \\ 1 & 0 & 1 & 1 & 0 \\ 1 & 0 & 1 & 0 & 1 \end{array} \right). \] Nous remarquons que tout vecteur de la forme \((\alpha+\beta,-\alpha,-\alpha, -\beta,-\beta)\) donne le vecteur nul lorsqu’il est multiplié par la matrice \(X\). Les valeurs \(\mu, \, a_i, \, b_i\) pour \(i=1\) ou 2 ne sont donc pas identifiables de manière unique. Le modèle est en fait sur-paramétré : nous avons 5 paramètres inconnus pour seulement 4 observations.
Definition 5.1 Le modèle est dit singulier quand la matrice \(X\) est non injective, c’est-à-dire s’il existe \(\theta\neq 0_k\) tel que \(X\theta =0_n\).
Nous rappelons que \(Ker(X)=\{u \in \mathbb{R}^k; \, Xu=0_{n}\}\) désigne le noyau de \(X\). Nous pouvons faire deux remarques :
- \(X\widehat{\theta}\) reste unique, puisque c’est la projection de \(Y\) sur \(Im(X)\).
- \(\widehat{\theta}\) ne peut être unique puisque si \(\widehat{\theta}\) est solution et si \(u \in Ker(X)\) alors \(\widehat{\theta} + u\) est encore solution.
Si \(X\) n’est pas régulière, alors la matrice \(X'X\) n’est pas inversible. Pour contourner ce problème, nous définissons alors un inverse généralisé de \((X'X)\).
Definition 5.2 Soit \(M\) une matrice. Alors la matrice \(M^-\) est une matrice inverse généralisée de \(M\) si \(MM^-M=M\).
Cette construction est toujours possible. En effet, \((X'X)\) définit une application bijective de \(Ker(X)^{\perp}\) sur lui-même. Il suffit donc simplement de négliger la partie contenue dans le noyau : on prend l’inverse sur \(Ker(X)^{\perp}\), complété arbitrairement sur \(Ker(X)\). La définition de \((X'X)^-\) est donc loin d’être unique ! Il est alors possible de généraliser les résultats du cas régulier.
Proposition 5.2 Si \((X'X)^-\) est une matrice inverse généralisée de \(X'X\) alors \(\widehat{\theta}=(X'X)^-X'Y\) est une solution des équations normales : \[(X'X)\widehat{\theta} = X'Y.\]
Proof. On commence par remarquer que \[\forall \omega\in\mathbb{R}^k,\, <X\omega,P_{[X]^\perp} Y> = <\omega, X'P_{[X]^\perp} Y> = 0\] donc \[ X'Y = X'P_{[X]}Y + X' P_{[X]^{\perp}} Y = X'P_{[X]}Y. \] Ainsi, \(\exists u\in\mathbb{R}^k,\ X'Y = X'Xu\). Finalement, \[ (X'X)\widehat\theta = (X'X) (X'X)^{-}X'Y = (X'X) (X'X)^{-}X'Xu = X'Xu = X'Y. \]
Remark. Cet estimateur n’est pas unique et dépend de la définition choisie pour \((X'X)^-\). Par contre, le vecteur \(X\widehat{\theta}\) reste unique, même si la matrice \(X\) est singulière. Ce vecteur correspond en effet à la projection orthogonale de \(Y\) sur \(Im(X)\).
En règle générale, nous préférons lever l’indétermination sur \(\hat\theta\) en fixant des contraintes, souvent afin de donner un sens plus intuitif à \(\theta\).
5.2.1 Contraintes d’identifiabilité
Proposition 5.3 Supposons la matrice \(X\) singulière de rang \(r<k\) de sorte qu’il y ait \(k-r\) paramètres redondants. Soit \(M\) une matrice à \(k-r\) lignes et \(k\) colonnes, supposée de rang \(k-r\) et telle que : \[Ker(M) \cap Ker(X)= \lbrace 0_{k} \rbrace.\] Alors,
- la matrice \((X'X+M'M)\) est inversible et son inverse est une matrice inverse généralisée de \(X'X\)
- le vecteur \(\widehat{\theta} = (X'X+M'M)^{-1}X'Y\) est l’unique solution du système \(\displaystyle \left\lbrace \begin{array}{c} X'X\alpha = X'Y \\ M\alpha= 0_{k-r}. \end{array} \right.\)
Exercise 5.1 L’objectif est de démontrer la proposition 5.3.
Pour montrer que \(X'X+M'M\) est inversible : montrez que la matrice \[ A=\left(\begin{array}{c} X\\ M \end{array}\right)\in\mathcal M_{n+k-r,k}(\mathbb{R}) \] est injective et donc \(A'A\) est inversible.
Considérez le problème de minimisation suivant : \[ g:\alpha\mapsto \|Y-X\alpha\|^2 + \|M\alpha\|^2 . \] Ecrivez \(g(\alpha)\) sous la forme \(g(\alpha)=\|\tilde{Y} - A \alpha\|^2\) avec \(\tilde {Y}\) à préciser. Déduisez-en que \(\widehat{\theta}\) est solution du système \(\displaystyle \left\lbrace \begin{array}{c} X'X\alpha = X'Y \\ M\alpha= 0_{k-r}. \end{array} \right.\)
Montrez l’unicité de la solution
Le choix de la contrainte n’est pas toujours évident. Par ailleurs, pour chaque contrainte \(M\), nous aurons un estimateur correspondant ce qui est parfois gênant.
Example 5.1 Prenons l’exemple de l’analyse de variance à un facteur avec effet différentiel : nous supposons pour simplifier que \(I=4\). Le modèle s’écrit donc de la façon suivante : \[Y_{i,j}=\mu + \alpha_i +\varepsilon_{ij} \mbox{ pour } i=1,\cdots,4 \mbox{ et } j=1.\] La matrice \(X\) associée au modèle est : \[ X= \left( \begin{array}{ccccc} 1 & 1 & 0 & 0 & 0 \\ 1 & 0 & 1 & 0 & 0 \\ 1 & 0 & 0 & 1 & 0 \\ 1 & 0 & 0 & 0 & 1 \end{array} \right).\] Il nous faut poser une contrainte (dite d’identifiabilité) sur le vecteur \(\theta\) au travers du choix d’une matrice à \(1\) ligne et \(k\) colonnes.
Si on considère \(M= (0 \ 1 \ 1 \ 1 \ 1)\), la contrainte correspondante est \[ M\theta = 0 \Leftrightarrow \alpha_1 + \alpha_2 + \alpha_3 + \alpha_4= 0.\] On impose donc que la somme des effets différentiels est nulle.
Si on considère \(M=(0,1,0,0,0)\), la contrainte correspondante est \(M\theta = \alpha_1=0\). On impose donc que la première modalité est la référence.
5.2.2 Fonctions estimables et contrastes
En présence d’une matrice singulière, il est donc toujours possible de construire un estimateur. Qu’en est-il des tests ? En particulier, ces contraintes sont-elles systématiquement nécessaires ?
La plupart des quantités que nous avons voulu tester sont des fonctions de \(\theta\) qui ne dépendent pas de la solution particulière des équations normales, c’est-à-dire du type de contraintes d’identifiabilité choisi. Ces fonctions sont appelées estimables car elles sont intrinsèques.
Definition 5.3 Une combinaison linéaire \(C\theta\) est dite fonction estimable (de paramètre \(\theta\)) si elle ne dépend pas du choix particulier d’une solution des équations normales. On caractérise ces fonctions comme étant celles qui s’écrivent \(C\theta=DX\theta\) où \(D\) est une matrice de plein rang.
Definition 5.4 On appelle contraste une fonction estimable \(C\theta\) telle que \(C \mathbb{1}=0\), où \(\un\) désigne le vecteur unité.
En analyse de variance, la plupart des combinaisons linéaires que l’on teste sont en fait des contrastes (cf chapitre 7). Dans l’exemple précédent, \(\alpha_1 - \alpha_2\) est un contraste.
5.3 Orthogonalité
5.3.1 Orthogonalité pour les modèles réguliers
L’orthogonalité est une notion qui peut notablement simplifier la résolution et la compréhension d’un modèle linéaire. Un modèle linéaire admet le plus souvent une décomposition naturelle des paramètres \(\theta\) (cf exemple ci-dessous) et conséquemment une décomposition de la matrice \(X\) associée au modèle. On va s’intéresser ici à l’orthogonalité éventuelle des différents espaces associés à cette décomposition (l’orthogonalité sera toujours comprise par la suite au sens d’orthogonalité liée au produit scalaire euclidien usuel). Le problème sera plus ou moins délicat suivant que le modèle est régulier ou non. En premier lieu, illustrons par deux exemples ce que l’on entend par décomposition des paramètres.
Example 5.2 Soit le modèle de régression linéaire multiple sur trois variables \(x^{(1)}, \, x^{(2)}\) et \(x^{(3)}\) : \[Y_i=\mu + \theta_1x_i^{(1)}+\theta_2 x^{(2)}_i +\theta_3 x^{(3)}_i + \varepsilon_i, i=1, \cdots, n>4.\] Le vecteur \(\theta\) comprend 4 coordonnées : \(\mu, \, \theta_1, \, \theta_2, \theta_3\) et la matrice \(X\) quatre colonnes. Assez naturellement ici, on peut considérer la décomposition, plus précisément on parlera par la suite de partition en quatre éléments. La partition de la matrice revient alors à l’écrire comme concaténation de 4 vecteurs colonnes. L’orthogonalité de la partition correspondra alors strictement à l’orthogonalité des 4 droites vectorielles : \([\mathbb{1}], \, [x^{(1)}], \, [x^{(2)}]\) et \([x^{(3)}]\).
Example 5.3 Soit le modèle de régression quadratique sur \(x^{(1)}\) et \(x^{(2)}\) : \[Y_i=\mu + \theta_1x^{(1)}_i+\theta_2 x^{(2)}_i +\gamma_1 \left(x^{(1)}_i\right)^2 +\gamma_2\left(x^{(2)}_i\right)^2+ \delta x^{(1)}_i x^{(2)}_i + \varepsilon_i, i=1, \cdots, n>6.\] Ici plutôt que de demander comme précédemment l’orthogonalité de chacun des régresseurs (ce qui serait beaucoup demander), on peut définir la partition naturelle correspondant à :
- la constante \(\mu\) ;
- les effets linéaires \(\theta_1, \, \theta_2\) ;
- les effets carrés \(\gamma_1, \, \gamma_2\) ;
- l’effet produit \(\delta\).
L’orthogonalité de la partition est alors définie comme l’orthogonalité des sous-espaces vectoriels : \([\mathbb{1}], \, [(x^{(1)}, \, x^{(2)})], \, \left[\left((x^{(1)})^2, \, (x^{(2)})^2\right)\right]\) et \([x^{(1)} x^{(2)}]\).
En conséquence, on voit bien, à partir de ces deux exemples, qu’il faudra parler de modèle avec partition orthogonale plutôt que de modèle orthogonal.
Formalisons ces exemples dans une définition.
Definition 5.5 Soit un modèle linéaire général régulier \(Y=X\theta + \varepsilon.\) Considérons une partition en \(m\) termes de \(X\) et de \(\theta\), soit \[Y= X_1\theta_1 + \cdots +X_m\theta_m +\varepsilon,\] où la matrice \(X_j\) est une matrice de taille \((n,k_j)\) et \(\theta_j \in \mathbb{R}^{k_j}\) avec \(k_j \in \{1,\cdots,k\}\) pour \(j=1,\cdots,m\) et avec \(\sum_{j=1}^{m} k_j=k\). On dit que cette partition est orthogonale si les sous-espaces vectoriels de \(\mathbb{R}^n\), \([X_1], \cdots, [X_m],\) sont orthogonaux.
Une conséquence simple de l’orthogonalité d’un modèle linéaire est que la matrice d’information \(X'X\) a une structure bloc diagonale, chaque bloc étant associé à chaque élément de la partition.
Le plus souvent, la partition du vecteur de paramètres \(\theta\) en différents effets vient
- en régression, des différentes variables ;
- en analyse de la variance, des décompositions en interactions.
L’orthogonalité donne aux modèles statistiques les deux propriétés suivantes :
Proposition 5.4 Soit un modèle linéaire régulier muni d’une partition orthogonale : \[Y= X_1\theta_1 + \cdots +X_m\theta_m +\varepsilon.\] Alors
- les estimateurs des moindres carrés des différents effets \(\widehat{\theta}_1, \dots, \widehat\theta_m\) sont non-corrélés et indépendants sous l’hypothèse gaussienne.
- pour \(l=1,\cdots, m\), l’expression de l’estimateur \(\widehat{\theta}_l\) ne dépend pas de la présence ou non des autres termes \(\theta_j\) dans le modèle.
L’orthogonalité apporte une simplification des calculs : elle permet d’obtenir facilement une expression explicite des estimateurs. Par ailleurs, elle donne une indépendance approximative entre les tests des différents effets. Les tests portant sur des effets orthogonaux ne sont liés que par l’estimation du \(\sigma^2\).
5.3.2 Orthogonalité pour les modèles non-réguliers
Lorsque le modèle est singulier, il est nécessaire de rajouter des contraintes. Il est alors raisonnable d’effectuer cette démarche en tenant compte de la partition, i.e. \(C_j \theta_j =0\) où \(X_j|_{Ker(C_j)}\) sont injectives.
Definition 5.6 Soit la partition suivante d’un modèle linéaire \[Y=X_1\theta_1+\cdots+X_m\theta_m +\varepsilon.\] Soit un système de contraintes \(C_1\theta_1=0, \cdots, C_m \theta_m=0\) qui rendent le modèle identifiable. On dit que ces contraintes rendent la partition orthogonale si les sous-espaces vectoriels \[V_j=\left\{X_j \theta_j; \, \theta_j \in Ker(C_j)\right\}, \, j=1, \cdots, m\] sont orthogonaux.
Cette notion est proche du cas régulier. Cependant, la notion d’orthogonalité dépend des contraintes choisies. L’idée sera en général de choisir des contraintes qui rendent le modèle orthogonal. On verra que cette définition prend tout son sens avec l’exemple incontournable du modèle d’analyse de la variance à deux facteurs croisés (cf chapitre 7).
5.4 En résumé
Dans ce chapitre, il est attendu que vous ayez compris
- la problématique de l’estimation des paramètres pour un modèle linéaire singulier
- l’intérêt d’avoir l’orthogonalité dans un modèle linéaire
Les résultats énoncés dans ce chapitre ne sont pas à connaitre. Il faudra savoir les mettre en application dans le cadre de l’ANOVA (voir Chapitre 7) et de l’ANCOVA (voir Chapitre 8).
References
Azaïs, Jean-Marc, and Jean-Marc Bardet. 2005. Le Modèle Linéaire Par L’exemple: Régression, Analyse de La Variance et Plans d’expériences Illustrés Avec R, Sas et Splus. Dunod.
Guyon, X. 2001. “Modele Linéaire et économétrie.” Ellipse, Paris.
Pinheiro, José, and Douglas Bates. 2006. Mixed-Effects Models in S and S-Plus. Springer Science & Business Media.