Importation, modélisation, scores
# Q.0 + 1 - Importation des données
#librairie
library(xlsx)
#charger
donnees <- read.xlsx("Low_Birth_Weight_Data.xlsx",sheetIndex = 1, header=TRUE, stringsAsFactors=TRUE,encoding="UTF-8")
#premieres lignes
print(head(donnees))
## LOW AGE LWT SMOKE HT UI FTV PTL
## 1 y 28 120 1 0 1 0 1
## 2 y 29 130 0 0 1 1 0
## 3 y 34 187 1 1 0 0 0
## 4 y 25 105 0 1 0 0 1
## 5 y 25 85 0 0 1 0 0
## 6 y 27 150 0 0 0 0 0
# description
str(donnees)
## 'data.frame': 189 obs. of 8 variables:
## $ LOW : Factor w/ 2 levels "n","y": 2 2 2 2 2 2 2 2 2 2 ...
## $ AGE : num 28 29 34 25 25 27 23 24 24 21 ...
## $ LWT : num 120 130 187 105 85 150 97 124 132 165 ...
## $ SMOKE: num 1 0 1 0 0 0 0 0 0 1 ...
## $ HT : num 0 0 1 1 0 0 0 0 1 1 ...
## $ UI : num 1 1 0 0 1 0 1 0 0 0 ...
## $ FTV : num 0 1 0 0 0 0 1 1 0 1 ...
## $ PTL : num 1 0 0 1 0 0 0 1 0 0 ...
# Q.2 + 3 - modélisation : LOW vs. l'ensemble des variables explicatives
modele <- glm(LOW ~ ., data = donnees, family = binomial)
# résultats détaillés - 3 coefficients (de variables) présentent une p-value < 0.05
print(summary(modele))
##
## Call:
## glm(formula = LOW ~ ., family = binomial, data = donnees)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.8223 -0.7914 -0.5983 0.9523 2.0628
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1.491472 1.115765 1.337 0.18131
## AGE -0.046676 0.036941 -1.264 0.20639
## LWT -0.013955 0.006766 -2.063 0.03915 *
## SMOKE 0.446009 0.355204 1.256 0.20924
## HT 1.832962 0.699416 2.621 0.00877 **
## UI 0.668032 0.464022 1.440 0.14996
## FTV -0.255064 0.364908 -0.699 0.48456
## PTL 1.336773 0.458652 2.915 0.00356 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 234.67 on 188 degrees of freedom
## Residual deviance: 202.77 on 181 degrees of freedom
## AIC: 218.77
##
## Number of Fisher Scoring iterations: 4
# Q.4 - probabilités d'affectation
scores <- predict(modele, newdata = donnees, type = "response")
print(summary(scores))
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.03686 0.18203 0.26000 0.31217 0.40345 0.88381
Courbe ROC
# Q.6 - coder la cibe en 0/1
y <- ifelse(donnees$LOW == "y",1,0)
# vérification
print(table(donnees$LOW,y))
## y
## 0 1
## n 130 0
## y 0 59
# *******************************
# Q.8 - création de la courbe ROC
# *******************************
# Q.7 - récupération du nombre de positifs
npos <- sum(y)
print(npos)
## [1] 59
# nombre de négatifs
nneg <- sum(1-y)
print(nneg)
## [1] 130
# Q.8a - créer un data frame avec y et scores
df <- data.frame(y,scores)
print(head(df))
## y scores
## 1 1 0.7232608
## 2 1 0.2203922
## 3 1 0.3950290
## 4 1 0.8838059
## 5 1 0.4517460
## 6 1 0.1344577
# Q.8b - trier le data frame avec les scores decroissants
dfSorted <- df[order(scores,decreasing=TRUE),]
# premières lignes
print(head(dfSorted))
## y scores
## 4 1 0.8838059
## 12 1 0.8563734
## 153 0 0.8099383
## 26 1 0.7915973
## 36 1 0.7892067
## 54 1 0.7705263
# dernières lignes
print(tail(dfSorted))
## y scores
## 189 0 0.07040073
## 148 0 0.06716471
## 98 0 0.05924674
## 170 0 0.05452911
## 127 0 0.04255867
## 82 0 0.03686464
# Q.8c - taux de faux positifs
tfp <- cumsum(1-dfSorted$y)/nneg
print(head(tfp))
## [1] 0.000000000 0.000000000 0.007692308 0.007692308 0.007692308 0.007692308
# Q.8d - taux de vrais positifs
tvp <- cumsum(dfSorted$y)/npos
print(head(tvp))
## [1] 0.01694915 0.03389831 0.03389831 0.05084746 0.06779661 0.08474576
# Q.8e - graphique courbe ROC + ajout de la diagonale
plot(tfp,tvp,main="Courbe ROC",type="l",col="blue")
abline(a=0,b=1)
# Q.9 - courbe ROC avec la librairie ROCR
# Q.9a - utilisation de la librairie ROCR
library(ROCR)
# Q.9b - construire un objet prediction avec les scores et la var. cible
pred <- ROCR::prediction(scores,donnees$LOW)
print(pred)
## A prediction instance
## with 189 data points
# propriétés - résultats détaillés
# attributes(pred)
# Q.9c - objet perfromance mesurée
graphs_roc <- ROCR::performance(pred,measure="tpr",x.measure="fpr")
print(graphs_roc)
## A performance instance
## 'False positive rate' vs. 'True positive rate' (alpha: 'Cutoff')
## with 183 data points
# Q.9d - graphique courbe ROC
ROCR::plot(graphs_roc,xlab="TFP",ylab="TFP",col="darkblue")
abline(a=0,b=1)
Calcul de l’AUC
# Q.10a -- rangs - tout simplement inverse du numéro pour la base triée
# selon les scores
rk <- seq(length(y),1,-1)
print(rk)
## [1] 189 188 187 186 185 184 183 182 181 180 179 178 177 176 175 174 173 172
## [19] 171 170 169 168 167 166 165 164 163 162 161 160 159 158 157 156 155 154
## [37] 153 152 151 150 149 148 147 146 145 144 143 142 141 140 139 138 137 136
## [55] 135 134 133 132 131 130 129 128 127 126 125 124 123 122 121 120 119 118
## [73] 117 116 115 114 113 112 111 110 109 108 107 106 105 104 103 102 101 100
## [91] 99 98 97 96 95 94 93 92 91 90 89 88 87 86 85 84 83 82
## [109] 81 80 79 78 77 76 75 74 73 72 71 70 69 68 67 66 65 64
## [127] 63 62 61 60 59 58 57 56 55 54 53 52 51 50 49 48 47 46
## [145] 45 44 43 42 41 40 39 38 37 36 35 34 33 32 31 30 29 28
## [163] 27 26 25 24 23 22 21 20 19 18 17 16 15 14 13 12 11 10
## [181] 9 8 7 6 5 4 3 2 1
# Q10.b - somme des rangs des positifs
# c.-à-d. des individus appartenant à la modalité cible
s_plus <- sum(rk * dfSorted$y)
print(s_plus)
## [1] 7492
# Q10.c - stat de Mann-Whitney
U_plus <- s_plus - 0.5*(npos*(npos+1))
print(U_plus)
## [1] 5722
# Q10.d - valeur de l'auc
auc <- U_plus/(npos*nneg)
print(auc)
## [1] 0.7460235
# Q.11 - avec la librairie ROCR - objet performance
auc_roc <- ROCR::performance(pred,measure="auc")
print(auc_roc)
## A performance instance
## 'Area under the ROC curve'
# pour obtenir la valeur de l'AUC
print(auc_roc@y.values)
## [[1]]
## [1] 0.7460235
Courbe ROC et AUC en Leave-One-Out
# Q.12 - vecteur préparatoire des résultats
scores_lvo <- rep(0,length(y))
# Q.13 - boucle de calcul leave-one-out
# tour à tour, chaque observation est traitée
# comme un individu supplémentaire
for (i in 1:length(y)){
#modélisation
temp_modele <- glm(LOW ~ ., data = donnees[-i,], family = binomial)
#prédiction
scores_lvo[i] <- predict(temp_modele, newdata = donnees[i,], type = "response")
}
#summary des scores
print(summary(scores_lvo))
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.03757 0.17929 0.25295 0.31212 0.40586 0.87359
# Q.14 - calcul de prediction au sens de ROCR
pred_lvo <- ROCR::prediction(scores_lvo,donnees$LOW)
#performance au sens de ROCR
graphs_roc_lvo <- ROCR::performance(pred_lvo,measure="tpr",x.measure="fpr")
#courbe ROC en LVO en vert
ROCR::plot(graphs_roc_lvo,xlab="TFP",ylab="TFP",col="green", main="Courbe ROC")
abline(a=0,b=1)
#rajouter la courbe en resubstitution
#qui est effectivement surévaluée
lines(tfp,tvp,col="blue")
#ajouter une légende
legend(0.6,0.4,legend=c("Resubstitution","Leave-one-out"),col=c("blue","green"),lty=1:2)
# Q.15 - calcul de l'AUC
# on obtient la valeur corrigée, plus réaliste
# plus conforme aux véritables performances
# du modèle dans la population
auc_roc_lvo <- ROCR::performance(pred_lvo,measure="auc")
print(auc_roc_lvo@y.values)
## [[1]]
## [1] 0.6827901