Modèle de régression linéaire simple, moindres carrés, tests et intervalles sur les coefficients, intervalles de confiance et de prédiction, coefficient de détermination, analyse des résidus, test F et analyse de la variance à un facteur.
Objectifs du chapitre
À la fin de ce chapitre, vous serez capable de:
formuler le modèle de régression linéaire simple Yi=β0+β1xi+εi à erreurs gaussiennes, en énoncer les quatre hypothèses et dire lesquelles servent à quoi;
calculer les estimateurs des moindres carrés, démontrer qu'ils sont sans biais, calculer leurs variances et estimer σ2 avec le bon nombre de degrés de liberté;
construire des intervalles de confiance et des tests de Student sur la pente et l'ordonnée à l'origine, et distinguer l'intervalle de confiance d'une droite moyenne de l'intervalle de prédiction d'une nouvelle observation;
décomposer la variabilité totale en SCT=SCE+SCR, interpréter le coefficient de détermination R2 et dresser la table d'analyse de la variance d'une régression;
lire un graphique des résidus et y reconnaître une courbure, une variance non constante ou un point aberrant;
comparer k moyennes par l'analyse de la variance à un facteur, et montrer que pour k=2 elle redonne exactement le test de Student à variance combinée du chapitre 12.
Du nuage de points au modèle statistique
Ce que les chapitres précédents ont déjà fait
Trois chapitres ont préparé celui-ci sans le dire. Le chapitre 9 a décrit une série double, l'essai d'usure: dix pièces prélevées d'heure en heure après l'affûtage de l'outil, dont l'écart du diamètre au nominal croît avec la durée d'usinage, avec un coefficient de corrélation r=0,948. Le chapitre 10 a ajusté à ce nuage la droite des moindres carrésy^=−14+4x (en micromètres, x en heures), et montré qu'elle est l'estimation du maximum de vraisemblance lorsque les erreurs sont gaussiennes. Le chapitre 12 a appris à décider: un test, un risque, une valeur p.
Il reste à assembler ces pièces. Une droite ajustée sur dix points n'est qu'une estimation: avec dix autres pièces, on aurait trouvé une autre pente. De combien la pente estimée peut-elle se tromper? La dérive est-elle réelle, ou la pente de 4 µm/h pourrait-elle être un artefact de l'aléa? Quel diamètre attendre de la pièce qui sortira à la huitième heure, et avec quelle marge? Ces questions sont celles des chapitres 11 et 12, posées à un modèle dont l'espérance dépend d'une variable.
Le modèle tient en quatre hypothèses, et chacune a son rôle:
linéarité: E[Yi]=β0+β1xi — sans elle, la droite estime autre chose que ce que l'on croit;
indépendance des erreurs — sans elle, les variances calculées plus bas sont fausses, en général trop optimistes;
homoscédasticité: toutes les erreurs ont la même variance σ2, quelle que soit xi;
normalité des erreurs — elle ne sert ni à calculer les estimateurs ni à établir leurs espérances et variances, mais elle donne les lois exactes de Student et de Fisher des tests et des intervalles.
Les xi ne sont pas aléatoires. Dans l'essai d'usure, c'est l'expérimentateur qui a choisi de prélever une pièce à la fin de chaque heure: les durées sont fixées, seuls les diamètres sont aléatoires. Lorsque les deux variables sont observées au hasard — la taille et le poids de personnes tirées au sort —, on raisonne conditionnellement aux valeurs de x observées, et tous les résultats de ce chapitre restent valables sous cette lecture. Mais l'asymétrie demeure: régresser Y sur x et x sur Y ne donne pas la même droite, parce que ce ne sont pas les mêmes écarts que l'on minimise.
Question 13.1
Laquelle de ces affirmations ne fait pas partie des hypothèses du modèle de régression linéaire simple (13.1)?
Les estimateurs des moindres carrés
Rappel et écriture linéaire
Le théorème 10.7 du chapitre 10 donne les estimateurs
que nous écrivons avec des majuscules, puisque ce sont des variables aléatoires: b0=−14 et b1=4 en sont les estimations sur l'essai d'usure. Une remarque simple ouvre tout le reste. Comme ∑i(xi−xˉ)=0, le terme en Yˉ disparaît du numérateur, et
B1=i=1∑nciYi,ci=Sxxxi−xˉ.(13.3)
La pente estimée est une combinaison linéaire des réponses, à coefficients connus, qui vérifient trois identités immédiates:
i∑ci=0,i∑cixi=1,i∑ci2=Sxx1.
(La deuxième vient de ∑i(xi−xˉ)xi=∑i(xi−xˉ)2=Sxx.) Le chapitre 7 sait tout calculer sur une combinaison linéaire de variables indépendantes.
Explorateur 13.1 · Une droite candidate contre la droite des moindres carrés
Déplacez l'ordonnée à l'origine et la pente de la droite en trait plein: les segments verticaux sont les écarts des dix pièces de l'essai d'usure à cette droite, et le premier readout additionne leurs carrés. Le minimum, 148, ne bouge jamais; l'excès se décompose exactement en un terme de niveau et un terme de pente, et il ne s'annule que lorsque la droite candidate se confond avec la droite tiretée.
Ordonnée à l'origine b0 (µm)−4,0
Pente b1 (µm/h)2,5
Somme des carrés S(b0, b1)
364,3
Minimum Smin (moindres carrés)
148,0
Terme de niveau n(ȳ − b0 − b1x̄)²
30,6
Terme de pente Sxx(b1 − 4)²
185,6
L'explorateur rappelle ce que «moindres carrés» veut dire. Partez de la droite proposée, de pente 2,5: la somme des carrés des écarts verticaux vaut 364,3. Rapprochez la pente de 4 et l'ordonnée à l'origine de −14, et regardez les deux termes d'excès s'effondrer l'un après l'autre; le minimum, 148, ne bouge jamais. La décomposition affichée est exactement celle de la démonstration par complétion des carrés du chapitre 10: un terme de niveau, nul dès que la droite passe par le point moyen, et un terme de pente, nul dès que la pente vaut Sxy/Sxx.
Espérances et variances
Démonstration. Par (13.3) et la linéarité de l'espérance, E[B1]=∑ici(β0+β1xi)=β0∑ici+β1∑icixi=β1. Les Yi étant non corrélées et de variance σ2, Var(B1)=∑ici2σ2=σ2/Sxx. Ensuite E[Yˉ]=β0+β1xˉ, donc E[B0]=β0+β1xˉ−β1xˉ=β0. La clé du reste est que Yˉ et B1 sont non corrélées: par bilinéarité de la covariance (chapitre 7),
Cov(Yˉ,B1)=i∑n1ciσ2=nσ2i∑ci=0.
Alors Var(B0)=Var(Yˉ)+xˉ2Var(B1)=σ2/n+xˉ2σ2/Sxx, et Cov(B0,B1)=Cov(Yˉ,B1)−xˉVar(B1)=−xˉσ2/Sxx. □
La formule Var(B1)=σ2/Sxx est un résultat de plan d'expérience. La précision de la pente ne dépend pas seulement du nombre de points, mais de leur étalement en x: dix pièces prélevées entre la première et la dixième heure donnent Sxx=82,5; dix pièces prélevées entre la quatrième et la sixième heure donneraient une pente bien plus incertaine. Pour mesurer une dérive, on la laisse se développer. L'exercice 13.5 pousse ce raisonnement jusqu'à ses limites.
Le signe de Cov(B0,B1) se comprend sur le nuage: si xˉ>0 et que la pente est surestimée, la droite, qui passe toujours par le point moyen, pivote autour de lui et coupe l'axe vertical trop bas.
Le meilleur estimateur linéaire
Les moindres carrés ont été choisis pour leur commodité. Ils sont aussi les meilleurs dans une classe précise, et sans hypothèse de normalité.
Démonstration.E[B~]=β0∑iai+β1∑iaixi vaut β1 pour tout (β0,β1) si et seulement si ∑iai=0 et ∑iaixi=1. Pour un tel estimateur, ∑iaici=(∑iaixi−xˉ∑iai)/Sxx=1/Sxx. L'inégalité de Cauchy–Schwarz dans Rn donne alors
d'où Var(B~)=σ2∑iai2≥σ2/Sxx=Var(B1), avec égalité si et seulement si (ai) est proportionnel à (ci), donc égal à (ci) à cause de la contrainte. □
Le résultat se généralise à β0 et à toute combinaison linéaire des coefficients; on le résume par l'acronyme anglais BLUE (best linear unbiased estimator). Il dit pourquoi les moindres carrés restent le choix par défaut même quand la normalité est douteuse — et il dit aussi sa limite: il ne compare les moindres carrés qu'aux estimateurs linéaires. Face à des erreurs à queues épaisses, un estimateur non linéaire et robuste (fondé sur les valeurs absolues, comme la médiane) peut faire beaucoup mieux; c'est la remarque du chapitre 10 sur les erreurs de Laplace.
Estimer la variance des erreurs
Les variances (13.4) et (13.5) contiennent σ2, inconnu. Il s'estime sur les résidus
ei=Yi−Y^i,Y^i=B0+B1xi,
écarts des observations à la droite estimée — à ne pas confondre avec les erreurs εi, écarts à la droite vraie, qui ne sont jamais observées. La somme de leurs carrés est la somme des carrés résiduelle, SCR=∑iei2; c'est le minimum Smin du chapitre 10.
Les résidus ne sont pas libres: les deux équations normales imposent ∑iei=0 et ∑ixiei=0. Des n résidus, n−2 seulement sont libres — on a dépensé deux degrés de liberté pour estimer deux coefficients —, et le bon diviseur est n−2.
Démonstration du point 1. Posons SYY=∑i(Yi−Yˉ)2. Comme ei=(Yi−Yˉ)−B1(xi−xˉ), le développement du carré et ∑i(xi−xˉ)(Yi−Yˉ)=B1Sxx donnent
SCR=SYY−2B1⋅B1Sxx+B12Sxx=SYY−B12Sxx.
Or Yi−Yˉ=β1(xi−xˉ)+(εi−εˉ), donc SYY=β12Sxx+2β1∑i(xi−xˉ)(εi−εˉ)+∑i(εi−εˉ)2. Le terme du milieu est d'espérance nulle, et le dernier a pour espérance (n−1)σ2 par le théorème 9.4 appliqué aux erreurs. Donc E[SYY]=β12Sxx+(n−1)σ2. Par ailleurs E[B12]=Var(B1)+β12=σ2/Sxx+β12, d'où E[B12Sxx]=σ2+β12Sxx. La différence vaut (n−2)σ2. □
Le point 2 est immédiat: B1 et B0 sont des combinaisons linéaires des Yi, qui sont des normales indépendantes, donc elles sont normales (stabilité de la loi normale, chapitre 7), et le théorème 13.1 donne leurs paramètres.
Le point 3 est admis, pour la même raison que le théorème 9.5: il demande une rotation orthogonale de Rn. L'idée, elle, est celle du chapitre 9. Le vecteur des ajustements (Y^1,…,Y^n) est la projection orthogonale du vecteur des observations sur le plan engendré par (1,…,1) et (x1,…,xn); les résidus sont la projection sur l'espace orthogonal à ce plan, de dimension n−2. Un vecteur gaussien de composantes indépendantes et de même variance a des projections indépendantes sur deux sous-espaces orthogonaux, et le carré de la norme de la seconde, divisé par σ2, est un khi-deux à autant de degrés de liberté que la dimension: n−2. Avec un seul paramètre de position, le chapitre 9 projetait sur une droite et laissait n−1 dimensions; ici, deux paramètres, un plan, n−2 dimensions.
Intervalles de confiance et tests sur les coefficients
Le pivot de Student
Démonstration. C'est le raisonnement du théorème 9.6, mot pour mot. Z=(B1−β1)/(σ/Sxx)∼N(0,1) par le point 2 du théorème 13.3, C=(n−2)S2/σ2∼χn−22 par le point 3, et Z et C sont indépendantes. Donc T1=Z/C/(n−2) suit tn−2 par définition (9.13). Même chose pour T0. □
Tout s'en déduit comme aux chapitres 11 et 12. L'intervalle de confiance de niveau 1−α pour la pente est
b1±tn−2;1−α/2Sxxs,(13.8)
et le test de H0:β1=0 — «la variable explicative n'a aucun effet linéaire sur la réponse» — rejette au seuil α lorsque
∣tobs∣=s/Sxx∣b1∣≥tn−2;1−α/2.(13.9)
Par la dualité du chapitre 12, ce test rejette exactement lorsque l'intervalle (13.8) exclut 0.
Question 13.2
Sur l'essai d'usure, quelle est la demi-largeur de l'intervalle de confiance à 99% pour la pente β1, en µm/h? On donne t8;0,995=3,3554.
Prévoir: intervalle de confiance et intervalle de prédiction
Deux questions de prévision se posent en un point x0, et elles n'ont pas la même réponse.
Quel est le diamètre moyen des pièces produites à la huitième heure? C'est l'estimation de μ0=β0+β1x0, un paramètre: la hauteur de la droite vraie en x0.
Quel diamètre aura la pièce qui sortira à la huitième heure du prochain poste? C'est la prévision d'une variable aléatoireY0=β0+β1x0+ε0, qui ajoute à la droite vraie l'erreur propre de cette pièce.
L'estimation est la même dans les deux cas, Y^0=B0+B1x0; l'incertitude ne l'est pas.
Démonstration. Écrivons Y^0=Yˉ+B1(x0−xˉ), ce qui découle de B0=Yˉ−B1xˉ. Comme Yˉ et B1 sont non corrélées (démonstration du théorème 13.1),
Var(Y^0)=nσ2+(x0−xˉ)2Sxxσ2.
Y^0 est normale d'espérance μ0, indépendante de S2; le raisonnement du théorème 13.4 donne (Y^0−μ0)/(S1/n+(x0−xˉ)2/Sxx)∼tn−2, d'où (13.10). Pour le point 2, Y0 ne dépend que de ε0, indépendante des n premières observations, donc de Y^0: Var(Y0−Y^0)=σ2+Var(Y^0), et Y0−Y^0 est normale centrée, indépendante de S2. Le même pivot de Student conclut. □
Deux lectures. Le «1» sous la racine change tout: quand n devient grand, l'intervalle de confiance se resserre vers zéro — on finit par connaître la droite exactement —, tandis que l'intervalle de prédiction tend vers ±z1−α/2σ: aucune quantité de données ne supprime la dispersion propre de la prochaine pièce. Les deux intervalles sont les plus étroits en x0=xˉ et s'élargissent comme une hyperbole quand on s'éloigne du centre des données: la droite est mieux connue au point moyen, où elle passe toujours, que loin de lui, où une petite erreur de pente se paie cher.
Figure 13.1. Essai d'usure (données fictives): les dix pièces (points), la droite des moindres carrés ŷ = −14 + 4x (trait plein), la bande de confiance à 95 % de la droite moyenne (aplat) et les limites de la bande de prédiction à 95 % d'une nouvelle pièce (tirets bleus). Les deux bandes sont les plus étroites au point moyen, x = 5,5 h, et s'évasent de part et d'autre. La zone grisée est le domaine observé, de 1 à 10 h; au-delà, jusqu'à 15 h, les bandes ne valent que si la dérive reste linéaire, ce que les données ne peuvent pas dire.
Question 13.3
On répète l'essai d'usure avec un nombre n de pièces de plus en plus grand, les durées restant réparties de la même façon entre 1 et 10 heures. Que deviennent, en x0=8 h, les demi-largeurs de l'intervalle de confiance (13.10) et de l'intervalle de prédiction (13.11)?
Décomposition de la variabilité et coefficient de détermination
Trois sommes de carrés
La réponse y varie d'une pièce à l'autre. Une partie de cette variation est expliquée par la droite — les pièces des dernières heures sont plus grosses parce que l'outil est plus usé —, le reste ne l'est pas. Cette intuition est une identité exacte.
Démonstration. Écrivons yi−yˉ=(y^i−yˉ)+ei et élevons au carré. Le double produit est nul, car y^i−yˉ=b1(xi−xˉ) et
i∑(y^i−yˉ)ei=b1(i∑xiei−xˉi∑ei)=0
par les deux équations normales. Ensuite SCE=∑ib12(xi−xˉ)2=b12Sxx=Sxy2/Sxx, et Sxy2/(SxxSyy)=r2 donne SCE=r2Syy=r2SCT. □
Géométriquement, (13.12) est le théorème de Pythagore dans Rn: le vecteur des écarts à la moyenne est la somme du vecteur des ajustements centrés et du vecteur des résidus, et ces deux vecteurs sont orthogonaux — c'est ce que disent les équations normales.
Sur l'essai d'usure, SCT=Syy=1468, SCE=42×82,5=1320 et SCR=148; on vérifie 1320+148=1468, et
R2=14681320=0,8992=0,94832.
La durée d'usinage explique 90% de la variabilité des écarts au nominal; les 10% restants sont la dispersion propre des pièces.
Le test F de la régression et la table d'analyse de la variance
La décomposition (13.12) conduit à un second test de H0:β1=0. Si la pente est nulle, la droite ne devrait rien expliquer de plus que le hasard: SCE devrait être du même ordre que la variance des erreurs. On compare donc SCE, qui a un degré de liberté, à SCR/(n−2), qui estime σ2.
Démonstration. Par le théorème 13.6, SCE=B12Sxx, donc F=B12Sxx/S2=T2. Sous H0, T∼tn−2 par le théorème 13.4, et le carré d'une variable tν est une variable F1,ν: si T=Z/C/ν, alors T2=(Z2/1)/(C/ν), avec Z2∼χ12 indépendante de C — c'est la relation F1,ν=tν2 du chapitre 9. Enfin {T2≥c2}={∣T∣≥c}, et F1,ν;1−α=tν;1−α/22. □
En régression simple, le test F n'apporte donc rien de neuf; il prend tout son sens en régression multiple, où il teste d'un coup la nullité de plusieurs pentes, et en analyse de la variance, où il compare plusieurs groupes. Mais il s'accompagne d'une présentation standard, la table d'analyse de la variance, que tous les logiciels impriment et qu'il faut savoir lire.
Source
Somme de carrés
Degrés de liberté
Carré moyen
F
Régression
SCE
1
CME=SCE/1
CME/CMR
Résiduelle
SCR
n−2
CMR=SCR/(n−2)=s2
Totale
SCT
n−1
Les degrés de liberté s'additionnent comme les sommes de carrés: (n−1)=1+(n−2). Le carré moyen d'une ligne est sa somme de carrés divisée par ses degrés de liberté, et le carré moyen résiduel est l'estimation s2 de σ2.
Question 13.4
Une régression linéaire simple sur n=12 points donne SCE=45 et SCR=15. Que vaut la statistique F du test de H0:β1=0?
Analyse des résidus
Toute l'inférence qui précède repose sur les quatre hypothèses du modèle, et aucune ne se vérifie sur R2, sur t ou sur F. Elles se vérifient sur les résidus, qui sont les seuls témoins observables des erreurs. Si le modèle est juste, les résidus doivent ressembler à un bruit sans structure, centré sur zéro, d'amplitude constante.
Le graphique de base porte les résidus ei en fonction de xi (ou des valeurs ajustées y^i, ce qui revient au même en régression simple). On y cherche trois défauts.
Une courbure: les résidus dessinent un U ou un arc. La relation n'est pas linéaire; la droite est trop haute au centre et trop basse aux bords, ou l'inverse. Remède: transformer une variable (logarithme, racine), ajouter un terme en x2, ou restreindre le domaine.
Un entonnoir: la dispersion des résidus croît (ou décroît) avec x. L'homoscédasticité est violée; les erreurs types calculées sont fausses, trop petites là où la dispersion est grande. Remède: transformer la réponse, ou pondérer les observations.
Un point isolé: un résidu très grand devant les autres. On le repère sur le résidu réduitei/s, qui dépasse rarement 2 en valeur absolue si les erreurs sont normales. Un tel point n'est pas à supprimer d'office — c'est peut-être la seule pièce qui dit quelque chose —, mais à examiner: erreur de saisie, pièce d'un autre lot, phénomène réel? Un point éloigné en x est encore plus dangereux, car il tire la droite à lui: c'est le quatrième panneau du quartet d'Anscombe.
Deux graphiques complètent le premier: les résidus dans l'ordre de prélèvement, qui révèlent une dépendance dans le temps (une suite de résidus de même signe), et la droite de Henry des résidus, qui juge la normalité comme au chapitre 12.
Figure 13.2. Trois graphiques des résidus en fonction de x. À gauche, l'essai d'usure: dix résidus sans structure visible, entre −6 et +6 µm, répartis des deux côtés de zéro sur toute la plage. Au centre, une relation courbe ajustée par une droite: les résidus dessinent un U, positifs aux bords et négatifs au centre. À droite, une erreur dont l'écart-type croît avec x: les résidus s'ouvrent en entonnoir. Les panneaux du centre et de droite sont des données simulées à partir d'un générateur ensemencé, construites pour montrer chacune un seul défaut; chaque panneau a sa propre échelle verticale.
Sur l'essai d'usure (figure 13.2, à gauche), les résidus réduits ei/s vont de −1,39 à +1,39: aucun point suspect. Les signes alternent sans motif — +,−,+,−,+,+,−,+,−,− —, sans courbure ni entonnoir apparents. Avec dix points, cette inspection n'a pas beaucoup de puissance, et il faut le dire: elle exclut un défaut grossier, elle ne certifie pas le modèle.
Analyse de la variance à un facteur
Le problème: plus de deux moyennes
Le chapitre 12 sait comparer deux moyennes. L'atelier a trois machines, une centrale de béton en compare trois, un essai agronomique quatre engrais. Faire tous les tests de Student deux à deux est une mauvaise idée: avec k groupes, il y a k(k−1)/2 comparaisons, et le chapitre 12 a chiffré ce que coûtent des tests multiples non corrigés. L'analyse de la variance (ANOVA) pose d'abord une seule question globale: les k moyennes sont-elles toutes égales?
Les hypothèses sont celles du test de Student à variance combinée, étendues à k groupes: indépendance, normalité, et même variance dans tous les groupes. L'alternative ne dit pas quelles moyennes diffèrent; c'est une question qui vient après.
La décomposition, encore
On note Yˉ⋅j la moyenne du groupe j et Yˉ la moyenne générale des N observations. Exactement comme en régression, l'écart d'une observation à la moyenne générale se coupe en deux: Yij−Yˉ=(Yˉ⋅j−Yˉ)+(Yij−Yˉ⋅j) — l'écart de son groupe à l'ensemble, et son écart à l'intérieur de son groupe.
Démonstration. En sommant le carré de la décomposition, le double produit vaut 2∑j(Yˉ⋅j−Yˉ)∑i(Yij−Yˉ⋅j)=0, car les écarts à la moyenne de chaque groupe somment à zéro. Pour SCR: dans le groupe j, ∑i(Yij−Yˉ⋅j)2=(nj−1)Sj2, d'espérance (nj−1)σ2 par le théorème 9.4; on somme, ∑j(nj−1)=N−k. Pour SCE, la forme de König–Huygens pondérée donne SCE=∑jnjYˉ⋅j2−NYˉ2. Or Yˉ⋅j a pour espérance μj et pour variance σ2/nj, et Yˉ a pour espérance μˉ et pour variance σ2/N; donc E[njYˉ⋅j2]=σ2+njμj2 et E[NYˉ2]=σ2+Nμˉ2. Il vient E[SCE]=kσ2+∑jnjμj2−σ2−Nμˉ2, et ∑jnjμj2−Nμˉ2=∑jnj(μj−μˉ)2 par König–Huygens une dernière fois. □
(13.15) contient toute l'idée de la méthode. Le carré moyen résiduelCMR=SCR/(N−k) estime σ2toujours, que H0 soit vraie ou non: c'est la variance combinée du chapitre 12, étendue à k groupes. Le carré moyen inter-groupesCME=SCE/(k−1) estime σ2si H0 est vraie, et une quantité plus grande sinon. Leur rapport est voisin de 1 sous H0 et grand sous H1: c'est d'une comparaison de deux variances que l'on tire une conclusion sur des moyennes — d'où le nom de la méthode.
Démonstration, avec une étape admise dans le cas général. Le dénominateur: SCR/σ2=∑j(nj−1)Sj2/σ2 est une somme de k khi-deux indépendants à nj−1 degrés de liberté (théorème 9.5, groupes indépendants), donc un χN−k2 par additivité. L'indépendance: SCE ne dépend que des moyennes de groupe Yˉ⋅j, et SCR que des variances Sj2; or dans chaque groupe moyenne et variance empiriques sont indépendantes (théorème 9.5), et les groupes sont indépendants entre eux, donc SCE et SCR sont indépendantes. Le numérateur: lorsque les groupes ont la même taillen, sous H0 les k moyennes Yˉ⋅j sont i.i.d. de loi N(μ,σ2/n), et SCE=n∑j(Yˉ⋅j−Yˉ)2 est (k−1) fois la variance empirique de ces k moyennes multipliée par n; le théorème 9.5 appliqué à elles donne SCE/σ2∼χk−12. Pour des tailles inégales, le même résultat demande une projection orthogonale pondérée et il est admis. Le quotient des deux khi-deux réduits suit alors Fk−1,N−k par définition. On rejette pour les grandes valeurs seulement, puisque (13.15) montre que H1 augmente l'espérance du numérateur sans toucher au dénominateur. □
La table d'analyse de la variance a la même forme qu'en régression, avec k−1 degrés de liberté pour la ligne «inter-groupes» et N−k pour la ligne «intra-groupes».
Figure 13.3. Résistance à la compression (MPa) de cinq éprouvettes de béton pour chacune de trois centrales (données fictives). Les points sont les mesures (deux mesures égales sont décalées latéralement), les traits pleins courts les moyennes de groupe, 32, 37 et 32 MPa, et la ligne tiretée la moyenne générale, 33,67 MPa. Sur la centrale 2, l'accolade de droite montre un écart intra-groupe (une mesure à sa moyenne de groupe) et l'accolade de gauche l'écart inter-groupes (la moyenne de groupe à la moyenne générale). L'analyse de la variance compare la taille typique de ces deux sortes d'écarts.
Deux groupes: l'ANOVA est le test de Student
Démonstration. Pour k=2, Yˉ=(n1Yˉ⋅1+n2Yˉ⋅2)/N, donc Yˉ⋅1−Yˉ=Nn2(Yˉ⋅1−Yˉ⋅2) et Yˉ⋅2−Yˉ=−Nn1(Yˉ⋅1−Yˉ⋅2). D'où
et n1n2/N=1/(1/n1+1/n2). Le dénominateur SCR/(N−2)=((n1−1)S12+(n2−1)S22)/(n1+n2−2) est la variance combinée Sp2. Le quotient est T2, et l'équivalence des deux tests suit comme au théorème 13.7. □
Question 13.5
Remettez dans l'ordre les étapes d'une analyse de la variance à un facteur.
Glissez les éléments pour les mettre dans le bon ordre
1.
Diviser chaque somme de carrés par ses degrés de liberté, k−1 et N−k
2.
Calculer SCE et SCR, et contrôler que leur somme vaut SCT
3.
En cas de rejet, comparer les paires de groupes avec une correction pour tests multiples
4.
Calculer les moyennes de groupe et la moyenne générale
5.
Comparer F=CME/CMR au quantile Fk−1,N−k;1−α
6.
Vérifier sur les données l'indépendance, la normalité approximative et des dispersions comparables dans les groupes
Ce qui vient après
Ce chapitre clôt le cours. La régression et l'analyse de la variance sont deux visages d'un même objet, le modèle linéaire: dans les deux cas, l'espérance de la réponse est une combinaison linéaire de paramètres inconnus, on l'estime par les moindres carrés, et l'on teste par des rapports de sommes de carrés. La régression multiple, où la réponse dépend de plusieurs variables explicatives, s'écrit en notation matricielle Y=Xβ+ε et réutilise mot pour mot ce chapitre: β^=(X⊤X)−1X⊤Y, n−p degrés de liberté pour p paramètres, un test F pour un groupe de coefficients. L'analyse de la variance à plusieurs facteurs ajoute la notion d'interaction et celle de plan d'expérience. Et l'économétrie affronte la question que ce chapitre a seulement signalée: comment estimer un effet causal quand x n'a pas été fixé par l'expérimentateur.
Synthèse
Le modèle de régression linéaire simpleYi=β0+β1xi+εi suppose des xi connus, des erreurs indépendantes, de même variance σ2, et — pour les lois exactes — normales. Les estimateurs des moindres carrés B1=SxY/Sxx et B0=Yˉ−B1xˉ sont des combinaisons linéaires des Yi, sans biais, de variances σ2/Sxx et σ2(1/n+xˉ2/Sxx), et les meilleurs des estimateurs linéaires sans biais (Gauss–Markov).
La variance des erreurs s'estime par S2=SCR/(n−2): deux coefficients estimés, deux degrés de liberté perdus. Les pivots (Bj−βj)/σBj suivent , d'où les intervalles et les tests. Sur l'essai d'usure: µm/h, intervalle à , — la dérive est réelle.
En un point x0, l'intervalle de confiance encadre la droite moyenne et se referme quand n croît; l'intervalle de prédiction encadre une nouvelle observation et garde la largeur de σ. Les deux s'élargissent loin de xˉ, et aucun ne protège d'une extrapolation hors du domaine observé.
SCT=SCE+SCR (Pythagore dans Rn), R2=SCE/SCT=r2, et le test F de la régression, , est le carré du test de Student sur la pente. ne valide ni la linéarité ni une prévision: ce sont les qui jugent le modèle — courbure, entonnoir, point isolé.
L'analyse de la variance à un facteur teste l'égalité de k moyennes par F=CME/CMR∼Fk−1,N−k: le carré moyen intra-groupes estime toujours σ2, le carré moyen inter-groupes ne l'estime que sous H0. Pour , et l'on retrouve le test de Student à variance combinée.
Un rejet global ne dit pas quelles moyennes diffèrent, et un non-rejet global peut coexister avec une paire «significative» testée seule: les comparaisons deux à deux se font après, avec une correction pour tests multiples. Les trois machines de l'atelier donnent F=2,97, p=0,056, alors que A et B seules donnaient p=0,022.
Série d'exercices du chapitre 13Exercice 1 sur 5
Question 13.6
Une régression linéaire simple donne R2=0,90. Quelle lecture est correcte?
Problème guidé 13.1 · Résistance du béton et rapport eau/ciment
Un laboratoire (données fictives) mesure la résistance à la compression yi (MPa) de huit éprouvettes gâchées avec des rapports eau/ciment xi de 0,40; 0,45; 0,50; 0,55; 0,60; 0,65; 0,70; 0,75. Les résistances valent 49; 46; 43; 42; 36; 35; 31; 28. On donne xˉ=0,575, yˉ=38,75, Sxx=0,105, Sxy=−6,3 et Syy=383,5, ainsi que t6;0,975=2,4469. On adopte le modèle de régression linéaire simple à erreurs normales.
1
1. La droite des moindres carrés
Question
Calculez la pente b1, en MPa par unité de rapport eau/ciment.
2. La dispersion autour de la droite
3. Le test de la pente
4. Prévoir une éprouvette
Exercices
Vous pouvez afficher le corrigé directement sous chaque énoncé après avoir cherché la solution.
Exercice 13.1 · Une régression à partir des sommes
Pour étudier le temps de séchage y (en heures) d'une peinture en fonction de l'humidité relative x (en %), on a fait n=10 essais (données fictives) aux humidités 40,45,…,85%. On donne
Soit la droite des moindres carrés y^i=b0+b1xi d'une série double, et ei=yi−y^i ses résidus.
Montrez que ∑iei=0 et ∑ixiei=0, et déduisez-en que .
Solution
1. Par (10.19), ei=(yi−yˉ)−b1(xi−xˉ). En sommant, . Ensuite
Exercice 13.3 · Prévoir quand réaffûter, et le danger d'extrapoler
On reprend l'essai d'usure: y^=−14+4x µm, s=4,3012 µm, n=10, xˉ=5,5 h, Sxx=82,5, t8;0,975=2,3060 et t8;0,95=1,8595.
Donnez l'intervalle de confiance à 95% de l'écart moyen et l'intervalle de prédiction à 95% de l'écart d'une pièce, à la douzième heure.
Le chef d'atelier veut réaffûter avant que le diamètre moyen n'atteigne +40 µm. Quelle durée la droite estimée suggère-t-elle?
Plus prudent, il veut qu'une pièce isolée ne dépasse +40 µm qu'avec une probabilité d'au plus 5%. On utilise la borne supérieure de prédiction unilatéraley^0+t8;0,95s1+1/n+(x0−xˉ)2/Sxx. Calculez-la pour et h, et concluez.
Solution
1. À x0=12: y^0=34 µm et 1/10+6,52/82,5=0,61212. L'intervalle de confiance vaut
Exercice 13.4 · Trois huiles de coupe
Un atelier compare la rugosité Ra (en µm) obtenue avec trois huiles de coupe (données fictives):
Huile 1
1,2
1,5
1,3
1,4
Huile 2
1,6
1,8
1,5
1,9
1,7
Huile 3
1,4
1,3
1,6
1,5
1,2
1,4
Calculez les moyennes de groupe, la moyenne générale, SCE, SCR et SCT, et vérifiez la décomposition.
Dressez la table d'analyse de la variance et concluez aux seuils de 5% et de 1% (F2,12;0,95=3,885, F2,12;0,99=6,927).
On ne compare plus que les huiles 1 et 3. Calculez la statistique de Student à variance combinée et la statistique F de l'analyse de la variance réduite à ces deux groupes, et vérifiez la relation entre elles.
Solution
1. Moyennes: yˉ1=5,4/4=1,35, yˉ2=8,5/5=1,70, , et la moyenne générale vaut µm.
Exercice 13.5 · Où placer les points d'un essai?
On veut estimer une pente dans le modèle (13.1) avec n=10 observations, les durées xi pouvant être choisies librement entre 1 et 10 heures.
Comparez Var(B1) pour le plan de l'essai d'usure (xi=1,2,…,10) et pour le plan «aux extrémités» (cinq pièces à 1 h, cinq à 10 h).
Quel est le défaut grave du second plan, que la variance ne montre pas?
Sur le plan de l'essai d'usure, on considère l'estimateur «des extrémités» B~=(Y10−Y1)/9. Montrez qu'il est sans biais, calculez sa variance, comparez-la à celle de B1, et reliez le résultat au théorème de Gauss–Markov. Que vaut sur les données?
En quel x0 la variance de Y^0 est-elle minimale? Interprétez.
Solution
1. Pour xi=1,…,10: Sxx=82,5 et Var(B1)=σ2/82,5=0,01212σ2. Pour cinq points à et cinq à : , chaque écart vaut , et , soit fois moins. Parmi tous les plans à dix points dans , celui-ci maximise , donc minimise la variance de la pente.
Références
Morgenthaler, S., Introduction à la statistique, Presses polytechniques et universitaires romandes, Lausanne — chapitres sur la régression linéaire et l'analyse de la variance, dans l'esprit d'un cours d'ingénieur de l'EPFL.
Saporta, G., Probabilités, analyse des données et statistique, Technip, Paris — chapitres sur le modèle linéaire, la régression simple et multiple, et l'analyse de la variance.
Wackerly, D., Mendenhall, W. et Scheaffer, R., Mathematical Statistics with Applications, Cengage — chapitre 11 (modèles linéaires et moindres carrés) et chapitre 13 (analyse de la variance), avec les démonstrations des lois d'échantillonnage.
Rice, J. A., Mathematical Statistics and Data Analysis, Duxbury — chapitre sur l'analyse de la variance et chapitre sur les moindres carrés linéaires, y compris le point de vue matriciel.
Draper, N. R. et Smith, H., Applied Regression Analysis, Wiley — la référence classique pour l'analyse des résidus et le diagnostic d'un modèle de régression.
Galton, F., «Regression Towards Mediocrity in Hereditary Stature», Journal of the Anthropological Institute of Great Britain and Ireland, vol. 15, 1886, p. 246–263 — l'article qui a donné son nom à la régression.
Anscombe, F. J., «Graphs in Statistical Analysis», The American Statistician, vol. 27, nº 1, 1973, p. 17–21 — les quatre nuages de même droite et de même R2, à relire avant toute régression.
indépendantes
N(0,σ2)
x
explicative
Y
expliquée
x↦β0+β1x
droite de régression
x
B0∼N(β0,σ2(1/n+xˉ2/Sxx))
(n−2)S2/σ2∼χn−22
S2
(B0,B1)
18
22
26
−1
b1=4
∑iei=0
∑ixiei=3−8+6−24+25+6−28+48−18−10=0
5%
trop petits
7
21
1
10
H0:β1=0
r
n=10
∣r∣≥0,632
5%
n=100
∣r∣≥0,197
0,2
4%
1−α
μ0=β0+β1x0
ε0∼N(0,σ2)
n
Var(Y0−Y^0)=σ2(1+n1+Sxx(x0−xˉ)2)
conditionnels au modèle
H0
F≥F1,n−2;1−α
exactement
α
SCR=36
1%
p
0,00075
p
0,00065
0,0167
5
tn−2
b1=4
95%
[2,91;5,09]
tobs=8,45
SCE/(SCR/(n−2))
R2
résidus
k=2
F=T2
R2
s
Testez H0:β1=0 au seuil de 5% et donnez l'intervalle de confiance à 95% de la pente (t8;0,975=2,3060). Traduisez le résultat pour une hausse de 10 points d'humidité.
Donc b1=89,5/2062,5=0,043394 h par point d'humidité et b0=3,08−0,043394×62,5=0,3679 h.
2.SCE=Sxy2/Sxx=8010,25/2062,5=3,88376, SCR=3,956−3,88376=0,07224, donc R2=3,88376/3,956=0,9817, s2=0,07224/8=0,009030 h² et s=0,0950 h, soit environ 6 minutes.
3. L'erreur type de la pente vaut s/Sxx=0,09503/45,415=0,0020924, donc
tobs=0,00209240,043394=20,74>2,3060:
on rejette H0 très nettement. L'intervalle à 95% vaut 0,043394±2,3060×0,0020924=[0,03857;0,04822] h par point. Une hausse de 10 points d'humidité allonge le séchage de 0,39 à 0,48 h, soit de 23 à 29 minutes environ.
∑iy^iei=0
Montrez que la moyenne des valeurs ajustées y^i est yˉ.
Démontrez la décomposition SCT=SCE+SCR et l'égalité R2=r2.
Montrez que R2=0 si et seulement si b1=0, et interprétez la droite des moindres carrés dans ce cas.
∑iei=0−b1×0=0
i∑(xi−xˉ)ei=Sxy−b1Sxx=Sxy−Sxy=0,
et ∑ixiei=∑i(xi−xˉ)ei+xˉ∑iei=0. Enfin ∑iy^iei=b0∑iei+b1∑ixiei=0. Ce sont les deux équations normales: le vecteur des résidus est orthogonal au vecteur (1,…,1) et au vecteur des xi, donc à toute combinaison des deux.
2.∑iy^i=∑i(yi−ei)=∑iyi, donc la moyenne des y^i est yˉ.
3.yi−yˉ=(y^i−yˉ)+ei. En élevant au carré et en sommant, le double produit 2∑i(y^i−yˉ)ei=2∑iy^iei−2yˉ∑iei=0 par la question 1, d'où SCT=SCE+SCR. Comme y^i−yˉ=b1(xi−xˉ), SCE=b12Sxx=Sxy2/Sxx, et
R2=SCTSCE=SxxSyySxy2=r2.
4.R2=0 équivaut à SCE=b12Sxx=0, donc à b1=0 puisque Sxx>0. La droite des moindres carrés est alors horizontale, y^=yˉ: on retrouve le modèle constant, dont l'estimateur des moindres carrés est la moyenne (chapitre 10). Connaître x n'améliore en rien la prévision linéaire de y — ce qui, rappelons-le, n'exclut pas un lien non linéaire.
x0=11
x0=11,1
Quelle réserve faut-il formuler sur toutes ces réponses?
3. À x0=11: y^0=30, 1+0,1+5,52/82,5=1,46667, et la borne vaut 30+1,8595×4,3012×1,21106=30+9,69=39,69 µm, sous 40. À x0=11,1: y^0=30,4, 1,48012, et la borne vaut 30,4+9,73=40,13 µm, au-dessus. Le seuil est atteint vers 11,1 h: pour protéger chaque pièce et non la seule moyenne, il faut réaffûter deux heures et demie plus tôt que ne le suggérait la droite. C'est l'écart entre σ/n et σ du chapitre 9, transposé à la régression.
4. Toutes ces réponses sont des extrapolations: l'essai s'arrête à 10 h, et le modèle linéaire n'a été vérifié qu'entre 1 et 10 h. Une usure d'outil s'accélère souvent en fin de vie, ce qui rendrait les durées calculées trop optimistes. La seule réponse honnête est de prolonger l'essai au-delà de 10 heures avant de fixer la règle de réaffûtage — ou d'adopter, en attendant, la durée la plus prudente, 11 h, en la présentant comme provisoire.
Les sommes de carrés intra-groupes valent 0,05 (huile 1: écarts −0,15; 0,15; −0,05; 0,05), 0,10 (huile 2) et 0,10 (huile 3), d'où SCR=0,25. Directement, SCT=0,59733=0,34733+0,25.
2.
Source
Somme de carrés
ddl
Carré moyen
F
Inter-huiles
0,34733
2
0,17367
8,336
Intra-huiles
0,25
12
0,020833
Totale
0,59733
14
8,336>6,927: on rejette l'égalité des trois rugosités moyennes au seuil de 1% (valeur p: 0,0054). Les moyennes suggèrent que c'est l'huile 2 qui se distingue, avec une rugosité plus élevée d'environ 0,3 µm; on le vérifierait par des comparaisons deux à deux avec correction de Bonferroni.
3. Pour les huiles 1 et 3 seules: sp2=(0,05+0,10)/(4+6−2)=0,01875 et
L'analyse de la variance à deux groupes donne SCE=0,006, SCR=0,15 et F=0,006/(0,15/8)=0,32=(−0,5657)2: c'est le théorème 13.10. Aucune différence entre les huiles 1 et 3 (p=0,59).
b~
1
10
xˉ=5,5
±4,5
Sxx=10×4,52=202,5
Var(B1)=σ2/202,5=0,00494σ2
2,45
[1,10]
Sxx
2. Avec deux valeurs de x seulement, toute relation passe par deux points: on ne peut plus voir une courbure, ni tester la linéarité. Le plan optimal pour la variance est aveugle à l'hypothèse dont dépend tout le reste. En pratique, on garde des points intermédiaires — au moins un au centre — pour pouvoir lire les résidus.
3.E[B~]=((β0+10β1)−(β0+β1))/9=β1: il est sans biais. Il est linéaire en Y (a10=1/9, a1=−1/9, les autres nuls), et
Var(B~)=812σ2=0,02469σ2>82,5σ2=0,01212σ2,
soit 2,04 fois la variance de B1. C'est ce qu'annonce Gauss–Markov: parmi les estimateurs linéaires sans biais, celui des moindres carrés a la plus petite variance; B~ jette l'information des huit points intérieurs. Sur les données, b~=(25−(−7))/9=3,556 µm/h, contre b1=4.
4.Var(Y^0)=σ2(1/n+(x0−xˉ)2/Sxx) est minimale en x0=xˉ, où elle vaut σ2/n — la variance d'une moyenne de n observations. La droite est la mieux connue au point moyen, par lequel elle passe toujours; plus on s'en éloigne, plus l'incertitude sur la pente s'ajoute. C'est la forme en hyperbole des bandes de la figure 13.1.