#Chapitre 4 : solution des exercices
#Complement en ligne
#Exercice 4.4
#q1
CancerSein <- read.csv("https://tinyurl.com/y3l6sh59")
#q2
layout(c(1,2))
plot(Survie~Traitement, data=CancerSein)
plot(Age~Traitement,data=CancerSein)
layout(1)
tapply(CancerSein$Age,CancerSein$Traitement,mean)
## A B C
## 49.00000 49.23889 42.86667
aov1 <- aov(Survie~Traitement,data=CancerSein)
shapiro.test(residuals(aov1))
##
## Shapiro-Wilk normality test
##
## data: residuals(aov1)
## W = 0.9784, p-value = 0.3434
library(car)
## Loading required package: carData
bartlett.test(residuals(aov1),CancerSein$Traitement)
##
## Bartlett test of homogeneity of variances
##
## data: residuals(aov1) and CancerSein$Traitement
## Bartlett's K-squared = 0.12452, df = 2, p-value = 0.9396
summary(aov1)
## Df Sum Sq Mean Sq F value Pr(>F)
## Traitement 2 55.13 27.567 6.895 0.00204 **
## Residuals 59 235.89 3.998
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#Exercice 4.5
#q1
SidaChat <- read.csv("https://tinyurl.com/yxe6yxem")
#q2
layout(matrix(c(1,1,2,3),byrow=T,nrow=2))
plot(LnT4~Jours, data=SidaChat)
plot(LnT4~Sexe, data=SidaChat)
plot(Jours~Sexe, data=SidaChat)
layout(1)
#Exercice 4.6
#q1
chal <- read.csv("https://tinyurl.com/yyb3cztf")
#q2
plot(Defaillance~Temperature,data=chal)
cdplot(as.factor(Defaillance)~Temperature, data=chal,
ylab="Defaillance")
#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)
## Loading required package: MASS
hnp(chal.glm, sim = 99, conf = 0.95)
## Binomial model
summary(chal.glm)
##
## Call:
## glm(formula = Defaillance ~ Temperature, family = "binomial",
## data = chal)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.2125 -0.8253 -0.4706 0.5907 2.0512
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 10.87535 5.70291 1.907 0.0565 .
## Temperature -0.17132 0.08344 -2.053 0.0400 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 28.975 on 23 degrees of freedom
## Residual deviance: 23.030 on 22 degrees of freedom
## AIC: 27.03
##
## Number of Fisher Scoring iterations: 4
anova(chal.glm,test = "Chisq")
## Analysis of Deviance Table
##
## Model: binomial, link: logit
##
## Response: Defaillance
##
## Terms added sequentially (first to last)
##
##
## Df Deviance Resid. Df Resid. Dev Pr(>Chi)
## NULL 23 28.975
## Temperature 1 5.9441 22 23.030 0.01477 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
ggiraphExtra::ggPredict(chal.glm)
#q4
confint(chal.glm, parm="Temperature")
## Waiting for profiling to be done...
## 2.5 % 97.5 %
## -0.3779407 -0.0305637
exp(coef(chal.glm))["Temperature"]
## Temperature
## 0.8425515
exp(confint(chal.glm, parm="Temperature"))
## Waiting for profiling to be done...
## 2.5 % 97.5 %
## 0.6852711 0.9698986
#q5
predict(chal.glm,newdata = list(Temperature=31),type = "response")
## 1
## 0.9961828
#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: lattice
## Loading required package: survival
## Loading required package: Formula
## Loading required package: ggplot2
##
## Attaching package: 'Hmisc'
## 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.7
#q1
read.csv("https://tinyurl.com/y3shxcsd")
## Calorie Vitamine Poids
## 1 1 1 84
## 2 1 2 62
## 3 2 1 87
## 4 2 2 103
## 5 1 1 66
## 6 1 2 84
## 7 2 1 92
## 8 2 2 107
## 9 1 1 82
## 10 1 2 73
## 11 2 1 77
## 12 2 2 95
## 13 1 1 62
## 14 1 2 75
## 15 2 1 88
## 16 2 2 96
## 17 1 1 66
## 18 1 2 59
## 19 2 1 89
## 20 2 2 90
## 21 1 1 56
## 22 1 2 74
## 23 2 1 101
## 24 2 2 116
## 25 1 1 79
## 26 1 2 74
## 27 2 1 95
## 28 2 2 112
## 29 1 1 89
## 30 1 2 74
## 31 2 1 91
## 32 2 2 92
#Exercice 4.8
#q1
read.csv("https://tinyurl.com/y4w2qv9t")
## Melangeur Casseur Resistance
## 1 M1 C1 5280
## 2 M1 C1 5520
## 3 M1 C1 4760
## 4 M1 C1 5800
## 5 M1 C2 4340
## 6 M1 C2 4400
## 7 M1 C2 5020
## 8 M1 C2 6200
## 9 M1 C3 4160
## 10 M1 C3 5180
## 11 M1 C3 5320
## 12 M1 C3 4600
## 13 M2 C1 4420
## 14 M2 C1 5280
## 15 M2 C1 5580
## 16 M2 C1 4900
## 17 M2 C2 5340
## 18 M2 C2 4880
## 19 M2 C2 4960
## 20 M2 C2 6200
## 21 M2 C3 4180
## 22 M2 C3 4800
## 23 M2 C3 4600
## 24 M2 C3 4480
## 25 M3 C1 5360
## 26 M3 C1 6160
## 27 M3 C1 5680
## 28 M3 C1 5500
## 29 M3 C2 5720
## 30 M3 C2 4760
## 31 M3 C2 5620
## 32 M3 C2 5560
## 33 M3 C3 4460
## 34 M3 C3 4930
## 35 M3 C3 4680
## 36 M3 C3 5600
#Exercice 4.9
#q1
read.csv("https://tinyurl.com/y4deakfd")
## Total N.morts Niveau.de.dose Sexe
## 1 20 1 1 M
## 2 20 4 2 M
## 3 20 9 4 M
## 4 20 13 8 M
## 5 20 18 16 M
## 6 20 20 32 M
## 7 20 0 1 F
## 8 20 2 2 F
## 9 20 6 4 F
## 10 20 10 8 F
## 11 20 12 16 F
## 12 20 16 32 F
#Exercice 4.10
#q1
poly <- read.csv("https://tinyurl.com/yyhhcw37")
#q1
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
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")
