#Chapitre 4
#page 248
#Si le package n'est pas installe, enlever le commentaire
#puis executer la commande ci-dessous.
#if(!require("BioStatR")){install.packages("BioStatR")}
library(BioStatR)
str(Mesures)
## 'data.frame': 252 obs. of 3 variables:
## $ masse : num 28.6 20.6 29.2 32 24.5 29 28.9 18.2 7.9 15.5 ...
## $ taille: num 19.1 14.8 19.7 21.1 19.4 19.5 18.9 14.6 10.2 14.6 ...
## $ espece: Factor w/ 4 levels "bignone","glycine blanche",..: 2 2 2 2 2 2 2 2 2 2 ...
summary(Mesures)
## masse taille espece
## Min. : 1.00 Min. : 4.80 bignone :70
## 1st Qu.: 4.50 1st Qu.:11.00 glycine blanche :54
## Median : 8.40 Median :13.20 glycine violette:56
## Mean :11.13 Mean :13.37 laurier rose :72
## 3rd Qu.:14.60 3rd Qu.:15.30
## Max. :49.20 Max. :27.00
#page 249
#Si le package n'est pas installe, enlever le commentaire
#puis executer la commande ci-dessous.
#if(!require("ISLR")){install.packages("ISLR")}
library(ISLR)
summary(Hitters[,17:20])
## Assists Errors Salary NewLeague
## Min. : 0.0 Min. : 0.00 Min. : 67.5 A:176
## 1st Qu.: 7.0 1st Qu.: 3.00 1st Qu.: 190.0 N:146
## Median : 39.5 Median : 6.00 Median : 425.0
## Mean :106.9 Mean : 8.04 Mean : 535.9
## 3rd Qu.:166.0 3rd Qu.:11.00 3rd Qu.: 750.0
## Max. :492.0 Max. :32.00 Max. :2460.0
## NA's :59
Mes.B = subset(Mesures,Mesures$espece=="bignone")
model<-lm(masse~taille,data=Mes.B)
summary(model)
##
## Call:
## lm(formula = masse ~ taille, data = Mes.B)
##
## Residuals:
## Min 1Q Median 3Q Max
## -4.137 -1.551 -0.455 1.529 10.699
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -10.01645 1.06941 -9.366 7.39e-14 ***
## taille 1.53572 0.09002 17.059 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 2.463 on 68 degrees of freedom
## Multiple R-squared: 0.8106, Adjusted R-squared: 0.8078
## F-statistic: 291 on 1 and 68 DF, p-value: < 2.2e-16
Mes.B$masse
## [1] 10.9 6.6 22.5 33.7 20.6 16.6 14.2 13.8 14.0 8.7 14.2 10.6 10.9 3.3
## [15] 9.7 9.3 17.2 10.1 9.0 7.1 7.1 1.5 4.1 8.0 7.4 7.2 6.9 2.9
## [29] 2.4 10.7 13.8 10.9 10.3 8.8 9.0 8.2 9.6 9.0 5.3 1.5 6.7 2.9
## [43] 2.9 3.5 3.4 4.9 4.7 4.7 5.2 2.1 2.2 1.4 2.7 1.0 2.5 5.5
## [57] 2.7 6.7 7.3 2.9 3.8 7.6 3.6 3.0 5.8 5.3 3.2 4.4 3.4 2.9
#page 250
round(fitted(model), 2)
## 111 112 113 114 115 116 117 118 119 120 121 122
## 9.64 6.11 17.63 23.00 17.32 15.02 16.71 14.09 16.55 10.56 11.33 12.41
## 123 124 125 126 127 128 129 130 131 132 133 134
## 11.79 3.65 9.95 8.72 15.32 12.56 10.56 7.80 9.79 -0.03 4.57 8.11
## 135 136 137 138 139 140 141 142 143 144 145 146
## 10.87 9.79 7.80 4.42 4.57 11.48 10.72 8.57 7.80 10.56 6.72 8.72
## 147 148 149 150 151 152 153 154 155 156 157 158
## 10.56 9.18 6.11 0.73 10.26 4.27 2.12 5.80 1.50 3.04 5.49 2.73
## 159 160 161 162 163 164 165 166 167 168 169 170
## 6.57 2.73 0.89 -0.19 -0.03 -2.65 1.35 4.27 3.19 5.19 11.33 3.34
## 171 172 173 174 175 176 177 178 179 180
## 4.27 11.02 2.73 2.42 7.18 8.41 7.34 0.27 6.26 3.65
round(residuals(model), 2)
## 111 112 113 114 115 116 117 118 119 120 121 122
## 1.26 0.49 4.87 10.70 3.28 1.58 -2.51 -0.29 -2.55 -1.86 2.87 -1.81
## 123 124 125 126 127 128 129 130 131 132 133 134
## -0.89 -0.35 -0.25 0.58 1.88 -2.46 -1.56 -0.70 -2.69 1.53 -0.47 -0.11
## 135 136 137 138 139 140 141 142 143 144 145 146
## -3.47 -2.59 -0.90 -1.52 -2.17 -0.78 3.08 2.33 2.50 -1.76 2.28 -0.52
## 147 148 149 150 151 152 153 154 155 156 157 158
## -0.96 -0.18 -0.81 0.77 -3.56 -1.37 0.78 -2.30 1.90 1.86 -0.79 1.97
## 159 160 161 162 163 164 165 166 167 168 169 170
## -1.37 -0.63 1.31 1.59 2.73 3.65 1.15 1.23 -0.49 1.51 -4.03 -0.44
## 171 172 173 174 175 176 177 178 179 180
## -0.47 -3.42 0.87 0.58 -1.38 -3.11 -4.14 4.13 -2.86 -0.75
#Si le package n'est pas installe, enlever le commentaire
#puis executer la commande ci-dessous.
#if(!require("broom")){install.packages("broom")}
library(broom)
(model.diag.metrics<-augment(model))
## # A tibble: 70 x 10
## .rownames masse taille .fitted .se.fit .resid .hat .sigma .cooksd
## <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 111 10.9 12.8 9.64 0.320 1.26 0.0168 2.48 2.27e-3
## 2 112 6.6 10.5 6.11 0.306 0.491 0.0154 2.48 3.16e-4
## 3 113 22.5 18 17.6 0.661 4.87 0.0721 2.40 1.64e-1
## 4 114 33.7 21.5 23.0 0.954 10.7 0.150 2.04 1.96e+0
## 5 115 20.6 17.8 17.3 0.645 3.28 0.0687 2.45 7.02e-2
## 6 116 16.6 16.3 15.0 0.529 1.58 0.0461 2.47 1.05e-2
## 7 117 14.2 17.4 16.7 0.614 -2.51 0.0620 2.46 3.65e-2
## 8 118 13.8 15.7 14.1 0.485 -0.294 0.0388 2.48 2.99e-4
## 9 119 14 17.3 16.6 0.606 -2.55 0.0605 2.46 3.67e-2
## 10 120 8.7 13.4 10.6 0.344 -1.86 0.0195 2.47 5.80e-3
## # … with 60 more rows, and 1 more variable: .std.resid <dbl>
#page 252
tidy(model)
## # A tibble: 2 x 5
## term estimate std.error statistic p.value
## <chr> <dbl> <dbl> <dbl> <dbl>
## 1 (Intercept) -10.0 1.07 -9.37 7.39e-14
## 2 taille 1.54 0.0900 17.1 2.89e-26
glance(model)
## # A tibble: 1 x 11
## r.squared adj.r.squared sigma statistic p.value df logLik AIC BIC
## <dbl> <dbl> <dbl> <dbl> <dbl> <int> <dbl> <dbl> <dbl>
## 1 0.811 0.808 2.46 291. 2.89e-26 2 -161. 329. 336.
## # … with 2 more variables: deviance <dbl>, df.residual <int>
with(Mes.B,plot(taille,masse,pch=20,xlab="Taille",ylab="Masse"))
abline(model, lwd=2, col="red")
with(Mes.B,lines(lowess(taille, masse), lty=2, lwd=2, col="blue"))
with(Mes.B,lines(smooth.spline(taille,masse),col="orange",lwd=2,
lty=4))
with(Mes.B,segments(taille, masse, taille, fitted(model), lty=2,
col="red"))
legend("topleft",lty=c(1,2,4),legend = c("lm","lowess",
"smooth.spline"),lwd=2,col=c("red","blue","orange"))
#page 253
#Si le package n'est pas installe, enlever le commentaire
#puis executer la commande ci-dessous.
#if(!require("ggiraphExtra")){install.packages("ggiraphExtra")}
library(ggiraphExtra)
ggPredict(model)
shapiro.test(residuals(model))
##
## Shapiro-Wilk normality test
##
## data: residuals(model)
## W = 0.92886, p-value = 0.0006601
with(Mes.B,plot(fitted(model),masse-fitted(model),pch=20,xlab=
"Valeurs pr\u00e9dites",main="R\u00e9sidus en fonction de la valeur pr\u00e9dite",
ylab="R\u00e9sidus = valeurs observ\u00e9es - valeurs pr\u00e9dites"))
abline(h=0)
with(Mes.B,lines(lowess(fitted(model),masse-fitted(model)),col=
"red",lwd=2))
#page 254
library(MASS)
model.r <- lqs(masse ~ taille, data=Mes.B)
summary(model.r)
## Length Class Mode
## crit 1 -none- numeric
## sing 1 -none- character
## coefficients 2 -none- numeric
## bestone 2 -none- numeric
## fitted.values 70 -none- numeric
## residuals 70 -none- numeric
## scale 2 -none- numeric
## terms 3 terms call
## call 3 -none- call
## xlevels 0 -none- list
## model 2 data.frame list
#page 255
coefficients(model.r)
## (Intercept) taille
## -7.070691 1.217949
model.r$bestone
## [1] 13 52
coef(lm(masse ~ taille, data=Mes.B[model.r$bestone,]))
## (Intercept) taille
## -6.394872 1.217949
with(Mes.B,1-sum(residuals(model.r)^2)/sum((masse-mean(masse))^2))
## [1] 0.760885
#page 256
with(Mes.B,1-sum(residuals(model)^2)/sum((masse-mean(masse))^2))
## [1] 0.8105863
plot(masse ~ taille, data=Mes.B, xlab="Taille", ylab="Masse")
abline(model.r, lty=1)
abline(model, lty=2)
with(Mes.B,points(taille[model.r$bestone],masse[model.r$bestone],
pch=19,col="red"))
abline(lm(masse ~ taille, data=Mes.B[model.r$bestone,]), lty=3)
legend("topleft", legend=c("R\u00e9sistante : moindres carr\u00e9s tronqu\u00e9s",
"Moindres carr\u00e9s ordinaires",
"Droite ajust\u00e9e entre les deux observations"), lty=1:3)
shapiro.test(residuals(model.r))
##
## Shapiro-Wilk normality test
##
## data: residuals(model.r)
## W = 0.81713, p-value = 7.025e-08
#page 257
model2<-lm(masse~taille+I(taille^2),data=Mes.B)
summary(model2)
##
## Call:
## lm(formula = masse ~ taille + I(taille^2), data = Mes.B)
##
## Residuals:
## Min 1Q Median 3Q Max
## -4.3005 -1.0546 -0.1031 1.0367 4.1643
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 5.48783 2.17254 2.526 0.01391 *
## taille -1.21792 0.36571 -3.330 0.00141 **
## I(taille^2) 0.11298 0.01476 7.656 1e-10 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.812 on 67 degrees of freedom
## Multiple R-squared: 0.899, Adjusted R-squared: 0.896
## F-statistic: 298.1 on 2 and 67 DF, p-value: < 2.2e-16
plot(masse~taille,data=Mes.B,xlab="Taille",ylab="Masse",pch=20)
with(Mes.B,lines(lowess(taille,masse),lty=3,lwd=2,col="blue"))
with(Mes.B,points(sort(taille),fitted(model2)[order(taille)]
,col="red"))
with(Mes.B,lines(sort(taille),fitted(model2)[order(taille)],lty=1,
col="red",lwd=2))
with(Mes.B,segments(taille,masse,taille,fitted(model2),lty=2,
col="red"))
legend("topleft",legend=c("Moindres carr\u00e9s ordinaires",
"Ajustement local"),lty=c(1,3),lwd=2,col=c("red","blue"))
#page 258
shapiro.test(residuals(model2))
##
## Shapiro-Wilk normality test
##
## data: residuals(model2)
## W = 0.98448, p-value = 0.5394
#page 259
confint(model2)
## 2.5 % 97.5 %
## (Intercept) 1.15143209 9.8242312
## taille -1.94787700 -0.4879672
## I(taille^2) 0.08352244 0.1424286
anova(model2)
## Analysis of Variance Table
##
## Response: masse
## Df Sum Sq Mean Sq F value Pr(>F)
## taille 1 1765.83 1765.83 537.576 < 2.2e-16 ***
## I(taille^2) 1 192.55 192.55 58.618 1.003e-10 ***
## Residuals 67 220.08 3.28
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
library(ModStatR)
my.confidence.region(model2, which=1)
my.confidence.region(model2, which=2)
my.confidence.region(model2, which=3)
#page 260
set.seed(314)
model2.r<-lqs(masse~taille+I(taille^2),data=Mes.B)
rbind(coefficients(model2),coefficients(model2.r))
## (Intercept) taille I(taille^2)
## [1,] 5.487832 -1.217922 0.1129755
## [2,] 3.849245 -1.102565 0.1162078
#page 261
shapiro.test(residuals(model2.r))
##
## Shapiro-Wilk normality test
##
## data: residuals(model2.r)
## W = 0.96603, p-value = 0.05436
plot(masse~taille,data=Mes.B,xlab="Taille",ylab="Masse",pch=20)
with(Mes.B,lines(lowess(taille, masse), lty=3, lwd=2, col="blue"))
with(Mes.B,points(sort(taille),fitted(model2.r)[order(taille)],
col="green"))
with(Mes.B,lines(sort(taille),fitted(model2.r)[order(taille)],
lty=2,col="green",lwd=2))
with(Mes.B,segments(taille,masse,taille,fitted(model2.r),lty=2,
col="green"))
legend("topleft",legend=c("R\u00e9sistante : moindres carr\u00e9s trim\u00e9s",
"Ajustement local"),lty=c(2,3),lwd=2,col=c("green","blue"))
#page 261
plot(masse~taille,data=Mes.B,xlab="Taille",ylab="Masse",pch=20)
with(Mes.B,lines(lowess(taille,masse),lty=3,lwd=2,col="blue"))
with(Mes.B,points(sort(taille),fitted(model2)[order(taille)],
col="red"))
with(Mes.B,lines(sort(taille),fitted(model2)[order(taille)],
col="red",lwd=2,lty=1))
with(Mes.B,points(sort(taille),fitted(model2.r)[order(taille)],
col="green"))
with(Mes.B,lines(sort(taille),fitted(model2.r)[order(taille)],
col="green",lwd=2,lty=2))
legend("topleft",legend=c("Moindres carr\u00e9s ordinaires",
"R\u00e9sistante : moindres carr\u00e9s trim\u00e9s","Ajustement local"),lty=1:3
,lwd=2,col=c("red","green","blue"))
#page 262
data("mtcars")
View(mtcars)
help(mtcars)
#page 263
head(mtcars,n=10)
## mpg cyl disp hp drat wt qsec vs am gear carb
## Mazda RX4 21.0 6 160.0 110 3.90 2.620 16.46 0 1 4 4
## Mazda RX4 Wag 21.0 6 160.0 110 3.90 2.875 17.02 0 1 4 4
## Datsun 710 22.8 4 108.0 93 3.85 2.320 18.61 1 1 4 1
## Hornet 4 Drive 21.4 6 258.0 110 3.08 3.215 19.44 1 0 3 1
## Hornet Sportabout 18.7 8 360.0 175 3.15 3.440 17.02 0 0 3 2
## Valiant 18.1 6 225.0 105 2.76 3.460 20.22 1 0 3 1
## Duster 360 14.3 8 360.0 245 3.21 3.570 15.84 0 0 3 4
## Merc 240D 24.4 4 146.7 62 3.69 3.190 20.00 1 0 4 2
## Merc 230 22.8 4 140.8 95 3.92 3.150 22.90 1 0 4 2
## Merc 280 19.2 6 167.6 123 3.92 3.440 18.30 1 0 4 4
?mtcars
str(mtcars)
## 'data.frame': 32 obs. of 11 variables:
## $ mpg : num 21 21 22.8 21.4 18.7 18.1 14.3 24.4 22.8 19.2 ...
## $ cyl : num 6 6 4 6 8 6 8 4 4 6 ...
## $ disp: num 160 160 108 258 360 ...
## $ hp : num 110 110 93 110 175 105 245 62 95 123 ...
## $ drat: num 3.9 3.9 3.85 3.08 3.15 2.76 3.21 3.69 3.92 3.92 ...
## $ wt : num 2.62 2.88 2.32 3.21 3.44 ...
## $ qsec: num 16.5 17 18.6 19.4 17 ...
## $ vs : num 0 0 1 1 0 1 0 1 1 1 ...
## $ am : num 1 1 1 0 0 0 0 0 0 0 ...
## $ gear: num 4 4 4 3 3 3 3 4 4 4 ...
## $ carb: num 4 4 1 1 2 1 4 2 2 4 ...
#page 264
mtcars2 <- within(mtcars, {
vs <- factor(vs, labels = c("V", "S"))
am <- factor(am, labels = c("automatic", "manual"))
cyl <- ordered(cyl)
gear <- ordered(gear)
carb <- ordered(carb)
})
summary(mtcars2)
## mpg cyl disp hp drat
## Min. :10.40 4:11 Min. : 71.1 Min. : 52.0 Min. :2.760
## 1st Qu.:15.43 6: 7 1st Qu.:120.8 1st Qu.: 96.5 1st Qu.:3.080
## Median :19.20 8:14 Median :196.3 Median :123.0 Median :3.695
## Mean :20.09 Mean :230.7 Mean :146.7 Mean :3.597
## 3rd Qu.:22.80 3rd Qu.:326.0 3rd Qu.:180.0 3rd Qu.:3.920
## Max. :33.90 Max. :472.0 Max. :335.0 Max. :4.930
## wt qsec vs am gear carb
## Min. :1.513 Min. :14.50 V:18 automatic:19 3:15 1: 7
## 1st Qu.:2.581 1st Qu.:16.89 S:14 manual :13 4:12 2:10
## Median :3.325 Median :17.71 5: 5 3: 3
## Mean :3.217 Mean :17.85 4:10
## 3rd Qu.:3.610 3rd Qu.:18.90 6: 1
## Max. :5.424 Max. :22.90 8: 1
subsetmtcars<-mtcars[,c(1,3,4,5,6,7)]
#page 265
library(MVN)
## Registered S3 method overwritten by 'GGally':
## method from
## +.gg ggplot2
## sROC 0.1-2 loaded
set.seed(1133)
result1 = mvn(data = subsetmtcars,
mvnTest = "mardia",
univariateTest = "SW", univariatePlot = "histogram",
multivariatePlot = "qq",
multivariateOutlierMethod = "adj",
showOutliers = TRUE, showNewData = TRUE)
result1$multivariateNormality
## Test Statistic p value Result
## 1 Mardia Skewness 77.6519454787057 0.0293424460277968 NO
## 2 Mardia Kurtosis 0.241701622075929 0.80901137258975 YES
## 3 MVN <NA> <NA> NO
result1$multivariateOutliers
## Observation Mahalanobis Distance Outlier
## Maserati Bora Maserati Bora 107.537 TRUE
## Ford Pantera L Ford Pantera L 85.115 TRUE
## Camaro Z28 Camaro Z28 35.230 TRUE
## Duster 360 Duster 360 30.679 TRUE
## Merc 240D Merc 240D 20.804 TRUE
## Merc 230 Merc 230 16.813 TRUE
## Fiat 128 Fiat 128 16.123 TRUE
#page 266
library(corrplot); set.seed(1133)
## corrplot 0.84 loaded
permmtcars <- perm.cor.mtest(subsetmtcars)
permmtcars$p<.05/choose(ncol(subsetmtcars),2)
## [,1] [,2] [,3] [,4] [,5] [,6]
## [1,] TRUE TRUE TRUE TRUE TRUE FALSE
## [2,] TRUE TRUE TRUE TRUE TRUE FALSE
## [3,] TRUE TRUE TRUE FALSE TRUE TRUE
## [4,] TRUE TRUE FALSE TRUE TRUE FALSE
## [5,] TRUE TRUE TRUE TRUE TRUE FALSE
## [6,] FALSE FALSE TRUE FALSE FALSE TRUE
corrplot(permmtcars$cor,p.mat=permmtcars$p,pch.col="white",insig=
"label_sig",sig.level=.05/choose(ncol(subsetmtcars),2))
#page 267
fit.mtcars<-lm(mpg~.,data=subsetmtcars)
fit.mtcars
##
## Call:
## lm(formula = mpg ~ ., data = subsetmtcars)
##
## Coefficients:
## (Intercept) disp hp drat wt
## 16.53357 0.00872 -0.02060 2.01577 -4.38546
## qsec
## 0.64015
shapiro.test(residuals(fit.mtcars))
##
## Shapiro-Wilk normality test
##
## data: residuals(fit.mtcars)
## W = 0.92929, p-value = 0.03743
fit2.mtcars<-lm(mpg~.+disp:hp,data=subsetmtcars)
round(coef(fit2.mtcars), 3)
## (Intercept) disp hp drat wt qsec
## 48.542 -0.052 -0.101 -1.241 -3.826 0.156
## disp:hp
## 0.000
shapiro.test(residuals(fit2.mtcars))
##
## Shapiro-Wilk normality test
##
## data: residuals(fit2.mtcars)
## W = 0.96295, p-value = 0.33
covratio(fit2.mtcars)
## Mazda RX4 Mazda RX4 Wag Datsun 710
## 1.4219125 1.5291914 0.7999250
## Hornet 4 Drive Hornet Sportabout Valiant
## 1.3236334 1.3714173 1.3215512
## Duster 360 Merc 240D Merc 230
## 1.0349281 1.5335287 3.2039954
## Merc 280 Merc 280C Merc 450SE
## 1.6039647 1.5770987 1.2297063
## Merc 450SL Merc 450SLC Cadillac Fleetwood
## 1.2891756 1.4466979 1.4781420
## Lincoln Continental Chrysler Imperial Fiat 128
## 1.7258001 0.8602740 0.6893901
## Honda Civic Toyota Corolla Toyota Corona
## 1.8711733 0.7327578 0.6084108
## Dodge Challenger AMC Javelin Camaro Z28
## 1.3597375 1.1814885 1.3136352
## Pontiac Firebird Fiat X1-9 Porsche 914-2
## 0.6392936 1.0386684 1.6124398
## Lotus Europa Ford Pantera L Ferrari Dino
## 0.9815851 2.0231153 1.6305812
## Maserati Bora Volvo 142E
## 2.6989138 1.3923920
dffits(fit2.mtcars)
## Mazda RX4 Mazda RX4 Wag Datsun 710
## -0.24499825 -0.06142428 -0.43740362
## Hornet 4 Drive Hornet Sportabout Valiant
## 0.42129269 0.43178720 -0.50472680
## Duster 360 Merc 240D Merc 230
## -0.61997624 0.12104459 0.33521636
## Merc 280 Merc 280C Merc 450SE
## 0.15377223 -0.20569571 0.35176504
## Merc 450SL Merc 450SLC Cadillac Fleetwood
## 0.20551285 -0.08194123 -0.45924043
## Lincoln Continental Chrysler Imperial Fiat 128
## -0.35087731 1.14031935 0.83684620
## Honda Civic Toyota Corolla Toyota Corona
## -0.28216079 0.88220835 -0.69922754
## Dodge Challenger AMC Javelin Camaro Z28
## -0.42399835 -0.36768988 -0.30281726
## Pontiac Firebird Fiat X1-9 Porsche 914-2
## 1.01527533 -0.64558078 0.14852889
## Lotus Europa Ford Pantera L Ferrari Dino
## 0.77778997 -0.48300413 0.41626973
## Maserati Bora Volvo 142E
## 0.69740095 -0.17128938
dfbetas(fit2.mtcars)
## (Intercept) disp hp drat
## Mazda RX4 -0.086088964 -0.008702825 0.037413221 -0.02182775
## Mazda RX4 Wag -0.013905727 -0.001244996 0.004697082 -0.01138912
## Datsun 710 -0.199000746 0.213768560 0.161143042 0.20081906
## Hornet 4 Drive -0.147833016 0.258699076 0.172982801 0.02940048
## Hornet Sportabout -0.110321619 0.289163636 0.125480071 0.08556140
## Valiant -0.100639457 0.076873454 0.026170932 0.30250979
## Duster 360 -0.192659793 0.169152653 0.130640793 0.26778134
## Merc 240D 0.041392073 -0.056394703 -0.064825546 -0.03027678
## Merc 230 -0.270944867 0.136040401 0.200444094 0.15199460
## Merc 280 -0.052943995 0.032781470 0.057422653 0.07861283
## Merc 280C 0.112678571 -0.068519304 -0.108977999 -0.12312406
## Merc 450SE 0.086316949 -0.084937691 0.022680395 -0.10647032
## Merc 450SL 0.011437015 0.003299762 0.064900462 -0.06344620
## Merc 450SLC 0.010068878 -0.007832338 -0.035470818 0.01731092
## Cadillac Fleetwood -0.050955216 0.074799187 0.197571813 0.03020046
## Lincoln Continental -0.075817382 0.132907459 0.177850197 0.04794989
## Chrysler Imperial 0.208147068 -0.505859889 -0.534896799 -0.07738021
## Fiat 128 0.409561167 -0.569954842 -0.553165508 -0.36941505
## Honda Civic 0.040714693 -0.036475491 0.057714677 -0.13168997
## Toyota Corolla 0.127010905 -0.314149654 -0.332797689 -0.19440346
## Toyota Corona 0.118445210 0.053317944 -0.114126131 0.14605461
## Dodge Challenger -0.208970200 -0.020400961 0.087875946 0.21152215
## AMC Javelin 0.047753232 -0.232191461 -0.114468764 -0.04700651
## Camaro Z28 0.060958292 -0.043185128 -0.051390620 -0.12193544
## Pontiac Firebird -0.211668988 0.639113539 0.155131235 0.23752923
## Fiat X1-9 -0.376862918 0.406080280 0.450567492 0.31954438
## Porsche 914-2 0.006220057 0.032183416 -0.010442674 0.06683811
## Lotus Europa 0.386612706 -0.176780826 -0.144536638 -0.40222001
## Ford Pantera L 0.219251092 -0.224726647 -0.178023752 -0.29473442
## Ferrari Dino 0.101253314 -0.042161872 0.102496298 -0.05311804
## Maserati Bora 0.024161864 -0.232079735 0.131016572 -0.14427475
## Volvo 142E 0.068270224 -0.030326014 -0.065962875 -0.08792691
## wt qsec disp:hp
## Mazda RX4 -0.0796827179 0.174052558 0.033599180
## Mazda RX4 Wag -0.0322453171 0.039296803 0.012111456
## Datsun 710 0.0051321367 0.070643004 -0.186838984
## Hornet 4 Drive -0.2276727214 0.205746109 -0.209213731
## Hornet Sportabout -0.2952023145 0.112955401 -0.174087055
## Valiant 0.0604536870 -0.112223907 -0.051120210
## Duster 360 0.2720698950 -0.017987111 -0.273636379
## Merc 240D 0.0670394034 -0.032995419 0.047302428
## Merc 230 -0.0771680415 0.302027348 -0.141540078
## Merc 280 0.0968644547 -0.014904655 -0.075741481
## Merc 280C -0.0992211299 -0.034097436 0.120668871
## Merc 450SE 0.2118821070 -0.113911976 -0.024114438
## Merc 450SL 0.0230840124 0.006973465 -0.047258618
## Merc 450SLC -0.0042643070 -0.019571187 0.024649961
## Cadillac Fleetwood -0.1154118693 0.038514240 -0.158237807
## Lincoln Continental -0.1689816365 0.073284335 -0.162627257
## Chrysler Imperial 0.6386748656 -0.251218092 0.547582294
## Fiat 128 0.0765837714 -0.131784778 0.605010812
## Honda Civic 0.0569160272 -0.004434124 -0.025174681
## Toyota Corolla -0.3027869038 0.222517856 0.453658820
## Toyota Corona 0.2607683913 -0.414392281 -0.041631752
## Dodge Challenger 0.0763635356 0.158292138 -0.008322643
## AMC Javelin 0.1327519934 -0.006044884 0.193724740
## Camaro Z28 -0.0003569387 0.014465490 0.024557457
## Pontiac Firebird -0.5211737987 0.144242532 -0.347249954
## Fiat X1-9 0.0476906991 0.152970258 -0.471354122
## Porsche 914-2 0.0174679399 -0.066414748 -0.028382664
## Lotus Europa -0.3741524726 -0.145848581 0.204137330
## Ford Pantera L 0.2041053791 -0.113085953 0.139074107
## Ferrari Dino 0.1513489108 -0.188847641 -0.100036080
## Maserati Bora -0.0367873163 0.094199143 0.137913963
## Volvo 142E -0.0639020279 -0.016923245 0.069537042
library(car)
## Loading required package: carData
vif(fit2.mtcars)
## disp hp drat wt qsec disp:hp
## 48.793223 26.669457 4.462626 7.224276 3.719492 76.606107
#page 268
library(perturb)
colldiag(fit2.mtcars)
## Condition
## Index Variance Decomposition Proportions
## intercept disp hp drat wt qsec
## 1 1.000 0.000 0.001 0.001 0.000 0.000 0.000
## 2 4.393 0.000 0.027 0.021 0.007 0.001 0.002
## 3 9.946 0.000 0.053 0.330 0.012 0.037 0.003
## 4 21.905 0.001 0.772 0.182 0.102 0.392 0.007
## 5 32.009 0.032 0.012 0.142 0.606 0.449 0.159
## 6 71.166 0.966 0.136 0.324 0.272 0.120 0.829
plot(fit2.mtcars)
influence.measures(fit2.mtcars)
## Influence measures of
## lm(formula = mpg ~ . + disp:hp, data = subsetmtcars) :
##
## dfb.1_ dfb.disp dfb.hp dfb.drat dfb.wt dfb.qsec
## Mazda RX4 -0.08609 -0.00870 0.0374 -0.0218 -0.079683 0.17405
## Mazda RX4 Wag -0.01391 -0.00124 0.0047 -0.0114 -0.032245 0.03930
## Datsun 710 -0.19900 0.21377 0.1611 0.2008 0.005132 0.07064
## Hornet 4 Drive -0.14783 0.25870 0.1730 0.0294 -0.227673 0.20575
## Hornet Sportabout -0.11032 0.28916 0.1255 0.0856 -0.295202 0.11296
## Valiant -0.10064 0.07687 0.0262 0.3025 0.060454 -0.11222
## Duster 360 -0.19266 0.16915 0.1306 0.2678 0.272070 -0.01799
## Merc 240D 0.04139 -0.05639 -0.0648 -0.0303 0.067039 -0.03300
## Merc 230 -0.27094 0.13604 0.2004 0.1520 -0.077168 0.30203
## Merc 280 -0.05294 0.03278 0.0574 0.0786 0.096864 -0.01490
## Merc 280C 0.11268 -0.06852 -0.1090 -0.1231 -0.099221 -0.03410
## Merc 450SE 0.08632 -0.08494 0.0227 -0.1065 0.211882 -0.11391
## Merc 450SL 0.01144 0.00330 0.0649 -0.0634 0.023084 0.00697
## Merc 450SLC 0.01007 -0.00783 -0.0355 0.0173 -0.004264 -0.01957
## Cadillac Fleetwood -0.05096 0.07480 0.1976 0.0302 -0.115412 0.03851
## Lincoln Continental -0.07582 0.13291 0.1779 0.0479 -0.168982 0.07328
## Chrysler Imperial 0.20815 -0.50586 -0.5349 -0.0774 0.638675 -0.25122
## Fiat 128 0.40956 -0.56995 -0.5532 -0.3694 0.076584 -0.13178
## Honda Civic 0.04071 -0.03648 0.0577 -0.1317 0.056916 -0.00443
## Toyota Corolla 0.12701 -0.31415 -0.3328 -0.1944 -0.302787 0.22252
## Toyota Corona 0.11845 0.05332 -0.1141 0.1461 0.260768 -0.41439
## Dodge Challenger -0.20897 -0.02040 0.0879 0.2115 0.076364 0.15829
## AMC Javelin 0.04775 -0.23219 -0.1145 -0.0470 0.132752 -0.00604
## Camaro Z28 0.06096 -0.04319 -0.0514 -0.1219 -0.000357 0.01447
## Pontiac Firebird -0.21167 0.63911 0.1551 0.2375 -0.521174 0.14424
## Fiat X1-9 -0.37686 0.40608 0.4506 0.3195 0.047691 0.15297
## Porsche 914-2 0.00622 0.03218 -0.0104 0.0668 0.017468 -0.06641
## Lotus Europa 0.38661 -0.17678 -0.1445 -0.4022 -0.374152 -0.14585
## Ford Pantera L 0.21925 -0.22473 -0.1780 -0.2947 0.204105 -0.11309
## Ferrari Dino 0.10125 -0.04216 0.1025 -0.0531 0.151349 -0.18885
## Maserati Bora 0.02416 -0.23208 0.1310 -0.1443 -0.036787 0.09420
## Volvo 142E 0.06827 -0.03033 -0.0660 -0.0879 -0.063902 -0.01692
## dfb.dsp. dffit cov.r cook.d hat inf
## Mazda RX4 0.03360 -0.2450 1.422 0.008808 0.1512
## Mazda RX4 Wag 0.01211 -0.0614 1.529 0.000561 0.1358
## Datsun 710 -0.18684 -0.4374 0.800 0.026150 0.0824
## Hornet 4 Drive -0.20921 0.4213 1.324 0.025607 0.1905
## Hornet Sportabout -0.17409 0.4318 1.371 0.026949 0.2084
## Valiant -0.05112 -0.5047 1.322 0.036543 0.2212
## Duster 360 -0.27364 -0.6200 1.035 0.053543 0.1900
## Merc 240D 0.04730 0.1210 1.534 0.002173 0.1525
## Merc 230 -0.14154 0.3352 3.204 0.016668 0.5938 *
## Merc 280 -0.07574 0.1538 1.604 0.003504 0.1938
## Merc 280C 0.12067 -0.2057 1.577 0.006251 0.1974
## Merc 450SE -0.02411 0.3518 1.230 0.017830 0.1361
## Merc 450SL -0.04726 0.2055 1.289 0.006174 0.0890
## Merc 450SLC 0.02465 -0.0819 1.447 0.000997 0.0968
## Cadillac Fleetwood -0.15824 -0.4592 1.478 0.030577 0.2498
## Lincoln Continental -0.16263 -0.3509 1.726 0.018096 0.2926
## Chrysler Imperial 0.54758 1.1403 0.860 0.172497 0.3079
## Fiat 128 0.60501 0.8368 0.689 0.092174 0.1826
## Honda Civic -0.02517 -0.2822 1.871 0.011765 0.3225 *
## Toyota Corolla 0.45366 0.8822 0.733 0.102915 0.2055
## Toyota Corona -0.04163 -0.6992 0.608 0.063810 0.1269
## Dodge Challenger -0.00832 -0.4240 1.360 0.025983 0.2020
## AMC Javelin 0.19372 -0.3677 1.181 0.019389 0.1302
## Camaro Z28 0.02456 -0.3028 1.314 0.013331 0.1394
## Pontiac Firebird -0.34725 1.0153 0.639 0.133282 0.2216
## Fiat X1-9 -0.47135 -0.6456 1.039 0.057987 0.1998
## Porsche 914-2 -0.02838 0.1485 1.612 0.003271 0.1961
## Lotus Europa 0.20414 0.7778 0.982 0.083029 0.2303
## Ford Pantera L 0.13907 -0.4830 2.023 0.034226 0.4045 *
## Ferrari Dino -0.10004 0.4163 1.631 0.025320 0.2817
## Maserati Bora 0.13791 0.6974 2.699 0.071236 0.5588 *
## Volvo 142E 0.06954 -0.1713 1.392 0.004323 0.1087
influencePlot(fit2.mtcars)
## StudRes Hat CookD
## Merc 230 0.2772254 0.5938467 0.01666835
## Chrysler Imperial 1.7095034 0.3079354 0.17249684
## Toyota Corona -1.8342492 0.1268803 0.06381047
## Pontiac Firebird 1.9028828 0.2215907 0.13328179
## Maserati Bora 0.6197073 0.5587835 0.07123631
influenceIndexPlot(fit2.mtcars)
summary(fit2.mtcars)
##
## Call:
## lm(formula = mpg ~ . + disp:hp, data = subsetmtcars)
##
## Residuals:
## Min 1Q Median 3Q Max
## -3.6810 -1.5658 -0.4375 1.4929 3.5889
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 48.5415672 14.5103848 3.345 0.00260 **
## disp -0.0517714 0.0227449 -2.276 0.03166 *
## hp -0.1010258 0.0303969 -3.324 0.00274 **
## drat -1.2406837 1.5944591 -0.778 0.44380
## wt -3.8259532 1.1085781 -3.451 0.00199 **
## qsec 0.1563953 0.4355541 0.359 0.72256
## disp:hp 0.0003046 0.0001033 2.949 0.00682 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 2.247 on 25 degrees of freedom
## Multiple R-squared: 0.8879, Adjusted R-squared: 0.861
## F-statistic: 33.01 on 6 and 25 DF, p-value: 1.035e-10
anova(fit2.mtcars)
## Analysis of Variance Table
##
## Response: mpg
## Df Sum Sq Mean Sq F value Pr(>F)
## disp 1 808.89 808.89 160.2156 2.269e-12 ***
## hp 1 33.67 33.67 6.6680 0.0160616 *
## drat 1 30.15 30.15 5.9713 0.0219401 *
## wt 1 70.51 70.51 13.9655 0.0009702 ***
## qsec 1 12.71 12.71 2.5171 0.1251833
## disp:hp 1 43.91 43.91 8.6973 0.0068211 **
## Residuals 25 126.22 5.05
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#Si le package n'est pas installe, enlever le commentaire
#puis executer la commande ci-dessous.
#if(!require("ISLR")){install.packages("ISLR")}
library(ISLR)
data(Hitters)
#Si le package n'est pas installe, enlever le commentaire
#puis executer la commande ci-dessous.
#if(!require("mice")){install.packages("mice")}
library(mice)
## Loading required package: lattice
## Registered S3 methods overwritten by 'lme4':
## method from
## cooks.distance.influence.merMod car
## influence.merMod car
## dfbeta.influence.merMod car
## dfbetas.influence.merMod car
##
## Attaching package: 'mice'
## The following objects are masked from 'package:base':
##
## cbind, rbind
md.pattern(Hitters)
## AtBat Hits HmRun Runs RBI Walks Years CAtBat CHits CHmRun CRuns CRBI
## 263 1 1 1 1 1 1 1 1 1 1 1 1
## 59 1 1 1 1 1 1 1 1 1 1 1 1
## 0 0 0 0 0 0 0 0 0 0 0 0
## CWalks League Division PutOuts Assists Errors NewLeague Salary
## 263 1 1 1 1 1 1 1 1 0
## 59 1 1 1 1 1 1 1 0 1
## 0 0 0 0 0 0 0 59 59
#page 269
md.pairs(Hitters)
## $rr
## AtBat Hits HmRun Runs RBI Walks Years CAtBat CHits CHmRun CRuns
## AtBat 322 322 322 322 322 322 322 322 322 322 322
## Hits 322 322 322 322 322 322 322 322 322 322 322
## HmRun 322 322 322 322 322 322 322 322 322 322 322
## Runs 322 322 322 322 322 322 322 322 322 322 322
## RBI 322 322 322 322 322 322 322 322 322 322 322
## Walks 322 322 322 322 322 322 322 322 322 322 322
## Years 322 322 322 322 322 322 322 322 322 322 322
## CAtBat 322 322 322 322 322 322 322 322 322 322 322
## CHits 322 322 322 322 322 322 322 322 322 322 322
## CHmRun 322 322 322 322 322 322 322 322 322 322 322
## CRuns 322 322 322 322 322 322 322 322 322 322 322
## CRBI 322 322 322 322 322 322 322 322 322 322 322
## CWalks 322 322 322 322 322 322 322 322 322 322 322
## League 322 322 322 322 322 322 322 322 322 322 322
## Division 322 322 322 322 322 322 322 322 322 322 322
## PutOuts 322 322 322 322 322 322 322 322 322 322 322
## Assists 322 322 322 322 322 322 322 322 322 322 322
## Errors 322 322 322 322 322 322 322 322 322 322 322
## Salary 263 263 263 263 263 263 263 263 263 263 263
## NewLeague 322 322 322 322 322 322 322 322 322 322 322
## CRBI CWalks League Division PutOuts Assists Errors Salary
## AtBat 322 322 322 322 322 322 322 263
## Hits 322 322 322 322 322 322 322 263
## HmRun 322 322 322 322 322 322 322 263
## Runs 322 322 322 322 322 322 322 263
## RBI 322 322 322 322 322 322 322 263
## Walks 322 322 322 322 322 322 322 263
## Years 322 322 322 322 322 322 322 263
## CAtBat 322 322 322 322 322 322 322 263
## CHits 322 322 322 322 322 322 322 263
## CHmRun 322 322 322 322 322 322 322 263
## CRuns 322 322 322 322 322 322 322 263
## CRBI 322 322 322 322 322 322 322 263
## CWalks 322 322 322 322 322 322 322 263
## League 322 322 322 322 322 322 322 263
## Division 322 322 322 322 322 322 322 263
## PutOuts 322 322 322 322 322 322 322 263
## Assists 322 322 322 322 322 322 322 263
## Errors 322 322 322 322 322 322 322 263
## Salary 263 263 263 263 263 263 263 263
## NewLeague 322 322 322 322 322 322 322 263
## NewLeague
## AtBat 322
## Hits 322
## HmRun 322
## Runs 322
## RBI 322
## Walks 322
## Years 322
## CAtBat 322
## CHits 322
## CHmRun 322
## CRuns 322
## CRBI 322
## CWalks 322
## League 322
## Division 322
## PutOuts 322
## Assists 322
## Errors 322
## Salary 263
## NewLeague 322
##
## $rm
## AtBat Hits HmRun Runs RBI Walks Years CAtBat CHits CHmRun CRuns
## AtBat 0 0 0 0 0 0 0 0 0 0 0
## Hits 0 0 0 0 0 0 0 0 0 0 0
## HmRun 0 0 0 0 0 0 0 0 0 0 0
## Runs 0 0 0 0 0 0 0 0 0 0 0
## RBI 0 0 0 0 0 0 0 0 0 0 0
## Walks 0 0 0 0 0 0 0 0 0 0 0
## Years 0 0 0 0 0 0 0 0 0 0 0
## CAtBat 0 0 0 0 0 0 0 0 0 0 0
## CHits 0 0 0 0 0 0 0 0 0 0 0
## CHmRun 0 0 0 0 0 0 0 0 0 0 0
## CRuns 0 0 0 0 0 0 0 0 0 0 0
## CRBI 0 0 0 0 0 0 0 0 0 0 0
## CWalks 0 0 0 0 0 0 0 0 0 0 0
## League 0 0 0 0 0 0 0 0 0 0 0
## Division 0 0 0 0 0 0 0 0 0 0 0
## PutOuts 0 0 0 0 0 0 0 0 0 0 0
## Assists 0 0 0 0 0 0 0 0 0 0 0
## Errors 0 0 0 0 0 0 0 0 0 0 0
## Salary 0 0 0 0 0 0 0 0 0 0 0
## NewLeague 0 0 0 0 0 0 0 0 0 0 0
## CRBI CWalks League Division PutOuts Assists Errors Salary
## AtBat 0 0 0 0 0 0 0 59
## Hits 0 0 0 0 0 0 0 59
## HmRun 0 0 0 0 0 0 0 59
## Runs 0 0 0 0 0 0 0 59
## RBI 0 0 0 0 0 0 0 59
## Walks 0 0 0 0 0 0 0 59
## Years 0 0 0 0 0 0 0 59
## CAtBat 0 0 0 0 0 0 0 59
## CHits 0 0 0 0 0 0 0 59
## CHmRun 0 0 0 0 0 0 0 59
## CRuns 0 0 0 0 0 0 0 59
## CRBI 0 0 0 0 0 0 0 59
## CWalks 0 0 0 0 0 0 0 59
## League 0 0 0 0 0 0 0 59
## Division 0 0 0 0 0 0 0 59
## PutOuts 0 0 0 0 0 0 0 59
## Assists 0 0 0 0 0 0 0 59
## Errors 0 0 0 0 0 0 0 59
## Salary 0 0 0 0 0 0 0 0
## NewLeague 0 0 0 0 0 0 0 59
## NewLeague
## AtBat 0
## Hits 0
## HmRun 0
## Runs 0
## RBI 0
## Walks 0
## Years 0
## CAtBat 0
## CHits 0
## CHmRun 0
## CRuns 0
## CRBI 0
## CWalks 0
## League 0
## Division 0
## PutOuts 0
## Assists 0
## Errors 0
## Salary 0
## NewLeague 0
##
## $mr
## AtBat Hits HmRun Runs RBI Walks Years CAtBat CHits CHmRun CRuns
## AtBat 0 0 0 0 0 0 0 0 0 0 0
## Hits 0 0 0 0 0 0 0 0 0 0 0
## HmRun 0 0 0 0 0 0 0 0 0 0 0
## Runs 0 0 0 0 0 0 0 0 0 0 0
## RBI 0 0 0 0 0 0 0 0 0 0 0
## Walks 0 0 0 0 0 0 0 0 0 0 0
## Years 0 0 0 0 0 0 0 0 0 0 0
## CAtBat 0 0 0 0 0 0 0 0 0 0 0
## CHits 0 0 0 0 0 0 0 0 0 0 0
## CHmRun 0 0 0 0 0 0 0 0 0 0 0
## CRuns 0 0 0 0 0 0 0 0 0 0 0
## CRBI 0 0 0 0 0 0 0 0 0 0 0
## CWalks 0 0 0 0 0 0 0 0 0 0 0
## League 0 0 0 0 0 0 0 0 0 0 0
## Division 0 0 0 0 0 0 0 0 0 0 0
## PutOuts 0 0 0 0 0 0 0 0 0 0 0
## Assists 0 0 0 0 0 0 0 0 0 0 0
## Errors 0 0 0 0 0 0 0 0 0 0 0
## Salary 59 59 59 59 59 59 59 59 59 59 59
## NewLeague 0 0 0 0 0 0 0 0 0 0 0
## CRBI CWalks League Division PutOuts Assists Errors Salary
## AtBat 0 0 0 0 0 0 0 0
## Hits 0 0 0 0 0 0 0 0
## HmRun 0 0 0 0 0 0 0 0
## Runs 0 0 0 0 0 0 0 0
## RBI 0 0 0 0 0 0 0 0
## Walks 0 0 0 0 0 0 0 0
## Years 0 0 0 0 0 0 0 0
## CAtBat 0 0 0 0 0 0 0 0
## CHits 0 0 0 0 0 0 0 0
## CHmRun 0 0 0 0 0 0 0 0
## CRuns 0 0 0 0 0 0 0 0
## CRBI 0 0 0 0 0 0 0 0
## CWalks 0 0 0 0 0 0 0 0
## League 0 0 0 0 0 0 0 0
## Division 0 0 0 0 0 0 0 0
## PutOuts 0 0 0 0 0 0 0 0
## Assists 0 0 0 0 0 0 0 0
## Errors 0 0 0 0 0 0 0 0
## Salary 59 59 59 59 59 59 59 0
## NewLeague 0 0 0 0 0 0 0 0
## NewLeague
## AtBat 0
## Hits 0
## HmRun 0
## Runs 0
## RBI 0
## Walks 0
## Years 0
## CAtBat 0
## CHits 0
## CHmRun 0
## CRuns 0
## CRBI 0
## CWalks 0
## League 0
## Division 0
## PutOuts 0
## Assists 0
## Errors 0
## Salary 59
## NewLeague 0
##
## $mm
## AtBat Hits HmRun Runs RBI Walks Years CAtBat CHits CHmRun CRuns
## AtBat 0 0 0 0 0 0 0 0 0 0 0
## Hits 0 0 0 0 0 0 0 0 0 0 0
## HmRun 0 0 0 0 0 0 0 0 0 0 0
## Runs 0 0 0 0 0 0 0 0 0 0 0
## RBI 0 0 0 0 0 0 0 0 0 0 0
## Walks 0 0 0 0 0 0 0 0 0 0 0
## Years 0 0 0 0 0 0 0 0 0 0 0
## CAtBat 0 0 0 0 0 0 0 0 0 0 0
## CHits 0 0 0 0 0 0 0 0 0 0 0
## CHmRun 0 0 0 0 0 0 0 0 0 0 0
## CRuns 0 0 0 0 0 0 0 0 0 0 0
## CRBI 0 0 0 0 0 0 0 0 0 0 0
## CWalks 0 0 0 0 0 0 0 0 0 0 0
## League 0 0 0 0 0 0 0 0 0 0 0
## Division 0 0 0 0 0 0 0 0 0 0 0
## PutOuts 0 0 0 0 0 0 0 0 0 0 0
## Assists 0 0 0 0 0 0 0 0 0 0 0
## Errors 0 0 0 0 0 0 0 0 0 0 0
## Salary 0 0 0 0 0 0 0 0 0 0 0
## NewLeague 0 0 0 0 0 0 0 0 0 0 0
## CRBI CWalks League Division PutOuts Assists Errors Salary
## AtBat 0 0 0 0 0 0 0 0
## Hits 0 0 0 0 0 0 0 0
## HmRun 0 0 0 0 0 0 0 0
## Runs 0 0 0 0 0 0 0 0
## RBI 0 0 0 0 0 0 0 0
## Walks 0 0 0 0 0 0 0 0
## Years 0 0 0 0 0 0 0 0
## CAtBat 0 0 0 0 0 0 0 0
## CHits 0 0 0 0 0 0 0 0
## CHmRun 0 0 0 0 0 0 0 0
## CRuns 0 0 0 0 0 0 0 0
## CRBI 0 0 0 0 0 0 0 0
## CWalks 0 0 0 0 0 0 0 0
## League 0 0 0 0 0 0 0 0
## Division 0 0 0 0 0 0 0 0
## PutOuts 0 0 0 0 0 0 0 0
## Assists 0 0 0 0 0 0 0 0
## Errors 0 0 0 0 0 0 0 0
## Salary 0 0 0 0 0 0 0 59
## NewLeague 0 0 0 0 0 0 0 0
## NewLeague
## AtBat 0
## Hits 0
## HmRun 0
## Runs 0
## RBI 0
## Walks 0
## Years 0
## CAtBat 0
## CHits 0
## CHmRun 0
## CRuns 0
## CRBI 0
## CWalks 0
## League 0
## Division 0
## PutOuts 0
## Assists 0
## Errors 0
## Salary 0
## NewLeague 0
library(dplyr)
##
## Attaching package: 'dplyr'
## The following object is masked from 'package:car':
##
## recode
## The following object is masked from 'package:MASS':
##
## select
## The following objects are masked from 'package:stats':
##
## filter, lag
## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, union
#Si le package n'est pas installe, enlever le commentaire
#puis executer la commande ci-dessous.
#if(!require("finalfit")){install.packages("finalfit")}
library(finalfit)
Hitters %>%
missing_plot()
explanatory = setdiff(colnames(Hitters),"Salary")
dependent = "Salary"
Hitters %>%
missing_pattern(dependent, explanatory)
## AtBat Hits HmRun Runs RBI Walks Years CAtBat CHits CHmRun CRuns CRBI
## 263 1 1 1 1 1 1 1 1 1 1 1 1
## 59 1 1 1 1 1 1 1 1 1 1 1 1
## 0 0 0 0 0 0 0 0 0 0 0 0
## CWalks League Division PutOuts Assists Errors NewLeague Salary
## 263 1 1 1 1 1 1 1 1 0
## 59 1 1 1 1 1 1 1 0 1
## 0 0 0 0 0 0 0 59 59
#page 270
#Si le package n'est pas installe, enlever le commentaire
#puis executer la commande ci-dessous.
#if(!require("naniar")){install.packages("naniar")}
library(naniar);library(ggplot2)
Hitters %>%
bind_shadow() %>%
ggplot(aes(x = Hits,
fill = Salary_NA)) +
geom_density(alpha = 0.5)
#page 271
gg_miss_var(Hitters)
try(gg_miss_upset(Hitters))
## Error : upset plots for missing data requre at least two variables to have missing data, only one variable, 'Salary' has missing values.
## Backtrace:
## █
## 1. ├─rmarkdown::render(...)
## 2. │ └─knitr::knit(...)
## 3. │ └─knitr:::process_file(text, output)
## 4. │ ├─base::withCallingHandlers(...)
## 5. │ ├─knitr:::process_group(group)
## 6. │ └─knitr:::process_group.block(group)
## 7. │ └─knitr:::call_block(x)
## 8. │ └─knitr:::block_exec(params)
## 9. │ ├─knitr:::in_dir(...)
## 10. │ └─knitr:::evaluate(...)
## 11. │ └─evaluate::evaluate(...)
## 12. │ └─evaluate:::evaluate_call(...)
## 13. │ ├─evaluate:::timing_fn(...)
## 14. │ ├─base:::handle(...)
## 15. │ ├─base::withCallingHandlers(...)
## 16. │ ├─base::withVisible(eval(expr, envir, enclos))
## 17. │ └─base::eval(expr, envir, enclos)
## 18. │ └─base::eval(expr, envir, enclos)
## 19. ├─base::try(gg_miss_upset(Hitters))
## 20. │ └─base::tryCatch(...)
## 21. │ └─base:::tryCatchList(expr, classes, parentenv, handlers)
## 22. │ └─base:::tryCatchOne(expr, names, parentenv, handlers[[1L]])
## 23. │ └─base:::doTryCatch(return(expr), name, parentenv, handler)
## 24. └─naniar::gg_miss_upset(Hitters)
#page 272
#Exercice 4.1
data(anscombe)
str(anscombe)
## 'data.frame': 11 obs. of 8 variables:
## $ x1: num 10 8 13 9 11 14 6 4 12 7 ...
## $ x2: num 10 8 13 9 11 14 6 4 12 7 ...
## $ x3: num 10 8 13 9 11 14 6 4 12 7 ...
## $ x4: num 8 8 8 8 8 8 8 19 8 8 ...
## $ y1: num 8.04 6.95 7.58 8.81 8.33 ...
## $ y2: num 9.14 8.14 8.74 8.77 9.26 8.1 6.13 3.1 9.13 7.26 ...
## $ y3: num 7.46 6.77 12.74 7.11 7.81 ...
## $ y4: num 6.58 5.76 7.71 8.84 8.47 7.04 5.25 12.5 5.56 7.91 ...
#page 273
#Exercice 4.2
Hitters = na.omit(Hitters)
#q4
#Si le package n'est pas installe, enlever le commentaire
#puis executer la commande ci-dessous.
#if(!require(leaps)){install.packages("leaps")}
library(leaps)
regfit.full = regsubsets(Salary ~ ., data = Hitters)
summary(regfit.full)
## Subset selection object
## Call: regsubsets.formula(Salary ~ ., data = Hitters)
## 19 Variables (and intercept)
## Forced in Forced out
## AtBat FALSE FALSE
## Hits FALSE FALSE
## HmRun FALSE FALSE
## Runs FALSE FALSE
## RBI FALSE FALSE
## Walks FALSE FALSE
## Years FALSE FALSE
## CAtBat FALSE FALSE
## CHits FALSE FALSE
## CHmRun FALSE FALSE
## CRuns FALSE FALSE
## CRBI FALSE FALSE
## CWalks FALSE FALSE
## LeagueN FALSE FALSE
## DivisionW FALSE FALSE
## PutOuts FALSE FALSE
## Assists FALSE FALSE
## Errors FALSE FALSE
## NewLeagueN FALSE FALSE
## 1 subsets of each size up to 8
## Selection Algorithm: exhaustive
## AtBat Hits HmRun Runs RBI Walks Years CAtBat CHits CHmRun CRuns
## 1 ( 1 ) " " " " " " " " " " " " " " " " " " " " " "
## 2 ( 1 ) " " "*" " " " " " " " " " " " " " " " " " "
## 3 ( 1 ) " " "*" " " " " " " " " " " " " " " " " " "
## 4 ( 1 ) " " "*" " " " " " " " " " " " " " " " " " "
## 5 ( 1 ) "*" "*" " " " " " " " " " " " " " " " " " "
## 6 ( 1 ) "*" "*" " " " " " " "*" " " " " " " " " " "
## 7 ( 1 ) " " "*" " " " " " " "*" " " "*" "*" "*" " "
## 8 ( 1 ) "*" "*" " " " " " " "*" " " " " " " "*" "*"
## CRBI CWalks LeagueN DivisionW PutOuts Assists Errors NewLeagueN
## 1 ( 1 ) "*" " " " " " " " " " " " " " "
## 2 ( 1 ) "*" " " " " " " " " " " " " " "
## 3 ( 1 ) "*" " " " " " " "*" " " " " " "
## 4 ( 1 ) "*" " " " " "*" "*" " " " " " "
## 5 ( 1 ) "*" " " " " "*" "*" " " " " " "
## 6 ( 1 ) "*" " " " " "*" "*" " " " " " "
## 7 ( 1 ) " " " " " " "*" "*" " " " " " "
## 8 ( 1 ) " " "*" " " "*" "*" " " " " " "
regfit19.full = regsubsets(Salary ~ ., data = Hitters, nvmax = 19)
reg.summary = summary(regfit19.full)
names(reg.summary)
## [1] "which" "rsq" "rss" "adjr2" "cp" "bic" "outmat" "obj"
plot(reg.summary$cp, xlab = "Number of Variables", ylab = "Cp")
which.min(reg.summary$cp)
## [1] 10
points(10, reg.summary$cp[10], pch = 20, col = "red")
plot(regfit19.full, scale = "Cp")
coef(regfit19.full, 10)
## (Intercept) AtBat Hits Walks CAtBat
## 162.5354420 -2.1686501 6.9180175 5.7732246 -0.1300798
## CRuns CRBI CWalks DivisionW PutOuts
## 1.4082490 0.7743122 -0.8308264 -112.3800575 0.2973726
## Assists
## 0.2831680
#page 274
#q5
regfit.fwd = regsubsets(Salary ~ ., data = Hitters, nvmax = 19,
method = "forward")
summary(regfit.fwd)
## Subset selection object
## Call: regsubsets.formula(Salary ~ ., data = Hitters, nvmax = 19, method = "forward")
## 19 Variables (and intercept)
## Forced in Forced out
## AtBat FALSE FALSE
## Hits FALSE FALSE
## HmRun FALSE FALSE
## Runs FALSE FALSE
## RBI FALSE FALSE
## Walks FALSE FALSE
## Years FALSE FALSE
## CAtBat FALSE FALSE
## CHits FALSE FALSE
## CHmRun FALSE FALSE
## CRuns FALSE FALSE
## CRBI FALSE FALSE
## CWalks FALSE FALSE
## LeagueN FALSE FALSE
## DivisionW FALSE FALSE
## PutOuts FALSE FALSE
## Assists FALSE FALSE
## Errors FALSE FALSE
## NewLeagueN FALSE FALSE
## 1 subsets of each size up to 19
## Selection Algorithm: forward
## AtBat Hits HmRun Runs RBI Walks Years CAtBat CHits CHmRun CRuns
## 1 ( 1 ) " " " " " " " " " " " " " " " " " " " " " "
## 2 ( 1 ) " " "*" " " " " " " " " " " " " " " " " " "
## 3 ( 1 ) " " "*" " " " " " " " " " " " " " " " " " "
## 4 ( 1 ) " " "*" " " " " " " " " " " " " " " " " " "
## 5 ( 1 ) "*" "*" " " " " " " " " " " " " " " " " " "
## 6 ( 1 ) "*" "*" " " " " " " "*" " " " " " " " " " "
## 7 ( 1 ) "*" "*" " " " " " " "*" " " " " " " " " " "
## 8 ( 1 ) "*" "*" " " " " " " "*" " " " " " " " " "*"
## 9 ( 1 ) "*" "*" " " " " " " "*" " " "*" " " " " "*"
## 10 ( 1 ) "*" "*" " " " " " " "*" " " "*" " " " " "*"
## 11 ( 1 ) "*" "*" " " " " " " "*" " " "*" " " " " "*"
## 12 ( 1 ) "*" "*" " " "*" " " "*" " " "*" " " " " "*"
## 13 ( 1 ) "*" "*" " " "*" " " "*" " " "*" " " " " "*"
## 14 ( 1 ) "*" "*" "*" "*" " " "*" " " "*" " " " " "*"
## 15 ( 1 ) "*" "*" "*" "*" " " "*" " " "*" "*" " " "*"
## 16 ( 1 ) "*" "*" "*" "*" "*" "*" " " "*" "*" " " "*"
## 17 ( 1 ) "*" "*" "*" "*" "*" "*" " " "*" "*" " " "*"
## 18 ( 1 ) "*" "*" "*" "*" "*" "*" "*" "*" "*" " " "*"
## 19 ( 1 ) "*" "*" "*" "*" "*" "*" "*" "*" "*" "*" "*"
## CRBI CWalks LeagueN DivisionW PutOuts Assists Errors NewLeagueN
## 1 ( 1 ) "*" " " " " " " " " " " " " " "
## 2 ( 1 ) "*" " " " " " " " " " " " " " "
## 3 ( 1 ) "*" " " " " " " "*" " " " " " "
## 4 ( 1 ) "*" " " " " "*" "*" " " " " " "
## 5 ( 1 ) "*" " " " " "*" "*" " " " " " "
## 6 ( 1 ) "*" " " " " "*" "*" " " " " " "
## 7 ( 1 ) "*" "*" " " "*" "*" " " " " " "
## 8 ( 1 ) "*" "*" " " "*" "*" " " " " " "
## 9 ( 1 ) "*" "*" " " "*" "*" " " " " " "
## 10 ( 1 ) "*" "*" " " "*" "*" "*" " " " "
## 11 ( 1 ) "*" "*" "*" "*" "*" "*" " " " "
## 12 ( 1 ) "*" "*" "*" "*" "*" "*" " " " "
## 13 ( 1 ) "*" "*" "*" "*" "*" "*" "*" " "
## 14 ( 1 ) "*" "*" "*" "*" "*" "*" "*" " "
## 15 ( 1 ) "*" "*" "*" "*" "*" "*" "*" " "
## 16 ( 1 ) "*" "*" "*" "*" "*" "*" "*" " "
## 17 ( 1 ) "*" "*" "*" "*" "*" "*" "*" "*"
## 18 ( 1 ) "*" "*" "*" "*" "*" "*" "*" "*"
## 19 ( 1 ) "*" "*" "*" "*" "*" "*" "*" "*"
plot(regfit.fwd, scale = "Cp")
#q6
set.seed(314)
train = sample(seq(263), 180, replace = FALSE)
regfit.fwd = regsubsets(Salary ~ ., data = Hitters[train, ],
nvmax = 19, method = "forward", nbest=9)
n.vars=NULL
val.errors=NULL
x.test = model.matrix(Salary ~ ., data = Hitters[-train, ])
for (i in 1:length(summary(regfit.fwd)$rss)) {
coefi = coef(regfit.fwd, id = i)
pred = x.test[, names(coefi)] %*% coefi
n.vars = c(n.vars,length(names(coefi))-1)
val.errors = c(val.errors,mean((Hitters$Salary[-train]-pred)^2))
}
#page 275
#q6 (suite)
mod.ordre=as.numeric(unlist(sapply(table(n.vars),
function(x){1:x})))
plot(n.vars,sqrt(val.errors), ylab ="Root MSE",pch =as.character(
mod.ordre),col=n.vars)
lines(n.vars[mod.ordre==1],sqrt(val.errors)[mod.ordre==1],
lwd=2,lty=2)
#page 276
#q7
predict.regsubsets = function(object, newdata, id, ...) {
form = as.formula(object$call[[2]])
mat = model.matrix(form, newdata)
coefi = coef(object, id = id)
mat[, names(coefi)] %*% coefi
}
#q8
library(glmnet)
## Loading required package: Matrix
## Loading required package: foreach
## Loaded glmnet 2.0-18
x=model.matrix(Salary~.-1,data=Hitters)
y=Hitters$Salary
lasso_model_cv=cv.glmnet(x,y)
plot(lasso_model_cv)
n_best_m=which(lasso_model_cv$lambda==lasso_model_cv$lambda.min)
lasso_model_cv$glmnet.fit$beta[,n_best_m]
## AtBat Hits HmRun Runs RBI
## -1.547343e+00 5.660897e+00 0.000000e+00 0.000000e+00 0.000000e+00
## Walks Years CAtBat CHits CHmRun
## 4.729691e+00 -9.595837e+00 0.000000e+00 0.000000e+00 5.108207e-01
## CRuns CRBI CWalks LeagueA LeagueN
## 6.594856e-01 3.927505e-01 -5.291586e-01 -3.206508e+01 3.285872e-14
## DivisionW PutOuts Assists Errors NewLeagueN
## -1.192990e+02 2.724045e-01 1.732025e-01 -2.058508e+00 0.000000e+00
#page 277
#Exercice 4.3
#q1
CancerSein <- read.csv("https://tinyurl.com/y3l6sh59")
#page 279
#Exercice 4.4
#q1
SidaChat <- read.csv("https://tinyurl.com/yxe6yxem")
#page 282
#Exercice 4.5
#q1
Vitamines <- read.csv("https://tinyurl.com/y3shxcsd")
#page 283
#Exercice 4.6
#q1
Beton <- read.csv("https://tinyurl.com/y4w2qv9t")
#page 284
#Exercice 4.7
#q1
chal <- read.csv("https://tinyurl.com/yyb3cztf")
#q2
cdplot(as.factor(Defaillance)~Temperature, data=chal,
ylab="Defaillance")
#page 285
#q3
chal.glm <- glm(Defaillance~Temperature,data=chal,
family="binomial")
#Si le package n'est pas installe, enlever le commentaire
#puis executer la commande ci-dessous.
#if(!require(hnp)){install.packages("hnp")}
library(hnp)
hnp(chal.glm, sim = 99, conf = 0.95)
## Binomial model
#page 286
#q6
#Si le package n'est pas installe, enlever le commentaire
#puis executer la commande ci-dessous.
#if(!require(rms)){install.packages("rms")}
library(rms)
## Loading required package: Hmisc
## Loading required package: survival
## Loading required package: Formula
##
## Attaching package: 'Hmisc'
## The following objects are masked from 'package:dplyr':
##
## src, summarize
## The following objects are masked from 'package:base':
##
## format.pval, units
## Loading required package: SparseM
##
## Attaching package: 'SparseM'
## The following object is masked from 'package:base':
##
## backsolve
##
## Attaching package: 'rms'
## The following objects are masked from 'package:car':
##
## Predict, vif
chal.lrm <- lrm(Defaillance~Temperature, data=chal, x=TRUE, y=TRUE)
print(chal.lrm)
## Logistic Regression Model
##
## lrm(formula = Defaillance ~ Temperature, data = chal, x = TRUE,
## y = TRUE)
##
## Model Likelihood Discrimination Rank Discrim.
## Ratio Test Indexes Indexes
## Obs 24 LR chi2 5.94 R2 0.313 C 0.723
## 0 17 d.f. 1 g 1.426 Dxy 0.445
## 1 7 Pr(> chi2) 0.0148 gr 4.164 gamma 0.461
## max |deriv| 4e-06 gp 0.247 tau-a 0.192
## Brier 0.157
##
## Coef S.E. Wald Z Pr(>|Z|)
## Intercept 10.8753 5.7031 1.91 0.0565
## Temperature -0.1713 0.0834 -2.05 0.0401
##
residuals(chal.lrm, "gof")
## Sum of squared errors Expected value|H0 SD
## 3.7572772 3.7345914 0.1360759
## Z P
## 0.1667141 0.8675950
#Exercice 4.8
#q1
Cypermethrine <- read.csv("https://tinyurl.com/y4deakfd")
#page 288
#Exercice 4.9
#q1
poly <- read.csv("https://tinyurl.com/yyhhcw37")
#q2
poly_glm1 <- glm(nombre~traitement+age, family=poisson(),data=poly)
library(hnp)
hnp(poly_glm1, sim = 99, conf = 0.95)
## Poisson model
summary(poly_glm1)
##
## Call:
## glm(formula = nombre ~ traitement + age, family = poisson(),
## data = poly)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -4.2212 -3.0536 -0.1802 1.4459 5.8301
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 3.169941 0.168210 18.84 < 2e-16 ***
## traitementplacebo 1.359083 0.117643 11.55 < 2e-16 ***
## age -0.038830 0.005955 -6.52 7.02e-11 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for poisson family taken to be 1)
##
## Null deviance: 378.66 on 19 degrees of freedom
## Residual deviance: 179.54 on 17 degrees of freedom
## AIC: 273.88
##
## Number of Fisher Scoring iterations: 5
confint(poly_glm1)
## Waiting for profiling to be done...
## 2.5 % 97.5 %
## (Intercept) 2.83903843 3.498933
## traitementplacebo 1.13376052 1.595497
## age -0.05074339 -0.027393
poly_glm2 <- glm(nombre~traitement+age, family=quasipoisson(),
data=poly)
library(hnp)
hnp(poly_glm2, sim = 99, conf = 0.95)
## Quasi-Poisson model
summary(poly_glm2)
##
## Call:
## glm(formula = nombre ~ traitement + age, family = quasipoisson(),
## data = poly)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -4.2212 -3.0536 -0.1802 1.4459 5.8301
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 3.16994 0.55095 5.754 2.34e-05 ***
## traitementplacebo 1.35908 0.38533 3.527 0.00259 **
## age -0.03883 0.01951 -1.991 0.06284 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for quasipoisson family taken to be 10.72805)
##
## Null deviance: 378.66 on 19 degrees of freedom
## Residual deviance: 179.54 on 17 degrees of freedom
## AIC: NA
##
## Number of Fisher Scoring iterations: 5
confint(poly_glm2)
## Waiting for profiling to be done...
## 2.5 % 97.5 %
## (Intercept) 2.06974676 4.247700305
## traitementplacebo 0.65393034 2.184666229
## age -0.07969808 -0.003027068
poly_glm3 <- glm.nb(nombre~traitement+age,data=poly)
library(hnp)
hnp(poly_glm3, sim = 99, conf = 0.95)
## Negative binomial model (using MASS package)
summary(poly_glm3)
##
## Call:
## glm.nb(formula = nombre ~ traitement + age, data = poly, init.theta = 1.719491,
## link = log)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.83270 -1.13898 -0.08851 0.33637 1.89785
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 3.15791 0.55753 5.664 1.48e-08 ***
## traitementplacebo 1.36812 0.36903 3.707 0.000209 ***
## age -0.03856 0.02095 -1.840 0.065751 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for Negative Binomial(1.7195) family taken to be 1)
##
## Null deviance: 36.734 on 19 degrees of freedom
## Residual deviance: 22.002 on 17 degrees of freedom
## AIC: 164.88
##
## Number of Fisher Scoring iterations: 1
##
##
## Theta: 1.719
## Std. Err.: 0.607
##
## 2 x log-likelihood: -156.880
confint(poly_glm3)
## Waiting for profiling to be done...
## 2.5 % 97.5 %
## (Intercept) 2.0665705 4.231950890
## traitementplacebo 0.6336661 2.096075160
## age -0.0776086 0.004319903
#page 289
with(poly,plot(nombre~age,type="n",ylab="Nombre de polypes",
xlab="\u00c2ge"))
with(poly,points(age[traitement=="placebo"],fitted(poly_glm2)[
traitement=="placebo"],pch="P",col="red"))
xv1<-seq(0,50,.05)
yv1<-predict(poly_glm2,list(traitement=as.factor(rep("placebo",
length(xv1))),age=xv1))
lines(xv1,exp(yv1),col="red")
with(poly,points(age[traitement=="placebo"],
nombre[traitement=="placebo"],pch="p"))
with(poly,points(age[traitement=="medicament"],fitted(
poly_glm2)[traitement=="medicament"],pch="D",col="blue"))
with(poly,points(age[traitement=="medicament"],
nombre[traitement=="medicament"],pch="d"))
yv2<-predict(poly_glm2,list(traitement=as.factor(rep("medicament",
length(xv1))),age=xv1))
lines(xv1,exp(yv2),col="blue")
with(poly,plot(nombre~age,type="n",ylab="Nombre de polypes",
xlab="\u00c2ge"))
with(poly,points(age[traitement=="placebo"],fitted(poly_glm3)[
traitement=="placebo"],pch="P",col="red"))
xv1<-seq(0,50,.05)
yv1<-predict(poly_glm3,list(traitement=as.factor(rep("placebo",
length(xv1))),age=xv1))
lines(xv1,exp(yv1),col="red")
with(poly,points(age[traitement=="placebo"],
nombre[traitement=="placebo"],pch="p"))
with(poly,points(age[traitement=="medicament"],fitted(
poly_glm3)[traitement=="medicament"],pch="D",col="blue"))
with(poly,points(age[traitement=="medicament"],
nombre[traitement=="medicament"],pch="d"))
yv2<-predict(poly_glm3,list(traitement=as.factor(rep("medicament",
length(xv1))),age=xv1))
lines(xv1,exp(yv2),col="blue")
