Ce document est un guide pour réaliser un score avec R à travers un cas pratique
Il ne revient pas sur les aspects mathématiques de la construction d'un score, mais insite plutôt sur les commandes en R et sur l'exploitation des résultats.
Il se subdivise en quatre parties :
1 - Description de l'objectif de l'étude et des données
2 - Préparation des données et premières analyses
3 - Construction et validation du score
4 - Interprétation des résultats
Dans une ultime partie, nous donnons quelques recommandations de lecture et les principaux documents qui nous ont permis de construire cette fiche.
Le scoring est une technique de hiérarchisation des données qui permet d'évaluer par une note ou un score la probabilité qu'un individu réponde à une sollicitation ou appartienne à la cible recherchée.
Le score est en général obtenu à partir des données quantitatives et qualitatives disponibles sur l'individu (données socio-demo, comportement d'achat, réponses précédentes, …) auxquelles sont appliquées un modèle de scoring.
De manière générale, la technique de modélisation utilisée est la régression logistique. elle fait partie des techniques d'apprentissage supervisée, c'est à dire que l'on souhaite en général expliquer l'appartenance à une catégorie à partir de descripteurs receuillis sur un échantillon de population dans le but de généraliser l'apprentissage
Quelques exemples d'applications :
On s'intéresse pour notre étude de cas à une base de données qui contient des données 462 patients pour lesquels on souhaite prédire l'exposition à un infarctus.
La base est disponible ici
C'est d'ailleurs un site où des données sont disponibles afin de s'entraîner sur la modélisation en employant diverses techniques.
setwd("D:/AID/Dossiers Internes/Tutoriels/Logistic Regression")
getwd()
## [1] "D:/AID/Dossiers Internes/Tutoriels/Logistic Regression"
HeartDisease <- read.csv(file = "saheartdisease.csv", header = TRUE, sep = ",")
head(HeartDisease)
## row.names sbp tobacco ldl adiposity famhist typea obesity alcohol age
## 1 1 160 12.00 5.73 23.11 Present 49 25.30 97.20 52
## 2 2 144 0.01 4.41 28.61 Absent 55 28.87 2.06 63
## 3 3 118 0.08 3.48 32.28 Present 52 29.14 3.81 46
## 4 4 170 7.50 6.41 38.03 Present 51 31.99 24.26 58
## 5 5 134 13.60 3.50 27.78 Present 60 25.99 57.34 49
## 6 6 132 6.20 6.47 36.21 Present 62 30.77 14.14 45
## chd
## 1 1
## 2 1
## 3 0
## 4 1
## 5 1
## 6 0
str(HeartDisease)
## 'data.frame': 462 obs. of 11 variables:
## $ row.names: int 1 2 3 4 5 6 7 8 9 10 ...
## $ sbp : int 160 144 118 170 134 132 142 114 114 132 ...
## $ tobacco : num 12 0.01 0.08 7.5 13.6 6.2 4.05 4.08 0 0 ...
## $ ldl : num 5.73 4.41 3.48 6.41 3.5 6.47 3.38 4.59 3.83 5.8 ...
## $ adiposity: num 23.1 28.6 32.3 38 27.8 ...
## $ famhist : Factor w/ 2 levels "Absent","Present": 2 1 2 2 2 2 1 2 2 2 ...
## $ typea : int 49 55 52 51 60 62 59 62 49 69 ...
## $ obesity : num 25.3 28.9 29.1 32 26 ...
## $ alcohol : num 97.2 2.06 3.81 24.26 57.34 ...
## $ age : int 52 63 46 58 49 45 38 58 29 53 ...
## $ chd : int 1 1 0 1 1 0 0 1 0 1 ...
Je ne fais aucun commentaires sur les commandes R que j'utilise, sauf quelques exceptions, je renvoie les débutants à quelques notes de cours pour le logiciels R, ici ou encore ici ou enfin là qui en plus de présenter les codes R est une bonne introduction à la modélisation en statistique.
Pour notre modélisation, on s'intéressera à cette dernière variable.
On souhaite déterminer la probabilité d'avoir une maladie du coeur compte tenu des descripteurs relevés sur les 462 patients chd = f(obesite, age, antécédents familiaux, etc… )
Cet exemple illustre bien l'utilisation que l'on peut faire de la régression logistique, la formulation du problème est en général la même, cela peut être pour une compagnie d'assurance, de déterminer les facteurs à risque afin de provisionner (tarification), déterminer les profils de clients appétents à une nouvelle offre commerciale, et en macroéconomie aussi, cette approche est utilisée pour quantifier le risque pays.
Comme toute bonne démarche de modélisation, la construction d'un bon score se fait par une succession d'étapes plus ou moins fondamentales en fonction des praticiens. Néanmoins tous s'accordent plus ou moins à respecter les suivantes :
Nous suivrons ces étapes et résoudrons le problème posé par notre étude de cas.
attach(HeartDisease)
names(HeartDisease)
## [1] "row.names" "sbp" "tobacco" "ldl" "adiposity"
## [6] "famhist" "typea" "obesity" "alcohol" "age"
## [11] "chd"
base = subset(HeartDisease, select = c(chd, age, tobacco, sbp, famhist, alcohol,
obesity, ldl, adiposity, typea))
head(base)
## chd age tobacco sbp famhist alcohol obesity ldl adiposity typea
## 1 1 52 12.00 160 Present 97.20 25.30 5.73 23.11 49
## 2 1 63 0.01 144 Absent 2.06 28.87 4.41 28.61 55
## 3 0 46 0.08 118 Present 3.81 29.14 3.48 32.28 52
## 4 1 58 7.50 170 Present 24.26 31.99 6.41 38.03 51
## 5 1 49 13.60 134 Present 57.34 25.99 3.50 27.78 60
## 6 0 45 6.20 132 Present 14.14 30.77 6.47 36.21 62
tail(base)
## chd age tobacco sbp famhist alcohol obesity ldl adiposity typea
## 457 0 57 0.4 170 Present 2.06 33.10 4.11 42.06 56
## 458 0 58 0.4 214 Absent 0.00 28.45 5.98 31.72 64
## 459 1 52 4.2 182 Absent 18.72 28.61 4.41 32.10 52
## 460 0 55 3.0 108 Absent 26.64 20.09 1.59 15.23 40
## 461 0 40 5.4 118 Absent 23.97 27.35 11.61 30.79 64
## 462 1 46 0.0 132 Present 0.00 14.70 4.82 33.41 62
str(base)
## 'data.frame': 462 obs. of 10 variables:
## $ chd : int 1 1 0 1 1 0 0 1 0 1 ...
## $ age : int 52 63 46 58 49 45 38 58 29 53 ...
## $ tobacco : num 12 0.01 0.08 7.5 13.6 6.2 4.05 4.08 0 0 ...
## $ sbp : int 160 144 118 170 134 132 142 114 114 132 ...
## $ famhist : Factor w/ 2 levels "Absent","Present": 2 1 2 2 2 2 1 2 2 2 ...
## $ alcohol : num 97.2 2.06 3.81 24.26 57.34 ...
## $ obesity : num 25.3 28.9 29.1 32 26 ...
## $ ldl : num 5.73 4.41 3.48 6.41 3.5 6.47 3.38 4.59 3.83 5.8 ...
## $ adiposity: num 23.1 28.6 32.3 38 27.8 ...
## $ typea : int 49 55 52 51 60 62 59 62 49 69 ...
summary(base)
## chd age tobacco sbp
## Min. :0.000 Min. :15.0 Min. : 0.000 Min. :101
## 1st Qu.:0.000 1st Qu.:31.0 1st Qu.: 0.052 1st Qu.:124
## Median :0.000 Median :45.0 Median : 2.000 Median :134
## Mean :0.346 Mean :42.8 Mean : 3.636 Mean :138
## 3rd Qu.:1.000 3rd Qu.:55.0 3rd Qu.: 5.500 3rd Qu.:148
## Max. :1.000 Max. :64.0 Max. :31.200 Max. :218
## famhist alcohol obesity ldl
## Absent :270 Min. : 0.00 Min. :14.7 Min. : 0.98
## Present:192 1st Qu.: 0.51 1st Qu.:23.0 1st Qu.: 3.28
## Median : 7.51 Median :25.8 Median : 4.34
## Mean : 17.04 Mean :26.0 Mean : 4.74
## 3rd Qu.: 23.89 3rd Qu.:28.5 3rd Qu.: 5.79
## Max. :147.19 Max. :46.6 Max. :15.33
## adiposity typea
## Min. : 6.74 Min. :13.0
## 1st Qu.:19.77 1st Qu.:47.0
## Median :26.11 Median :53.0
## Mean :25.41 Mean :53.1
## 3rd Qu.:31.23 3rd Qu.:60.0
## Max. :42.49 Max. :78.0
On remarque que la variabe cible chd a été traitée comme une variable numérique : On rémédie à cela en la transformant en type facteur
base$chd <- factor(base$chd)
summary(base)
## chd age tobacco sbp famhist
## 0:302 Min. :15.0 Min. : 0.000 Min. :101 Absent :270
## 1:160 1st Qu.:31.0 1st Qu.: 0.052 1st Qu.:124 Present:192
## Median :45.0 Median : 2.000 Median :134
## Mean :42.8 Mean : 3.636 Mean :138
## 3rd Qu.:55.0 3rd Qu.: 5.500 3rd Qu.:148
## Max. :64.0 Max. :31.200 Max. :218
## alcohol obesity ldl adiposity
## Min. : 0.00 Min. :14.7 Min. : 0.98 Min. : 6.74
## 1st Qu.: 0.51 1st Qu.:23.0 1st Qu.: 3.28 1st Qu.:19.77
## Median : 7.51 Median :25.8 Median : 4.34 Median :26.11
## Mean : 17.04 Mean :26.0 Mean : 4.74 Mean :25.41
## 3rd Qu.: 23.89 3rd Qu.:28.5 3rd Qu.: 5.79 3rd Qu.:31.23
## Max. :147.19 Max. :46.6 Max. :15.33 Max. :42.49
## typea
## Min. :13.0
## 1st Qu.:47.0
## Median :53.0
## Mean :53.1
## 3rd Qu.:60.0
## Max. :78.0
On va essayer d'identifier les corrélations éventuelles entre les descripteurs
Mais avant cela : Nous allons supprimer toutes les valeurs manquantes possibles du jeu de données
base <- na.omit(base)
Analysons les distributions
attach(base)
## The following object(s) are masked from 'HeartDisease':
##
## adiposity, age, alcohol, chd, famhist, ldl, obesity, sbp,
## tobacco, typea
par(mfrow = c(3, 2))
hist(x = age, col = "lightblue", main = "Age", xlab = "", ylab = "")
hist(x = sbp, col = "orange", main = "pression sanguine systolique", xlab = "",
ylab = "")
hist(x = alcohol, col = "red", main = "Consommation cour. d'alcool", xlab = "",
ylab = "")
hist(x = tobacco, col = "violet", main = "Qté de tabac cumulée", xlab = "",
ylab = "")
hist(x = ldl, col = "green2", main = "taux de cholesterol", xlab = "", ylab = "")
hist(x = obesity, col = "slategray", main = "Obésité", xlab = "", ylab = "")
Attention, lorsqu'on mène une analyse graphique, c'est dans le but de détecter d'éventuelles colinérités, ou au moins d'en avoir quelques idées. Les variables conso d'alcool et Qté de tabac semblent être distribuées de la même façon, de même que le taux de cholestérol et l'obésité.
Un autre outil d'analyse est aussi de réaliser des nuages de points pour toutes les vairiables. On peut éventuellement colorer les points selon la variable cible
par(mfrow = c(1, 1))
pairs(base, col = base$chd)
De plus jolis graphiques (en cas de présentation client par exemple) sont possibles à travers le package ggplot
par exemple :
par(mfrow = c(1, 2))
require(ggplot2)
## Loading required package: ggplot2
qplot(alcohol, tobacco, data = base, facets = chd ~ .)
Les valeurs aberrantes dépendent de la distribution et surtout de l'objet de l'étude. Dans la littérature, le traitement des valeurs manquantes ou aberrantes fait parfois l'objet d'interminables discussions dont le praticien devrait tenir compte. Mais, pris par le temps du projet, il est plus simple en général de décider ce qui est aberrant en connaissance du métier et de la problématique. Regardons encore une fois les statistiques de base
summary(base)
## chd age tobacco sbp famhist
## 0:302 Min. :15.0 Min. : 0.000 Min. :101 Absent :270
## 1:160 1st Qu.:31.0 1st Qu.: 0.052 1st Qu.:124 Present:192
## Median :45.0 Median : 2.000 Median :134
## Mean :42.8 Mean : 3.636 Mean :138
## 3rd Qu.:55.0 3rd Qu.: 5.500 3rd Qu.:148
## Max. :64.0 Max. :31.200 Max. :218
## alcohol obesity ldl adiposity
## Min. : 0.00 Min. :14.7 Min. : 0.98 Min. : 6.74
## 1st Qu.: 0.51 1st Qu.:23.0 1st Qu.: 3.28 1st Qu.:19.77
## Median : 7.51 Median :25.8 Median : 4.34 Median :26.11
## Mean : 17.04 Mean :26.0 Mean : 4.74 Mean :25.41
## 3rd Qu.: 23.89 3rd Qu.:28.5 3rd Qu.: 5.79 3rd Qu.:31.23
## Max. :147.19 Max. :46.6 Max. :15.33 Max. :42.49
## typea
## Min. :13.0
## 1st Qu.:47.0
## Median :53.0
## Mean :53.1
## 3rd Qu.:60.0
## Max. :78.0
la distribution de la consommation du tabac est très étalée, celle de l'alcool aussi. Les autres distributions semblent plutot cohérentes.
Pour l'instant, on ne fait rien sur ces valeurs que l'on considère à priori comme aberrante compte tenu de la distribution.
C'est une question fréquente en analyse exploratoire. La discrétisation des variables continues.
Dans le jeu de données, les variables :
Doit-on discrétiser les variables continues? Oui, la plupart du temps. Mais Comment? en lien avec la variable cible? Par connaissance métier? Distribution basée sur les quantiles?
Aucune réponse définitive. D'un point de vue général le choix de la méthode dépend en général du problème, du temps que l'on souhaite y consacrer. Toujours garder en mémoire : Aucun découpage n'est bon à priori et en fonction des résultats pratiques, ne pas hésiter à revenir sur son découpage.
la variable age est en général la plus simple à découper. les problèmes cardiaques n'affectent pas de la même façon selon les âges, comme indiqué ici
On fait le choix de discrétiser les variables age et tobacco en faisant une distinction entre petits, moyens et gros fumeur
detach(HeartDisease)
attach(base)
## The following object(s) are masked from 'base (position 4)':
##
## adiposity, age, alcohol, chd, famhist, ldl, obesity, sbp,
## tobacco, typea
Breaksage = c(0, 15, 25, 45, max(age))
age.d = cut(age, breaks = Breaksage, include.lowest = TRUE)
summary(age.d)
## [0,15] (15,25] (25,45] (45,64]
## 3 68 170 221
Breakstob = c(0, 0.05, 2, max(tobacco))
tobacco.d = cut(tobacco, breaks = Breakstob, include.lowest = TRUE)
summary(tobacco.d)
## [0,0.05] (0.05,2] (2,31.2]
## 116 118 228
Un découpage n'est pertinent que lorsque les modalités de la variable cible y sont correctement représentées
base2 <- cbind(base, age.d, tobacco.d)
detach(base)
attach(base2)
## The following object(s) are masked _by_ '.GlobalEnv':
##
## age.d, tobacco.d
## The following object(s) are masked from 'base':
##
## adiposity, age, alcohol, chd, famhist, ldl, obesity, sbp,
## tobacco, typea
xtabs(~chd + age.d, data = base2)
## age.d
## chd [0,15] (15,25] (25,45] (45,64]
## 0 3 65 123 111
## 1 0 3 47 110
xtabs(~chd + famhist, data = base2)
## famhist
## chd Absent Present
## 0 206 96
## 1 64 96
xtabs(~chd + tobacco.d, data = base2)
## tobacco.d
## chd [0,0.05] (0.05,2] (2,31.2]
## 0 100 83 119
## 1 16 35 109
La catégorie de moins de 15 ans n'est pas du tout représentative dans l'échantillon. les moins de 25ans non plus.
La moitié des personnes de plus de 45 ans sont atteintes d'un problème cardiaque.
On peut penser à une forme d'hérédité dans les problèmes cardiaques au vu des croisements
Il est clair que le tabagisme a une vraie influence sur les problèmes cardiaques, car on retrouve une proportion non négligeable de personnes atteintes quelque soit leur quantité de tabac consommé.
Ces premières analyses nous indiquent que :
base3 <- subset(x = base2, subset = (age > 15))
detach(base2)
attach(base3)
## The following object(s) are masked _by_ '.GlobalEnv':
##
## age.d, tobacco.d
## The following object(s) are masked from 'base':
##
## adiposity, age, alcohol, chd, famhist, ldl, obesity, sbp,
## tobacco, typea
# On enlève dans la base de travail les variables age et tobacco
# continues, au profit des discrétisés
base3 <- subset(base3, select = -c(age, tobacco))
Nous maintenant échantilloner notre base afin d'avoir un échantillon d'appretnissage et un échantillon test
# Tirage aléaoire et sans remise des 65% des individus de l'échantillon On
# initialise le tirage aléatoire afin de retomber sur nos pieds à chaque
# fois
set.seed(111)
d = sort(sample(nrow(base3), nrow(base3) * 0.65))
# Echantillon d'apprentissage
appren <- base3[d, ]
# Echantillon de test
test <- base3[-d, ]
summary(appren)
## chd sbp famhist alcohol obesity
## 0:183 Min. :101 Absent :165 Min. : 0.00 Min. :14.7
## 1:115 1st Qu.:124 Present:133 1st Qu.: 0.70 1st Qu.:23.2
## Median :134 Median : 8.33 Median :26.0
## Mean :140 Mean : 18.77 Mean :26.2
## 3rd Qu.:152 3rd Qu.: 24.84 3rd Qu.:28.7
## Max. :218 Max. :147.19 Max. :46.6
## ldl adiposity typea age.d
## Min. : 0.98 Min. : 7.12 Min. :13.0 [0,15] : 0
## 1st Qu.: 3.26 1st Qu.:21.07 1st Qu.:47.0 (15,25]: 38
## Median : 4.37 Median :26.46 Median :53.5 (25,45]:106
## Mean : 4.73 Mean :25.79 Mean :53.3 (45,64]:154
## 3rd Qu.: 5.81 3rd Qu.:31.68 3rd Qu.:60.8
## Max. :14.16 Max. :42.06 Max. :75.0
## tobacco.d
## [0,0.05]: 69
## (0.05,2]: 83
## (2,31.2]:146
##
##
##
summary(test)
## chd sbp famhist alcohol obesity
## 0:116 Min. :103 Absent :102 Min. : 0.0 Min. :17.8
## 1: 45 1st Qu.:124 Present: 59 1st Qu.: 0.0 1st Qu.:22.9
## Median :134 Median : 4.8 Median :25.6
## Mean :136 Mean :14.2 Mean :25.9
## 3rd Qu.:144 3rd Qu.:18.7 3rd Qu.:28.3
## Max. :198 Max. :92.6 Max. :45.7
## ldl adiposity typea age.d tobacco.d
## Min. : 1.43 Min. : 6.74 Min. :29.0 [0,15] : 0 [0,0.05]:44
## 1st Qu.: 3.30 1st Qu.:18.96 1st Qu.:46.0 (15,25]:30 (0.05,2]:35
## Median : 4.24 Median :25.84 Median :53.0 (25,45]:64 (2,31.2]:82
## Mean : 4.81 Mean :24.89 Mean :52.8 (45,64]:67
## 3rd Qu.: 5.76 3rd Qu.:30.28 3rd Qu.:59.0
## Max. :15.33 Max. :42.49 Max. :78.0
Avant de se lancer dans la construction des modèles, décrivons très rapidement de quoi est constitué notre échantillon d'apprentissge
attach(appren)
## The following object(s) are masked _by_ '.GlobalEnv':
##
## age.d, tobacco.d
## The following object(s) are masked from 'base3':
##
## adiposity, age.d, alcohol, chd, famhist, ldl, obesity, sbp,
## tobacco.d, typea
## The following object(s) are masked from 'base':
##
## adiposity, alcohol, chd, famhist, ldl, obesity, sbp, typea
str(appren)
## 'data.frame': 298 obs. of 10 variables:
## $ chd : Factor w/ 2 levels "0","1": 2 1 2 2 2 2 2 1 1 1 ...
## $ sbp : int 160 118 170 134 114 206 134 118 112 117 ...
## $ famhist : Factor w/ 2 levels "Absent","Present": 2 2 2 2 2 1 2 1 2 2 ...
## $ alcohol : num 97.2 3.81 24.26 57.34 6.72 ...
## $ obesity : num 25.3 29.1 32 26 23.1 ...
## $ ldl : num 5.73 3.48 6.41 3.5 4.59 2.95 4.44 1.88 2.29 2.44 ...
## $ adiposity: num 23.1 32.3 38 27.8 14.6 ...
## $ typea : int 49 52 51 60 62 72 65 59 54 35 ...
## $ age.d : Factor w/ 4 levels "[0,15]","(15,25]",..: 4 4 4 4 4 4 3 2 4 4 ...
## $ tobacco.d: Factor w/ 3 levels "[0,0.05]","(0.05,2]",..: 3 2 3 3 3 3 3 1 3 2 ...
On souhaite modéliser chd = f(variables explicatives)
Dans la vraie vie, on dispose de beaucoup plus de variables et il faut donc avoir recours à un critère qui sélectionne les variables à introduire dans le modèle. Nous le faisons ici à titre purement illustratif, on pourrait s'en passer dans cet exemple. Nous laissons aussi de côté les analyses des interactions entre variables qui constituent une étape à part entière de l'analyse car la lecture des résultats ne se fait pas de la même façon
Quelles variables introduire?
Pourquoi sélectionner? : Pour plusieurs raisons, voir par exemple le cours de Ricco Rakotomalala ici
Mais si on ne devait retenir qu'une, ce serait la suivante : ** Un modèle avec peu de variables sera plus facilement généralisable en terme de robustesse (principe du rasoir d'Occam)**
Plusieurs façons de sélectionner :
# modèle trivial réduit à la constante
str_constant <- "~ 1"
# modèle complet incluant toutes les explicatives potentielles
str_all <- "~sbp+famhist+alcohol+obesity+ldl+adiposity+typea+age.d+tobacco.d"
require(MASS)
## Loading required package: MASS
modele <- glm(chd ~ 1, data = appren, family = binomial)
modele.forward <- stepAIC(modele, scope = list(lower = str_constant, upper = str_all),
trace = TRUE, data = appren, direction = "forward")
## Start: AIC=399.5
## chd ~ 1
##
## Df Deviance AIC
## + tobacco.d 2 362 368
## + age.d 2 364 370
## + famhist 1 373 377
## + ldl 1 377 381
## + adiposity 1 384 388
## + sbp 1 386 390
## + typea 1 394 398
## + alcohol 1 395 399
## <none> 397 399
## + obesity 1 397 401
##
## Step: AIC=367.7
## chd ~ tobacco.d
##
## Df Deviance AIC
## + famhist 1 345 353
## + age.d 2 348 358
## + ldl 1 351 359
## + sbp 1 355 363
## + adiposity 1 356 364
## + typea 1 356 364
## <none> 362 368
## + obesity 1 362 370
## + alcohol 1 362 370
##
## Step: AIC=352.7
## chd ~ tobacco.d + famhist
##
## Df Deviance AIC
## + age.d 2 333 345
## + ldl 1 337 347
## + sbp 1 339 349
## + typea 1 340 350
## + adiposity 1 340 350
## <none> 345 353
## + obesity 1 345 355
## + alcohol 1 345 355
##
## Step: AIC=344.9
## chd ~ tobacco.d + famhist + age.d
##
## Df Deviance AIC
## + typea 1 327 341
## + ldl 1 327 341
## + sbp 1 330 344
## <none> 333 345
## + obesity 1 332 346
## + adiposity 1 333 347
## + alcohol 1 333 347
##
## Step: AIC=340.8
## chd ~ tobacco.d + famhist + age.d + typea
##
## Df Deviance AIC
## + ldl 1 322 338
## + sbp 1 324 340
## <none> 327 341
## + obesity 1 326 342
## + adiposity 1 327 343
## + alcohol 1 327 343
##
## Step: AIC=337.8
## chd ~ tobacco.d + famhist + age.d + typea + ldl
##
## Df Deviance AIC
## + obesity 1 318 336
## + sbp 1 319 337
## <none> 322 338
## + adiposity 1 322 340
## + alcohol 1 322 340
##
## Step: AIC=336.3
## chd ~ tobacco.d + famhist + age.d + typea + ldl + obesity
##
## Df Deviance AIC
## + sbp 1 315 335
## <none> 318 336
## + adiposity 1 317 337
## + alcohol 1 318 338
##
## Step: AIC=335.2
## chd ~ tobacco.d + famhist + age.d + typea + ldl + obesity + sbp
##
## Df Deviance AIC
## <none> 315 335
## + adiposity 1 314 336
## + alcohol 1 315 337
# affichage du modèle final
summary(modele.forward)
##
## Call:
## glm(formula = chd ~ tobacco.d + famhist + age.d + typea + ldl +
## obesity + sbp, family = binomial, data = appren)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.837 -0.866 -0.397 0.910 2.775
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -6.08323 1.52932 -3.98 7e-05 ***
## tobacco.d(0.05,2] 0.62751 0.46053 1.36 0.17301
## tobacco.d(2,31.2] 1.21942 0.42779 2.85 0.00437 **
## famhistPresent 0.98254 0.27518 3.57 0.00036 ***
## age.d(25,45] 1.55762 0.81618 1.91 0.05634 .
## age.d(45,64] 2.04607 0.82368 2.48 0.01299 *
## typea 0.03712 0.01463 2.54 0.01117 *
## ldl 0.20200 0.07751 2.61 0.00916 **
## obesity -0.07447 0.03704 -2.01 0.04437 *
## sbp 0.01124 0.00645 1.74 0.08118 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 397.46 on 297 degrees of freedom
## Residual deviance: 315.18 on 288 degrees of freedom
## AIC: 335.2
##
## Number of Fisher Scoring iterations: 5
Analysons l'algorithmique faite par R.
Dans le modèle de base,
AIC = 367
Il rajoute la consommation de tabac en variable explicative et calcule l'AIC par ordre croissant des variables introduites: AIC= 367; c'est le meilleur modèle à ce stade
On continue :
Il teste tous les modèles à deux variables sans remettre en cause la consommation de tabac puisque déjà sélectionné :il retient l'antécédant familial
etc…
Le modèle qu'il propose au final
Call:
glm(formula = chd ~ tobacco.d + famhist + age.d + typea + ldl + obesity + sbp, family = binomial, data = appren) retient toutes variables à l'exeption de la consommation d'alcool sans doute très correlé avec la consommation du tabac comme on l'avait décelé dans l'analyse descriptive au début de notre propos.
modele <- glm(chd ~ 1, data = appren, family = binomial)
modele.stepwise <- stepAIC(modele, scope = list(lower = str_constant, upper = str_all),
trace = TRUE, data = appren, direction = "both")
## Start: AIC=399.5
## chd ~ 1
##
## Df Deviance AIC
## + tobacco.d 2 362 368
## + age.d 2 364 370
## + famhist 1 373 377
## + ldl 1 377 381
## + adiposity 1 384 388
## + sbp 1 386 390
## + typea 1 394 398
## + alcohol 1 395 399
## <none> 397 399
## + obesity 1 397 401
##
## Step: AIC=367.7
## chd ~ tobacco.d
##
## Df Deviance AIC
## + famhist 1 345 353
## + age.d 2 348 358
## + ldl 1 351 359
## + sbp 1 355 363
## + adiposity 1 356 364
## + typea 1 356 364
## <none> 362 368
## + obesity 1 362 370
## + alcohol 1 362 370
## - tobacco.d 2 397 399
##
## Step: AIC=352.7
## chd ~ tobacco.d + famhist
##
## Df Deviance AIC
## + age.d 2 333 345
## + ldl 1 337 347
## + sbp 1 339 349
## + typea 1 340 350
## + adiposity 1 340 350
## <none> 345 353
## + obesity 1 345 355
## + alcohol 1 345 355
## - famhist 1 362 368
## - tobacco.d 2 373 377
##
## Step: AIC=344.9
## chd ~ tobacco.d + famhist + age.d
##
## Df Deviance AIC
## + typea 1 327 341
## + ldl 1 327 341
## + sbp 1 330 344
## <none> 333 345
## + obesity 1 332 346
## + adiposity 1 333 347
## + alcohol 1 333 347
## - age.d 2 345 353
## - tobacco.d 2 346 354
## - famhist 1 348 358
##
## Step: AIC=340.8
## chd ~ tobacco.d + famhist + age.d + typea
##
## Df Deviance AIC
## + ldl 1 322 338
## + sbp 1 324 340
## <none> 327 341
## + obesity 1 326 342
## + adiposity 1 327 343
## + alcohol 1 327 343
## - typea 1 333 345
## - age.d 2 340 350
## - tobacco.d 2 341 351
## - famhist 1 341 353
##
## Step: AIC=337.8
## chd ~ tobacco.d + famhist + age.d + typea + ldl
##
## Df Deviance AIC
## + obesity 1 318 336
## + sbp 1 319 337
## <none> 322 338
## + adiposity 1 322 340
## + alcohol 1 322 340
## - ldl 1 327 341
## - typea 1 327 341
## - age.d 2 332 344
## - tobacco.d 2 333 345
## - famhist 1 335 349
##
## Step: AIC=336.3
## chd ~ tobacco.d + famhist + age.d + typea + ldl + obesity
##
## Df Deviance AIC
## + sbp 1 315 335
## <none> 318 336
## + adiposity 1 317 337
## - obesity 1 322 338
## + alcohol 1 318 338
## - typea 1 325 341
## - ldl 1 326 342
## - tobacco.d 2 329 343
## - age.d 2 331 345
## - famhist 1 331 347
##
## Step: AIC=335.2
## chd ~ tobacco.d + famhist + age.d + typea + ldl + obesity + sbp
##
## Df Deviance AIC
## <none> 315 335
## + adiposity 1 314 336
## - sbp 1 318 336
## + alcohol 1 315 337
## - obesity 1 319 337
## - typea 1 322 340
## - age.d 2 324 340
## - ldl 1 322 340
## - tobacco.d 2 325 341
## - famhist 1 328 346
# affichage du modèle final
summary(modele.stepwise)
##
## Call:
## glm(formula = chd ~ tobacco.d + famhist + age.d + typea + ldl +
## obesity + sbp, family = binomial, data = appren)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.837 -0.866 -0.397 0.910 2.775
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -6.08323 1.52932 -3.98 7e-05 ***
## tobacco.d(0.05,2] 0.62751 0.46053 1.36 0.17301
## tobacco.d(2,31.2] 1.21942 0.42779 2.85 0.00437 **
## famhistPresent 0.98254 0.27518 3.57 0.00036 ***
## age.d(25,45] 1.55762 0.81618 1.91 0.05634 .
## age.d(45,64] 2.04607 0.82368 2.48 0.01299 *
## typea 0.03712 0.01463 2.54 0.01117 *
## ldl 0.20200 0.07751 2.61 0.00916 **
## obesity -0.07447 0.03704 -2.01 0.04437 *
## sbp 0.01124 0.00645 1.74 0.08118 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 397.46 on 297 degrees of freedom
## Residual deviance: 315.18 on 288 degrees of freedom
## AIC: 335.2
##
## Number of Fisher Scoring iterations: 5
Expliciter la démarche comme on l'a fait pour Backward serait un peu compliqué, alors, regardons la sélection au final faite par l'algoritme :
glm(formula = chd ~ tobacco.d + famhist + age.d + typea + ldl + obesity + sbp, family = binomial, data = appren) Soient, les mêmes variables que dans la sélection forward.
Avant de continuer, je propose de construire nous même notre fonction pour réaliser une régression logistique sous R. Ceci uniquement dans le bit d'alléger l'écriture de nos modèles.La fonction utilisée pour la régression logistique est simplement la fonction glm de R.
logit = function(formula, lien = "logit", data = NULL) {
glm(formula, family = binomial(link = lien), data)
}
On introduit toutes les variables dans notre modèle sauf la conso d'alcool :
m.logit <- logit(chd ~ sbp + famhist + obesity + ldl + adiposity + typea + age.d +
tobacco.d, data = appren)
# résultats du modèle
summary(m.logit)
##
## Call:
## glm(formula = formula, family = binomial(link = lien), data = data)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.807 -0.865 -0.389 0.899 2.769
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -5.73451 1.57212 -3.65 0.00026 ***
## sbp 0.01032 0.00649 1.59 0.11222
## famhistPresent 0.97583 0.27557 3.54 0.00040 ***
## obesity -0.11060 0.05263 -2.10 0.03561 *
## ldl 0.18719 0.07936 2.36 0.01834 *
## adiposity 0.03376 0.03401 0.99 0.32095
## typea 0.03858 0.01476 2.61 0.00897 **
## age.d(25,45] 1.41550 0.82662 1.71 0.08682 .
## age.d(45,64] 1.78113 0.86283 2.06 0.03899 *
## tobacco.d(0.05,2] 0.66933 0.46699 1.43 0.15177
## tobacco.d(2,31.2] 1.25793 0.43401 2.90 0.00375 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 397.46 on 297 degrees of freedom
## Residual deviance: 314.19 on 287 degrees of freedom
## AIC: 336.2
##
## Number of Fisher Scoring iterations: 5
Par défaut, R ne donne pas les odds-ratios. Mais, il se calcule aisément en se rappelant qu'un odd-ratio n'est rien de plus que l'estiamtion des coefs dela régression. On peut donc obtenir les Odd-ratio et leurs intervalles de confiance en faisant
exp(cbind(OR = coef(m.logit), confint(m.logit)))
## Waiting for profiling to be done...
## OR 2.5 % 97.5 %
## (Intercept) 0.003232 0.000128 0.06374
## sbp 1.010369 0.997704 1.02357
## famhistPresent 2.653379 1.552120 4.58203
## obesity 0.895300 0.804319 0.98983
## ldl 1.205851 1.036067 1.41574
## adiposity 1.034335 0.968189 1.10680
## typea 1.039331 1.010404 1.07081
## age.d(25,45] 4.118525 0.973345 28.85054
## age.d(45,64] 5.936559 1.289700 43.73757
## tobacco.d(0.05,2] 1.952935 0.798525 5.05115
## tobacco.d(2,31.2] 3.518125 1.546387 8.59172
Note :
summary(m.logit)
##
## Call:
## glm(formula = formula, family = binomial(link = lien), data = data)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.807 -0.865 -0.389 0.899 2.769
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -5.73451 1.57212 -3.65 0.00026 ***
## sbp 0.01032 0.00649 1.59 0.11222
## famhistPresent 0.97583 0.27557 3.54 0.00040 ***
## obesity -0.11060 0.05263 -2.10 0.03561 *
## ldl 0.18719 0.07936 2.36 0.01834 *
## adiposity 0.03376 0.03401 0.99 0.32095
## typea 0.03858 0.01476 2.61 0.00897 **
## age.d(25,45] 1.41550 0.82662 1.71 0.08682 .
## age.d(45,64] 1.78113 0.86283 2.06 0.03899 *
## tobacco.d(0.05,2] 0.66933 0.46699 1.43 0.15177
## tobacco.d(2,31.2] 1.25793 0.43401 2.90 0.00375 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 397.46 on 297 degrees of freedom
## Residual deviance: 314.19 on 287 degrees of freedom
## AIC: 336.2
##
## Number of Fisher Scoring iterations: 5
# On peur aussi afficher les différents résultats disponibles dans l'objet
# m.logit par
attributes(m.logit)
## $names
## [1] "coefficients" "residuals" "fitted.values"
## [4] "effects" "R" "rank"
## [7] "qr" "family" "linear.predictors"
## [10] "deviance" "aic" "null.deviance"
## [13] "iter" "weights" "prior.weights"
## [16] "df.residual" "df.null" "y"
## [19] "converged" "boundary" "model"
## [22] "call" "formula" "terms"
## [25] "data" "offset" "control"
## [28] "method" "contrasts" "xlevels"
##
## $class
## [1] "glm" "lm"
Après avoir obtenu un modèle, il faut diagnostiquer la régression afin de valider ou non le modèle.
L'analyse des résidus est de ce point de vue très importante. Il est important de noter qu'en régression logistique, on s'intéresse la plupart du temps aux résidus de déviance. Ils prennent généralement les valeurs qui oscillent entre -2 et 2. On construit généralement un index plot pour détecter les valeurs aberrantes (en dehors des lignes)
par(mfrow = c(1, 1))
plot(rstudent(m.logit), type = "p", cex = 0.5, ylab = "Résidus studentisés ",
col = "springgreen2", ylim = c(-3, 3))
abline(h = c(-2, 2), col = "red")
On s'intéresse aussi en général à la déviance du modèle : les test de rapport des vraisemblances et le calcul de la p_value à l'écart de degré de liberté entre le modèle réduit à la cste et le modèle retenu donne la significativité globale du modèle.
(chi2 <- with(m.logit, null.deviance - deviance))
## [1] 83.28
(ddl <- with(m.logit, df.null - df.residual))
## [1] 10
(pvalue <- pchisq(chi2, ddl, lower.tail = F))
## [1] 1.141e-13
Nous allons tenter de valider maintenant sur l'échantillon test que nous avons précédemment défini;
Voici les étapes que nous allons suivre pour valider notre modèle
Sur l'échantillon d'apprentissage et sur l'échantillon test :
Mais avant, pour prévoir avec R, une fonction * predict()* quelque soit le modèle développé et le package utilisé
Soit donc :
appren.p <- cbind(appren, predict(m.logit, newdata = appren, type = "link",
se = TRUE))
head(appren.p)
## chd sbp famhist alcohol obesity ldl adiposity typea age.d tobacco.d
## 1 1 160 Present 97.20 25.30 5.73 23.11 49 (45,64] (2,31.2]
## 3 0 118 Present 3.81 29.14 3.48 32.28 52 (45,64] (0.05,2]
## 4 1 170 Present 24.26 31.99 6.41 38.03 51 (45,64] (2,31.2]
## 5 1 134 Present 57.34 25.99 3.50 27.78 60 (45,64] (2,31.2]
## 8 1 114 Present 6.72 23.11 4.59 14.60 62 (45,64] (2,31.2]
## 11 1 206 Absent 56.06 26.81 2.95 32.27 72 (45,64] (2,31.2]
## fit se.fit residual.scale
## 1 0.8758 0.3215 1
## 3 -0.5666 0.4244 1
## 4 0.9472 0.3324 1
## 5 0.6958 0.3125 1
## 8 0.6443 0.4848 1
## 11 0.8836 0.5821 1
tail(appren.p)
## chd sbp famhist alcohol obesity ldl adiposity typea age.d
## 456 1 128 Absent 47.42 23.96 2.83 26.48 48 (25,45]
## 457 0 170 Present 2.06 33.10 4.11 42.06 56 (45,64]
## 458 0 214 Absent 0.00 28.45 5.98 31.72 64 (45,64]
## 460 0 108 Absent 26.64 20.09 1.59 15.23 40 (45,64]
## 461 0 118 Absent 23.97 27.35 11.61 30.79 64 (25,45]
## 462 1 132 Present 0.00 14.70 4.82 33.41 62 (45,64]
## tobacco.d fit se.fit residual.scale
## 456 (2,31.2] -1.1152 0.3976 1
## 457 (0.05,2] 0.1342 0.4953 1
## 458 (0.05,2] 0.4361 0.5743 1
## 460 (2,31.2] -1.4484 0.5031 1
## 461 (2,31.2] 0.8129 0.6354 1
## 462 [0,0.05] 1.1802 0.8901 1
On peut utiliser les types : link ou response, selon que l'on veut la probabilité ou non.
En tout état de cause on peut toujours obtenir les probabilités avec les fonctions sigmoides correspondance à l'inverse des fit par la fonction logistique.
Soit,
appren.p <- within(appren.p, {
PredictedProb <- plogis(fit)
LL <- plogis(fit - (1.96 * se.fit))
UL <- plogis(fit + (1.96 * se.fit))
})
tail(appren.p)
## chd sbp famhist alcohol obesity ldl adiposity typea age.d
## 456 1 128 Absent 47.42 23.96 2.83 26.48 48 (25,45]
## 457 0 170 Present 2.06 33.10 4.11 42.06 56 (45,64]
## 458 0 214 Absent 0.00 28.45 5.98 31.72 64 (45,64]
## 460 0 108 Absent 26.64 20.09 1.59 15.23 40 (45,64]
## 461 0 118 Absent 23.97 27.35 11.61 30.79 64 (25,45]
## 462 1 132 Present 0.00 14.70 4.82 33.41 62 (45,64]
## tobacco.d fit se.fit residual.scale UL LL PredictedProb
## 456 (2,31.2] -1.1152 0.3976 1 0.4168 0.13074 0.2469
## 457 (0.05,2] 0.1342 0.4953 1 0.7512 0.30227 0.5335
## 458 (0.05,2] 0.4361 0.5743 1 0.8266 0.33413 0.6073
## 460 (2,31.2] -1.4484 0.5031 1 0.3865 0.08057 0.1902
## 461 (2,31.2] 0.8129 0.6354 1 0.8868 0.39354 0.6927
## 462 [0,0.05] 1.1802 0.8901 1 0.9491 0.36252 0.7650
appren.p <- cbind(appren.p, pred.chd = factor(ifelse(appren.p$PredictedProb >
0.5, 1, 0)))
head(appren.p)
## chd sbp famhist alcohol obesity ldl adiposity typea age.d tobacco.d
## 1 1 160 Present 97.20 25.30 5.73 23.11 49 (45,64] (2,31.2]
## 3 0 118 Present 3.81 29.14 3.48 32.28 52 (45,64] (0.05,2]
## 4 1 170 Present 24.26 31.99 6.41 38.03 51 (45,64] (2,31.2]
## 5 1 134 Present 57.34 25.99 3.50 27.78 60 (45,64] (2,31.2]
## 8 1 114 Present 6.72 23.11 4.59 14.60 62 (45,64] (2,31.2]
## 11 1 206 Absent 56.06 26.81 2.95 32.27 72 (45,64] (2,31.2]
## fit se.fit residual.scale UL LL PredictedProb pred.chd
## 1 0.8758 0.3215 1 0.8185 0.5611 0.7059 1
## 3 -0.5666 0.4244 1 0.5659 0.1981 0.3620 0
## 4 0.9472 0.3324 1 0.8318 0.5734 0.7205 1
## 5 0.6958 0.3125 1 0.7872 0.5208 0.6673 1
## 8 0.6443 0.4848 1 0.8312 0.4241 0.6557 1
## 11 0.8836 0.5821 1 0.8834 0.4360 0.7076 1
# Matrice de confusion
(m.confusion <- as.matrix(table(appren.p$pred.chd, appren.p$chd)))
##
## 0 1
## 0 151 46
## 1 32 69
m.confusion <- unclass(m.confusion)
# Taux d'erreur
Tx_err <- function(y, ypred) {
mc <- table(y, ypred)
error <- (mc[1, 2] + mc[2, 1])/sum(mc)
print(error)
}
Tx_err(appren.p$pred.chd, appren.p$chd)
## [1] 0.2617
On réalise la même opération sur l'échantillon test :
test.p <- cbind(test, predict(m.logit, newdata = test, type = "response", se = TRUE))
test.p <- cbind(test.p, pred.chd <- factor(ifelse(test.p$fit > 0.5, 1, 0)))
(m.confusiontest <- as.matrix(table(test.p$pred.chd, test.p$chd)))
##
## 0 1
## 0 95 23
## 1 21 22
m.confusiontest <- unclass(m.confusiontest)
# calcul du taux d'erreur sur l'échantillon test
Tx_err(test.p$pred.chd, test.p$chd)
## [1] 0.2733
Cette courbe, ou plutôt l'aire sous elle, représente la sensibilité/spécifité du modèle. Un modèle est bon si les positifs (les 1) ont été prédit positifs et les 0 ont été prévus 0.
Généralement, on s'intéresse à la fois à la forme de la courbe et à l'aire sous elle : 1–> Modèle idéal, 0.5 –> Modèle aléatoire;
Principe de la courbe ROC : si le test donne un résultat numérique avec un seuil t tel que la prédiction est positive si x > t, et la prédiction est négative si x < t, alors au fur et à mesure que t augmente :
La courbe ROC représente l'évolution de la sensibilité (taux de vrais positifs) en fonction de 1 - spécificité (taux de faux positifs) quand on fait varier le seuil t.
C'est une courbe croissante entre le point (0,0) et le point (1, 1) et en principe au-dessus de la première bissectrice.Une prédiction random donnerait la première bissectrice. Meilleure est la prédiction, plus la courbe est au-dessus la première bissectrice. Une prédiction idéale est l'horizontale y=1 sur ]0,1] et le point (0,0). L'aire sous la courbe ROC (AUC, Area Under the Curve) donne un indicateur de la qualité de la prédiction (1 pour une prédiction idéale, 0.5 pour une prédiction random). Pour une bonne compréhension de cette aire et de son implémentation sous R, je renvoie ici
Une librairie est disponible pour calculer rapidement cette aire
require(ROCR)
## Loading required package: ROCR
## Loading required package: gplots
## Loading required package: gtools
## Attaching package: 'gtools'
## The following object(s) are masked _by_ '.GlobalEnv':
##
## logit
## Loading required package: gdata
## gdata: Unable to locate valid perl interpreter gdata: gdata: read.xls()
## will be unable to read Excel XLS and XLSX files gdata: unless the 'perl='
## argument is used to specify the location gdata: of a valid perl
## intrpreter. gdata: gdata: (To avoid display of this message in the future,
## please gdata: ensure perl is installed and available on the executable
## gdata: search path.)
## gdata: Unable to load perl libaries needed by read.xls() gdata: to support
## 'XLX' (Excel 97-2004) files.
## ```
## gdata: Run the function 'installXLSXsupport()' gdata: to automatically
## download and install the perl gdata: libaries needed to support Excel XLS
## and XLSX formats.
## Attaching package: 'gdata'
## The following object(s) are masked from 'package:stats':
##
## nobs
## The following object(s) are masked from 'package:utils':
##
## object.size
## Loading required package: caTools
## Loading required package: grid
## Loading required package: KernSmooth
## KernSmooth 2.23 loaded Copyright M. P. Wand 1997-2009
## Attaching package: 'gplots'
## The following object(s) are masked from 'package:stats':
##
## lowess
Pred = prediction(appren.p$PredictedProb, appren.p$chd)
Perf = performance(Pred, "tpr", "fpr")
plot(Perf, colorize = TRUE, main = "ROC apprentissage")
Pour avoir l'aire sous la courbe, on utilisera plutôt
perf <- performance(Pred, "auc")
perf@y.values[[1]]
## [1] 0.7946
On va mettre les deux courbes ROC côte à côte (apprenissage et tests)
Predtest = prediction(test.p$fit, test.p$chd)
Perftest = performance(Predtest, "tpr", "fpr")
perftest <- performance(Predtest, "auc")
perftest@y.values[[1]]
## [1] 0.7573
par(mfrow = c(1, 2))
plot(Perf, colorize = TRUE, main = "ROC apprentissage - AUC= 0.8")
plot(Perftest, colorize = TRUE, main = "ROC Test - AUC =0.75 ")
On peut dans un souci de présentation être intéressé par voir où les probabilités sont les plus élevées par rapport à une variable par exemple :
ggplot(appren.p, aes(x = age.d, y = PredictedProb)) + geom_ribbon(aes(ymin = LL,
ymax = UL, fill = age.d), alpha = 0.15) + geom_line(aes(colour = age.d),
size = 1)
A faire
A faire