library("GLMsData")
library("splines")
library("MASS")
library("statmod")
data(downs)
head(downs)
## Age Births DS
## 1 17.0 13555 16
## 2 18.5 13675 15
## 3 19.5 18752 16
## 4 20.5 22005 22
## 5 21.5 23896 16
## 6 22.5 24667 12
n <- dim(downs)[1]
Statistiques descriptives
summary(downs)
## Age Births DS
## Min. :17.00 Min. : 249 Min. : 7.0
## 1st Qu.:24.75 1st Qu.: 4179 1st Qu.:13.0
## Median :32.00 Median :11485 Median :16.0
## Mean :31.99 Mean :11829 Mean :17.3
## 3rd Qu.:39.25 3rd Qu.:19090 3rd Qu.:21.5
## Max. :47.10 Max. :24807 Max. :33.0
downs$Taux <- downs$DS/downs$Births
downs$Taux.1000 <- 1000*downs$Taux
summary(downs$Taux.1000)
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.4865 0.8549 1.3099 6.2579 6.4561 36.9128
plot(downs$Age, downs$Births, type = "h", lwd = 3,
xlab = "Age de la mere", ylab = "Nombre de naissances")
plot(downs$Age, downs$DS, pch = 19,
xlab = "Age de la mere", ylab = "Nombre de naissances avec syndrome de Down")
plot(downs$Age, downs$Taux.1000, pch = 19,
xlab = "Age de la mere", ylab = "Taux pour 1000 naissances")
plot(downs$Age, log(downs$Taux), pch = 19,
xlab = "Age de la mere", ylab = "log(DS/Births)")
Définition du modèle
modele.lineaire <- glm(DS ~ offset(log(Births)) + Age,
family = poisson, data = downs)
summary(modele.lineaire)
##
## Call:
## glm(formula = DS ~ offset(log(Births)) + Age, family = poisson,
## data = downs)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -10.553981 0.214179 -49.28 <2e-16 ***
## Age 0.137191 0.006463 21.23 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for poisson family taken to be 1)
##
## Null deviance: 622.56 on 29 degrees of freedom
## Residual deviance: 183.03 on 28 degrees of freedom
## AIC: 326.1
##
## Number of Fisher Scoring iterations: 5
modele.quadratique <- glm(DS ~ offset(log(Births)) + Age + I(Age^2),
family = poisson, data = downs)
summary(modele.quadratique)
##
## Call:
## glm(formula = DS ~ offset(log(Births)) + Age + I(Age^2), family = poisson,
## data = downs)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -2.3425778 0.6587995 -3.556 0.000377 ***
## Age -0.4071225 0.0432067 -9.423 < 2e-16 ***
## I(Age^2) 0.0084991 0.0006705 12.677 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for poisson family taken to be 1)
##
## Null deviance: 622.56 on 29 degrees of freedom
## Residual deviance: 46.05 on 27 degrees of freedom
## AIC: 191.12
##
## Number of Fisher Scoring iterations: 4
modele.logit <- glm(cbind(DS, Births - DS) ~ Age,
family = binomial(link = "logit"), data = downs)
cbind(coef(modele.lineaire), coef(modele.logit))
## [,1] [,2]
## (Intercept) -10.5539807 -10.5634589
## Age 0.1371908 0.1375859
Modélisation de l’âge avec des splines
modele.ns.2 <- glm(DS ~ offset(log(Births)) + ns(Age, df = 2), family = poisson, data = downs)
modele.ns.3 <- glm(DS ~ offset(log(Births)) + ns(Age, df = 3), family = poisson, data = downs)
modele.ns.4 <- glm(DS ~ offset(log(Births)) + ns(Age, df = 4), family = poisson, data = downs)
modele.ns.5 <- glm(DS ~ offset(log(Births)) + ns(Age, df = 5), family = poisson, data = downs)
modele.ns.6 <- glm(DS ~ offset(log(Births)) + ns(Age, df = 6), family = poisson, data = downs)
AIC.modeles <- c(extractAIC(modele.lineaire, k = 2)[2],
extractAIC(modele.quadratique, k = 2)[2],
extractAIC(modele.ns.2, k = 2)[2],
extractAIC(modele.ns.3, k = 2)[2],
extractAIC(modele.ns.4, k = 2)[2],
extractAIC(modele.ns.5, k = 2)[2],
extractAIC(modele.ns.6, k = 2)[2])
BIC.modeles <- c(extractAIC(modele.lineaire, k = log(n))[2],
extractAIC(modele.quadratique, k = log(n))[2],
extractAIC(modele.ns.2, k = log(n))[2],
extractAIC(modele.ns.3, k = log(n))[2],
extractAIC(modele.ns.4, k = log(n))[2],
extractAIC(modele.ns.5, k = log(n))[2],
extractAIC(modele.ns.6, k = log(n))[2])
data.frame(modele = c("lineaire", "quadratique", "ns df=2", "ns df=3",
"ns df=4", "ns df=5", "ns df=6"),
AIC = round(AIC.modeles, 2), BIC = round(BIC.modeles, 2))
## modele AIC BIC
## 1 lineaire 326.10 328.90
## 2 quadratique 191.12 195.32
## 3 ns df=2 187.60 191.80
## 4 ns df=3 188.73 194.33
## 5 ns df=4 188.96 195.97
## 6 ns df=5 185.38 193.78
## 7 ns df=6 185.78 195.59
modele.bs <- glm(DS ~ offset(log(Births)) + bs(Age, knots = c(25, 35)),
family = poisson, data = downs)
extractAIC(modele.bs, k = 2)[2]
## [1] 184.9074
extractAIC(modele.bs, k = log(n))[2]
## [1] 193.3145
modele <- modele.ns.5
summary(modele)
##
## Call:
## glm(formula = DS ~ offset(log(Births)) + ns(Age, df = 5), family = poisson,
## data = downs)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -6.6986 0.1901 -35.228 < 2e-16 ***
## ns(Age, df = 5)1 -0.1620 0.2511 -0.645 0.519
## ns(Age, df = 5)2 -0.1007 0.3118 -0.323 0.747
## ns(Age, df = 5)3 2.6314 0.2234 11.780 < 2e-16 ***
## ns(Age, df = 5)4 2.4182 0.4845 4.991 6.01e-07 ***
## ns(Age, df = 5)5 3.9408 0.2483 15.871 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for poisson family taken to be 1)
##
## Null deviance: 622.564 on 29 degrees of freedom
## Residual deviance: 34.307 on 24 degrees of freedom
## AIC: 185.38
##
## Number of Fisher Scoring iterations: 4
p <- length(coef(modele)) - 1
stat.test.D <- deviance(modele.lineaire) - deviance(modele)
df.test <- modele.lineaire$df.residual - modele$df.residual
stat.test.D
## [1] 148.7215
1 - pchisq(stat.test.D, df = df.test)
## [1] 0
Interprétation du modèle
axe.age <- seq(min(downs$Age), max(downs$Age), 0.1)
nouvelles.donnees <- predict(modele, newdata = data.frame(Age = axe.age, Births = 1000),
se.fit = TRUE)
eta.chapeau <- nouvelles.donnees$fit
se.eta.chapeau <- nouvelles.donnees$se.fit
taux.chapeau <- exp(eta.chapeau)
borne.inf <- exp(eta.chapeau - 1.96*se.eta.chapeau)
borne.sup <- exp(eta.chapeau + 1.96*se.eta.chapeau)
plot(downs$Age, downs$Taux.1000, pch = 19,
xlab = "Age de la mere", ylab = "Taux pour 1000 naissances")
lines(axe.age, taux.chapeau, col = "red")
lines(axe.age, borne.inf, lty = 2, col = "blue")
lines(axe.age, borne.sup, lty = 2, col = "blue")
plot(downs$Age, log(downs$Taux.1000), pch = 19,
xlab = "Age de la mere", ylab = "log(Taux pour 1000 naissances)")
lines(axe.age, eta.chapeau, col = "red")
lines(axe.age, eta.chapeau - 1.96*se.eta.chapeau, lty = 2, col = "blue")
lines(axe.age, eta.chapeau + 1.96*se.eta.chapeau, lty = 2, col = "blue")
ages <- c(20, 25, 30, 35, 40, 45)
prediction.ages <- predict(modele, newdata = data.frame(Age = ages, Births = 1000),
se.fit = TRUE)
round(data.frame(Age = ages,
Taux.1000 = exp(prediction.ages$fit),
B.Inf = exp(prediction.ages$fit - 1.96*prediction.ages$se.fit),
B.Sup = exp(prediction.ages$fit + 1.96*prediction.ages$se.fit)), 3)
## Age Taux.1000 B.Inf B.Sup
## 1 20 0.851 0.707 1.025
## 2 25 0.723 0.603 0.866
## 3 30 1.001 0.806 1.241
## 4 35 1.685 1.342 2.115
## 5 40 8.000 6.764 9.462
## 6 45 27.552 21.380 35.506
taux.ages <- exp(prediction.ages$fit)
taux.ages/taux.ages[1]
## 1 2 3 4 5 6
## 1.000000 0.848949 1.175182 1.979064 9.396580 32.360783
Validation du modèle
Adéquation globale
Deviance <- deviance(modele)
Pearson <- sum(resid(modele, type = "pearson")^2)
df <- modele$df.residual
(pvaleur.Deviance <- 1 - pchisq(Deviance, df = df))
## [1] 0.07935752
(pvaleur.Pearson <- 1 - pchisq(Pearson, df = df))
## [1] 0.05725585
Résidus
residus.deviance <- resid(modele)
residus.Pearson <- resid(modele, type = "pearson")
residus.quantile <- qresid(modele)
residus.standardises <- rstandard(modele)
residus.studentises <- rstudent(modele)
plot(residus.deviance, xlab = "Numero d'observation", ylab = "Residus deviance")
abline(h = 0, col = "red")
abline(h = 3, lty = 2, col = "green")
abline(h = -3, lty = 2, col = "green")
plot(residus.Pearson, xlab = "Numero d'observation", ylab = "Residus Pearson")
abline(h = 0, col = "red")
abline(h = 3, lty = 2, col = "green")
abline(h = -3, lty = 2, col = "green")
scatter.smooth(residus.standardises ~ downs$Age,
xlab = "Age de la mere", ylab = "Residus standardises")
abline(h = 0, col = "red")
plot(fitted(modele), residus.standardises,
xlab = "Valeurs predites", ylab = "Residus standardises")
abline(h = 0, col = "red")
qqnorm(residus.quantile, xlab = "Quantiles theoriques",
ylab = "Quantiles empiriques", main = "")
qqline(residus.quantile)
plot(fitted(modele), downs$DS, xlab = "Valeurs predites de DS",
ylab = "Valeurs observees de DS")
abline(0, 1, col = "red")
Observations aberrantes et influentes
which(abs(as.vector(residus.standardises)) > 3)
## integer(0)
which(abs(as.vector(residus.studentises)) > 3)
## integer(0)
h <- hatvalues(modele)
(indices.h <- which(as.vector(h) > 3*(p+1)/n))
## [1] 1
d1 <- cooks.distance(modele)
max(d1)
## [1] 0.3979281
qf(0.5, p+1, n-(p+1))
## [1] 0.9168687
(indices.d1 <- which(as.vector(d1) > qf(0.5, p+1, n-(p+1))))
## integer(0)
d2 <- dffits(modele)
max(abs(d2))
## [1] 1.384648
3*sqrt(p/(n-p))
## [1] 1.341641
(indices.d2 <- which(as.vector(abs(d2)) > 3*sqrt(p/(n-p))))
## [1] 30
d3 <- dfbetas(modele)
max(abs(d3))
## [1] 1.262701
3/sqrt(n)
## [1] 0.5477226
d3.matrice <- abs(d3) > 3/sqrt(n)
(indices.d3 <- which(as.vector(rowSums(d3.matrice)) > 0))
## [1] 11 13 19 28 30
d4 <- covratio(modele)
max(abs(1-d4))
## [1] 2.21738
3*(p+1)/n
## [1] 0.6
(indices.d4 <- which(as.vector(abs(1-d4)) > 3*(p+1)/n))
## [1] 1 2 30
intersect(indices.d2, indices.d3)
## [1] 30
plot(d1, type = "h", xlab = "Numero d'observation", ylab = "Distance de Cook")
abline(h = qf(0.5, p+1, n-(p+1)), lty = 2, col = "red")
plot(h, xlab = "Numero d'observation", ylab = "Leviers")
abline(h = 3*(p+1)/n, lty = 2, col = "red")
indices <- intersect(indices.d2, indices.d3)
modele.sans <- glm(DS ~ offset(log(Births)) + ns(Age, df = 5),
family = poisson, data = downs[-indices, ])
data.frame(coef.avec = coef(modele), coef.sans = coef(modele.sans))
## coef.avec coef.sans
## (Intercept) -6.6985589 -6.6948587
## ns(Age, df = 5)1 -0.1619720 -0.2180973
## ns(Age, df = 5)2 -0.1006517 -0.1061848
## ns(Age, df = 5)3 2.6313948 2.0608604
## ns(Age, df = 5)4 2.4181918 2.2075907
## ns(Age, df = 5)5 3.9408316 3.8623114
data.frame(s.e.avec = sqrt(diag(vcov(modele))),
s.e.sans = sqrt(diag(vcov(modele.sans))))
## s.e.avec s.e.sans
## (Intercept) 0.1901494 0.1923201
## ns(Age, df = 5)1 0.2510948 0.2555288
## ns(Age, df = 5)2 0.3117636 0.3214309
## ns(Age, df = 5)3 0.2233863 0.2262517
## ns(Age, df = 5)4 0.4845314 0.4819683
## ns(Age, df = 5)5 0.2483096 0.2068306