Régression linéaire multiple et ANOVA à 2 facteurs

16 Pages·2016·0.26 MB·French

by  Charlotte Baey

TP 4 - Régression linéaire multiple et ANOVA à 2
facteurs
Charlotte Baey
2 décembre 2016
Plan de la séance
• fin du TP précédent sur la régression linéaire multiple, notamment concernant la sélection de modèles
• ANOVA à 2 facteurs
• résumé des principales notions vu pendant le cours et des fonctions R associées
• séance de questions
Exercice 1 : fin du TP sur la régression linéaire multiple
Finir les questions non traitées.
Exercice 2 : ANOVA à 2 facteurs
Les données étudiées se trouve dans le package MASS et s’intitulent crabs. L’objet de ce TP est d’étudier la
variable RW.
1. Préciserlesvariables‘’expliquées”,lesvariablesexplicativesainsiqueleurtype(qualitative,quantitative)
pour l’ensemble du jeu de données. Dans la suite, notre variable d’intérêt sera la variable RW.
library(MASS) # chargement de la librairie contenant le jeu de données
data() # donne la liste des jeux de données disponibles
help(crabs)
attach(crabs) # pour éviter de devoir rajouter crabs$ à chaque appel des variables de la table
str(crabs)
## 'data.frame': 200 obs. of 8 variables:
## $ sp : Factor w/ 2 levels "B","O": 1 1 1 1 1 1 1 1 1 1 ...
## $ sex : Factor w/ 2 levels "F","M": 2 2 2 2 2 2 2 2 2 2 ...
## $ index: int 1 2 3 4 5 6 7 8 9 10 ...
## $ FL : num 8.1 8.8 9.2 9.6 9.8 10.8 11.1 11.6 11.8 11.8 ...
## $ RW : num 6.7 7.7 7.8 7.9 8 9 9.9 9.1 9.6 10.5 ...
## $ CL : num 16.1 18.1 19 20.1 20.3 23 23.8 24.5 24.2 25.2 ...
## $ CW : num 19 20.8 22.4 23.1 23 26.5 27.1 28.4 27.8 29.3 ...
## $ BD : num 7 7.4 7.7 8.2 8.2 9.8 9.8 10.4 9.7 10.3 ...
2. Faire une étude succinte de statistique descriptive de la variable RW (calculs de moyennes, variances,
boîtes à moustaches), globale et en fonction des facteurs sp, sex et de l’interaction sp/sex. On pourra
utiliser la fonction tapply.
summary(RW) # statistiques descriptives de la variable RW tous facteurs confondus
1
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 6.50 11.00 12.80 12.74 14.30 20.20
tapply(RW,sp,mean) # moyenne selon les valeurs du facteur sp
## B O
## 11.928 13.549
tapply(RW,sex,mean) # moyenne selon les valeurs du facteur sex
## F M
## 13.487 11.990
tapply(RW,sp:sex,mean) # moyenne selon les valeurs de l'interaction sp:sex
## B:F B:M O:F O:M
## 12.138 11.718 14.836 12.262
tapply(RW,sp,var) # variance selon les valeurs du facteur sp
tapply(RW,sex,var) # variance selon les valeurs du facteur sex
tapply(RW,sp:sex,var) # variance selon les valeurs du facteur sp:sex
## B O
## 5.195168 6.788787
## F M
## 7.511445 4.667778
## B:F B:M O:F O:M
## 5.947302 4.459057 5.515004 4.820771
par(mfrow=c(1,4))
boxplot(RW)
boxplot(RW~sp)
boxplot(RW~sex)
boxplot(RW~sp:sex)
0 0 0 0
2 2 2 2
6 6 6 6
1 1 1 1
2 2 2 2
1 1 1 1
8 8 8 8
6 6 6 6
B O F M B.F O.F B.M O.M
3. Etudier graphiquement l’intéraction entre les facteurs sexe et espèce à l’aide de la commande
interaction.plot.
par(mfrow=c(1,2))
interaction.plot(sp,sex,RW)
interaction.plot(sex,sp,RW)
2
sex sp
W 0 F W 0 O
. .
4 4
R 1 M R 1 B

f f
o o
0 0

n . n .
3 3
a a
1 1
e e
m m
0 0
. .
2 2
1 1
B O F M
sp sex
4. Ajuster deux modèles d’analyse de la variance à deux facteurs (le modèle complet, c’est-à-dire avec
interaction, et le modèle additif).
# Modele avec interaction RW_ijk=mu+delta_ij+eps_ijk
modeleComp=lm(RW ~ sp + sex + sp:sex)
summary(modeleComp)
# on constate que l'effet espèce est significatif, ainsi que l'effet de l'interaction,
# mais pas l'effet du sexe
##
## Call:
## lm(formula = RW ~ sp + sex + sp:sex)
##
## Residuals:
## Min 1Q Median 3Q Max
## -5.638 -1.280 0.000 1.688 5.364
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 12.1380 0.3220 37.691 < 2e-16 ***
## spO 2.6980 0.4554 5.924 1.39e-08 ***
## sexM -0.4200 0.4554 -0.922 0.357561
## spO:sexM -2.1540 0.6441 -3.344 0.000988 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 2.277 on 196 degrees of freedom
## Multiple R-squared: 0.2287, Adjusted R-squared: 0.2169
## F-statistic: 19.38 on 3 and 196 DF, p-value: 4.813e-11
# Modele constant
modeleCst=lm(RW~1)
summary(modeleCst)
##
## Call:
## lm(formula = RW ~ 1)
3
##
## Residuals:
## Min 1Q Median 3Q Max
## -6.2385 -1.7385 0.0615 1.5615 7.4615
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 12.739 0.182 70.01 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 2.573 on 199 degrees of freedom
# Modele additif
modeleAdd=lm(RW~sp+sex)
summary(modeleAdd)
##
## Call:
## lm(formula = RW ~ sp + sex)
##
## Residuals:
## Min 1Q Median 3Q Max
## -6.1765 -1.5195 0.0025 1.6213 5.9025
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 12.6765 0.2860 44.321 < 2e-16 ***
## spO 1.6210 0.3303 4.908 1.92e-06 ***
## sexM -1.4970 0.3303 -4.533 1.01e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 2.335 on 197 degrees of freedom
## Multiple R-squared: 0.1847, Adjusted R-squared: 0.1765
## F-statistic: 22.32 on 2 and 197 DF, p-value: 1.834e-09
5. Tester l’effet de l’intéraction et de chacun des deux facteurs. Quel modèle préférez-vous ?
anova(modeleComp) # Test l'hypothèse nulle d'absence d'effet de chaque facteur du modèle,
## Analysis of Variance Table
##
## Response: RW
## Df Sum Sq Mean Sq F value Pr(>F)
## sp 1 131.38 131.382 25.336 1.087e-06 ***
## sex 1 112.05 112.050 21.608 6.133e-06 ***
## sp:sex 1 58.00 57.996 11.184 0.0009884 ***
## Residuals 196 1016.36 5.186
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
4
# donc ici test l'absence d'effet de l'espèce, du sexe, et de l'interaction espèce*sexe.
# Toutes les variables sont significatices au niveau 5% (p-valeurs inférieures à 5%)
anova(modeleCst,modeleComp) # on préfère le modèle 2 donc le modèle complet
## Analysis of Variance Table
##
## Model 1: RW ~ 1
## Model 2: RW ~ sp + sex + sp:sex
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 199 1317.8
## 2 196 1016.4 3 301.43 19.376 4.813e-11 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
anova(modeleAdd,modeleComp) # on préfère aussi le modèle complet
## Analysis of Variance Table
##
## Model 1: RW ~ sp + sex
## Model 2: RW ~ sp + sex + sp:sex
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 197 1074.4
## 2 196 1016.4 1 57.996 11.184 0.0009884 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Remarque : par défaut, les sorties de la fonction lm n’utilisent pas les mêmes contraintes que celles vues en
cours. Dans notre exemple d’étude d’une interaction de 2 facteurs, chacun à 2 niveaux, l’un des 4 niveaux
possibles est choisi comme niveau de référence, puis tous les autres effets sont calculés par rapport à celui-ci.
Si on reprend la sortie du modèle modeleComp, on obtient :
modeleComp$coefficients
## (Intercept) spO sexM spO:sexM
## 12.138 2.698 -0.420 -2.154
et si l’on compare cela à la moyenne obtenue dans chacun des quatre groupes formés par l’interaction des
deux facteurs :
tapply(RW,sp:sex,mean) # moyenne selon les valeurs de l'interaction sp:sex
## B:F B:M O:F O:M
## 12.138 11.718 14.836 12.262
on retrouve que la catégorie B:F a été prise comme référence (on retrouve sa moyenne dans l’intercept), puis
on obtient les moyennes dans chaque groupe de la façon suivante :
modeleComp$coefficients[1] # catégorie de référence B:F
## (Intercept)
## 12.138
5
# on rajoute l'effet "spO", c'est-à-dire que l'on passe de la catégorie B:F à la catégorie # O:F
modeleComp$coefficients[1] + modeleComp$coefficients[2]
## (Intercept)
## 14.836
# on rajoute l'effet "sexM", c'est-à-dire que l'on passe de la catégorie B:F à la catégorie
# B:M
modeleComp$coefficients[1] + modeleComp$coefficients[3]
## (Intercept)
## 11.718
# on rajoute l'effet "spO" , l'effet "sexM", et l'effet de l'interaction "sp0:sexM",
# c'est-à-dire que l'on passe de la catégorie B:F à la catégorie O:M
modeleComp$coefficients[1] + modeleComp$coefficients[2] + modeleComp$coefficients[3] +
modeleComp$coefficients[4]
## (Intercept)
## 12.262
Pour obtenir les mêmes contraintes que dans le cours, il faut utiliser l’option suivante :
modeleBis <- lm(RW ~ sp + sex + sp:sex, contrasts = list(sp=contr.sum, sex=contr.sum))
summary(modeleBis)
##
## Call:
## lm(formula = RW ~ sp + sex + sp:sex, contrasts = list(sp = contr.sum,
## sex = contr.sum))
##
## Residuals:
## Min 1Q Median 3Q Max
## -5.638 -1.280 0.000 1.688 5.364
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 12.7385 0.1610 79.111 < 2e-16 ***
## sp1 -0.8105 0.1610 -5.034 1.09e-06 ***
## sex1 0.7485 0.1610 4.648 6.13e-06 ***
## sp1:sex1 -0.5385 0.1610 -3.344 0.000988 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 2.277 on 196 degrees of freedom
## Multiple R-squared: 0.2287, Adjusted R-squared: 0.2169
## F-statistic: 19.38 on 3 and 196 DF, p-value: 4.813e-11
Dans ce cas, on obtient directement la moyenne dans chaque groupe en utilisant individuellement chaque
coefficient :
6
modeleBis$coefficients[1] # moyenne dans toute la population
## (Intercept)
## 12.7385
modeleBis$coefficients[1] + modeleBis$coefficients[2] # moyenne dans le groupe sp=B
## (Intercept)
## 11.928
modeleBis$coefficients[1] - modeleBis$coefficients[2] # moyenne dans le groupe sp=O
## (Intercept)
## 13.549
modeleBis$coefficients[1] + modeleBis$coefficients[3] # moyenne dans le groupe sex=F
## (Intercept)
## 13.487
modeleBis$coefficients[1] - modeleBis$coefficients[3] # moyenne dans le groupe sex=M
## (Intercept)
## 11.99
# moyenne dans le groupe sp=B et sex=F, en comptant l'interaction B:F
modeleBis$coefficients[1] + modeleBis$coefficients[2] + modeleBis$coefficients[3] +
modeleBis$coefficients[4]
## (Intercept)
## 12.138
# moyenne dans le groupe sp=B et sex=M, moins l'interaction B:M (signe moins devant le
# coefficient pour le sexe donc aussi devant l'interaction)
modeleBis$coefficients[1] + modeleBis$coefficients[2] - modeleBis$coefficients[3] -
modeleBis$coefficients[4]
## (Intercept)
## 11.718
# moyenne dans le groupe sp=O et sex=F, moins l'interaction O:F (signe moins devant le
# coefficient pour l'espèce donc aussi devant l'interaction)
modeleBis$coefficients[1] - modeleBis$coefficients[2] + modeleBis$coefficients[3] -
modeleBis$coefficients[4]
## (Intercept)
## 14.836
7
# moyenne dans le groupe sp=O et sex=M, plus l'interaction O:M (signe moins devant les
# deux coefficients pour le sexe et l'espèce donc signe plus devant l'interaction)
modeleBis$coefficients[1] - modeleBis$coefficients[2] - modeleBis$coefficients[3] +
modeleBis$coefficients[4]
## (Intercept)
## 12.262
6. Que donne le crière BIC ?
# Critere BIC
library(stats4)
BIC(modeleCst)
## [1] 955.2514
BIC(modeleComp)
## [1] 919.201
BIC(modeleAdd)
## [1] 925.0015
7. Vérifier à l’aide d’un test de Shapiro-Wilk que les variables observées peuvent être considérées de loi
gaussienne.
shapiro.test(RW[sp:sex=="B:F"])
##
## Shapiro-Wilk normality test
##
## data: RW[sp:sex == "B:F"]
## W = 0.98459, p-value = 0.7538
shapiro.test(RW[sp:sex=="O:F"])
##
## Shapiro-Wilk normality test
##
## data: RW[sp:sex == "O:F"]
## W = 0.98038, p-value = 0.5683
shapiro.test(RW[sp:sex=="B:M"])
##
## Shapiro-Wilk normality test
##
## data: RW[sp:sex == "B:M"]
## W = 0.97707, p-value = 0.4361
8
shapiro.test(RW[sp:sex=="O:M"])
##
## Shapiro-Wilk normality test
##
## data: RW[sp:sex == "O:M"]
## W = 0.98691, p-value = 0.8496
8. Vérifier l’hypothèse d’homoscédasticité (égalité des variances).
# Hypothèse nulle = égalité des variances
var.test(RW[sp:sex=="B:F"], RW[sp:sex=="O:F"]) # p-valeur > 5% -> on ne rejette pas
##
## F test to compare two variances
##
## data: RW[sp:sex == "B:F"] and RW[sp:sex == "O:F"]
## F = 1.0784, num df = 49, denom df = 49, p-value = 0.7927
## alternative hypothesis: true ratio of variances is not equal to 1
## 95 percent confidence interval:
## 0.6119582 1.9003191
## sample estimates:
## ratio of variances
## 1.078386
var.test(RW[sp:sex=="B:F"], RW[sp:sex=="B:M"]) # p-valeur > 5% -> on ne rejette pas
##
## F test to compare two variances
##
## data: RW[sp:sex == "B:F"] and RW[sp:sex == "B:M"]
## F = 1.3338, num df = 49, denom df = 49, p-value = 0.3168
## alternative hypothesis: true ratio of variances is not equal to 1
## 95 percent confidence interval:
## 0.7568757 2.3503326
## sample estimates:
## ratio of variances
## 1.333758
var.test(RW[sp:sex=="B:F"], RW[sp:sex=="O:M"]) # p-valeur > 5% -> on ne rejette pas
##
## F test to compare two variances
##
## data: RW[sp:sex == "B:F"] and RW[sp:sex == "O:M"]
## F = 1.2337, num df = 49, denom df = 49, p-value = 0.465
## alternative hypothesis: true ratio of variances is not equal to 1
## 95 percent confidence interval:
## 0.7000855 2.1739814
## sample estimates:
## ratio of variances
## 1.233683
9
var.test(RW[sp:sex=="B:M"], RW[sp:sex=="O:F"]) # p-valeur > 5% -> on ne rejette pas
##
## F test to compare two variances
##
## data: RW[sp:sex == "B:M"] and RW[sp:sex == "O:F"]
## F = 0.80853, num df = 49, denom df = 49, p-value = 0.4597
## alternative hypothesis: true ratio of variances is not equal to 1
## 95 percent confidence interval:
## 0.4588226 1.4247858
## sample estimates:
## ratio of variances
## 0.808532
var.test(RW[sp:sex=="B:M"], RW[sp:sex=="O:M"]) # p-valeur > 5% -> on ne rejette pas
##
## F test to compare two variances
##
## data: RW[sp:sex == "B:M"] and RW[sp:sex == "O:M"]
## F = 0.92497, num df = 49, denom df = 49, p-value = 0.786
## alternative hypothesis: true ratio of variances is not equal to 1
## 95 percent confidence interval:
## 0.524897 1.629967
## sample estimates:
## ratio of variances
## 0.9249676
var.test(RW[sp:sex=="O:F"], RW[sp:sex=="O:M"]) # p-valeur > 5% -> on ne rejette pas
##
## F test to compare two variances
##
## data: RW[sp:sex == "O:F"] and RW[sp:sex == "O:M"]
## F = 1.144, num df = 49, denom df = 49, p-value = 0.6396
## alternative hypothesis: true ratio of variances is not equal to 1
## 95 percent confidence interval:
## 0.6491976 2.0159589
## sample estimates:
## ratio of variances
## 1.144009
9. Etudier les résidus du modèle sélectionné ci-dessus.
residus=modeleComp$residuals
hist(residus,main="Histogramme des résidus",freq=F)
curve(dnorm(x,0,sqrt(var(residus))),add=T)
10