Analyse de la variance¶
Reprenons les données sur le poids de naissance. En effet, lors de l'étude, d'autres variables que celle que nous avons regardées au dernier TP ont été observées. Le fichier poids_naiss_2 contient les variables suivantes :
- le poids de naissance du bébé (en grammes) (BWT = birth weight),
- tabagimse durant la grossesse (Y=oui; N=non) (SMOKE)
- race de la mère (1=blanche; 2=noire; 3=autre) (RACE)
- antécédents d'hypertensions (Y=oui; N=non) (HT)
- présence d'irritabilité utérine (Y=oui; N=non) (UI)
naiss <- read.table('poids_naiss_2',header=T,sep=':')
head(naiss)
La fonction table calcule les effectifs d'un variable :
table(naiss$RACE)
table(naiss$SMOKE)
table(naiss$HT)
table(naiss$UI)
On voit bien que pour chaque variable la taille des groupes est variables, autrement dit le plan d'expérience est non équilibré en chaque variable.
Modèle à un facteur¶
Pour commencer nous considérons le modèle à un facteur où le poids de naissance est la variable d'intérêt et la race de la mère est le facteur.
Représentation des données¶
Afin de visualiser l'impact du facteur sur la variable à expliquer on peut tracer des boxplots :
boxplot(BWT~RACE,data=naiss, xlab='RACE')
On observe des légères différences entre les trois boxplots. Mais sont-elles significatives ?
Estimation des paramètres¶
Comparons d'abord les différentes façons de définir le modèle linéaire avec R. D'abord, en utilisant la fonction lm de la même façon que dans les TP précédents on obtient :
mod1 <- lm(BWT~RACE,data=naiss)
summary(mod1)
Dans mod1 la variable RACE est interpréter comme une variable quantitative. Plus précisément, on considère le modèle $$\text{BWT}_i = \beta_0 +\beta_1\text{RACE}_i+\varepsilon_i,$$ et les estimateurs dans le tableau coefficients sont les estiamteurs de $\beta_0$ et $\beta_1$. Mais comme la variable RACE est de type qualitative, ce modèle n'a pas de sens d'un point de vue d'interprétation !
Voyons ce qui se passe en la définissant comme une variable qualitative, c'est à dire en la déclarant comme une variable de type factor.
naiss_2 <- naiss
naiss_2$RACE <- as.factor(naiss_2$RACE)
mod2 <- lm(BWT~RACE,data=naiss_2)
summary(mod2)
Tout d'abord, on observe que le tableau Coefficients est composé de trois lignes (Intercept, RACE2 et RACE3) contrairement au modèle mod1, où il n'y a que deux (Intercept et RACE). Quel est alors le modèle utilisé ici ? En effet, il s'agit d'un modèle de type analyse de la variance. Notons les observations $\text{BWT}_{i,j}$ où l'indice $i$ indique le niveau du facteur (c'est-à-dire de RACE$=i$). De plus, on note $\mu_i$ les moyennes par groupe : $$\mu_i= \mathbb E[\text{BWT}_{i,j}],\qquad i=1,\dots,I\quad j=1,\dots,n_i.$$ Or, le modèle mod2 correspond au modèle donné par $$\text{BWT}_{ij} = \mu_1+ \delta_i +\varepsilon_{ij},$$ avec la contrainte $\delta_1=0$. Autrement dit, les $\delta_i$ pour $i=2,\dots,I$ vérifient $$\delta_i= \mu_i-\mu_1,$$ c'est-à-dire l'écart entre les moyennes $\mu_i$ et $\mu_1$. La colonne Estimate contient les estimateurs de $\mu_1$ (Intercept), de $\delta_2$ (RACE2) et de $\delta_3$ (RACE3).
Dans le modèle mod2, le niveau $i=1$ joue un rôle particulier : il est le niveau de référence, car les autres paramètres dépendent de lui. Dans une application où il y a p.ex. un groupe de contrôle auquel différents traitements sont comparés, il est naturel que ce groupe de contrôle soit le groupe de référence. Dans la définition du modèle par lm on peut préciser le niveau qui doit servir de référence par l'option base. Par exemple : pour choisir le deuxième groupe comme niveau de référence, on écrit :
mod2bis <- lm(BWT~C(RACE,base=2),data=naiss_fac)
ce qui donne lieu au modèle
$$\text{BWT}_{ij} = \mu_2+ \delta_i +\varepsilon_{ij},$$
avec la contrainte $\delta_2=0$.
Cette paramétrisation du modèle ne correspond à aucune de celles que nous avons vues en cours. Pour travailler avec le modèle donné par $$\text{BWT}_{ij} = m+ \alpha_i +\varepsilon_{ij},$$ où $m=\frac1n\sum_{i=1}^In_i\mu_i$ et sous la contrainte $\sum_{i=1}^In_i\alpha_i=0$, on écrit
n_i <- table(naiss_2$RACE)
I <- length(n_i)
contrainte <- rbind(diag(I-1),-n_i[-I]/n_i[I])
mod3 <- lm(BWT~C(RACE,contrainte),data=naiss_2)
summary(mod3)
Les valeurs de la colonne Estimate sont : l'estimateur de la moyenne globale $m$ par la moyenne empirique de toutes les observations $\bar X_{\cdot\cdot}$ (Intercept) et les estimateurs $\hat\alpha_1$ (C(RACE, contrainte)1) et $\hat\alpha_2$ (C(RACE, contrainte)2). La valeur de $\hat\alpha_3$ découle de la contrainte $\sum_{i=1}^In_i\alpha_i=0$. Ainsi, ici $\hat\alpha_3=(-n_i\hat\alpha_1-n_2\hat\alpha_2)/n_3$.
La matrice contrainte utilisée dans la défintion de mod3 précise la contrainte $\sum_{i=1}^In_i\alpha_i=0$. Au cas où le plan d'expérience est équilibré (i.e. $n_1=\dots=n_I$), on obtient le même modèle par l'instruction :
mod3bis <- lm(BWT~C(RACE,sum),data=naiss_fac)
Plus précisément, c'est la contrainte $\sum_{i=1}^I\alpha_i=0$.
D'ailleurs, quand le plan d'expérience n'est pas équilibré, on peut toujours utiliser la définition mod3bis. Mais il sera plus difficile d'interpréter les paramètres $\alpha_i$.
Test de significativité des coefficients¶
Quelque soit le modèle (mod1, mod2 ou mod3), les deux dernières colonnes t value et Pr(>|t|) du tableau Coefficients correspondent aux différents tests de nullité du coefficient. Pour une interprétation correcte de ces tests, il faut bien avoir en tête la paramétrisation utilisée. Ainsi, dans le modèle mod1 la première ligne du tableau (Intercept) correspond au test $H_0:\beta_0=0$, dans le modèle mod2 au test $H_0:\mu_1=0$ et dans le modèle mod3 au test de nullité de la moyenne globale $H_0:m=0$. Cela explique les différentes valeurs de la statistique de test t value dans la ligne Intercept.
De même, les autres lignes du tableau Coefficients correspondent
- dans le modèle mod1 : au test $H_0:\beta_1=0$ ;
- dans le modèle mod2 : au test $H_0:\delta_i=0$ pour $i=2,\dots,I$, c'est-à-dire au test d'égalité des moyennes $H_0:\mu_i=\mu_1$ pour $i=2,\dots,I$ ;
- dans le modèle mod3 : au test $H_0:\alpha_i=0$ pour $i=2,\dots,I$, c'est-à-dire au test $H_0:\mu_i=m$ pour $i=2,\dots,I$.
Tableau de l'analyse de la variance¶
Pour l'analyse de la variance, il faut impérativement utiliser un des modèles où la variable explicative est déclarée comme factor, à savoir mod2 ou mod3.
La fonction anova, appliquée à un objet renvoyé par lm, effectue le calcul de la table de la variance.
Remarquons qu'il existe aussi la fonction aov pour faire une analyse de la variance, mais elle est moins complète.
Pour l'analyse de la variance, la paramétrisation du modèle n'a pas d'importance, car dans les deux cas, les sous-espaces vectoriels sont les mêmes.
On obtient le tableau de l'analyse de la variance par l'instruction :
anova(mod2)
Ou par la commande :
anova(mod3)
Les colonnes du tableau sont
- Df = Degree of Freedom = nombre de degrés de liberté
- Sum Sq = Sum of Squares = somme des carrés des écarts
- Mean Sq = Mean sum of Squares = carré moyen des écarts
- F value = valeur de la statistique de test $\mathcal F$
- Pr(>F) = p-valeur associée au test
La première ligne (RACE) correspond au facteur : c'est ici où on trouve le résultat du test de l'absence d'effet dû au facteur, à savoir $$H_0:\mu_1=\dots=\mu_I.$$
La deuxième ligne (Residuals) donne les informations sur les résidus.
Dans cet exemple, on conclut qu'il y a effectivement un effet dû à la race de la mère sur le poids de naissance de son bébé.
Exercice 1¶
Considérer les autres variables dans le tableau des données (SMOKE, UI et HT).
- Faites des représentations graphiques appropriées pour ces variables et interpréter-les.
- En utilisant la fonction anova, faites les tests d'absence d'effet dû aux différentes variables (SMOKE, UI et HT). Quelles sont les conclusions de ces tests ?
boxplot(BWT~SMOKE*RACE, data=naiss_2, xlab='SMOKE/RACE')
Pour analyser l'impacet des facteurs, une autre représentations graphiques est préférable :
with(naiss_2,interaction.plot(RACE,SMOKE,BWT))
On voit bien que la courbe des moyennes de poids associées aux non-fumeurs est plus élevée que celle des fumeurs quelque soit la race de la mère. Il semble d'y avoir un vrai effet de tabagisme sur le poids de naissance.
On peut tracer le même graphique en permutant les deux facteurs :
with(naiss_2,interaction.plot(SMOKE,RACE,BWT))
L'effet dû à la race est un peu moins visible que l'effet du au tabagisme.
Estimation des paramètres¶
Similaire au modèle à un facteur, plusieurs paramétrisations du modèle à deux facteurs avec interaction entre les facteurs sont possibles.
Le modèle par défaut dans R est donné par $$\text{BWT}_{ijk} = \mu_{11}+ \delta_i + \nu_j + \tau_{ij} +\varepsilon_{ij},$$ avec les contraintes $\delta_1=0$, $\nu_1=0$, $\tau_{1j}=0, j=1,\dots,J$ et $\tau_{i1}=0,i=1,\dots,I$. On l'obtient par l'instruction suivante :
mod4 <- lm(BWT~SMOKE*RACE, data= naiss_2)
summary(mod4)
La colonne Estimate contient les estimateurs de :
- $\mu_{11}$ (Intercept)
- $\delta_2$ (SMOKEY) (correspond au niveau Y du facteur SMOKE)
- $\nu_2, \nu_3$ (RACE2, RACE3)
- $\tau_{22}, \tau_{23}$ (SMOKEY:RACE2,SMOKEY:RACE3)
Le modèle vu en cours pour un plan est équilibré, à savoir $$\text{BWT}_{ijk} = m+ \alpha_i + \beta_j + \gamma_{ij} +\varepsilon_{ij},\qquad i=1,\dots,I\quad j=1,\dots,J,$$ avec $m=\frac1{IJ}\sum_i\sum_j\mu_{ij}$ et sous les contraintes $$\sum_i\alpha_i=0, \sum_j\beta_j=0, \sum_i\gamma_{ij}=0, j=1,\dots,J\quad\text{ et }\quad \sum_j\gamma_{ij}=0,i=1,\dots,I$$ est obtenu par l'instruction
mod5 <- lm(BWT~C(SMOKE,sum)*C(RACE,sum), data= naiss_2)
summary(mod5)
Attention : dans notre exemple le plan d'expérience n'est pas équilibré, donc il faut éviter d'interpréter les coefficients estimés. En revanche, on peut utiliser ce modèle pour faire l'analyse de la variance.
Analyse de la variance¶
On obtient le tableau de l'analyse de la variance par la fonction anova. Là encore, le choix des contraintes sur les paramètres n'a pas d'importance, car les sous-espaces vectoriels intervenant dans l'analyse de la variance sont les mêmes quelques soient les contraintes d'identifiabilité. Ainsi, on obtient
anova(mod4)
anova(mod5)
L'interprétation est la même que dans le cas à un facteur. La dernière colonne donne les p-valeurs des différents tests, à savoir :
- le test de l'absence d'effet principal dû au facteur SMOKE
- le test de l'absence d'effet principal dû au facteur RACE
- le test de l'absence d'interaction entre les facteurs.
Nous voyons que les tests sur les facteurs sont tous les deux significatifs. En revanche, la p-valeur du test d'interaction est relativement élevé avec une valeur de $0.1$. Donc, dans ce cas on ne peut pas rejeter l'hypothèse nulle au niveau habituel de $0.05$ et un modèle sans interaction est envisageable vu les données.
On définit le modèle sans interaction par la commande :
mod6 <- lm(BWT~SMOKE+RACE, data= naiss_2)
anova(mod6)
Exercice 2¶
Pour les données du poids de naissance, considérons les trois modèles suivants à deux facteurs:
- les facteurs sont SMOKE et UI
- les facteurs sont RACE et HT
- les facteurs sont UI et HT
- Faites des représentations graphiques appropriées et interpréter-les.
- Effectuer les différents tests de l'analyse de la variance (avec et sans interaction). Quelles sont les conclusions de ces tests ? Dans le modèle avec les facteurs hypertension et irritabilité utérine le résultat ne correspond pas à ce qu'on attend. Expliquer et en déduire une condition nécessaire sur les données pour l'analyse de la variance.
- En bas de la sortie de summary(mod) on trouve toujours la valeur du coefficient de détermination $R^2$ et du coefficient ajusté $R^2_a$ associés au modèle considéré. Utiliser ces valeurs afin de comparer tous les modèles à deux variables entre eux et en choisir le meilleur.
Devoir à rendre¶
On reprend les données d'ozone.
Analyse de la variance¶
Le jeu de données ozone contient deux variables qualitatives. Analyser l'effet de ces deux variables sur la concentration d'ozone. On pourra faire des représenations graphiques et effectuer l'analyse de la variance à un et à deux facteurs. Il y t-il des niveaux de variables qui ne sont pas significativement différents ?
Sélection de variables¶
Le jeu de données ozone contient de nombreuses variables explicatives. Le but est d'appliquer des critères classiques de sélection de variables pour en choisir les plus pertinentes pour expliquer la concentration d'ozone. Dans cette partie on met de côté les variables qualitatives et on se concentre sur les variables quantitatives seulement.
Pour le choix de variables, on peut appliquer la fonction regsubsets du package leaps (à charger par l'instruction library(leaps)) qui fait une recherche exhaustive. Utiliser l'aide pour comprendre le fonctionnenment de cette fonction. Appliquer-la à vos données afin de déterminer le meilleur modèle selon
- le coefficient de détermination $R^2$,
- le coefficient de détermination ajusté $R^2_a$,
- le $C_p$ de Mallows,
- le critère BIC.
On pourra tracer des graphiques pour représenter les résultats par une plot de l'objet renvoyé par regsubsets.
Commenter vos résultats.
Modalités du devoir¶
Le devoir doit être préparé en binôme.
Il est dû au plus tard le 20 décembre avant minuit. Tout document rendu après cette date ne sera pas noté.
Le devoir est à déposer dans votre boîte de dépôt Sakai.
Vous devez indiquer le nom de chacun des deux auteurs sur les documents à rendre.
Le rapport est à rédiger en R Markdown.
Vous devez rendre un fichier html (ou pdf) et le fichier Rmd associé.