Dans ce TP, vous allez appliquer la régression logistique sur le cas d'étude portant sur les maladies cardio-vasculaires.
Importez les données
On charge tout d'abord la librairie ggplot2 :
library(ggplot2)On importe ensuite les données :
maladie <- read.table("maladie.txt",header=TRUE,sep=";",dec=".")
Visualisez le nuage de points
Pour étudier le fait d'être malade en fonction de l'âge, on peut visualiser le nuage de points :
ggplot()+
geom_point(data=maladie,aes(x=age,y=chd))
Il y a des 0 et des 1, mais il est ici difficile de dire si l'on est plus ou moins malade en fonction de l'âge.
On voit également qu'une régression linéaire sur un tel nuage de points n'aurait aucun sens, car elle nous donnerait des valeurs qui ne seraient quasiment jamais sur 0 ni 1.
Calculez les proportions de malades
On peut calculer des classes d'âge et les proportions de malades associées.
maladie$cl_age <- cut(maladie$age,breaks=seq(15-1e-10,65,by=10))
prop <- prop.table(table(maladie$cl_age,maladie$chd),1)
prop_chd <- data.frame(age=c(15,rep(seq(25,55,by=10),each=2),65),
prop_chd=rep(prop[,2],each=2))On peut représenter ces proportions :
ggplot()+
geom_point(data=maladie,aes(x=age,y=chd))+
geom_line(data=prop_chd,aes(x=age,y=prop_chd,colour="blue"))+
xlab("age")+
ylab("chd")+
scale_colour_manual(values="blue",label="Proportion de malades",name=" ")
On y voit une fonction en escalier avec une forme de "S".
Effectuez la régression logistique
Effectuons donc une régression logistique de CHD en fonction de l'âge :
reg_log1 <- glm(chd~age,family="binomial",data=maladie)
summary(reg_log1)Voici le résultat :
Call:
glm(formula = chd ~ age, family = "binomial", data = maladie)
Deviance Residuals:
Min 1Q Median 3Q Max
-1.4321 -0.9215 -0.5392 1.0952 2.2433
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -3.521710 0.416031 -8.465 < 2e-16 ***
age 0.064108 0.008532 7.513 5.76e-14 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
(Dispersion parameter for binomial family taken to be 1)
Null deviance: 596.11 on 461 degrees of freedom
Residual deviance: 525.56 on 460 degrees of freedom
AIC: 529.56
Number of Fisher Scoring iterations: 4On obtient les paramètres estimés : et . Enregistrons-les :
beta1 <- reg_log1$coefficients[1]
beta2 <- reg_log1$coefficients[2]Dans le but de tracer la courbe logistique entre les abscisses x=15 et x=65 , on définit une séquence de 15 à 65 par pas de 500, puis on la place dans la variable x. On calcule ensuite les ordonnées de la courbe, grâce à l'expression de la courbe en S :
Nous plaçons ces ordonnées dans la variable y. Enfin, avec x et y, nous créons un dataframe :
x <- seq(15,65,len=500)
y <- exp(beta1+beta2*x)/(1+exp(beta1+beta2*x))
reg_log <- data.frame(age=x,prop_chd=y)Si l'on souhaite superposer la fonction de lien obtenue par régression logistique sur le graphique précédent, nous obtenons ceci :
ggplot()+
geom_point(data=maladie,aes(x=age,y=chd))+
geom_line(data=prop_chd,aes(x=age,y=prop_chd,colour="blue"))+
geom_line(data=reg_log,aes(x=age,y=prop_chd,colour="red"))+
xlab("age")+
ylab("chd")+
scale_colour_manual(values=c("blue","red"),
label=list("Proportion de malades","Courbe logistique"),
name=" ")
La courbe rouge est celle qui est obtenue par régression logistique.
Si l'on avait voulu considérer l'ensemble des variables médicales (et non pas seulement l'âge comme jusqu'à présent), nous aurions écrit :
reg_log2 <- glm(chd~sbp+tobacco+ldl+adiposity+famhist+typea+obesity+alcohol+age,
family="binomial",data=maladie)
summary(reg_log2)Call:
glm(formula = chd ~ sbp + tobacco + ldl + adiposity + famhist +
typea + obesity + alcohol + age, family = "binomial", data = maladie)
Deviance Residuals:
Min 1Q Median 3Q Max
-1.7781 -0.8213 -0.4387 0.8889 2.5435
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -6.1507209 1.3082600 -4.701 2.58e-06 ***
sbp 0.0065040 0.0057304 1.135 0.256374
tobacco 0.0793764 0.0266028 2.984 0.002847 **
ldl 0.1739239 0.0596617 2.915 0.003555 **
adiposity 0.0185866 0.0292894 0.635 0.525700
famhistPresent 0.9253704 0.2278940 4.061 4.90e-05 ***
typea 0.0395950 0.0123202 3.214 0.001310 **
obesity -0.0629099 0.0442477 -1.422 0.155095
alcohol 0.0001217 0.0044832 0.027 0.978350
age 0.0452253 0.0121298 3.728 0.000193 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
(Dispersion parameter for binomial family taken to be 1)
Null deviance: 596.11 on 461 degrees of freedom
Residual deviance: 472.14 on 452 degrees of freedom
AIC: 492.14
Number of Fisher Scoring iterations: 5
Certaines des variables obtenues ont des p-valeurs qui sont inférieures au niveau de test de 5 %, ce qui nous indique qu'elles sont bien significatives. Certaines autres ne sont pas en dessous de ce seuil.
On peut donc passer sur une procédure de sélection en retirant les variables non significatives au fur et à mesure, mais nous pouvons aussi sélectionner automatiquement un modèle avec une commande telle que stepAIC , qui sélectionne de manière automatique un modèle en se basant sur le critère AIC.
library(MASS)
stepAIC(reg_log2)Et voilà, vous avez appliqué une régression logistique permettant de traiter des variables qualitatives binaires (avec deux modalités). Sachez qu'il existe d'autres méthodes de classification permettant de traiter des variables quantitatives avec davantage de modalités !
Plus qu'un quiz, et vous arriverez à la cinquième partie de ce cours, qui porte sur l'analyse de la variance, l'ANOVA. Je vous retrouve là-bas !
.
