#Chapitre 9
require(BioStatR)
## Le chargement a nécessité le package : BioStatR
#page 378
#Exercice 9.1
#1)
lauriers<-subset(Mesures,subset=(Mesures$espece=="laurier rose"))
plot(taille~masse,data=lauriers,pch=19)
#page 379
#3)
droite_lauriers<-lm(taille~masse,data=lauriers)
coef(droite_lauriers)
## (Intercept) masse
## 6.413523 1.700114
#4)
fitted(droite_lauriers)
## 181 182 183 184 185 186 187 188
## 14.74408 16.95423 13.21398 12.02390 14.57407 15.93416 14.06404 17.12424
## 189 190 191 192 193 194 195 196
## 13.55400 13.04397 16.27419 14.40406 16.61421 17.46427 14.91409 15.76415
## 197 198 199 200 201 202 203 204
## 14.40406 16.10418 12.53393 15.59414 15.42413 14.91409 14.06404 13.89403
## 205 206 207 208 209 210 211 212
## 14.57407 14.06404 11.85389 14.40406 13.21398 16.27419 15.76415 13.89403
## 213 214 215 216 217 218 219 220
## 12.36392 13.89403 13.72401 13.38399 15.42413 14.40406 15.42413 14.40406
## 221 222 223 224 225 226 227 228
## 14.74408 13.38399 14.23405 14.57407 12.19391 12.19391 16.27419 14.57407
## 229 230 231 232 233 234 235 236
## 13.04397 12.19391 14.06404 12.02390 12.02390 12.53393 12.36392 12.87396
## 237 238 239 240 241 242 243 244
## 11.85389 12.87396 15.42413 16.27419 14.23405 11.85389 13.72401 11.00383
## 245 246 247 248 249 250 251 252
## 10.83382 10.49380 10.83382 11.85389 17.29426 12.19391 12.19391 11.00383
#page 380
#5)
abline(coef(droite_lauriers),col="red",lwd=2)
#6)
predict(droite_lauriers,data.frame(masse=4.8))
## 1
## 14.57407
#fonctionne comme predict(droite_lauriers,list(masse=4.8))
#7)
residuals(droite_lauriers)[lauriers$masse==4.8]
## 185 205 224 228
## 0.52592789 0.62592789 -0.07407211 0.52592789
#page 381
#8)
mean(lauriers$taille)
## [1] 13.91528
6.413523+1.700114*mean(lauriers$masse)
## [1] 13.91528
coef(droite_lauriers)[1]+coef(droite_lauriers)[2]*mean(lauriers$masse)
## (Intercept)
## 13.91528
#9)
summary(droite_lauriers)
##
## Call:
## lm(formula = taille ~ masse, data = lauriers)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.1339 -0.8590 0.1310 0.9259 2.3060
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 6.4135 0.5924 10.83 <2e-16 ***
## masse 1.7001 0.1309 12.99 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.12 on 70 degrees of freedom
## Multiple R-squared: 0.7068, Adjusted R-squared: 0.7026
## F-statistic: 168.8 on 1 and 70 DF, p-value: < 2.2e-16
#page 382
#10)
anova(droite_lauriers)
## Analysis of Variance Table
##
## Response: taille
## Df Sum Sq Mean Sq F value Pr(>F)
## masse 1 211.57 211.573 168.76 < 2.2e-16 ***
## Residuals 70 87.76 1.254
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#11)
summary(droite_lauriers)
##
## Call:
## lm(formula = taille ~ masse, data = lauriers)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.1339 -0.8590 0.1310 0.9259 2.3060
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 6.4135 0.5924 10.83 <2e-16 ***
## masse 1.7001 0.1309 12.99 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.12 on 70 degrees of freedom
## Multiple R-squared: 0.7068, Adjusted R-squared: 0.7026
## F-statistic: 168.8 on 1 and 70 DF, p-value: < 2.2e-16
#page 383
#12)
residus<-residuals(droite_lauriers)
shapiro.test(residus)
##
## Shapiro-Wilk normality test
##
## data: residus
## W = 0.96531, p-value = 0.04439
#page 384
plot(lauriers$masse,residus)
pdf("residusmasse.pdf")
plot(lauriers$masse,residus)
dev.off()
## quartz_off_screen
## 2
#Les r\'esidus ont l'air corrects => homosc\'edasticit\'e des erreurs ok et
#absence d'effet syst\'ematique
#Approche par permutation valide
#13)
if(!("lmPerm" %in% rownames(installed.packages()))){install.packages("lmPerm")}
library(lmPerm)
lmp(taille~masse,lauriers)
## [1] "Settings: unique SS : numeric variables centered"
##
## Call:
## lmp(formula = taille ~ masse, data = lauriers)
##
## Coefficients:
## (Intercept) masse
## 13.92 1.70
#page 385
perm_laurier<-lmp(taille~masse,lauriers,center=FALSE)
## [1] "Settings: unique SS "
summary(perm_laurier)
##
## Call:
## lmp(formula = taille ~ masse, data = lauriers, center = FALSE)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.1339 -0.8590 0.1309 0.9259 2.3060
##
## Coefficients:
## Estimate Iter Pr(Prob)
## masse 1.7 5000 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.12 on 70 degrees of freedom
## Multiple R-Squared: 0.7068, Adjusted R-squared: 0.7026
## F-statistic: 168.8 on 1 and 70 DF, p-value: < 2.2e-16
#page 386
#14)
confint(droite_lauriers)
## 2.5 % 97.5 %
## (Intercept) 5.232099 7.594947
## masse 1.439098 1.961131
predict(droite_lauriers,list(masse=c(4.8)),interval="confidence")
## fit lwr upr
## 1 14.57407 14.29212 14.85602
predict(droite_lauriers,list(masse=c(4.8)),interval="prediction")
## fit lwr upr
## 1 14.57407 12.32318 16.82496
#page 387
#Exercice 9.2
#1)
bignones<-subset(Mesures5,subset=(Mesures5$espece=="bignone"))[,c(1,4)]
plot(masse~masse_sec,data=bignones,pch=19)
pdf("figure94.pdf")
plot(masse~masse_sec,data=bignones,pch=19)
dev.off()
## quartz_off_screen
## 2
#3)a)
droite_bignones<-lm(masse~masse_sec,data=bignones)
coef(droite_bignones)
## (Intercept) masse_sec
## -0.5391407 4.8851935
#page 388
residus<-residuals(droite_bignones)
plot(bignones$masse_sec,residus)
pdf("figure95.pdf")
plot(bignones$masse_sec,residus)
dev.off()
## quartz_off_screen
## 2
#Les r\'esidus n'ont l'air corrects car ils pr\'esentent une forme en trompette,
#ce qui remet en question de l'homosc\'edasticit\'e des erreurs. Nous proc\'ederons
#dans la suite \`a un test pour nous assurer que ce d\'efaut est significatif au
#seuil de \alpha=5%. Par contre les r\'esidus ont l'air r\'epartis al\'eatoirement
#au-dessus ou en-dessous de l'axe des abscisses. Vous notez \'egalement l'absence
#d'un effet syst\'ematique qui se traduirait par exemple par une forme de banane.
#L'hypoth\`ese d'ind\'ependance n'est pas remise en question.
#Malgr\'e l'inhomog\'en\'eit\'e des variances l'estimation de la pente et de l'ordonn\'ee
#\`a l'origine reste sans biais. Il sera, par contre, n\'ecessaire tenir compte de
#l'h\'et\'erosc\'edastict\'e des erreurs pour la mise en oeuvre des proc\'edures de test et
#la construction des intervalles de confiance.
#page 389
#4)
fitted(droite_bignones)
## 111 112 113 114 115 116 117 118
## 14.116440 4.834572 16.070517 31.703136 16.559037 20.467191 9.719766 9.719766
## 119 120 121 122 123 124 125 126
## 11.185324 10.208285 15.093478 7.765688 7.277169 1.903456 6.788650 8.742727
## 127 128 129 130 131 132 133 134
## 11.185324 10.696804 13.627920 6.788650 7.277169 0.437898 6.300130 10.208285
## 135 136 137 138 139 140 141 142
## 8.742727 4.834572 5.811611 3.369014 4.834572 10.696804 16.559037 7.765688
## 143 144 145 146 147 148 149 150
## 9.231246 7.765688 11.673843 7.765688 5.323091 7.765688 3.857533 2.880495
## 151 152 153 154 155 156 157 158
## 4.834572 2.391975 4.346053 5.323091 4.834572 6.788650 5.323091 5.811611
## 159 160 161 162 163 164 165 166
## 3.857533 1.903456 1.903456 2.391975 1.414937 1.414937 3.369014 3.857533
## 167 168 169 170 171 172 173 174
## 3.857533 9.231246 4.834572 5.811611 4.346053 12.162362 6.788650 4.346053
## 175 176 177 178 179 180
## 9.231246 4.834572 6.788650 6.788650 6.300130 3.857533
#5)
plot(masse~masse_sec,data=bignones,pch=19)
abline(coef(droite_bignones),col="red",lwd=2)
pdf("figure96.pdf")
plot(masse~masse_sec,data=bignones,pch=19)
abline(coef(droite_bignones),col="red",lwd=2)
dev.off()
## quartz_off_screen
## 2
#6)
predict(droite_bignones,data.frame(masse_sec=2.5))
## 1
## 11.67384
plot(masse~masse_sec,data=bignones,pch=19)
abline(coef(droite_bignones),col="red",lwd=2)
points(2.5,predict(droite_bignones,data.frame(masse_sec=2.5)),pch=17,col="blue")
segments(2.5, bignones$masse[bignones$masse_sec==2.5],2.5,
predict(droite_bignones,data.frame(masse_sec=2.5)),lty=2,lwd=2)
pdf("figure96residusmasselinepoint.pdf")
plot(masse~masse_sec,data=bignones,pch=19)
abline(coef(droite_bignones),col="red",lwd=2)
points(2.5,predict(droite_bignones,data.frame(masse_sec=2.5)),pch=17,col="blue")
segments(2.5, bignones$masse[bignones$masse_sec==2.5],2.5,
predict(droite_bignones,data.frame(masse_sec=2.5)),lty=2,lwd=2)
dev.off()
## quartz_off_screen
## 2
#page 390
#7)
residuals(droite_bignones)[bignones$masse_sec==2.5]
## 145
## -2.673843
#8)
mean(bignones$masse)
## [1] 7.521429
-0.5391407+4.8851935*mean(bignones$masse_sec)
## [1] 7.521429
coef(droite_bignones)[1]+coef(droite_bignones)[2]*mean(bignones$masse_sec)
## (Intercept)
## 7.521429
#page 391
#9)
summary(droite_bignones)
##
## Call:
## lm(formula = masse ~ masse_sec, data = bignones)
##
## Residuals:
## Min 1Q Median 3Q Max
## -4.6279 -1.7444 -0.2961 1.4310 6.4295
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.5391 0.5657 -0.953 0.344
## masse_sec 4.8852 0.2912 16.774 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 2.497 on 68 degrees of freedom
## Multiple R-squared: 0.8054, Adjusted R-squared: 0.8025
## F-statistic: 281.4 on 1 and 68 DF, p-value: < 2.2e-16
#10)
anova(droite_bignones)
## Analysis of Variance Table
##
## Response: masse
## Df Sum Sq Mean Sq F value Pr(>F)
## masse_sec 1 1754.44 1754.44 281.36 < 2.2e-16 ***
## Residuals 68 424.01 6.24
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#page 392
#12) et 13)
residus<-residuals(droite_bignones)
shapiro.test(residus)
##
## Shapiro-Wilk normality test
##
## data: residus
## W = 0.98209, p-value = 0.4161
length(residus)
## [1] 70
#Les r\'esidus sont au nombre de 70 sup\'erieur ou \'egal \`a 30. Le test de normalit\'e
#est donc fiable. La $p$-valeur du test est strictement sup\'erieure \`a \alpha=5%,
#le test n'est pas significatif. Nous conservons, par d\'efaut, l'hypoth\`ese H0
#de normalit\'e des erreurs.
#page 393
#Le test de White est un cas particulier du test de Breusch-Pagan qui est
#disponible dans le biblioth\`eque lmtest
if(!("lmtest" %in% rownames(installed.packages()))){install.packages("lmtest")}
library(lmtest)
## Le chargement a nécessité le package : zoo
##
## Attachement du package : 'zoo'
## Les objets suivants sont masqués depuis 'package:base':
##
## as.Date, as.Date.numeric
bptest(droite_bignones, ~ masse_sec + I(masse_sec^2), data = bignones)
##
## studentized Breusch-Pagan test
##
## data: droite_bignones
## BP = 20.238, df = 2, p-value = 4.031e-05
## White test (Table 5.1, p. 113)
#bptest(cig_lm2, ~ income * price + I(income^2) + I(price^2), data = CigarettesB)
#Le test de White permet de s'int\'eresser aux deux hypoth\`eses :
#"H0 : les erreurs sont homosc\'edastiques"
#contre
#"H1 : les erreurs sont h\'et\'erosc\'edastiques".
#L'hypoth\`ese de normalit\'e des erreurs n'a \'et\'e remise en cause, le test de White
#est donc fiable. La $p$-valeur du test est inf\'erieure ou \'egale \`a \alpha=5%,
#le test est significatif. Nous rejetons l'hypoth\`ese H0 d'homosc\'edasticit\'e
#des erreurs et d\'ecidons que l'hypoth\`ese alternative d'h\'et\'erosc\'edasticit\'e
#des erreurs est vraie.
#Comme nous l'avions per\c cu graphiquement, les erreurs ne sont pas homosc\'edastiques,
#il faut tenir compte de cette inhomog\'en\'eit\'e des variances lors de l'estimation
#des param\`etres du mod\`ele puis de la mise en oeuvre des tests de student ou
#du test global de Fisher pour la r\'egression.
if(!("sandwich" %in% rownames(installed.packages()))){install.packages("sandwich")}
library(sandwich)
vcovHC(droite_bignones)
## (Intercept) masse_sec
## (Intercept) 0.2782291 -0.1714076
## masse_sec -0.1714076 0.1397390
#Estimation, tenant de l'inhomog\'en\'eit\'e des variances, de la matrice de
#variance-covariance des estimateurs \hat\beta_0 et \hat\beta_1.
coeftest(droite_bignones, df="inf", vcov=vcovHC)
##
## z test of coefficients:
##
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.53914 0.52747 -1.0221 0.3067
## masse_sec 4.88519 0.37382 13.0684 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#Tests de student des coefficient \beta_0 et \beta_1.
#page 394
waldtest(droite_bignones, vcov=vcovHC)
## Wald test
##
## Model 1: masse ~ masse_sec
## Model 2: masse ~ 1
## Res.Df Df F Pr(>F)
## 1 68
## 2 69 -1 170.78 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#Tests de Fihser global du mod\`ele de r\'egression lin\'eaire simple.
#Pour construire les intervalles de confiance autour des param\`etres,
#vous poouvez utiliser la biblioth\`eque hcci.
if(!("hcci" %in% rownames(installed.packages()))){install.packages("hcci")}
library(hcci)
?hcci
#L'aide de la biblioth\`eque HCCI vous apprend qu'il existe plusieurs proc\'edures
#permettant de tenir compte de l'h\'et\'erosc\'edasticit\'e. La fonction vcovHC utilise
#la m\'ethode HC3 par d\'efaut La fonction HC, la m\'ethode HC4 avec le param\`etre k=0.7
#par d\'efaut. Les m\'ethodes HC3, HC4 et HC5 sont recommend\'ees. En comparant leurs
#r\'esultats, vous constatez qu'elles aboutissent toutes aux m^emes conclusions
#au seuil de \alpha=5% : conservation, par d\'efaut, de "H0 : \beta_0=0" pour
#le test de l'ordonn\'ee \`a l'origine et d\'ecision que "H1 : \beta_1<>0" est vraie.
HC(droite_bignones,method=3)
## [,1] [,2]
## [1,] 0.2782291 -0.1714076
## [2,] -0.1714076 0.1397390
coeftest(droite_bignones, df="inf", vcov=HC(droite_bignones,method=3))
##
## z test of coefficients:
##
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.53914 0.52747 -1.0221 0.3067
## masse_sec 4.88519 0.37382 13.0684 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#page 395
vcovHC(droite_bignones,type="HC4")
## (Intercept) masse_sec
## (Intercept) 0.4035131 -0.2603917
## masse_sec -0.2603917 0.2022249
coeftest(droite_bignones, df="inf", vcov=vcovHC(droite_bignones,type="HC4"))
##
## z test of coefficients:
##
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.53914 0.63523 -0.8487 0.396
## masse_sec 4.88519 0.44969 10.8634 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
vcovHC(droite_bignones,type="HC4m")
## (Intercept) masse_sec
## (Intercept) 0.3020561 -0.1891167
## masse_sec -0.1891167 0.1526165
coeftest(droite_bignones, df="inf", vcov=vcovHC(droite_bignones,type="HC4m"))
##
## z test of coefficients:
##
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.53914 0.54960 -0.981 0.3266
## masse_sec 4.88519 0.39066 12.505 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#page 396
HC(droite_bignones,method=4,k=0.7)
## [,1] [,2]
## [1,] 0.4035131 -0.2603917
## [2,] -0.2603917 0.2022249
coeftest(droite_bignones, df="inf", vcov=HC(droite_bignones,method=4,k=0.7))
##
## z test of coefficients:
##
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.53914 0.63523 -0.8487 0.396
## masse_sec 4.88519 0.44969 10.8634 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
vcovHC(droite_bignones,type="HC5")
## (Intercept) masse_sec
## (Intercept) 0.4141695 -0.2662869
## masse_sec -0.2662869 0.2047638
coeftest(droite_bignones, df="inf", vcov=vcovHC(droite_bignones,type="HC5"))
##
## z test of coefficients:
##
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.53914 0.64356 -0.8377 0.4022
## masse_sec 4.88519 0.45251 10.7958 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
HC(droite_bignones,method=5)
## [,1] [,2]
## [1,] 0.4141695 -0.2662869
## [2,] -0.2662869 0.2047638
#page 397
coeftest(droite_bignones, df="inf", vcov=HC(droite_bignones,method=5))
##
## z test of coefficients:
##
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.53914 0.64356 -0.8377 0.4022
## masse_sec 4.88519 0.45251 10.7958 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#Passons \`a la construction d'intervalles de confiance sur les param\`etres
#\beta_0 et \beta_1 de la r\'egression lin\'eaire simple. Nous devons passer par
#cette \'etape de r\'e\'ecriture du mod\`ele pour pouvoir utiliser les fonctions Pboot
#et Tbootde la biblioth\`eque hcci.
y = bignones$masse
x = bignones$masse_sec
model = lm(y ~ x)
#Il est possible de "fixer" le point de d\'epart du g\'en\'erateur al\'eatoir
#pour avoir des r\'esultats reproductibles \`a l'aide de la fonction set.seed
set.seed(123456)
#Commencez par utiliser une technique de bootstrap simple.
#Bootstrap percentile simple.
Pboot(model, significance = 0.05, double = FALSE, J=1000, K = 100,
distribution = "rademacher")
## $beta
## [1] -0.5391407 4.8851935
##
## $ci_lower_simple
## [1] -1.492404 4.213572
##
## $ci_upper_simple
## [1] 0.4458551 5.5789370
##
## $ci_lower_double
## logical(0)
##
## $ci_upper_double
## logical(0)
#page 398
#Bootstrap t simple.
Tboot(model, significance = 0.05, double = FALSE, J=1000, K = 100,
distribution = "rademacher")
## $beta
## [1] -0.5391407 4.8851935
##
## $ci_lower_simple
## [1] -1.737704 4.001968
##
## $ci_upper_simple
## [1] 0.7309306 5.7797771
##
## $ci_lower_double
## logical(0)
##
## $ci_upper_double
## logical(0)
##
## $J
## [1] 1000
##
## $K
## [1] 100
#Utilisez maintenant une technique de bootstrap double.
#Bootstrap percentile double.
Pboot(model, significance = 0.05, double = TRUE, J=1000, K = 100,
distribution = "rademacher")
## $beta
## [1] -0.5391407 4.8851935
##
## $ci_lower_simple
## [1] -1.523889 4.220719
##
## $ci_upper_simple
## [1] 0.4474727 5.5874831
##
## $ci_lower_double
## [1] -1.657840 3.994579
##
## $ci_upper_double
## [1] 0.5401815 5.9500533
#Bootstrap t double.
Tboot(model, significance = 0.05, double = TRUE, J=1000, K = 100,
distribution = "rademacher")
## $beta
## [1] -0.5391407 4.8851935
##
## $ci_lower_simple
## [1] -1.730011 3.960551
##
## $ci_upper_simple
## [1] 0.743166 5.781661
##
## $ci_lower_double
## [1] -2.253324 3.854114
##
## $ci_upper_double
## [1] 1.248126 5.904192
##
## $J
## [1] 1000
##
## $K
## [1] 100
#Le mod\`ele \'etant h\'et\'erosc\'edastique, la construction d'intervalles de pr\'ediction
#n'est pas fiable