library(ggplot2)
library(dplyr)
library("ggfortify")
library("questionr")

Exemple 1

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.

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)
attach(maladie)
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(CHD, levels = c(0,1), labels = c("No","Yes"))
maladie$AGRP = factor(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 <fctr> 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  <fctr> 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
maladie.logit = glm(CHD~AGE, family=binomial(link="logit"),data=maladie)
coef(maladie.logit)
## (Intercept)         AGE 
##  -5.3094534   0.1109211
#Resumé complet 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

Une manière de tester 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
maladie.logit.2 = glm(CHD~I(AGE^2), family=binomial(link="logit"))
summary(maladie.logit.2)
## 
## Call:
## glm(formula = CHD ~ I(AGE^2), family = binomial(link = "logit"))
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -2.1019  -0.8149  -0.5149   0.8191   2.1253  
## 
## Coefficients:
##               Estimate Std. Error z value Pr(>|z|)    
## (Intercept) -2.9117530  0.6184278  -4.708 2.50e-06 ***
## I(AGE^2)     0.0012218  0.0002636   4.635 3.57e-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.39  on 98  degrees of freedom
## AIC: 111.39
## 
## Number of Fisher Scoring iterations: 4
# Motivation pour la prediction et courbe du ROC

n = nrow(maladie); 
p=0.75; ntrain=floor(p*n) ; ntest = n-ntrain
ind=sample(n); ind_train=ind[1:ntrain] ; ind_test=ind[(ntrain+1):n]
maladie.train=maladie[ind_train,] ; maladie.test=maladie[ind_test,]

#score
score.glm  =  predict.glm(maladie.logit, maladie.test, 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.test[,4],pred.glm)
##      pred.glm
##       No Yes
##   No   9   2
##   Yes  3  11
mean(maladie.test[,4]!=pred.glm)
## [1] 0.2
#Proportion de vrais positifs (au seuil fixé 0.5)
sum( score.glm >= 0.5 & maladie.test[,4]=="Yes")/sum(maladie.test[,4]=="Yes")
## [1] 0.7857143
#Proportion de vrais negatifs
sum( score.glm < 0.5 & maladie.test[,4]=="No")/sum(maladie.test[,4]=="No")
## [1] 0.8181818
#Proportion de faux poistifs
sum( score.glm >= 0.5 & maladie.test[,1]=="No")/sum(maladie.test[,4]=="No")
## [1] 0

Exemple 2

Nous traitons un problème de défaut bancaire. Nous cherchons à déterminer quels clients seront en défaut sur leur dette de carte de crédit (ici = 1 si le client fait défaut sur sa dette). La variable est la variable réponse.

Nous disposons d’un échantillon de taille \(10000\) et 3 variables explicativesLes variables explicatives sont les suivantes :

library(ISLR)
data(Default); attach(Default)
glimpse(Default)
## Observations: 10,000
## Variables: 4
## $ default <fctr> No, No, No, No, No, No, No, No, No, No, No, No, No, N...
## $ student <fctr> No, Yes, No, No, No, Yes, No, Yes, No, No, Yes, Yes, ...
## $ balance <dbl> 729.5265, 817.1804, 1073.5492, 529.2506, 785.6559, 919...
## $ income  <dbl> 44361.625, 12106.135, 31767.139, 35704.494, 38463.496,...
#résumé des données (un peu de stats descriptive)
summary(Default)
##  default    student       balance           income     
##  No :9667   No :7056   Min.   :   0.0   Min.   :  772  
##  Yes: 333   Yes:2944   1st Qu.: 481.7   1st Qu.:21340  
##                        Median : 823.6   Median :34553  
##                        Mean   : 835.4   Mean   :33517  
##                        3rd Qu.:1166.3   3rd Qu.:43808  
##                        Max.   :2654.3   Max.   :73554
qplot(data = Default, default)

qplot(data = Default, balance)

qplot(data = Default, income)

qplot(data = Default, student)

Default.logit=glm(default~balance + income + student,family=binomial(link="logit"),data=Default)
summary(Default.logit)
## 
## Call:
## glm(formula = default ~ balance + income + student, family = binomial(link = "logit"), 
##     data = Default)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -2.4691  -0.1418  -0.0557  -0.0203   3.7383  
## 
## Coefficients:
##               Estimate Std. Error z value Pr(>|z|)    
## (Intercept) -1.087e+01  4.923e-01 -22.080  < 2e-16 ***
## balance      5.737e-03  2.319e-04  24.738  < 2e-16 ***
## income       3.033e-06  8.203e-06   0.370  0.71152    
## studentYes  -6.468e-01  2.363e-01  -2.738  0.00619 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 2920.6  on 9999  degrees of freedom
## Residual deviance: 1571.5  on 9996  degrees of freedom
## AIC: 1579.5
## 
## Number of Fisher Scoring iterations: 8

Selection de modèles par selection de variables

# Sélection par AIC et backward
step(Default.logit,direction="backward")
## Start:  AIC=1579.54
## default ~ balance + income + student
## 
##           Df Deviance    AIC
## - income   1   1571.7 1577.7
## <none>         1571.5 1579.5
## - student  1   1579.0 1585.0
## - balance  1   2907.5 2913.5
## 
## Step:  AIC=1577.68
## default ~ balance + student
## 
##           Df Deviance    AIC
## <none>         1571.7 1577.7
## - student  1   1596.5 1600.5
## - balance  1   2908.7 2912.7
## 
## Call:  glm(formula = default ~ balance + student, family = binomial(link = "logit"), 
##     data = Default)
## 
## Coefficients:
## (Intercept)      balance   studentYes  
##  -10.749496     0.005738    -0.714878  
## 
## Degrees of Freedom: 9999 Total (i.e. Null);  9997 Residual
## Null Deviance:       2921 
## Residual Deviance: 1572  AIC: 1578
Default.logit.choisi=glm(default~balance + student,family=binomial(link="logit"),data=Default)
summary(Default.logit.choisi)
## 
## Call:
## glm(formula = default ~ balance + student, family = binomial(link = "logit"), 
##     data = Default)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -2.4578  -0.1422  -0.0559  -0.0203   3.7435  
## 
## Coefficients:
##               Estimate Std. Error z value Pr(>|z|)    
## (Intercept) -1.075e+01  3.692e-01 -29.116  < 2e-16 ***
## balance      5.738e-03  2.318e-04  24.750  < 2e-16 ***
## studentYes  -7.149e-01  1.475e-01  -4.846 1.26e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 2920.6  on 9999  degrees of freedom
## Residual deviance: 1571.7  on 9997  degrees of freedom
## AIC: 1577.7
## 
## Number of Fisher Scoring iterations: 8
odds.ratio(Default.logit.choisi)
##                OR 2.5 % 97.5 %         p    
## (Intercept) 0.000 0.000  0.000 < 2.2e-16 ***
## balance     1.006 1.005  1.006 < 2.2e-16 ***
## studentYes  0.489 0.365  0.651  1.26e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#score
score.glm  =  predict.glm(Default.logit.choisi, data=Default, 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(Default$default,pred.glm)
##      pred.glm
##         No  Yes
##   No  9628   39
##   Yes  228  105
mean(Default$default!=pred.glm)
## [1] 0.0267
#Estimation de l'erreurr par validation croisée
glm_KFold_CV <- function (D,name,K) {
formula=formula(paste(name,"~ ."))
n=nrow(D)
Error_glm=dim(K)
ind_Fold=sample(rep(1:K,ceiling(n/K)),n)
for (i in 1:K){
  Dtrain=D[which(ind_Fold!=i),]; Dtest=D[which(ind_Fold==i),]
  fit.glm=glm(formula,data=Dtrain,family=binomial)
  pred.glm  = as.numeric( predict(fit.glm, Dtest, type="response") >= 0.5 )
  Error_glm[i]=sum((as.numeric(Dtest[,name])-1)!=pred.glm)/nrow(Dtest)
}
mean(Error_glm)
}
VC.glm=glm_KFold_CV(Default,"default",5)
VC.glm
## [1] 0.0269
# Motivation pour la prediction et courbe du ROC

n = nrow(Default); 
p=0.75; ntrain=floor(p*n) ; ntest = n-ntrain
ind=sample(n); ind_train=ind[1:ntrain] ; ind_test=ind[(ntrain+1):n]
Default.train=Default[ind_train,] ; Default.test=Default[ind_test,]

fit.glm=glm(default~.,data=Default.train,family=binomial)

#score
score.glm  =  predict.glm(fit.glm, Default.test, 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(Default.test[,1],pred.glm)
##      pred.glm
##         No  Yes
##   No  2390   13
##   Yes   73   24
mean(Default.test[,1]!=pred.glm)
## [1] 0.0344
#Proportion de vrais positifs (au seuil fixé 0.5)
sum( score.glm >= 0.5 & Default.test[,1]=="Yes")/sum(Default.test[,1]=="Yes")
## [1] 0.2474227
#Proportion de vrais negatifs
sum( score.glm < 0.5 & Default.test[,1]=="No")/sum(Default.test[,1]=="No")
## [1] 0.9945901
#Proportion de faux poistifs
sum( score.glm >= 0.5 & Default.test[,1]=="No")/sum(Default.test[,1]=="No")
## [1] 0.005409904
#Taux de vrais positifs (TPR) et taux faux positifs (FPR) pour difeferents valeurs de s

#seuil s
s.glm=seq(0,1,.01)
#des fois c'est mieux avec les quantiles car plus des precision dans les rangs des scores
#s.glm=quantile(score.glm,probs=seq(0,1,.01))

#Initialisation :
absc.glm=numeric(length(s.glm));ordo.glm=numeric(length(s.glm))

for (i in 1:length(s.glm)){
  ordo.glm[i]=sum( score.glm>=s.glm[i] & Default.test[,1]=="Yes")/sum(Default.test[,1]=="Yes")
  absc.glm[i]=sum( score.glm>=s.glm[i] & Default.test[,1]=="No")/sum(Default.test[,1]=="No")
}
ROC = data.frame(FPR=absc.glm, TPR=ordo.glm)
#Courbe Roc : taux de vrais positifs (TPR) en fonction du taux de faux positifs (FPR)
ggplot(ROC,aes(x=FPR,y=TPR)) + geom_path(aes(),color="red")  +
   geom_segment(aes(x = 0, y = 0, xend = 1, yend = 1), linetype = 2, color = "black")

Exemple 3

Nous disposons de données recuilles à Hewlett-Packards Labs, qui classe 4601 e-mails comme spam ou non spam. Il y a 57 variables indiquant la fréquence de certains mots et de caractères dans l’e-mail. Ces données sont disponnibles sous la librarie kernlab du logiciel R. Le donneur des e-mails c’est George Forman de HP Labs.

library(kernlab)
data(spam);attach(spam)
names(spam)
##  [1] "make"              "address"           "all"              
##  [4] "num3d"             "our"               "over"             
##  [7] "remove"            "internet"          "order"            
## [10] "mail"              "receive"           "will"             
## [13] "people"            "report"            "addresses"        
## [16] "free"              "business"          "email"            
## [19] "you"               "credit"            "your"             
## [22] "font"              "num000"            "money"            
## [25] "hp"                "hpl"               "george"           
## [28] "num650"            "lab"               "labs"             
## [31] "telnet"            "num857"            "data"             
## [34] "num415"            "num85"             "technology"       
## [37] "num1999"           "parts"             "pm"               
## [40] "direct"            "cs"                "meeting"          
## [43] "original"          "project"           "re"               
## [46] "edu"               "table"             "conference"       
## [49] "charSemicolon"     "charRoundbracket"  "charSquarebracket"
## [52] "charExclamation"   "charDollar"        "charHash"         
## [55] "capitalAve"        "capitalLong"       "capitalTotal"     
## [58] "type"
fit.glm=glm(type~.,data=spam,family=binomial)
#summary(fit.glm)

# Decoupage (2 sous échantillons)
n = nrow(spam); d= ncol(spam)
p=0.75; ntrain=floor(p*n) ; ntest = n-ntrain
ind=sample(n); ind_train=ind[1:ntrain] ; ind_test=ind[(ntrain+1):n]
spam.train=spam[ind_train,] ; spam.test=spam[ind_test,]

# Ajustement du modèle logistique
fit.glm=glm(type~.,data=spam.train,family=binomial)
#score : estimation de p(x) = P(Y=spam|X=x) à l'aidele modèle logistique
score.glm  =  predict.glm(fit.glm, spam.test, 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("nonspam","spam"))
table(spam.test[,d],pred.glm)
##          pred.glm
##           nonspam spam
##   nonspam     670   32
##   spam         46  403
mean(spam.test[,d]!=pred.glm)
## [1] 0.06776716
# Courbe du ROC

#Proportion de vrais positifs (au seuil fixé 0.5)
sum( score.glm >= 0.5 & spam.test[,d]=="spam")/sum(spam.test[,d]=="spam")
## [1] 0.8975501
#Proportion de vrais negatifs
sum( score.glm < 0.5 & spam.test[,d]=="nonspam")/sum(spam.test[,d]=="nonspam")
## [1] 0.954416
#Proportion de faux poistifs
sum( score.glm >= 0.5 & spam.test[,d]=="nonspam")/sum(spam.test[,d]=="nonspam")
## [1] 0.04558405
#Taux de vrais positifs (TPR) et taux faux positifs (FPR) pour difeferents valeurs de s

#seuil s
s.glm=seq(0,1,.01)
#des fois c'est mieux avec les quantiles car plus des precision dans les rangs des scores
#s.glm=quantile(score.glm,probs=seq(0,1,.01))
#Initialisation :
absc.glm=numeric(length(s.glm));ordo.glm=numeric(length(s.glm))

for (i in 1:length(s.glm)){
  ordo.glm[i]=sum( score.glm>=s.glm[i] & spam.test[,d]=="spam")/sum(spam.test[,d]=="spam")
  absc.glm[i]=sum( score.glm>=s.glm[i] & spam.test[,d]=="nonspam")/sum(spam.test[,d]=="nonspam")
}


ROC = data.frame(FPR=absc.glm, TPR=ordo.glm)
#Courbe Roc : Taux de vrais positifs (TPR) en fonction de taux de faux positifs (FPR)
ggplot(ROC,aes(x=FPR,y=TPR)) + geom_path(aes(),color="red")  +
   geom_segment(aes(x = 0, y = 0, xend = 1, yend = 1), linetype = 2, color = "black")