On a relevé l’âge et la présence (1) ou l’absence (0) d’une maladie cardiovasculaire chez 100 individus. Les données sont stockées dans le fichier : sur une ligne donnée, la variable fournit l’âge d’un individu tandis que la variable prend la valeur 1 en cas de présence d’une maladie cardiovasculaire chez cet individu et la valeur 0 sinon. Les variables et donnent respectivement le numéro d’un individu et sa classe d’âge. Ces données sont proposées par Hosmer et Lemeshow

L’objectif est d’expliquer la présence/absence d’une maladie cardio vasculaire (notée aussi CHD), par l’âge des patients.

Lecture, analyse descriptive des données et visualisation

On acquière les données et on en visualise les premières lignes grâce aux commandes

#Données 
maladie= read.table('maladie_cardiovasculaire.txt',header=TRUE)

library(tidyverse)
glimpse(maladie)
## Observations: 100
## Variables: 4
## $ ID   <int> 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17...
## $ AGRP <int> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 2, 2, 2, ...
## $ AGE  <int> 20, 23, 24, 25, 25, 26, 26, 28, 28, 29, 30, 30, 30, 30, 3...
## $ CHD  <int> 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, ...
#Transformation
maladie$CHD = factor(maladie$CHD, levels = c(0,1), labels = c("No","Yes"))
maladie$AGRP = factor(maladie$AGRP)
glimpse(maladie)
## Observations: 100
## Variables: 4
## $ ID   <int> 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17...
## $ AGRP <fct> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 2, 2, 2, ...
## $ AGE  <int> 20, 23, 24, 25, 25, 26, 26, 28, 28, 29, 30, 30, 30, 30, 3...
## $ CHD  <fct> No, No, No, No, Yes, No, No, No, No, No, No, No, No, No, ...
summary(maladie)
##        ID              AGRP         AGE         CHD    
##  Min.   :  1.00   7      :17   Min.   :20.00   No :57  
##  1st Qu.: 25.75   2      :15   1st Qu.:34.75   Yes:43  
##  Median : 50.50   4      :15   Median :44.00           
##  Mean   : 50.50   5      :13   Mean   :44.38           
##  3rd Qu.: 75.25   3      :12   3rd Qu.:55.00           
##  Max.   :100.00   1      :10   Max.   :69.00           
##                   (Other):18

Dans ce fichier de données nous avons une seule variable explicative l’âge. Afin de visualiser la relation entre la présence ou l’absence d’une maladie coronarienne en fonctionde l’âge nous allons représenter les données sur un graphe. Que peut-on en remarquer ?

#Plot (avec la fonction qplot)
qplot(data = maladie, CHD)

qplot(data=maladie, AGE,CHD)

Pourquoi la régression linéaire de la variable CHD en fonction de la variable AGE n’a pas de sens ?

Pour représenter la relation fonctionnelle entre la probabilité de maladie coronarienne et l’âge nous allons définir des catégories d’âge et calculer le pourcentage de maladie coronarienne dans chacune de ces catégories

group.AGRP <- group_by(maladie,AGRP)
maladie.summary <- summarise(group.AGRP,AGE = (max(AGE)+min(AGE))/2, p = sum(CHD=="Yes")/n())
maladie.summary
## # A tibble: 8 x 3
##   AGRP    AGE     p
##   <fct> <dbl> <dbl>
## 1 1      24.5 0.100
## 2 2      32.0 0.133
## 3 3      37.0 0.250
## 4 4      42.0 0.333
## 5 5      47.0 0.462
## 6 6      52.0 0.625
## 7 7      57.0 0.765
## 8 8      64.5 0.800
ggplot(maladie.summary,aes(x=AGE,y=p))+ geom_point(aes(col="p")) +
  geom_point(data=maladie,aes(x = AGE,y=as.numeric(CHD=="Yes"),col="CHD")) + 
  scale_color_discrete(name = "Legend") +
   geom_segment(data = data.frame(x = c(19,29,34,39,44,49,54,59,69), y = rep(0,9)), 
                aes(x = x, y = y, xend = x, yend = y+1), linetype = 2, color = "black")

Modèle Logistique

Modéliser les modalités de CHD (présence ou absence) en termes de proportion par rapport à AGE.

maladie.logit = glm(CHD~AGE, family=binomial(link="logit"),data=maladie)
coef(maladie.logit)
## (Intercept)         AGE 
##  -5.3094534   0.1109211
#Resumé de la régression logistique
summary(maladie.logit)
## 
## Call:
## glm(formula = CHD ~ AGE, family = binomial(link = "logit"), data = maladie)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -1.9718  -0.8456  -0.4576   0.8253   2.2859  
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept) -5.30945    1.13365  -4.683 2.82e-06 ***
## AGE          0.11092    0.02406   4.610 4.02e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 136.66  on 99  degrees of freedom
## Residual deviance: 107.35  on 98  degrees of freedom
## AIC: 111.35
## 
## Number of Fisher Scoring iterations: 4
beta0=coef(maladie.logit)[1]
beta1=coef(maladie.logit)[2]
abscisse1 = seq(20,70,length=100)
maladie.abscisse1 = data.frame(abscisse1 = abscisse1)
ggplot(maladie.summary,aes(x=AGE,y=p)) + geom_point(aes(color = "p")) + 
  geom_line(data = maladie.abscisse1, aes(x = abscisse1, 
                                          y = plogis(beta0 + beta1*abscisse1), 
                                          color = "logist")) + xlab("AGE")+
  ylab("") + scale_color_discrete(name = "Legend")

Une manière de mesurer la qualité d’un modèle est le calcul d’une matrice de confusion, c’est-à-dire le tableau croisé des valeurs observées et celles des valeurs prédites en appliquant le modèle aux données d’origine.

#score
score.glm  =  predict.glm(maladie.logit, maladie, type="response")
#Matrice de confusion (au seuil fixé 0.5)
pred.glm  = as.numeric(score.glm >=0.5)
pred.glm = factor(pred.glm, levels = c(0,1), labels = c("No","Yes"))
table(maladie$CHD,pred.glm)
##      pred.glm
##       No Yes
##   No  45  12
##   Yes 14  29
mean(maladie$CHD!=pred.glm)
## [1] 0.26

ESt-ce que cette erreur empirique (calculé sur les données déjà observées) est un bon estimateur de la moyenne de l’erreur de prédiction pour une nouvelle donnée ?

Pourquoi la validation croisée est-elle une meilleur méthode ? (Expliquer au passage la méthode)