1 Introduzione

Con diagnosi intendiamo un processo in cui un insieme di osservazioni relative ad un paziente (sintomi, esami strumentali quali livello di glicemia, globuli rossi etc…) sono valutate per determinare con un certo grado di confidenza la patologia che lo affligge allo scopo di prescrivere la terapia più adatta.

Spesso per determinare con un elevato grado di sicurezza la patologia è necessario effettuare indagini strumentali costose e/o invasive per il paziente. Per cui è di interesse verificare se è possibile diagnosticare la patologia mediante esami alternativi poco costosi.

Per valutare l’efficacia diagnostica dei metodi alternativi è necessario effettuare uno studio per esaminare le osservazioni (indagini) alternative su una popolazione di pazienti su cui è già stata diagnosticata la patologia con sicurezza.

L’analisi delle osservazioni per trarre le opportune conclusioni si deve effettuare con i metodi della statistica. Si deve tenere presente che i metodi della statistica consentono di fare molto di più di quanto accennato sopra. La statistica è stata grandemente sviluppata da scienziati quali Fisher e Pearson ed altri nella prima metà del 1900. Noi non entreremo nel dettaglio dei metodi statistici se non con uno sguardo molto superficiale per porre le basi minime per la comprensione delle basi dell’intelligenza artificiale.

Negli ultimi 70 anni con l’avvento dei computer si è sviluppato un ramo della matematica denominato machine learning che ha molto in comune con i metodi statistici ma il cui obiettivo e le metodologie sono leggermente diversi.

Con il perfezionamento delle tecniche teoriche e matematiche e con il crescere della potenza di calcolo dei computer, tale settore si è sviluppato enormemente diventando l’attuale intelligenza artificiale. L’obiettivo di questo mini-corso è di porre le basi per la comprensione di questo settore.

1.1 Il problema fondamentale del machine learning

In termini molto semplici uno dei task principali del machine learning è di assegnare un certo individuo (e.g. un paziente) ad una certa classe (e.g. malato/sano). L’individuo è caratterizzato mediante un insieme valori numerici (per esempio i risultati delle indagini strumentali di cui sopra). Tali valori numerici sono detti features. La macchina cioè il computer può imparare ad assegnare un individuo alla classe giusta se gli vengono presentati un certo numero di esempi cioè di features insieme alle corrispondenti label (etichette, classe di appartenenza) corrette.

Con il linguaggio matematico se \(\mathbf{x}_k = [f_1,f_2,...f_N]\) sono \(N\) features che rappresentano l’individuo \(k-\)simo e se \(\omega_k \in {\omega_1, \ldots, \omega_K}\) è la sua classe tra le \(K\) possibili, la macchina impara ad assegnare \(\mathbf{x}_k\) alla classe \(\omega_k\). L’apprendimento avviene presentando un gran numero di esemplari.

Facciamo un esempio. In un albergo vi sono due congressi: uno di giocatori di basket ed uno di giocatori di ping-ping. Il portiere deve indicare ad ogni cliente la sala del suo congresso. Dopo che gli sono state presentati un certo numero di giocatori dei due tipi comincia a capire che gli individui alti e pesanti sono giocatori di basket mentre quelli piccoli e leggeri sono giocatori di ping-pong. Il portiere ha imparato ad assegnare la classe (basket o ping-pong) ad un giocatore usando le features (peso, altezza) dopo avere osservato un certo numero di individui (training).

Dalla figura emergono alcune considerazioni:

  • il trend delle due popolazioni sembra essere quello descritto (alti e pesanti i giocatori di basket, bassi e leggeri quelli di ping-pong): sembra esistere una linea grossolana di demarcazione, tuttavia alcuni individui si trovano border-line ed altri addirittura si trovano mescolalta con individui dell’altra popolazione

  • non potendo osservare tutti gli individui come facciamo a determinare il migliore algoritmo per separare le due popolazioni ?

  • una volta determinato il migliore algoritmo possibile con i dati a nostra disposizione (training-set) possiamo valutare la probabilità di errore di classificazione su un individuo nuovo mai visto prima ?

Di come sia possibile dare una risposta a queste domande ci occuperemo in questa lezione.

1.2 I software per l’intelligenza artificiale

Eviteremo i dettagli matematici del problema e ci concentreremo sin dall’inizio su come gestire questi problemi nella pratica. Pertanto sin dall’inizio useremo un software statistico: R.

Esistono software commerciali ed open source per fare analisi statistiche ed intelligenza artificiale. Alcuni software si basano su interfacce grafiche, altri su comandi da scrivere sulla tastiera (riga di comando). In generale, anche se l’interfaccia grafica risulta apparentemente semplice per l’utente novizio, in realtà la riga di comando alla lunga è più versatile e potente.

I principali software open-source (licenza d’uso gratuita) sono:

  • R:
    • orientato alla analisi statistica sin dall’inizio;
    • facile da usare
    • non richiede conoscenze spinte di programmazione, ma se si programma si possono fare cose molto sofisticate
    • strutturato con packages sviluppati dalla comunità degli statistici che contengono funzioni per fare quasi qualsiasi cosa
    • è integrato con Keras e TensorFlow che sono sviluppati in Python (vedi dopo)
    • grafici eccellenti
    • è basato su comandi scritti da tastiera piuttosto che su interfaccia grafica
  • Python
    • richiede conoscenze spinte di programmazione
    • strutturato in packages sviluppati dalla comunità
    • alcuni packages come Keras e TensorFlow sono diventati lo standard per la intelligenza artificiale
    • è basato su comandi scritti da tastiera piuttosto che su interfaccia grafica

Alcuni software commerciali

Su WIKIPEDIA c’è una lista pressocchè completa di software commerciali e free per la statistica. Un confronto tra le possibilità offerte dai vari software si trova su questa pagina

In questo mini-corso abbiamo scelto R perchè :

  • è molto versatile
  • anche se richiede di scrivere comandi da tastiera l’apprendimento è molto veloce e le possibilità sono pressocchè infinite
  • è free

1.3 Bibliografia

Sin da ora voglio richiamare lo studente sulla necessità di approfondire autonomamente gli argomenti trattai in questa lezione che assolutamente non può essere esaustiva data la vastità della questione trattata.

Esistono numerosi testi di introduzione alla statistica, ad R e alle tecniche di classificazione. Ci limitiamo a suggerirne qualcuno in cui lo studente potrà approfondire gli argomenti visti in questo mini-corso.

  • Kabacoff, R in Action, Manning Publications

    • eccellente introduzione ad R applicato alla statistica
    • richiede poche conoscenze teoriche di base

  • James, Witten, Hastie, Tibshirani, An Introduction to Statistical Learning, https://www.statlearning.com/

  • Hastie, Tibshirani, The Elements of Statistical Learning, https://hastie.su.domains/Papers/ESLII.pdf

  • Venable and Ripley, Modern Applied Statistics with S, Springer, 2002

    • applicazioni di statistica in S (un predecessore di R) e R
    • è orientata all’uso del software e richiede poche conoscenze teoriche
    • è più avanzato rispetto ad “R in Action”

  • Kuhn, Johnson, Applied Predictive Modelling

    • approccio pratico al machine learning

  • Bishop , Pattern recognition and Machine Learning, Springer, 2006

    • una introduzione completa al pattern recognition
    • richiede conscenze teoriche avanzate

  • Bishop, Neural networks and pattern recognition, Oxford press, 1995

    • introduzione completa alle reti neurali
    • conoscenze teoriche avanzate

Un elenco di libri su R è disponile al seguente link

Uno dei pacchetti più importanti per la classificazione è il pacchetto Classification And REgression Training CARET

2 Tipi di Classificatori

2.1 Ottimizzazione

Nei classificatori è importante il concetto di ottimizzazione o di minimizzazione di qualche criterio rispetto a qualche parametro: per esempio la distanza di una retta rispetto ad alcuni punti prefissati.

In genere, come si vede nella figura seguente il minimo può essere globale o locale. Per cui spesso il punto di partenza nella ricerca del minimo ci può far trovare un minimo locale o globale.

Il punto di partenza in genere viene scelta in maniera casuale. Questo è il motivo per cui risultati della classificazione possono essere diversi tra vari computer. A meno che non si utilizza un meccanismo di impostazione della casualità (i computer non possono generare numeri veramente casuali, ma solo pseudo-casuali con particolari algoritmi molto softisticati).

2.2 Dati di esempio

Per vedere come funzionano i vari tipi di classificatori li proveremo sui medesimi dati di esempio generati artificialmente.

Due classi, \(N=100\) individui per classe. Il numero di features in genere si indica con \(p\). In questo caso \(p=2\) le features siano \(x_1\) e \(x_2\). Partiamo da distribuzioni normali con medie differenti.

N<-100 # per ciascuna classe
set.seed(12345)# in modo che i risultati siano ripetibili su altri computer

x1 <- rnorm(N)
x2 <- rnorm(N)

dati <- data.frame(x1,x2,cl=rep(0,N))

x1 <- rnorm(N,mean = 2)
x2 <- rnorm(N,mean = 2)

dati <- rbind(dati,cbind(x1,x2,cl=rep(1,N)))
dati$cl <- factor(dati$cl)

plot(x2 ~ x1,
     data=dati,
     col=c('green','red')[dati$cl],
     pch=20, cex=2, asp=1)

2.3 Linear Discriminant Analysis

Alcuni classificatori molto semplici sono quelli lineari che tracciano una retta (oppure un piano nel caso di 3 features, oppure un iper-piano nel caso di \(p\) > 3 features) di separazione tra le classi.

LDA (Linear Discriminant Analysis) è un tipo di classificatore lineare.

Il criterio con cui viene scelta la retta di separazione è basato sulla proiezione dei dati su una retta tale che le medie delle popolazioni siano più distanti possibile. [il dettaglio matematico è complesso e non vi entriamo]

library(MASS)

modLDA <- lda(cl ~ ., dati)

previsioni<-predict(modLDA,newdata = dati[,c('x1','x2')])$class

table(predizioni=previsioni,classe=dati$cl)
##           classe
## predizioni  0  1
##          0 86  6
##          1 14 94
# tracciamo la curva di separazione
# cioè il luogo dei punti dove la probabilità di predizione è 0.5
x1l<-seq(-3,5,.1)
x2l<-seq(-3,5,.1)

ndati<-expand.grid(x1=x1l,x2=x2l)

previsioni <- predict(modLDA,newdata = ndati)$posterior[,1]

plot(x2 ~ x1,data=dati,
     col=c('green','red')[dati$cl],pch=20, cex=2,
  main='LDA', asp=1)
contour(x=x1l,y=x2l,
        z=matrix(previsioni,nrow=length(x1l),byrow = FALSE),
        levels=0.5,add=TRUE)

2.4 Support Vector Machine

library(e1071)

mod <- svm(cl ~ ., data=dati,probability=TRUE)

previsioni<-predict(mod, newdata = dati[,c('x1','x2')])

table(predizioni=previsioni,classe=dati$cl)
##           classe
## predizioni  0  1
##          0 94  6
##          1  6 94
previsioni <- predict(mod,newdata = ndati)

plot(x2 ~ x1,data=dati,
     col=c('green','red')[dati$cl],pch=20, cex=2,
  main='SVM', asp=1)
contour(x=x1l,y=x2l,
        z=matrix(previsioni,nrow=length(x1l),byrow = FALSE),
        levels=0.5,add=TRUE,drawlabels = FALSE)

2.5 Neural Networks

library(nnet)
set.seed(1111)
mod <- nnet(cl ~ ., data=dati, size=3, trace=FALSE)

previsioni<-predict(mod, newdata = dati[,c('x1','x2')],type='class') 

table(predizioni=previsioni,classe=dati$cl)
##           classe
## predizioni  0  1
##          0 94 10
##          1  6 90
previsioni <- predict(mod,newdata = ndati)

plot(x2 ~ x1,data=dati,
     col=c('green','red')[dati$cl],pch=20, cex=2,
  main='NEURAL-NET', asp=1)
contour(x=x1l,y=x2l,
        z=matrix(previsioni,nrow=length(x1l),byrow = FALSE),
        levels=0.5,add=TRUE,drawlabels = FALSE)

library(NeuralNetTools)
plotnet(mod)

2.6 K-nearest neighbours

Un altro approccio intuitivo alla classificazione consiste nell’assegnare a ciascun individuo la classe dell’individuo (o degli individui) più vicini

library(caret)

mod<-knn3(x=dati[,c('x1','x2')],
          y=factor(dati$cl),
          k=1)

previsioni<-predict(mod,newdata = dati[,c('x1','x2')],type='class')

table(predizioni=previsioni,classe=dati$cl)
##           classe
## predizioni   0   1
##          0 100   0
##          1   0 100
previsioni<-predict(mod,newdata = data.frame(ndati))[,1]

plot(x2 ~ x1,data=dati,
     col=c('green','red')[dati$cl],pch=20, cex=2,
  main='K-NN', asp=1)

contour(x=x1l,y=x2l,
        z=matrix(previsioni,nrow=length(x1l),byrow = FALSE),
        levels=0.5,add=TRUE,drawlabels = FALSE)

Si vede che la separazione è mlto frastagliata ed in generale non sembra plausibile (tra l’altro l’errore di predizione è 0).

Se consideriamo un numero \(k\) di vicini e assegnamo la classe in base alla classe maggioritaria tra i \(k\) avremo una separazione più plausibile

mod<-knn3(x=dati[,c('x1','x2')],
          y=factor(dati$cl),
          k=5) # cambiamo il numero di 'vicini'

previsioni<-predict(mod,newdata = dati[,c('x1','x2')],type='class')

table(predizioni=previsioni,classe=dati$cl)
##           classe
## predizioni  0  1
##          0 96  6
##          1  4 94
previsioni<-predict(mod,newdata = data.frame(ndati))[,1]

plot(x2 ~ x1,data=dati,
     col=c('green','red')[dati$cl],pch=20, cex=2,
  main='K-NN k=5', asp=1)

contour(x=x1l,y=x2l,
        z=matrix(previsioni,nrow=length(x1l),byrow = FALSE),
        levels=0.5,add=TRUE,drawlabels = FALSE)

2.7 CART

Classification And Regression Trees, sono delle tecniche di ripartizione dicotomica ricorsiva.

Sono molto semplici da interpretare pechè prevedono sono delle domande del tipo : \(x > soglia\) ?

library(rpart)

mod <- rpart(cl ~ ., data=dati, method = 'class')

previsioni<-predict(mod, newdata = dati[,c('x1','x2')],type='class')

table(predizioni=previsioni,classe=dati$cl)
##           classe
## predizioni  0  1
##          0 93  6
##          1  7 94
previsioni <- predict(mod,newdata = ndati)[,2]

plot(x2 ~ x1, data=dati, col=c('green','red')[dati$cl],
  pch=20, cex=2,main='CART', asp=1)

contour(x=x1l,y=x2l,
        z=matrix(previsioni,nrow=length(x1l),byrow = FALSE),
        levels=0.5,add=TRUE,drawlabels = FALSE)

library(rpart.plot)
rpart.plot(mod, type=1, extra=1)

2.8 General Linear Model (Logistic Regression)

mod <- glm(cl ~ ., data=dati,family='binomial')

previsioni<-factor(as.numeric(predict(mod, newdata = dati[,c('x1','x2')],type='response') >0.5))


table(predizioni=previsioni,classe=dati$cl)
##           classe
## predizioni  0  1
##          0 92  6
##          1  8 94
previsioni <- predict(mod,newdata = ndati,type='response')

plot(x2 ~ x1, data=dati, col=c('green','red')[dati$cl],
  pch=20, cex=2,main='GLM', asp=1)

contour(x=x1l,y=x2l,
        z=matrix(previsioni,nrow=length(x1l),byrow = FALSE),
        levels=0.5,add=TRUE,drawlabels = FALSE)

3 Misure di Performance

Diamo alcune definizioni.

  • training-set: L’apprendimento coinvolge un certo numero di individui di esempio.

  • validation-set: per alcuni tipi di classificatori è necessario ottimizzare alcuni parametri per ottimizzare l’addetramento, a tale scopo si usa il validation-set.

  • test-set: Per misurare le prestazioni si usa un set di individui che la macchina non ha mai visto.

  • over-fit: Una macchina può diventare molto brava a riconoscere gli individui del training-set ma commettere errori quando cerca di classificare nuovi individui mai visti prima.

  • error-rate: Quando la macchina è addestrata, si cerca di misurare la capacità di generalizzare il riconoscimento. Si deve misurare il numero di errori commessi su un insieme di dati mai visto durane l’addestramento.

Allo scopo di calcolare le prestazioni si premettono le seguenti ulteriori definizioni relative al caso in cui le classi sono solo due.

La nomenclatura deriva dal caso di un test diagnostico in cui se una persona è POSITIVA al test allora è MALATA.

Un tale test può avere degli errori come si vede nella tabella seguente: per esempio se un soggetto SANO risulta POSITIVO allora è un FALSO POSITIVO cioè FP; un soggetto MALATO ma che risulta NEGATIVO al test è un FALSO NEGATIVO cioè FN.

Un soggetto MALATO che risulta POSITIVO al test è un TRUE POSITIVE (TP), mentre un soggetto SANO NEGATIVO è un TRUE NEGATIVE.

MALATI SANI
POSITIVI TP FP
NEGATIVI FN TN

Con queste definizioni preliminari le principali misure di prestazione sono:

  • \(\textrm{Accuracy} = \frac{Corretti \ riconoscimenti}{Tutti} = \frac{TP+TN }{TP+FP+TN+FN}\)

  • \(\textrm{Sensitivity} = \frac{Veri \ positivi}{Totale \ Malati} = \frac{TP}{TP+FN}\)

  • \(\textrm{Specificity} = \frac{Veri \ negativi}{Totale \ Sani} = \frac{TN}{TN+FP}\)

La nomenclatura del settore dipende anche dal contesto: fuori dal campo medico per esempio la Sensitivity viene chiamata Recall. Sulla pagina web del package CARET si trovano varie definizioni.

Nel caso multi-classe la nomenclatura è simile (TP, FP, Accuracy, etc…) ma si aggiunge la Matrice di Confusione (di cui vedremo esempi nel seguito) che generalizza i concetti precedenti. Sul sito di WIKIPEDIA c’è un buon riassunto di tutte le misure di prestazione.

3.1 Esempi di valutazione prestazioni

Suddividiamo l’insieme complessivo dei dati in un training e un test set. Addestriamo sul training e verifichiamo gli errori sul test. Effettuiamo l’intera procedura sul data set di esempio.

library(MASS)
library(e1071)
library(caret)
library(nnet)
library(rpart)

N<-100
set.seed(123456)
dati <- data.frame(x1=c(rnorm(N),rnorm(N,mean=1)),
  x2=c(rnorm(N),rnorm(N,mean=1)),
  cl=factor(rep(c(1,2),each=N)))

set.seed(123456)
indiciTrain <- sample(1:(2*N),size=0.75*2*N)

training <- dati[indiciTrain,]
test <- dati[-indiciTrain,]

plot(x2 ~ x1, data=training, col=c(1:2)[training$cl],pch=20)
points(x2 ~ x1, data=test,col=c(3:4)[test$cl],pch=22)
legend(x='topleft',legend=c('train','test'),pch=c(20,22))

Addestramo vari classificatori

modLDA <- lda(cl ~ ., data = training)
prevLDA <- predict(modLDA, newdata = test)$class
confusionMatrix(data = prevLDA,reference = test$cl,positive='1')
## Confusion Matrix and Statistics
## 
##           Reference
## Prediction  1  2
##          1 17  4
##          2  6 23
##                                           
##                Accuracy : 0.8             
##                  95% CI : (0.6628, 0.8997)
##     No Information Rate : 0.54            
##     P-Value [Acc > NIR] : 0.0001186       
##                                           
##                   Kappa : 0.5948          
##                                           
##  Mcnemar's Test P-Value : 0.7518296       
##                                           
##             Sensitivity : 0.7391          
##             Specificity : 0.8519          
##          Pos Pred Value : 0.8095          
##          Neg Pred Value : 0.7931          
##              Prevalence : 0.4600          
##          Detection Rate : 0.3400          
##    Detection Prevalence : 0.4200          
##       Balanced Accuracy : 0.7955          
##                                           
##        'Positive' Class : 1               
## 
modSVM <- svm(cl ~., data=training)
prevSVM<-predict(modSVM,newdata = test)
confusionMatrix(data = prevSVM,reference = test$cl,positive='1')
## Confusion Matrix and Statistics
## 
##           Reference
## Prediction  1  2
##          1 18  8
##          2  5 19
##                                           
##                Accuracy : 0.74            
##                  95% CI : (0.5966, 0.8537)
##     No Information Rate : 0.54            
##     P-Value [Acc > NIR] : 0.002962        
##                                           
##                   Kappa : 0.4817          
##                                           
##  Mcnemar's Test P-Value : 0.579100        
##                                           
##             Sensitivity : 0.7826          
##             Specificity : 0.7037          
##          Pos Pred Value : 0.6923          
##          Neg Pred Value : 0.7917          
##              Prevalence : 0.4600          
##          Detection Rate : 0.3600          
##    Detection Prevalence : 0.5200          
##       Balanced Accuracy : 0.7432          
##                                           
##        'Positive' Class : 1               
## 
modCART <- rpart(cl ~., data=training)
prevCART<- factor(as.numeric(predict(modCART,newdata = test)[,2]>.5)+1)
confusionMatrix(data = prevCART,reference = test$cl,positive='1')
## Confusion Matrix and Statistics
## 
##           Reference
## Prediction  1  2
##          1 16  7
##          2  7 20
##                                           
##                Accuracy : 0.72            
##                  95% CI : (0.5751, 0.8377)
##     No Information Rate : 0.54            
##     P-Value [Acc > NIR] : 0.007101        
##                                           
##                   Kappa : 0.4364          
##                                           
##  Mcnemar's Test P-Value : 1.000000        
##                                           
##             Sensitivity : 0.6957          
##             Specificity : 0.7407          
##          Pos Pred Value : 0.6957          
##          Neg Pred Value : 0.7407          
##              Prevalence : 0.4600          
##          Detection Rate : 0.3200          
##    Detection Prevalence : 0.4600          
##       Balanced Accuracy : 0.7182          
##                                           
##        'Positive' Class : 1               
## 
modNNET <- nnet(cl ~., data=training,size=3,trace=FALSE)
prevNNET<- factor(as.numeric(predict(modNNET,newdata = test)>.5)+1)
confusionMatrix(data = prevNNET,reference = test$cl,positive='1')
## Confusion Matrix and Statistics
## 
##           Reference
## Prediction  1  2
##          1 15  5
##          2  8 22
##                                           
##                Accuracy : 0.74            
##                  95% CI : (0.5966, 0.8537)
##     No Information Rate : 0.54            
##     P-Value [Acc > NIR] : 0.002962        
##                                           
##                   Kappa : 0.4715          
##                                           
##  Mcnemar's Test P-Value : 0.579100        
##                                           
##             Sensitivity : 0.6522          
##             Specificity : 0.8148          
##          Pos Pred Value : 0.7500          
##          Neg Pred Value : 0.7333          
##              Prevalence : 0.4600          
##          Detection Rate : 0.3000          
##    Detection Prevalence : 0.4000          
##       Balanced Accuracy : 0.7335          
##                                           
##        'Positive' Class : 1               
##

3.2 ROC analysis

Nel caso di due classi c’è una modalità grafica di rappresentare le prestazioni. La Receiver Operating Characteristics è molto usata per le prestazioni di un problem a due classi (binario).

La curva migliore deve tendere verso lo spigolo in alto a sinistra. l’area sotto la curva (AUC) tende ad 1.

library(pROC)
previsioni<-predict(modNNET, newdata = test)
plot(myroc <- roc(test$cl, as.numeric(previsioni)))
text(.5,.1,paste('AUC',format(myroc$auc, digits=2)) )

4 Metodi di resampling

Uno dei problemi fondamentali della statistica (e delle varie discipline collegate di machine learning e intelligenza artificiale) consiste nel determinare il valore di una o più grandezze che caratterizzano una popolazione [per esempio la altezza media della popolazione italiana].

In genere è impossibile esaminare tutta la popolazione per cui si usa un campione di individui estratti a caso dalla popolazione sulla base dei quali si stima il valore cercato (nell’esempio la altezza media).

Si pone il problema di valutare l’accuratezza della stima basata sul campione. I metodi di resampling servono per aver un idea della accuratezza.

Senza entrare in troppi dettagli matematici nella prossima slide esaminiamo il problema della stima della media di una popolazione.

  • consideriamo una variabile aleatoria (popolazione) \(X\)
  • sia \(\mu_X\) la media e \(\sigma^2_X\) la varianza della popolazione
  • estraiamo un campione di \(N\) elementi \(\{ x_1, x_2, \ldots, x_N \}\)
  • la media del campione \[\begin{equation} \bar{x}=(1/N) \sum_i x_i \end{equation}\] è una variabile aleatoria
  • la expectation (valore atteso o medio) della media del campione è: \[\begin{equation} E[\bar{x}] = \mu_X \end{equation}\]
  • lo standard error è la deviazione standard della media del campione \[\begin{equation} se(\bar{x}) = \sigma_{\bar{x}} = \sigma_X / \sqrt{N} \end{equation}\]

4.1 Metodo bootstrap

La media campionaria è rappresentativa della media della popolazione ?

 Npopolazione<-10000
 set.seed(123456)
 popolazione <- rnorm(Npopolazione, mean=1, sd=1) # popolazione
 Nsamp <- 10 # dimensione campione
 campione <- popolazione[sample(1:Npopolazione,size=Nsamp)]
 mc <- mean(campione) # media del campione
 mpv<-mean(popolazione)
hist(popolazione,30,main='distribuzione popolazione')
rug(campione,col='darkgreen',lwd=1.5)
abline(v=mc,col='red',lwd=2)
abline(v=mpv,col='blue',lwd=2)
legend(x='topleft',legend=c('media del campione','media vera popolazione','campione'),pch=15,col=c('red','blue','darkgreen'))

# se avessimo estratto altri campioni ?
Ncampioni<-200
m <- matrix(NA,nrow=Nsamp,Ncampioni)
for(k in 1:Ncampioni){
  m[,k]<-  popolazione[sample(1:Npopolazione,size=Nsamp)]
}
m<-apply(m,MARGIN=2,mean)# distribuzione media campionaria
mm<-mean(m) # expectaion della media campione = stima media popolazione
ms<-sd(m) # stima dello SE = sigma/sqrt(Nsamp)
hist(m,30, xlim=c(-1,3),
  main='distribuzione media campionaria')
abline(v=mc,col='red',lwd=2)
abline(v=mpv, col='blue',lwd=2)
abline(v=mm,col='cyan',lwd=2)
legend(x='topleft',legend=c('media del campione','media popolazione','media della distribuzione media campionaria'),pch=15,col=c('red','blue','cyan'))

La media del campione quindi può essere polarizzata (si parla di bias).

Usando il bootstrap posso cercare di stimare l’accuratezza, cioè lo standard error = deviazione standard della distribuzione della media campionaria. Infatti la distribuzione della media campionaria è simile (dal punto di vista della dispersione, non del valor medio!) alla distribuzione della media campionaria con boostrap.

library(boot)
res<-boot(campione,statistic=function(a,ind){mean(a[ind])},R=300)
layout(matrix(1:2,1,2))
hist(m,30,main='distribuzione media campionaria',sub='partendo dalla popolazione',freq = FALSE)
hist(res$t,30,main='distribuzione media campionaria',sub='usando il bootstrap sul campione',freq = FALSE)

4.2 Metodi per valutare l’errore di classificazione

Attualmente i metodi usati per valutare l’errore di classificazione sono basti su idee simili a quelle del bootstrap.

  • Leave One Group Out (LOGCV)
  • Repeated LGOCV
  • Leave One Out (LOOCV)
  • Cross-Validation (CV)
  • Repeated Cross-Validation (repeated-CV)

Le varie modalità sono semplificate usando il package CARET.

4.2.1 Leave Group Out (LGOCV)

L’insieme dei dati disponibili si divide in due parti: training e test. Solitamente il training è il 75% dei dati, il restante è il test.

library(caret)

myctrl <- trainControl(method='LGOCV',p=0.75, number = 1)

modLDA<-train(x = dati[,-3], y=dati[,3],
  method='lda',
  trControl = myctrl)

modSVM<-train(x = dati[,-3], y=dati[,3],
  method='svmLinear2',
  trControl = myctrl)

modNNET<-train(x = dati[,-3], y=dati[,3],
  method='nnet',
  trControl = myctrl, trace=FALSE)

modKNN<-train(x = dati[,-3], y=dati[,3],
  method='knn',
  trControl = myctrl)



print(modLDA$results)
##   parameter Accuracy Kappa AccuracySD KappaSD
## 1      none     0.64  0.28         NA      NA
print(modSVM$results)
##   cost Accuracy Kappa AccuracySD KappaSD
## 1 0.25     0.76  0.52         NA      NA
## 2 0.50     0.76  0.52         NA      NA
## 3 1.00     0.76  0.52         NA      NA
print(modNNET$results)
##   size decay Accuracy Kappa AccuracySD KappaSD
## 1    1 0e+00     0.74  0.48         NA      NA
## 2    1 1e-04     0.74  0.48         NA      NA
## 3    1 1e-01     0.70  0.40         NA      NA
## 4    3 0e+00     0.72  0.44         NA      NA
## 5    3 1e-04     0.72  0.44         NA      NA
## 6    3 1e-01     0.74  0.48         NA      NA
## 7    5 0e+00     0.70  0.40         NA      NA
## 8    5 1e-04     0.76  0.52         NA      NA
## 9    5 1e-01     0.72  0.44         NA      NA
print(modKNN$results)
##   k Accuracy Kappa AccuracySD KappaSD
## 1 5     0.72  0.44         NA      NA
## 2 7     0.70  0.40         NA      NA
## 3 9     0.66  0.32         NA      NA

4.2.2 Repeated LGOCV

La divisione in due parti viene ripetuta molte volte, ogni volta con una selezione casuale del training e del tes set. La ripetizione ha lo scopo di valutare la accuratezza (segundo l’idea del bootstrap)

myctrl <- trainControl(method='LGOCV',p=0.75, number = 50)

modLDA<-train(x = dati[,-3], y=dati[,3],
  method='lda',
  trControl = myctrl)

modSVM<-train(x = dati[,-3], y=dati[,3],
  method='svmLinear2',
  trControl = myctrl)

modNNET<-train(x = dati[,-3], y=dati[,3],
  method='nnet',
  trControl = myctrl, trace=FALSE)

modKNN<-train(x = dati[,-3], y=dati[,3],
  method='knn',
  trControl = myctrl)



print(modLDA$results)
##   parameter Accuracy  Kappa AccuracySD    KappaSD
## 1      none   0.7832 0.5664 0.04487807 0.08975613
print(modSVM$results)
##   cost Accuracy  Kappa AccuracySD   KappaSD
## 1 0.25   0.7756 0.5512 0.05485268 0.1097054
## 2 0.50   0.7772 0.5544 0.05642550 0.1128510
## 3 1.00   0.7752 0.5504 0.05489285 0.1097857
print(modNNET$results)
##   size decay Accuracy  Kappa AccuracySD    KappaSD
## 1    1 0e+00   0.7684 0.5368 0.05223026 0.10446052
## 2    1 1e-04   0.7684 0.5368 0.05223026 0.10446052
## 3    1 1e-01   0.7732 0.5464 0.04958769 0.09917538
## 4    3 0e+00   0.7388 0.4776 0.04821889 0.09643778
## 5    3 1e-04   0.7336 0.4672 0.05382360 0.10764719
## 6    3 1e-01   0.7660 0.5320 0.05162660 0.10325321
## 7    5 0e+00   0.7228 0.4456 0.05584381 0.11168761
## 8    5 1e-04   0.7184 0.4368 0.05448685 0.10897369
## 9    5 1e-01   0.7688 0.5376 0.05081539 0.10163078
print(modKNN$results)
##   k Accuracy  Kappa AccuracySD    KappaSD
## 1 5   0.7316 0.4632 0.04917109 0.09834218
## 2 7   0.7396 0.4792 0.05735532 0.11471064
## 3 9   0.7396 0.4792 0.04526385 0.09052770

4.2.2.1 Come varia l’errore e la varianza al variare del numero di ripetizioni ?

Prendiamo il classificatore LDA come esempio e vediamo come varia la performance: si vede che in pratica dopo 100-200 ripetizioni le performance si assestano su un certo valore sia in termini di accuratezza sia in termini di variabilità. Nel grafico abbiamo rappresentato in nero la variabilità come 1 deviazione standard che corrisponde a circa il 65% dei casi, volendo dare una rappresentazione più affidabile si dovrebbe usare 2 deviazioni standard (in rosso)

ripetizioni<-c(1,10,20,30,40,50,100,200,300,400,500)
risultati <- rep(NA,length(ripetizioni))
variabilita <- risultati
conta <- 1
for (k in ripetizioni){
myctrl <- trainControl(method='LGOCV',p=0.75, number = k)

mod<-train(x = dati[,-3], y=dati[,3],
  method='lda',
  trControl = myctrl)
  
risultati[conta] <- mod$results$Accuracy
variabilita[conta] <- mod$results$AccuracySD
conta<-conta+1
}

plot(ripetizioni,risultati,type='l',ylim=c(.5,1))
arrows(ripetizioni,risultati,ripetizioni,risultati+variabilita, angle=90,length=.1)
arrows(ripetizioni,risultati,ripetizioni,risultati-variabilita, angle=90,length=.1)
arrows(ripetizioni,risultati,ripetizioni,risultati+2*variabilita, angle=90,length=.1,col='red')
arrows(ripetizioni,risultati,ripetizioni,risultati-2*variabilita, angle=90,length=.1,col='red')

4.2.3 Leave One Out (LOOCV)

Per addestrare la macchia si usa tutto il training set tranne 1 elemento, l’errore si valuta sul singolo elemento rimasto fuori; la preocedura si ripete per ogni elemento.

Se il training set è molto numeroso il LOOCV può impiegare molto tempo. L’accuratezza valutata con LOOCV sembra essere meno polarizzata (dato che si usa quasi tutto il traning set) ma possiede una ampia variabilità (essendo le stime tutte correlate tra loro la loro media ha una ampia deviazione standard).

myctrl <- trainControl(method='LOOCV')

modLDA<-train(x = dati[,-3], y=dati[,3],
  method='lda',
  trControl = myctrl)

modSVM<-train(x = dati[,-3], y=dati[,3],
  method='svmLinear2',
  trControl = myctrl)

modNNET<-train(x = dati[,-3], y=dati[,3],
  method='nnet',
  trControl = myctrl, trace=FALSE)

modKNN<-train(x = dati[,-3], y=dati[,3],
  method='knn',
  trControl = myctrl)



print(modLDA$results)
##   parameter Accuracy Kappa
## 1      none    0.785  0.57
print(modSVM$results)
##   cost Accuracy Kappa
## 1 0.25    0.770  0.54
## 2 0.50    0.755  0.51
## 3 1.00    0.755  0.51
print(modNNET$results)
##   size decay Accuracy Kappa
## 1    1 0e+00    0.755  0.51
## 2    1 1e-04    0.755  0.51
## 3    1 1e-01    0.760  0.52
## 4    3 0e+00    0.745  0.49
## 5    3 1e-04    0.740  0.48
## 6    3 1e-01    0.760  0.52
## 7    5 0e+00    0.745  0.49
## 8    5 1e-04    0.750  0.50
## 9    5 1e-01    0.765  0.53
print(modKNN$results)
##   k Accuracy Kappa
## 1 5    0.745  0.49
## 2 7    0.755  0.51
## 3 9    0.740  0.48

4.2.4 Cross Validation (CV)

Il training-set ed il test-set sono scelti all’interno del data-set. In genere per addestrare sufficientemente un algoritmo è necessaria una grande training-set, per cui i dati del test-set sono pochi e potrebbe essere difficoltoso valutare le prestazioni con affidabilità.

Per superare questo limite è stata inventata la tecnica detta Cross-Validation in cui il data-set complessivo viene suddiviso in tanti sottoinsiemi detti fold. Poniamo che siano stati ottenuti K fold, la procedura di addestramento-valutazione funziona nel modo seguente: la macchina viene addestrata sui primi K-1 fold e l’ultimo viene usato per valutare le prestazioni (sostanzialmente contare gli errori). Successivamente si usano i primi K-2 più l’ultimo fold per addestrare e il (K-1)-simo per valutare le prestazioni. La procedura va avanti finchè non si sono ottenuti K addestramenti e K valutazioni. L’insieme delle K valutazioni viene usato per stimare delle misure di prestazioni che saranno più affidabili che se la procedura fosse stata effettuata solo una volta. Nella figura si cerca di illustrare questo metodo.

Quella che abbiamo descritto è chiamata K-fold Cross-Validation.

Ulteriori informazioni possono essere trovate su WIKIPEDIA.

Nfold <- 5

myctrl <- trainControl(method='cv', number=Nfold)

modLDA<-train(x = dati[,-3], y=dati[,3],
  method='lda',
  trControl = myctrl)

modSVM<-train(x = dati[,-3], y=dati[,3],
  method='svmLinear2',
  trControl = myctrl)

modNNET<-train(x = dati[,-3], y=dati[,3],
  method='nnet',
  trControl = myctrl, trace=FALSE)

modKNN<-train(x = dati[,-3], y=dati[,3],
  method='knn',
  trControl = myctrl)



print(modLDA$results)
##   parameter Accuracy Kappa AccuracySD   KappaSD
## 1      none     0.78  0.56  0.1279648 0.2559297
print(modSVM$results)
##   cost Accuracy Kappa AccuracySD    KappaSD
## 1 0.25    0.755  0.51 0.04107919 0.08215838
## 2 0.50    0.750  0.50 0.04330127 0.08660254
## 3 1.00    0.750  0.50 0.04330127 0.08660254
print(modNNET$results)
##   size decay Accuracy Kappa AccuracySD    KappaSD
## 1    1 0e+00    0.775  0.55 0.05590170 0.11180340
## 2    1 1e-04    0.775  0.55 0.05590170 0.11180340
## 3    1 1e-01    0.770  0.54 0.05419871 0.10839742
## 4    3 0e+00    0.715  0.43 0.02850439 0.05700877
## 5    3 1e-04    0.750  0.50 0.05863020 0.11726039
## 6    3 1e-01    0.770  0.54 0.04107919 0.08215838
## 7    5 0e+00    0.710  0.42 0.03354102 0.06708204
## 8    5 1e-04    0.735  0.47 0.04873397 0.09746794
## 9    5 1e-01    0.760  0.52 0.03354102 0.06708204
print(modKNN$results)
##   k Accuracy Kappa AccuracySD   KappaSD
## 1 5    0.740  0.48 0.08767839 0.1753568
## 2 7    0.730  0.46 0.10062306 0.2012461
## 3 9    0.745  0.49 0.06471090 0.1294218

4.2.5 Repeated - CV

Nfold <- 5

myctrl <- trainControl(method='repeatedcv', 
  number=Nfold,
  repeats = 10)

modLDA<-train(x = dati[,-3], y=dati[,3],
  method='lda',
  trControl = myctrl)

modSVM<-train(x = dati[,-3], y=dati[,3],
  method='svmLinear2',
  trControl = myctrl)

modNNET<-train(x = dati[,-3], y=dati[,3],
  method='nnet',
  trControl = myctrl, trace=FALSE)

modKNN<-train(x = dati[,-3], y=dati[,3],
  method='knn',
  trControl = myctrl)



print(modLDA$results)
##   parameter Accuracy Kappa AccuracySD    KappaSD
## 1      none   0.7715 0.543  0.0479184 0.09583681
print(modSVM$results)
##   cost Accuracy Kappa AccuracySD   KappaSD
## 1 0.25   0.7655 0.531 0.07773746 0.1554749
## 2 0.50   0.7645 0.529 0.07285889 0.1457178
## 3 1.00   0.7645 0.529 0.07268362 0.1453672
print(modNNET$results)
##   size decay Accuracy Kappa AccuracySD    KappaSD
## 1    1 0e+00   0.7655 0.531 0.04680888 0.09361776
## 2    1 1e-04   0.7655 0.531 0.04680888 0.09361776
## 3    1 1e-01   0.7730 0.546 0.04513020 0.09026039
## 4    3 0e+00   0.7360 0.472 0.06168551 0.12337102
## 5    3 1e-04   0.7425 0.485 0.04903404 0.09806807
## 6    3 1e-01   0.7715 0.543 0.04226073 0.08452146
## 7    5 0e+00   0.7255 0.451 0.06728897 0.13457795
## 8    5 1e-04   0.7315 0.463 0.05950690 0.11901380
## 9    5 1e-01   0.7675 0.535 0.04552696 0.09105392
print(modKNN$results)
##   k Accuracy Kappa AccuracySD    KappaSD
## 1 5   0.7410 0.482 0.05434208 0.10868415
## 2 7   0.7490 0.498 0.04682795 0.09365591
## 3 9   0.7435 0.487 0.05170807 0.10341614

4.2.6 Vediamo i vari metodi su un unico classificatore LDA

myctrl <- trainControl(method='LGOCV', p=0.75,number=1)
modLDA1<-train(x = dati[,-3], y=dati[,3],
  method='lda',
  trControl = myctrl)

myctrl <- trainControl(method='LGOCV', p=0.75, number=50)
modLDA2<-train(x = dati[,-3], y=dati[,3],
  method='lda',
  trControl = myctrl)

myctrl <- trainControl(method='LOOCV')
modLDA3<-train(x = dati[,-3], y=dati[,3],
  method='lda',
  trControl = myctrl)

myctrl <- trainControl(method='cv',number=5)
modLDA4<-train(x = dati[,-3], y=dati[,3],
  method='lda',
  trControl = myctrl)

myctrl <- trainControl(method='repeatedcv',number=5,repeats=50)
modLDA5<-train(x = dati[,-3], y=dati[,3],
  method='lda',
  trControl = myctrl)


print(modLDA1$results)
##   parameter Accuracy Kappa AccuracySD KappaSD
## 1      none     0.78  0.56         NA      NA
print(modLDA2$results)
##   parameter Accuracy  Kappa AccuracySD  KappaSD
## 1      none   0.7688 0.5376  0.0531705 0.106341
print(modLDA3$results)
##   parameter Accuracy Kappa
## 1      none    0.785  0.57
print(modLDA4$results)
##   parameter Accuracy Kappa AccuracySD  KappaSD
## 1      none    0.775  0.55   0.121192 0.242384
print(modLDA5$results)
##   parameter Accuracy  Kappa AccuracySD   KappaSD
## 1      none   0.7718 0.5436 0.06160308 0.1232062

5 Altre problematiche

5.1 Sbilanciamento del data-set

Un altro problema consiste nel numero elevato di esemplari di alcune classi rispetto ad altre.

Per esempio se considero una patologia molto rara potrei avere 1 caso su 1000. Usando come data-set 1000 casi di cui solo 1 in cui la patologia risulta presente, il classificatore otterrebbe una accuratezza del 99.9 % seplicemente assegnando la classe ‘no-patologia’ a tutti gli esemplari anche a quello con la patologia !! In questo modo le feature non verrebbero proprio prese in considerazione.

Il data-set dovrebbe pertanto essere sempre bilanciato: avere un numero di esemplari grosso modo uguale per le varie classi.

6 Trade-off complessità algoritmi vs generalizzabilità

In genere più un algoritmo è complesso (con molti parametri) più riesce ad adattarsi ai dati di training: questo non è sempre positivo poichè può generare over-fitting. Questo compromesso tra complessità e generalizzabilità si può illustrare bene con un esempio di regressione.

La regressione consiste nel prevedere un valore numerico partendo da feature (la classificazione consiste nel prevedere la classe di un individuo che è un dato categorico, cioè presente in un numero finito di posisbilità)

6.1 Esempio di regressione: dati simulati

x <- seq(0,2*pi,.1)
y <- sin(x) # modello generativo

set.seed(1234)
indici <- sample(1:length(x),size=8)
xn <- x[indici]
noise <- rnorm(length(xn),mean = 0,sd = .1)
yn <- y[indici] + noise # misure 

set.seed(56789)
indicit <- sample(1:length(x),size=5)
xt <- x[indicit]
yt <- y[indicit] + rnorm(length(indicit), mean = 0, sd= .1)

plot(y ~ x,type='l',col='green')
points(yn ~ xn,type='p',pch=1,col='red')
points(yt ~ xt,type='p',pch=1,col='blue')
legend(x='topright',legend=c('ideale','training','test'),
       pch=15,col=c('green','red','blue'))

6.2 Fit con modelli di complessità crescente

6.2.1 Fit con modello lineare: ordine 0

mod0 <- lm(yn ~ 1)
plot(y ~ x,type='l',col='green')
points(yn ~ xn,type='p',pch=1)
abline(mod0,col='red',lty=3)
points(yn ~ xn,type='p',pch=1,col='red')
points(yt ~ xt,type='p',pch=1,col='blue')
legend(x='topright',legend=c('ideale','training','test'),
       pch=15,col=c('green','red','blue'))

6.2.2 Fit con modello lineare: ordine 1

mod1 <- lm(yn ~ 1 + xn)
plot(y ~ x,type='l',col='green')
points(yn ~ xn,type='p',pch=1)
abline(mod1,col='red',lty=3)
points(yn ~ xn,type='p',pch=1,col='red')
points(yt ~ xt,type='p',pch=1,col='blue')
legend(x='topright',legend=c('ideale','training','test'),
       pch=15,col=c('green','red','blue'))

6.2.3 Fit con modello lineare: ordine 2

mod2 <- lm(yn ~ 1+xn+I(xn^2))
plot(y ~ x,type='l',col='green')
points(yn ~ xn,type='p',pch=1)
lines(x,predict(mod2,newdata = data.frame(xn=x)),col='red',lty=3)
points(yn ~ xn,type='p',pch=1,col='red')
points(yt ~ xt,type='p',pch=1,col='blue')
legend(x='topright',legend=c('ideale','training','test'),
       pch=15,col=c('green','red','blue'))

6.2.4 Fit con modello lineare: ordine 3

6.2.5 Fit con modello lineare ordine 7

6.3 Errore di approssimazione

Al variare dell’ordine del modello l’errore di approssimazione ai dati diminuisce e tende a 0.

errori<-rep(NA,length(xn))
erroriTest <- errori
mod <- list()
f <- 'yn ~ -1 '
for(k in c(1:8)){
  f<-paste(f,  '+I(xn^',format(k-1),')',sep='')
  mod[[k]]<-lm(formula(f))
  errori[k]<- sqrt(mean(residuals(mod[[k]])^2))
  erroriTest[k]<- sqrt(mean((yt-predict(mod[[k]],newdata=data.frame(xn=xt)))^2))
  
}

L’errore sul training-set tende a 0 se faccio aumentare la complessità del modello.

plot(c(1:8),errori,type='b',ylim=c(0,2))
lines(c(1:8),erroriTest,type='b',col='red')
legend(x='top',legend=c('errore training','errore test'),pch=15,col=c('black','red'))

7 Explainability

Nel caso di molte feature è desiderabile cercare di capire come mai un classificatore prende una certa decisione, cioè a quali feature da maggiore importanza.

Il numero di feature in ambito clinico è molto aumentato: si pensi alla genetica e alla radiomica, per cui questo problema è divenuto particolarmente importante.

Per capire alcune idee fondamentali, generiamo dati con alcune feature che sono collegate alla classe mentre altre feature sono casuali.

N<-100
set.seed(123456)
dati <- data.frame(
  x1=c(rnorm(N),rnorm(N,mean=2)),
  x2=c(rnorm(N),rnorm(N,mean=2)),
  x3=rnorm(2*N,sd=2), # le feature x3 e x4 non sono collegate alla classe
  x4=rnorm(2*N,sd=2),
  cl=factor(rep(c(1,2),each=N))  
  )

set.seed(123456)
indiciTrain <- sample(1:(2*N),size=0.75*2*N)

training <- dati[indiciTrain,]
test <- dati[-indiciTrain,]


pairs(dati[,-5],col=c('red','blue')[dati$cl],pch=20,asp=1)

Nel grafico si vede che il plot di x3 vs x4 contiene punti rossi e blue sovrapposti, mentre le variabili x1 e x2 che sono quelle collegate alla classe danno dei grafici in cui gli agglomerati rossi e blue sono separati. [Questa è comunque una situazione ideale nella raltà le cose non sono mai così nette !]

Vediamo solo i traning

pairs(training[,-5],col=c('red','blue')[training$cl],pch=20,asp=1)

e solo i test

pairs(test[,-5],col=c('red','blue')[test$cl],pch=20,asp=1)

Ora vorremmo cercare di capire un classificatore quali feature considera per prendere la sua decisione.

Il package DALEX contiene la funzione model_parts che consente di esaminare le prestazioni del classificatore quando una certa feature viene permutata cioè i valori della feature sui vari individui sono scambiati tra loro: l’idea è che se una feature non è collegata alla classe dell’individuo, permutando i valori tra i vari individui questa feature non dovrebbe comunque avere influenza sulla classificazione; se una feature invece ha importanza, permutando i valori si avrà una perdita di prestazione del classificatore.

library(caret)
mod<-train( x=training[,-5],
  y=training[,5],
  method='svmLinear2',
  trControl = trainControl(method='repeatedcv',number=5,repeats=50),
  probability=TRUE)

library(DALEX)

explainer<-explain(mod,data=test[,-5],y=as.numeric(test[,5])-1,label='SVM')
## Preparation of a new explainer is initiated
##   -> model label       :  SVM 
##   -> data              :  50  rows  4  cols 
##   -> target variable   :  50  values 
##   -> predict function  :  yhat.train  will be used (  default  )
##   -> predicted values  :  No value for predict function target column. (  default  )
##   -> model_info        :  package caret , ver. 6.0.93 , task classification (  default  ) 
##   -> predicted values  :  numerical, min =  0.003830008 , mean =  0.5554298 , max =  0.9979057  
##   -> residual function :  difference between y and yhat (  default  )
##   -> residuals         :  numerical, min =  -0.6997434 , mean =  -0.01542976 , max =  0.6956231  
##   A new explainer has been created!
expl_svm<-model_parts(explainer,type='variable_importance')
perf_svm <- model_performance(explainer,label='SVM')
plot(expl_svm)

Vediamo se la cosa dipende dal classificatore:

library(caret)
mod<-train( x=training[,-5],
  y=training[,5],
  method='lda',
  trControl = trainControl(method='repeatedcv',number=5,repeats=50))

library(DALEX)

explainer<-explain(mod,data=test[,-5],y=as.numeric(test[,5])-1,label='LDA')
## Preparation of a new explainer is initiated
##   -> model label       :  LDA 
##   -> data              :  50  rows  4  cols 
##   -> target variable   :  50  values 
##   -> predict function  :  yhat.train  will be used (  default  )
##   -> predicted values  :  No value for predict function target column. (  default  )
##   -> model_info        :  package caret , ver. 6.0.93 , task classification (  default  ) 
##   -> predicted values  :  numerical, min =  0.000905956 , mean =  0.5492297 , max =  0.9987487  
##   -> residual function :  difference between y and yhat (  default  )
##   -> residuals         :  numerical, min =  -0.7005501 , mean =  -0.009229666 , max =  0.5778278  
##   A new explainer has been created!
expl_lda <- model_parts(explainer,type='variable_importance')
perf_lda<-model_performance(explainer,label='LDA')

plot(expl_lda,expl_svm)

Vediamo se le performance dei classificaori differiscono:

plot(perf_svm,perf_lda,geom='roc')

8 Feature selection

Quando ci sono molte feature è opportuno selezionare le feature prima della classificazione, infatti: - le feature inappropriate possono confondere il classificatore - lo spazio di ricerca delle soluzioni ha un numero eccessivo di dimensioni (curse fo dimensionality) che può portare a soluzioni non adeguate

Vediamo un esempio concreto. La radiomica consente di estrarre una notevole quantità di feature da una immagine radiologica. Uno degli obiettivi della radiogenomica è di individuare asociazioni tra radiologia e genetica.

Vediamo un certo gene se è associato alla radiomica.

d<-read.csv('dati-esempio/polmone-tot.csv')

m<-as.matrix(d[,-c(1,109:114)])
rownames(m)<-d$ID
colnames(m)<-names(d)[-c(1,109:114)]
m<-scale(m)# ogni feat è normalizza sui pz
kgene<-114
heatmap(m,scale='none',
  Rowv=TRUE,Colv=TRUE,mar=c(15,3),
  RowSideColors = c('red','blue')[d[,kgene]],
  col=cm.colors(9, alpha = 1, rev = FALSE),
  main=names(d)[kgene])

# scala le variabili in d
d[,-c(1,109:114)] <- scale(d[,-c(1,109:114)])

8.1 Selezione del miglior sottoinsieme

Vediamo quale è il miglior sottoinsieme di feture radiomiche cioè quello che ha la più alta accuratezza rispetto al gene individuato. Dovremmo esaminare tutte le singole, coppie, triple, etcc.

library(MASS)
nomiF<-names(d)[grep('original',names(d))]
Nf <- length(grep('original',names(d)))
indiciF<-grep('original',names(d))
# singole
perS<-rep(NA,Nf)# performance
for (nf in 1:Nf){
  x<-d[,c(indiciF[nf],kgene)]
  names(x)[2]<-'classe'
  mod<-lda(classe ~., x)
  perS[nf]<- sum(predict(mod)$class == x$classe) # accuratezza=corrette predizioni
}
# migliore feature singola
nomiF[which(perS==max(perS))]
## [1] "original_glrlm_RunVariance"                        
## [2] "original_glszm_SizeZoneNonUniformity"              
## [3] "original_gldm_DependenceVariance"                  
## [4] "original_gldm_LargeDependenceHighGrayLevelEmphasis"
# coppie
perD<-rep(NA,choose(Nf,2))
coppie<-matrix(NA,choose(Nf,2),2)
conta<-1
for(k1 in 1:(Nf-1))
  for(k2 in (k1+1):Nf){
    if (k1!=k2){
    coppie[conta,1]<-k1
    coppie[conta,2]<-k2
    conta<-conta+1
      
    }
}

for(conta in 1:choose(Nf,2)){
  x<-d[,c(kgene,indiciF[coppie[conta,1]],indiciF[coppie[conta,2]])]
  names(x)[1]<-'classe'
  mod<-lda(classe ~., x)
  perD[conta]<- sum(predict(mod)$class == x$classe) # accuratezza=corrette  
  conta<-conta+1
  }


coppie[which(perD==max(perD)),]
## [1]  7 12
nomiF[coppie[which(perD==max(perD)),1]]
## [1] "original_shape_Maximum2DDiameterSlice"
nomiF[coppie[which(perD==max(perD)),2]]
## [1] "original_shape_SurfaceArea"

Il problema di questo approccio è che se abbiamo \(p\) feature abbiamo \(2^p\) possibili sottoinsiemi con 1, 2 3, \(\ldots\) , \(p\) feature, se \(p=10\) il numero delle possibilità è \(2^10 \approx 10^6\).

Quindi non è fattibile nella pratica: ci vorrebbero computer molto potenti impegnati per molte ore/giorni con solo 10 feature (a seconda dei classificatori usati) : nel caso visto sopra abbiamo 107 feature !!

Un approccio possibile usa l’idea del bootstrap, cioè seleziona a caso un certo numero di feature (non tutte le possibili coppie, triple etc…).

Nrep <- 50
errori <- matrix(NA,nrow=Nf,Nrep)

for(k in 1:Nf){
  for(ks in 1:Nrep){
    ind_samp <- sample(1:Nf,k)
    x <- data.frame(classe=d[,kgene], d[,indiciF[ind_samp]])
    mod<-lda(classe ~ ., data=x)
    errori[k,ks]<- sum(d[,kgene] != predict(mod)$class)
  }  
}

#levelplot(errori,xlab='numero feature',ylab='ripetizioni')
plot(apply(errori,1,mean), type='l',xlab='numero feature',ylab='numero errori medio')

Si vede che per circa 58-60 feature il numero di errori è minimo, quindi si può cercare tra tutte le combinazioni con 60 feature. Questo dimostra che vi sono una serie di variabili non informative.

Vediamo due approcci intermedi che non usano tutti i sottoinsiemi ma solo una parte: backward e forward selection.

library(klaR)

xx<-d[,-c(1,110:114)]
names(xx)[108]<-'classe'
modb<-stepclass(classe ~ . , data=xx, 
  method = 'lda', direction='backward', fold=5, output=FALSE)

modf<-stepclass(classe ~ . , data=xx, 
  method = 'lda', direction='forward', fold=5, output=FALSE)

plot(modb)

plot(modf)

9 Principal Components Analysis (PCA)

La riduzioen delle features in genere è possibile perchè i dati dipendono da un numero limtato di variabili latenti, per esempio se le variabili sono 3 ma due sono correlate in realtà le variabili indipendenti (latenti) sono solo 2. Questo principio è alla base della PCA.

Vediamo sempre su dati simulati con poche variabili.

x<-c(rnorm(100),rnorm(100,mean=1))
y<-c(rnorm(100,mean=1),rnorm(100,mean=0))
w<-c(rnorm(100,mean=.5),rnorm(100,mean=-0.5))
z<-x*3+rnorm(100,mean=0,sd=2)
cl<-factor(rep(c('a','b'),each=100))


d<-cbind(x,y,z,w)

pairs(d,col=c('red','blue')[as.numeric(cl)])

pca<-prcomp(scale(d))
plot(pca)

pairs(pca$x,col=c('red','blue')[as.numeric(cl)])

set.seed(1234)
ind_train<-sample(1:dim(d)[1],0.8*dim(d)[1])

library(MASS)

x1<-data.frame(x,y,z,w,cl)
x1<-x1[ind_train,]
x1t<-x1[-ind_train,]
mod1<-lda(cl ~ . , data=x1)

x2<-data.frame(pca$x[,1:3],cl)
x2<-x2[ind_train,]
x2t<-x2[-ind_train,]
mod2<-lda(cl ~ . , data=x2)

x3<-data.frame(pca$x[,1],pca$x[,4],cl)
x3<-x3[ind_train,]
x3t<-x3[-ind_train,]
mod3<-lda(cl ~ . , data=x3)

# Errori
# questi due modelli sono equivalenti
print(sum(x1t$cl != predict(mod1,newdata=x1t)$class))
## [1] 9
print(sum(x2t$cl != predict(mod2,newdata=x2t)$class))
## [1] 9
# questo commette più errori
print(sum(x3t$cl != predict(mod3,newdata=x3t)$class))
## [1] 12

Gli \(n\) valori più alti nel plot(pca) ci dicono quante sono le variabili latenti (indipendenti). PC1 e PC2 nel plot a coppie sono le variabili latenti scoperte in automatico. (bisogna conservare sempre le prime \(n\) PC )

Gli errori fatti dai modelli che contengono solo le PC1 e PC2 sono grosso modo gli stessi fatti dal modello originale.