Classificazione II (notebook 5)

Statistica II - 750AA

Classificazione basata su Bayes

In questo notebook R consideriamo i metodi di classificazione basati sulla formula di Bayes: Naive Bayes, LDA (Linear Discriminant Analysis), QDA (Quadratic Discriminant Analysis), regressione logistica, regressione softmax. Utilizzeremo come in altre occasioni dataset generati ma anche dataset reali.

Iniziamo con caricare delle librerie necessarie. Installate i pacchetti se non sono già disponibili.

# Carica le librerie necessarie
library(ggplot2)
library(caret) 
library(e1071)  # Per Naive Bayes semplificato
library(naivebayes) # Naive Bayes avanzato
library(MASS)   # Per LDA e QDA

Datasets

Consideriamo 3 datasets: uno generato, uno dal pacchetto caret relativo a dati di individui per accesso al credito in Germania, e infine i dati sugli infortuni sul lavoro in Toscana.

Dataset generato

Generiamo una tabella con \(4\) features: 2 di tipo discreto, 2 di tipo numerico e una etichetta binaria, e alleniamo un modello di Naive Bayes su di esso.

# features (si può usare anche rbinom)
x1 <- sample( c(TRUE, FALSE), size=200, replace=TRUE, prob=c(0.2, 0.8))
x2 <- sample(  c(TRUE, FALSE), size=50, replace=TRUE, prob=c(0.5, 0.5))
y1 <- sample(  c(TRUE, FALSE), size=200, replace=TRUE, prob=c(0.3, 0.7))
y2 <- sample(  c(TRUE, FALSE), size=50, replace=TRUE, prob=c(0.2, 0.8))
z1 <- rnorm(200, mean=0, sd=1)
z2 <- rnorm(50, mean=2, sd=2)
w1 <- rpois(200, lambda=2)
w2 <- rpois(50, lambda=3)

# labels
data_generated <- data.frame("x" = c(x1, x2), "y"=c(y1, y2), "z"=c(z1, z2), "w"=c(w1, w2), "label"=as.factor(c(rep(-1, 200), rep(1,50))))

head(data_generated)

Partizioniamo il dataset in train e test usando il comando createDataPartition() dal pacchetto caret. Questo campiona in modo casuale come sample() ma cerca di rispettare la anche distribuzione delle diverse classi.

train_index <- createDataPartition(data_generated$label, p = 0.8, list = FALSE)
train_data_generated <- data_generated[train_index, ]
test_data_generated <- data_generated[-train_index, ]

print(paste("Taglia del train set:",nrow(train_data_generated) ))

[1] "Taglia del train set: 200"

print(paste("Taglia del test set:", nrow(test_data_generated) ))

[1] "Taglia del test set: 50"

Dataset German Credit Data

Nel pacchetto caret vi sono diversi dataset precaricati. Usiamo il dataset GermanCredit per addestrare e testare i modelli di classificazione basati su Bayes.

data("GermanCredit")

# Visualizza le prime righe del dataset

head(GermanCredit)

NA

Suddividiamo il dataset in training e test.

set.seed(42)

train_index_gc <- createDataPartition(GermanCredit$Class, p = 0.7, list = FALSE)
train_set_gc <- GermanCredit[train_index_gc, ]
test_set_gc <- GermanCredit[-train_index_gc, ]

print(paste("Taglia del train set:",nrow(train_set_gc) ))

[1] "Taglia del train set: 700"

print(paste("Taglia del test set:", nrow(test_set_gc) ))

[1] "Taglia del test set: 300"

Dataset Infortuni sul Lavoro in Toscana

I dati sugli infortuni sul lavoro in Italia divisi per regione e con alta granularità sono disponibili online sul sito Inail Open Data. Abbiamo scaricato il dataset relativo ai dati mensili e lo abbiamo salvato in formato CSV nella cartella datasets. Carichiamo il dataset usando il pacchetto rio.

library(rio)

infortuni_toscana <- import("datasets/DatiConCadenzaMensileInfortuniToscana.csv")

head(infortuni_toscana)

NA

Notiamo che in tutto il dataset comprende 61818 righe.

Selezioniamo solo alcune colonne di features rilevanti.

infortuni_toscana_tidy <- infortuni_toscana[ c(5, 7, 8, 9, 10, 11,13, 15)]


head(infortuni_toscana_tidy)

Un obiettivo potrebbe essere di classificare il luogo di nascita (ITALIA/ESTERO) dell’infortunato/a in base alle altre caratteristiche. Trasformiamo quindi la colonna LuogoNascita in una feature binaria (NatItal) ed estraiamo anche in questo caso il dataset in training e test set.

set.seed(42)

infortuni_toscana_tidy$NatItal <- factor(infortuni_toscana_tidy$LuogoNascita == "ITAL")

# rimuoviamo la colonna LuogoNascita
infortuni_toscana_tidy$LuogoNascita <- NULL

# estraiamo train e test set usando la colonnna NatItal come label

train_index_it <- createDataPartition(infortuni_toscana_tidy$NatItal, p = 0.7, list = FALSE)

train_set_it <- infortuni_toscana_tidy[train_index_it, ]
test_set_it <- infortuni_toscana_tidy[-train_index_it, ]

print(paste("Taglia del train set:",nrow(train_set_it) ))

[1] "Taglia del train set: 43274"

print(paste("Taglia del test set:", nrow(test_set_it) ))

[1] "Taglia del test set: 18544"

Naive Bayes

Ricordiamo che il classificatore Naive Bayes si basa sulla formula di Bayes, assumendo l’indipendenza condizionale tra le features data la label. La previsione per una nuova osservazione \(x=(x_1, x_2, \ldots, x_d)\) viene effettuata calcolando la probabilità a posteriori per ciascuna label \(\ell\) e scegliendo la label con la probabilità massima.

Ci sono vari pacchetti che contengono funzioni per calcolare classificatori Naive Bayes in R: caret (consigliato anche se un po’ più lento), e1071 e naivebayes.

Dataset generato

Vediamoli in azione sul dataset generato. Il comando è nel pacchetto caret è train() specificando nelle opzioni il metodo method="naive_bayes". Bisogna anche specificare una formula come in altri modelli di regressione/classificazione in R (questo è il primo caso che vediamo).

# Allenamento del modello Naive Bayes con caret
nb_model_caret <- train(label ~ ., data = train_data_generated, method = "naive_bayes")

#print(nb_model_caret)

# valuta sul test set
nb_pred_caret <- predict(nb_model_caret, test_data_generated)

confusionMatrix(nb_pred_caret, test_data_generated$label)

Confusion Matrix and Statistics

          Reference
Prediction -1  1
        -1 38  3
        1   2  7
                                          
               Accuracy : 0.9             
                 95% CI : (0.7819, 0.9667)
    No Information Rate : 0.8             
    P-Value [Acc > NIR] : 0.04803         
                                          
                  Kappa : 0.6753          
                                          
 Mcnemar's Test P-Value : 1.00000         
                                          
            Sensitivity : 0.9500          
            Specificity : 0.7000          
         Pos Pred Value : 0.9268          
         Neg Pred Value : 0.7778          
             Prevalence : 0.8000          
         Detection Rate : 0.7600          
   Detection Prevalence : 0.8200          
      Balanced Accuracy : 0.8250          
                                          
       'Positive' Class : -1              
                                          

Confrontiamo con il pacchetto e1071.

# Allenamento del modello Naive Bayes con e1071

nb_model_e1071 <- naiveBayes(label ~ x+z+w, data = train_data_generated)

# la formula indica che vogliamo prevedere la colonna "label" in funzione di tutte le altre colonne (features). Se avessimo voluto usare solo alcune features ad esempio la colonna x, avremmo scritto ad esempio label ~ x

print(nb_model_e1071)

Naive Bayes Classifier for Discrete Predictors

Call:
naiveBayes.default(x = X, y = Y, laplace = laplace)

A-priori probabilities:
Y
 -1   1 
0.8 0.2 

Conditional probabilities:
    x
Y      FALSE    TRUE
  -1 0.76875 0.23125
  1  0.47500 0.52500

    z
Y          [,1]      [,2]
  -1 0.08540609 0.9381225
  1  2.66902610 1.7484799

    w
Y       [,1]     [,2]
  -1 1.89375 1.403486
  1  2.82500 1.583367

Facciamo ora le previsioni sul test set e calcoliamo la matrice di confusione e l’accuratezza. Il comando è predict() come per altri modelli in R. La sintassi richiede di specificare il modello allenato e il dataset su cui fare le previsioni.

# Previsione sul test set

nb_pred_e1071 <- predict(nb_model_e1071, test_data_generated)

# il vettore contiene le etichette classificate.
# verifichiamo l'accuratezza e altre metriche con la matrice di confusione

nb_conf_matrix_e1071 <- confusionMatrix(nb_pred_e1071, test_data_generated$label)

print(nb_conf_matrix_e1071)

Confusion Matrix and Statistics

          Reference
Prediction -1  1
        -1 38  3
        1   2  7
                                          
               Accuracy : 0.9             
                 95% CI : (0.7819, 0.9667)
    No Information Rate : 0.8             
    P-Value [Acc > NIR] : 0.04803         
                                          
                  Kappa : 0.6753          
                                          
 Mcnemar's Test P-Value : 1.00000         
                                          
            Sensitivity : 0.9500          
            Specificity : 0.7000          
         Pos Pred Value : 0.9268          
         Neg Pred Value : 0.7778          
             Prevalence : 0.8000          
         Detection Rate : 0.7600          
   Detection Prevalence : 0.8200          
      Balanced Accuracy : 0.8250          
                                          
       'Positive' Class : -1              
                                          

Possiamo anche plottare la curva ROC per valutare la performance del classificatore al variare della soglia di accetazione (o equivalentemente delle probabilità a priori). Usiamo il pacchetto pROC per calcolare e plottare la curva ROC.

library(pROC)

# Ci servono le probabilità predette per la classe positiva, non soltanto l'etichetta predetta. Per questo specifichiamo type="raw" e prendiamo la seconda colonna (corrispondente alla classe +1)

nb_pred_prob_e1071 <- predict(nb_model_e1071, test_data_generated, type = "raw")[, 2]

# Calcola la ROC curve ed esegui un plot usando ggroc (pacchetto pROC). Specifichiamo legacy.axes=TRUE per avere FPR nelle ascisse

roc_nb_e1071 <- roc(test_data_generated$label, nb_pred_prob_e1071)

ggroc(roc_nb_e1071, legacy.axes = TRUE) +
  ggtitle("ROC Curve per Naive Bayes dataset generato") +
  xlab("False Positive Rate") +
  ylab("True Positive Rate")

auc(roc_nb_e1071)

Area under the curve: 0.9475

Per confronto, applichiamo invece la funzione naive_bayes() dal pacchetto naivebayes.

# Allenamento del modello Naive Bayes con naivebayes

nb_model_nb <- naive_bayes(label ~ ., data = train_data_generated)

print(nb_model_nb)

========================= Naive Bayes ==========================

Call:
naive_bayes.formula(formula = label ~ ., data = train_data_generated)

---------------------------------------------------------------- 
 
Laplace smoothing: 0

---------------------------------------------------------------- 
 
A priori probabilities: 

 -1   1 
0.8 0.2 

---------------------------------------------------------------- 
 
Tables: 

---------------------------------------------------------------- 
:: x (Bernoulli) 
---------------------------------------------------------------- 
       
x            -1       1
  FALSE 0.76875 0.47500
  TRUE  0.23125 0.52500

---------------------------------------------------------------- 
:: y (Bernoulli) 
---------------------------------------------------------------- 
       
y            -1       1
  FALSE 0.63125 0.75000
  TRUE  0.36875 0.25000

---------------------------------------------------------------- 
:: z (Gaussian) 
---------------------------------------------------------------- 
      
z              -1          1
  mean 0.08540609 2.66902610
  sd   0.93812248 1.74847988

---------------------------------------------------------------- 
:: w (Gaussian) 
---------------------------------------------------------------- 
      
w            -1        1
  mean 1.893750 2.825000
  sd   1.403486 1.583367

----------------------------------------------------------------

nb_pred_nb <- predict(nb_model_nb, test_data_generated[1:2])

# Confrontiamo le due previsioni

table(nb_pred_nb, nb_pred_e1071)

          nb_pred_e1071
nb_pred_nb -1  1
        -1 41  9
        1   0  0

Applichiamo infine \(k\)-nn sullo stesso dataset per confronto. Usiamo sempre il pacchetto caret per coerenza.

# Allenamento del modello KNN


# usiamo il pacchetto caret e cross validation per scegliere il valore di k
knn_model_generated <- train(label ~ ., data = train_data_generated, method="knn", tuneLength=10)

print(knn_model_generated)

k-Nearest Neighbors 

200 samples
  4 predictor
  2 classes: '-1', '1' 

No pre-processing
Resampling: Bootstrapped (25 reps) 
Summary of sample sizes: 200, 200, 200, 200, 200, 200, ... 
Resampling results across tuning parameters:

  k   Accuracy   Kappa    
   5  0.8956391  0.6405813
   7  0.9066797  0.6720492
   9  0.9053593  0.6676041
  11  0.9059815  0.6669415
  13  0.9059195  0.6615438
  15  0.9088043  0.6686325
  17  0.9056190  0.6547086
  19  0.9024574  0.6414763
  21  0.9051536  0.6544913
  23  0.9045415  0.6489380

Accuracy was used to select the optimal model using the
 largest value.
The final value used for the model was k = 15.

plot(knn_model_generated)

# Previsione sul test set
knn_pred_generated <- predict(knn_model_generated, test_data_generated)

# Matrice di confusione e accuratezza
knn_conf_matrix_generated <- confusionMatrix(knn_pred_generated, test_data_generated$label)

print(knn_conf_matrix_generated)

Confusion Matrix and Statistics

          Reference
Prediction -1  1
        -1 39  5
        1   1  5
                                          
               Accuracy : 0.88            
                 95% CI : (0.7569, 0.9547)
    No Information Rate : 0.8             
    P-Value [Acc > NIR] : 0.1034          
                                          
                  Kappa : 0.5588          
                                          
 Mcnemar's Test P-Value : 0.2207          
                                          
            Sensitivity : 0.9750          
            Specificity : 0.5000          
         Pos Pred Value : 0.8864          
         Neg Pred Value : 0.8333          
             Prevalence : 0.8000          
         Detection Rate : 0.7800          
   Detection Prevalence : 0.8800          
      Balanced Accuracy : 0.7375          
                                          
       'Positive' Class : -1              
                                          

Dataset German Credit Data

Applichiamo ora Naive Bayes sul dataset German Credit Data. Visualizziamo prima la correlazione tra le features usando una heatmap.

library(ggcorrplot)

# Calcola la matrice di correlazione relativa al training set per gli individui con una buona/cattiva classificazione di credito

cor_matrix_gc <- cor(train_set_gc[train_set_gc$Class=="Good", sapply(train_set_gc, is.numeric)])
# Plot della heatmap
ggcorrplot(cor_matrix_gc, 
           method = "square", 
           type = "full", 
           lab = FALSE, 
           title = "Matrice di Correlazione - German Credit Data", tl.cex=0)

Applichiamo Naive Bayes per determinare se un individuo ha un buon o cattivo credito (colonna Class).

# Allenamento del modello Naive Bayes

nb_model_gc <- naiveBayes(Class ~ ., data = train_set_gc)

# Previsione sul test set

nb_pred_gc <- predict(nb_model_gc, test_set_gc)

# Matrice di confusione e accuratezza

nb_conf_matrix_gc <- confusionMatrix(nb_pred_gc, test_set_gc$Class)

print(nb_conf_matrix_gc)

Confusion Matrix and Statistics

          Reference
Prediction Bad Good
      Bad   70   95
      Good  20  115
                                        
               Accuracy : 0.6167        
                 95% CI : (0.559, 0.672)
    No Information Rate : 0.7           
    P-Value [Acc > NIR] : 0.9992        
                                        
                  Kappa : 0.2628        
                                        
 Mcnemar's Test P-Value : 5.181e-12     
                                        
            Sensitivity : 0.7778        
            Specificity : 0.5476        
         Pos Pred Value : 0.4242        
         Neg Pred Value : 0.8519        
             Prevalence : 0.3000        
         Detection Rate : 0.2333        
   Detection Prevalence : 0.5500        
      Balanced Accuracy : 0.6627        
                                        
       'Positive' Class : Bad           
                                        

Plottiamo anche in questo caso la curva ROC.

# attenzione: se usate caret il comando per ottenere le probabilità predette è leggermente diverso "prob" invece di "raw"

nb_pred_prob_gc <- predict(nb_model_gc, test_set_gc, type = "raw")[, 2]

# Calcola la ROC curve ed esegui un plot usando ggroc (pacchetto pROC). Specifichiamo legacy.axes=TRUE per avere FPR nelle ascisse

roc_nb_e1071 <- roc(test_set_gc$Class, nb_pred_prob_gc, levels=c("Good", "Bad"))

ggroc(roc_nb_e1071, legacy.axes = TRUE) +
  ggtitle("ROC Curve per Naive Bayes dataset generato") +
  xlab("False Positive Rate") +
  ylab("True Positive Rate")

NA
NA

Ci chiediamo: quali delle features sono più importanti per la classificazione? Usiamo la funzione varImp() dal pacchetto caret per calcolare l’importanza delle variabili ed eventualmente effettuare una selezione delle features per aumentare l’intepretabilità del modello.

# Calcola l'importanza delle variabili

var_imp_gc <- varImp(nb_model_gc)

Error in UseMethod("varImp") : 
  no applicable method for 'varImp' applied to an object of class "naiveBayes"

Esercizio: Provate a rifare l’analisi di Naive Bayes usando solo le prime 5 features più importanti (in base alla classifica calcolata sopra). Come cambiano le performance del modello?

Dataset Infortuni sul Lavoro in Toscana

Applichiamo Naive Bayes sul dataset degli infortuni sul lavoro in Toscana per prevedere se l’infortunato è nato in Italia o all’estero (colonna NatItal).

# Allenamento del modello Naive Bayes

#nb_model_it <- train(NatItal ~ ., data =train_set_it, method="naive_bayes")

nb_model_it <- naiveBayes(NatItal ~ ., data = train_set_it)
# Previsione sul test set

nb_pred_it <- predict(nb_model_it, test_set_it)

Consideriamo la performance con la solita matrice di confusione.

# Matrice di confusione e accuratezza

nb_conf_matrix_it <- confusionMatrix(as.factor(nb_pred_it), as.factor(test_set_it$NatItal))

print(nb_conf_matrix_it)

Confusion Matrix and Statistics

          Reference
Prediction FALSE  TRUE
     FALSE   462   549
     TRUE   3642 13891
                                        
               Accuracy : 0.774         
                 95% CI : (0.7679, 0.78)
    No Information Rate : 0.7787        
    P-Value [Acc > NIR] : 0.9389        
                                        
                  Kappa : 0.1021        
                                        
 Mcnemar's Test P-Value : <2e-16        
                                        
            Sensitivity : 0.11257       
            Specificity : 0.96198       
         Pos Pred Value : 0.45697       
         Neg Pred Value : 0.79228       
             Prevalence : 0.22131       
         Detection Rate : 0.02491       
   Detection Prevalence : 0.05452       
      Balanced Accuracy : 0.53728       
                                        
       'Positive' Class : FALSE         
                                        

nb_pred_prob_it <- predict(nb_model_it, test_set_it, type = "raw")[, 2]

# Calcola la ROC curve ed esegui un plot usando ggroc (pacchetto pROC). Specifichiamo legacy.axes=TRUE per avere FPR nelle ascisse


roc_nb_it <- roc(test_set_it$NatItal, nb_pred_prob_it, levels=c("TRUE", "FALSE"))

ggroc(roc_nb_it, legacy.axes = TRUE) +
  ggtitle("ROC Curve per Naive Bayes dataset INAIL") +
  xlab("False Positive Rate") +
  ylab("True Positive Rate")

LDA e QDA

Applichiamo lda al dataset generato. Il comando per allenare il modello è lda() dal pacchetto MASS. Oppure il comando train() dal pacchetto caret specificando method="lda".`

# Allenamento del modello LDA

lda_model_generated <- lda(label ~ ., data = train_data_generated)
# Previsione sul test set
lda_pred_generated <- predict(lda_model_generated, test_data_generated)
# Matrice di confusione e accuratezza
lda_conf_matrix_generated <- confusionMatrix(lda_pred_generated$class, test_data_generated$label)
print(lda_conf_matrix_generated)

Confusion Matrix and Statistics

          Reference
Prediction -1  1
        -1 38  4
        1   2  6
                                          
               Accuracy : 0.88            
                 95% CI : (0.7569, 0.9547)
    No Information Rate : 0.8             
    P-Value [Acc > NIR] : 0.1034          
                                          
                  Kappa : 0.5946          
                                          
 Mcnemar's Test P-Value : 0.6831          
                                          
            Sensitivity : 0.9500          
            Specificity : 0.6000          
         Pos Pred Value : 0.9048          
         Neg Pred Value : 0.7500          
             Prevalence : 0.8000          
         Detection Rate : 0.7600          
   Detection Prevalence : 0.8400          
      Balanced Accuracy : 0.7750          
                                          
       'Positive' Class : -1              
                                          

Applichiamo anche QDA per confronto. Il comando è qda().

# Allenamento del modello QDA
qda_model_generated <- qda(label ~ ., data = train_data_generated)
# Previsione sul test set
qda_pred_generated <- predict(qda_model_generated, test_data_generated)
# Matrice di confusione e accuratezza
qda_conf_matrix_generated <- confusionMatrix(qda_pred_generated$class, test_data_generated$label)
print(qda_conf_matrix_generated)

Confusion Matrix and Statistics

          Reference
Prediction -1  1
        -1 38  3
        1   2  7
                                          
               Accuracy : 0.9             
                 95% CI : (0.7819, 0.9667)
    No Information Rate : 0.8             
    P-Value [Acc > NIR] : 0.04803         
                                          
                  Kappa : 0.6753          
                                          
 Mcnemar's Test P-Value : 1.00000         
                                          
            Sensitivity : 0.9500          
            Specificity : 0.7000          
         Pos Pred Value : 0.9268          
         Neg Pred Value : 0.7778          
             Prevalence : 0.8000          
         Detection Rate : 0.7600          
   Detection Prevalence : 0.8200          
      Balanced Accuracy : 0.8250          
                                          
       'Positive' Class : -1              
                                          

Confrontiamo le curve ROC per i due modelli (LDA e QDA).

lda_pred_prob_generated <- predict(lda_model_generated, test_data_generated)$posterior[, 2]
qda_pred_prob_generated <- predict(qda_model_generated, test_data_generated)$posterior[, 2]
# Calcola la ROC curve ed esegui un plot usando ggroc (pacchetto pROC). Specifichiamo legacy.axes=TRUE per avere FPR nelle ascisse
roc_lda_generated <- roc(test_data_generated$label, lda_pred_prob_generated)
roc_qda_generated <- roc(test_data_generated$label, qda_pred_prob_generated)
ggroc(roc_lda_generated, legacy.axes = TRUE, color="blue") +
  ggtitle("ROC Curve per LDA (blu) e QDA (verde) - Dataset Generato") +
  xlab("False Positive Rate") +
  ylab("True Positive Rate") +
  geom_line(data = ggroc(roc_qda_generated)$data, aes(x = 1 - specificity, y = sensitivity), color="green") +
  theme(legend.position = "bottom") +
  guides(color = guide_legend(override.aes = list(size = 4)))

Dataset German Credit Data

Applichiamo ora LDA sul dataset German Credit Data.

# Allenamento del modello LDA

# se applicato interamente al dataset abbiamo un errore perchè alcune colonne all'interno dei gruppi sono costanti. Questo crea dei problemi nella stima della varianza. Rimuoviamole.

#train_set_gc[c(25, 26, 30, 33, 36,  39, 42, 44)] <- NULL

gc_tidy <- train_set_gc[1:10]

lda_model_gc <- train(Class ~ ., data = gc_tidy, method="lda")
# Previsione sul test set
lda_pred_gc <- predict(lda_model_gc, test_set_gc[1:10])

# Matrice di confusione e accuratezza
lda_conf_matrix_gc <- confusionMatrix(lda_pred_gc, test_set_gc$Class)
print(lda_conf_matrix_gc)

Plottiamo anche in questo caso la curva ROC (sovrapponendola a quella di Naive Bayes per confronto).

lda_pred_prob_gc <- predict(lda_model_gc, test_set_gc, type="prob")[, 2]

# Calcola la ROC curve ed esegui un plot usando ggroc (pacchetto pROC). Specifichiamo legacy.axes=TRUE per avere FPR nelle ascisse

roc_lda_gc <- roc(test_set_gc$Class, lda_pred_prob_gc, levels=c("Good", "Bad"))

ggroc(roc_lda_gc, legacy.axes = TRUE, color="blue") +
  ggtitle("ROC Curve per Naive Bayes (rosso) e LDA (blu)- German Credit Data") +
  xlab("False Positive Rate") +
  ylab("True Positive Rate") +
  geom_line(data = ggroc(roc_nb_e1071)$data, aes(x = 1 - specificity, y = sensitivity), color="red") +
  theme(legend.position = "bottom") +
  guides(color = guide_legend(override.aes = list(size = 4)))

Vediamo che LDA migliora leggermente le performance rispetto a Naive Bayes in questo caso.

Dataset Infortuni sul Lavoro in Toscana

Applichiamo ora LDA e QDA sul dataset degli infortuni sul lavoro in Toscana per prevedere se l’infortunato è nato in Italia o all’estero (colonna NatItal). Usiamo prima LDA.

# Allenamento del modello LDA
lda_model_it <- lda(NatItal ~ ., data = train_set_it)
# Previsione sul test set
lda_pred_it <- predict(lda_model_it, test_set_it)
# Matrice di confusione e accuratezza
lda_conf_matrix_it <- confusionMatrix(lda_pred_it$class, test_set_it$NatItal)
print(lda_conf_matrix_it)

Confusion Matrix and Statistics

          Reference
Prediction FALSE  TRUE
     FALSE   233   232
     TRUE   3871 14208
                                          
               Accuracy : 0.7787          
                 95% CI : (0.7727, 0.7847)
    No Information Rate : 0.7787          
    P-Value [Acc > NIR] : 0.4971          
                                          
                  Kappa : 0.0596          
                                          
 Mcnemar's Test P-Value : <2e-16          
                                          
            Sensitivity : 0.05677         
            Specificity : 0.98393         
         Pos Pred Value : 0.50108         
         Neg Pred Value : 0.78588         
             Prevalence : 0.22131         
         Detection Rate : 0.01256         
   Detection Prevalence : 0.02508         
      Balanced Accuracy : 0.52035         
                                          
       'Positive' Class : FALSE           
                                          

Confrontiamo le previsioni dei due modelli (LDA e Naive Bayes).

# tabella di contingenza

table(lda_pred_it$class, nb_pred_it)

       nb_pred_it
        FALSE  TRUE
  FALSE   410    55
  TRUE    601 17478

Confrontiamo le curve ROC.

lda_pred_prob_it <- predict(lda_model_it, test_set_it)$posterior[, 2]
# Calcola la ROC curve ed esegui un plot usando ggroc (pacchetto pROC). Specifichiamo legacy.axes=TRUE per avere FPR nelle ascisse

roc_lda_it <- roc(test_set_it$NatItal, lda_pred_prob_it, levels=c("TRUE", "FALSE"))

ggroc(roc_lda_it, legacy.axes = TRUE, color="blue") +
  ggtitle("ROC Curve per Naive Bayes (rosso) e LDA (blu)- INAIL Toscana") +
  xlab("False Positive Rate") +
  ylab("True Positive Rate") +
  geom_line(data = ggroc(roc_nb_it)$data, aes(x = 1 - specificity, y = sensitivity), color="red") +
  theme(legend.position = "bottom") +
  guides(color = guide_legend(override.aes = list(size = 4)))

NA
NA

Non stupisce che il modello LDA abbia performance peggiori: non sta utilizzando le colonne discrete in modo ottimale.

Esercizio: Provate a rifare l’analisi di LDA usando solo le colonne numeriche del dataset e trasformando eventuali factor in numerici (come il genere). Come cambiano le performance del modello?

Regressione Logistica

Consideriamo il dataset generato, che contiene 4 features (tutte numeriche) e due sole labels. Il comando per applicare la regressione logistica in R è glm() specificando family=binomial. Possiamo in alternativa usare il pacchetto caret con method="multinom" (per la regressione logistica multinomiale, che funziona anche per il caso binario).

# Allenamento del modello di regressione logistica
log_model_generated <- glm(label ~ ., data = train_data_generated, family = binomial)

#con la funzione summary() otteniamo un riassunto del modello
summary(log_model_generated)

Call:
glm(formula = label ~ ., family = binomial, data = train_data_generated)

Coefficients:
            Estimate Std. Error z value Pr(>|z|)    
(Intercept)  -4.2609     0.7060  -6.036 1.58e-09 ***
xTRUE         1.4844     0.5812   2.554   0.0106 *  
yTRUE        -1.1254     0.6383  -1.763   0.0779 .  
z             1.6121     0.2958   5.451 5.01e-08 ***
w             0.3825     0.1850   2.067   0.0387 *  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 200.161  on 199  degrees of freedom
Residual deviance:  91.126  on 195  degrees of freedom
AIC: 101.13

Number of Fisher Scoring iterations: 6

# Previsione sul test set
log_pred_generated_prob <- predict(log_model_generated, test_data_generated, type = "response")

# Convertiamo le probabilità in etichette
log_pred_generated <- ifelse(log_pred_generated_prob > 0.5, 1, -1)

# Matrice di confusione e accuracy

log_conf_matrix_generated <- confusionMatrix(as.factor(log_pred_generated), test_data_generated$label)
print(log_conf_matrix_generated)

Confusion Matrix and Statistics

          Reference
Prediction -1  1
        -1 38  4
        1   2  6
                                          
               Accuracy : 0.88            
                 95% CI : (0.7569, 0.9547)
    No Information Rate : 0.8             
    P-Value [Acc > NIR] : 0.1034          
                                          
                  Kappa : 0.5946          
                                          
 Mcnemar's Test P-Value : 0.6831          
                                          
            Sensitivity : 0.9500          
            Specificity : 0.6000          
         Pos Pred Value : 0.9048          
         Neg Pred Value : 0.7500          
             Prevalence : 0.8000          
         Detection Rate : 0.7600          
   Detection Prevalence : 0.8400          
      Balanced Accuracy : 0.7750          
                                          
       'Positive' Class : -1              
                                          

Plottiamo ora le curve ROC per i tre metodi (Naive Bayes, LDA, Regressione Logistica) per confronto.

log_pred_prob_generated <- predict(log_model_generated, test_data_generated, type = "response")

# Calcola la ROC curve ed esegui un plot usando ggroc (pacchetto pROC). Specifichiamo legacy.axes=TRUE per avere FPR nelle ascisse
roc_log_generated <- roc(test_data_generated$label, log_pred_prob_generated)
ggroc(roc_log_generated, legacy.axes = TRUE, color="green") +
  ggtitle("ROC: NB (rosso), LDA (blu) e RL (verde) - Dataset Generato") +
  xlab("False Positive Rate") +
  ylab("True Positive Rate") +
  geom_line(data = ggroc(roc_nb_e1071)$data, aes(x = 1 - specificity, y = sensitivity), color="red") +
  geom_line(data = ggroc(roc_lda_generated)$data, aes(x = 1 - specificity, y = sensitivity), color="blue") +
  theme(legend.position = "bottom") +
  guides(color = guide_legend(override.aes = list(size = 4)))

Consideriamo infine un modello di classificazione nel caso multi-classe. Per semplicità aggiungiamo una terza classe al dataset generato.

# Aggiungiamo una terza classe al dataset generato
x3 <- sample( 0:1, size=80, replace=TRUE, prob=c(0.2, 0.8))
y3 <- sample( 0:1, size=80, replace=TRUE, prob=c(0.5, 0.5))
z3 <- rnorm(80, mean=0, sd=2)
w3 <- rpois(80, lambda=5)

data_generated_multi <- data.frame("x" = c(x1, x2, x3), "y"=c(y1, y2, y3), "z"=c(z1, z2, z3), "w"=c(w1, w2, w3), "label"=as.factor(c(rep(-1, 200), rep(1,50), rep(0,80))))

# Suddivisione in train e test set
train_index_multi <- createDataPartition(data_generated_multi$label, p = 0.7, list = FALSE)
train_data_generated_multi <- data_generated_multi[train_index_multi, ]
test_data_generated_multi <- data_generated_multi[-train_index_multi, ]

ggplot(data_generated_multi, aes(x =z, y = w, color = label)) +
  geom_point() +
  labs(title = "Dataset Fittizio Multi-classe", x = "feature 3", y = "feature 4")

Applichiamo prima una strategia euristica one-vs-all per la regressione logistica multi-classe.

# Allenamento del modello di regressione logistica one-vs-all
log_models_multi <- list()

classes <- levels(train_data_generated_multi$label)
for (cls in classes) {
  binary_labels <- ifelse(train_data_generated_multi$label == cls, 1, 0)
  log_models_multi[[cls]] <- glm(binary_labels ~ ., data = train_data_generated_multi, family = binomial)
}
# Previsione sul test set
log_pred_multi <- sapply(classes, function(cls) {
  predict(log_models_multi[[cls]], test_data_generated_multi, type = "response")
})
log_pred_multi_labels <- apply(log_pred_multi, 1, function(probs) {
  classes[which.max(probs)]
})
# Matrice di confusione e accuracy
log_conf_matrix_multi <- confusionMatrix(as.factor(log_pred_multi_labels), test_data_generated_multi$label)
print(log_conf_matrix_multi)

Confusion Matrix and Statistics

          Reference
Prediction -1  0  1
        -1 60  0  0
        0   0 24  0
        1   0  0 15

Overall Statistics
                                     
               Accuracy : 1          
                 95% CI : (0.9634, 1)
    No Information Rate : 0.6061     
    P-Value [Acc > NIR] : < 2.2e-16  
                                     
                  Kappa : 1          
                                     
 Mcnemar's Test P-Value : NA         

Statistics by Class:

                     Class: -1 Class: 0 Class: 1
Sensitivity             1.0000   1.0000   1.0000
Specificity             1.0000   1.0000   1.0000
Pos Pred Value          1.0000   1.0000   1.0000
Neg Pred Value          1.0000   1.0000   1.0000
Prevalence              0.6061   0.2424   0.1515
Detection Rate          0.6061   0.2424   0.1515
Detection Prevalence    0.6061   0.2424   0.1515
Balanced Accuracy       1.0000   1.0000   1.0000

Applichiamo ora una regressione multinomiale usando il pacchetto caret.

# Allenamento del modello di regressione logistica multi classe
log_model_multi <- train(label ~ ., data = train_data_generated_multi, trControl = trainControl(verboseIter = FALSE),  trace = FALSE, method = "multinom")

# Previsione sul test set
log_pred_multi_caret <- predict(log_model_multi, test_data_generated_multi)
# Matrice di confusione e accuracy
log_conf_matrix_multi_caret <- confusionMatrix(log_pred_multi_caret, test_data_generated_multi$label)
print(log_conf_matrix_multi_caret)

Confusion Matrix and Statistics

          Reference
Prediction -1  0  1
        -1 58  5  4
        0   1 18  2
        1   1  1  9

Overall Statistics
                                          
               Accuracy : 0.8586          
                 95% CI : (0.7741, 0.9205)
    No Information Rate : 0.6061          
    P-Value [Acc > NIR] : 3.558e-08       
                                          
                  Kappa : 0.7289          
                                          
 Mcnemar's Test P-Value : 0.187           

Statistics by Class:

                     Class: -1 Class: 0 Class: 1
Sensitivity             0.9667   0.7500  0.60000
Specificity             0.7692   0.9600  0.97619
Pos Pred Value          0.8657   0.8571  0.81818
Neg Pred Value          0.9375   0.9231  0.93182
Prevalence              0.6061   0.2424  0.15152
Detection Rate          0.5859   0.1818  0.09091
Detection Prevalence    0.6768   0.2121  0.11111
Balanced Accuracy       0.8679   0.8550  0.78810

Esercizio: Applicare la regressione logistica al dataset German Credit Data e confrontare le performance con Naive Bayes e LDA viste in precedenza.

Esercizi

1. Generare una tabella con 30 righe e 8 colonne con dati casuali. Aggiungere una quinta colonna contenente una classificazione binaria, generata anche essa casualmente. Calcolare la discrepanza tra la classificazione originaria e il risultato di una classificazione logistica. Se si ripete questo esperimento per 1000 volte, calcolare valore medio e deviazione standard della discrepanza.

2. Generare una tabella A con 170 righe e 7 colonne, in modo che l’ultima colonna contenga valori binari relativi all’appartenenza a una o l’altra di due classi. Generare una seconda tabella B con 5 righe e 6 colonne. Implementare una classificazione degli individui nella tabella B applicando alla tabella A la classificazione logistica e la analisi discriminante lineare. Confrontare i risultati dei due metodi.

3. Confrontare i metodi di regressione logistica e analisi discriminante lineare e quadratica sul dataset titanic usando l’autovalidazione, per stabilire il metodo migliore.

4. Determinare il migliore metodo di classificazione (tra quelli studiati nel corso) per il dataset Pima.tr (pacchetto MASS). Usare poi i modelli per classificare i nuovi dati contenuti nel dataset Pima.te e verificare se il miglior modello stabilito nella prima parte sia effettivamente il più efficace.

Nota sui blocchi

Potete creare un blocco iniziale e non visualizzare il codice di tutti i blocchi in un R notebook usande il comando: knitr::opts_chunk$set(warning = FALSE, message = FALSE, echo=FALSE)