#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")