Algoritmi diagnostici e intelligenza artificiale: hans on R

Parte I: Esplorazione dei dati

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 sicurezza la patologia è necessario effettuare indagini strumentali costose e/o invasive per il paziente. Per cui è di interesse verificare se è possibile rilevare la presenza di patologia mediante esami alternativi poco costosi.

Per valutare l’efficacia dei metodi alternativi è necessario preliminarmente effettuare uno studio per esaminare le medesime osservazioni su una popolazione di pazienti cui è stata diagnosticata la medesima patologia.

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 del machine learning.

In termini molto semplici uno dei task principali del machine learning è di assegnare un certo individuo ad una certa classe. 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 con le corrispondenti label (etichette, classe di appartenenza) corrette.

Con il linguaggio matematico se x = [f1,f2,…fN] sono N features che rappresentano l’individuo e se \(\omega_k \in {\omega_1, \ldots, \omega_K}\) è la sua classe tra le \(K\) possibili, la macchina impara ad assegnare x 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 giocato ri 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).

Dati di Esempio

In genere i dati raccolti consistono in una serie di valori misurati (decine o centinaia di variabili) o di risposte a questionari (anche qui possono essere decine di domande) che possono riguardare una popolazione molto numerosa di pazienti (centinaia o migliaia).

E’ molto difficile trarre conclusioni dall’esame visivo di una tabella composta da centinaia di colonne (le variabili) e migliaia di righe (i pazienti).

Vedremo nel seguito le tecniche statistiche disponibili in R per esaminare e riassumere i dati raccolti per trarre conclusioni.

La prima cosa da fare è acquisire dimestichezza con i dati che si vogliono analizzare.

Carichiamo un set di dati clinici di esempio riguardante pazienti che hanno effettuato una mammografia. I dati si possono scaricare al seguente indirizzo: dati-esempio.

library(readxl)
df <- read_excel(path='/Volumes/CORRENTE/COLLABORAZIONI/FORMISANO/MASTER/dati-esempio/dati-clinici-mammografia.xlsx')
df <- as.data.frame(df)
str(df)

## 'data.frame':    565 obs. of  20 variables:
##  $ codice2       : chr  NA NA NA NA ...
##  $ peso          : chr  NA NA "60" "73" ...
##  $ altezza       : chr  NA NA "1.78" "155" ...
##  $ BMI           : chr  NA NA "189370.03" "30.39" ...
##  $ eta           : chr  "60" "67" "47" "71" ...
##  $ CA            : chr  NA NA NA NA ...
##  $ eta_menarca   : chr  NA NA "15" "13" ...
##  $ estrogeni     : chr  "SI" NA "SI" "SI" ...
##  $ contraccettivi: chr  "SI" NA "NO" "NO" ...
##  $ nullipara     : chr  NA NA "NO" "NO" ...
##  $ gravidanza_30 : chr  NA NA "SI" "NO" ...
##  $ allattamento  : chr  NA NA "SI" "NO" ...
##  $ menopausa     : chr  NA NA "NO" "SI" ...
##  $ eta_menopausa : chr  NA NA NA "45" ...
##  $ alcol         : chr  NA NA "NO" "NO" ...
##  $ carne         : chr  NA NA "NO" "NO" ...
##  $ salumi        : chr  NA NA "SI" "NO" ...
##  $ verdura       : chr  NA NA "SI" "NO" ...
##  $ frutta        : chr  NA NA "SI" "NO" ...
##  $ densita       : chr  "A" "B" "C" "B" ...

Vedere i dati

Per avere un idea dei dati vediamo le prime 5 righe del dataset:

head(df,5)

##   codice2 peso altezza       BMI eta   CA eta_menarca estrogeni
## 1    <NA> <NA>    <NA>      <NA>  60 <NA>        <NA>        SI
## 2    <NA> <NA>    <NA>      <NA>  67 <NA>        <NA>      <NA>
## 3    <NA>   60    1.78 189370.03  47 <NA>          15        SI
## 4    <NA>   73     155     30.39  71 <NA>          13        SI
## 5    <NA>   65     164     24.17  40 <NA>          13        NO
##   contraccettivi nullipara gravidanza_30 allattamento menopausa
## 1             SI      <NA>          <NA>         <NA>      <NA>
## 2           <NA>      <NA>          <NA>         <NA>      <NA>
## 3             NO        NO            SI           SI        NO
## 4             NO        NO            NO           NO        SI
## 5             SI        NO            SI           NO        NO
##   eta_menopausa alcol carne salumi verdura frutta densita
## 1          <NA>  <NA>  <NA>   <NA>    <NA>   <NA>       A
## 2          <NA>  <NA>  <NA>   <NA>    <NA>   <NA>       B
## 3          <NA>    NO    NO     SI      SI     SI       C
## 4            45    NO    NO     NO      NO     NO       B
## 5          <NA>    NO    NO     NO      NO     NO       C

Tipi di dato

• Numerici (numeri reali, qualsiasi valore)
• Categorici (un numero predefinito di valori, fattori)

Convertiamo in fattori le variabili chr; ed in numeri le variabili numeriche.

df <- transform(df,
                codice2 = factor(codice2),
                peso = as.numeric(peso),
                altezza = as.numeric(altezza),
                BMI = as.numeric(BMI),
                eta = as.numeric(eta),
                CA = as.numeric(CA),
                eta_menarca = as.numeric(eta_menarca),
                estrogeni = factor(estrogeni),
                contraccettivi = factor(contraccettivi),
                nullipara = factor(nullipara),
                gravidanza_30 = factor(gravidanza_30),
                allattamento = factor(allattamento),
                menopausa = factor(menopausa),
                eta_menopausa = as.numeric(eta_menopausa),
                alcol = factor(alcol),
                carne = factor(carne),
                salumi = factor(salumi),
                frutta = factor(frutta),
                verdura = factor(verdura),
                densita = factor(densita)
                )

Vediamo i dati in forma sintetica

Vediamo un riassunto delle colonne:

summary(df)

##             codice2         peso           altezza      
##  000000000      :  1   Min.   :  9.00   Min.   :  1.57  
##  000111         :  1   1st Qu.: 60.00   1st Qu.:158.00  
##  01-64364_012722:  1   Median : 66.00   Median :160.00  
##  0101728_021121 :  1   Mean   : 68.77   Mean   :157.19  
##  0114468_111220 :  1   3rd Qu.: 75.00   3rd Qu.:166.00  
##  (Other)        :492   Max.   :166.00   Max.   :187.00  
##  NA's           : 68   NA's   :8        NA's   :10      
##       BMI                eta              CA         eta_menarca    
##  Min.   :    17.9   Min.   : 2.00   Min.   : 44.0   Min.   :  1.00  
##  1st Qu.:    22.7   1st Qu.:48.00   1st Qu.: 98.0   1st Qu.: 11.00  
##  Median :    25.4   Median :55.00   Median :102.0   Median : 12.00  
##  Mean   :  6619.2   Mean   :53.43   Mean   :101.6   Mean   : 12.43  
##  3rd Qu.:    29.4   3rd Qu.:64.00   3rd Qu.:106.0   3rd Qu.: 13.00  
##  Max.   :327160.5   Max.   :77.00   Max.   :150.0   Max.   :102.00  
##  NA's   :12         NA's   :504     NA's   :66      NA's   :13      
##  estrogeni  contraccettivi nullipara  gravidanza_30 allattamento
##  NO  :497   NO  :453       NO  :444   NO  :275      NO  :326    
##  SI  : 55   SI  :103       SI  :111   SI  :166      SI  :110    
##  NA's: 13   NA's:  9       NA's: 10   NA's:124      NA's:129    
##                                                                 
##                                                                 
##                                                                 
##                                                                 
##  menopausa  eta_menopausa     alcol      carne      salumi    verdura   
##  NO  :170   Min.   :  1.00   NO  :538   NO  :460   NO  :458   NO  :450  
##  SI  :384   1st Qu.: 45.00   SI  : 16   SI  : 94   SI  : 95   SI  :105  
##  NA's: 11   Median : 50.00   NA's: 11   NA's: 11   NA's: 12   NA's: 10  
##             Mean   : 47.42                                              
##             3rd Qu.: 53.00                                              
##             Max.   :535.00                                              
##             NA's   :191                                                 
##   frutta    densita   
##  NO  :436   A   : 59  
##  SI  :115   AA  :  3  
##  NA's: 14   B   :234  
##             C   :150  
##             D   : 94  
##             NA's: 25  
## 

Individuazione di errori

Una delle fasi più importanti nella esplorazione dei dati è acquisire dimestichezza con i valori raccolti: in quale range si trovano, quali sono le categorie possibili. E’ anche importante capire se i dati sono sporchi.

Vediamo alcuni grafici delle varie variabili numeriche, per individuare: valori troppo lontani dal resto della popolazione (outlier)

layout(matrix(1:8,4,2))
with(df, {
  plot(peso)
  plot(altezza)
  plot(BMI)
  plot(CA)
  plot(eta)
  plot(eta_menarca)
  plot(eta_menopausa)
})

Individazione di errori (2)

Nei dati reali spesso:

• i valori non sono inseriti [missing data] (e.g. perchè non disponibili): R gestisce il dato mancante con il valore NA
• sono inseriti in modo errato (e.g. unità di misura non corretta)

vediamo per esempio la variabile altezza:

df$altezza

##   [1]     NA     NA   1.78 155.00 164.00 162.00 150.00 168.00 163.00 160.00
##  [11] 165.00 158.00 158.00     NA   1.60   1.57 160.00   1.70   1.70   1.70
##  [21]   1.62   1.66   1.63   1.65 170.00 150.00     NA 158.00 170.00 165.00
##  [31] 165.00 148.00 154.00 170.00 158.00 155.00 165.00 153.00 163.00 156.00
##  [41] 173.00 168.00 170.00 170.00 170.00 160.00 158.00 165.00 153.00 153.00
##  [51] 150.00 165.00 156.00 162.00 158.00 152.00 156.00     NA 157.00 160.00
##  [61] 160.00 155.00 164.00 151.00 158.00 160.00 160.00 154.00 166.00 162.00
##  [71] 170.00 154.00 158.00 160.00 168.00 160.00 160.00 157.00 160.00 154.00
##  [81] 154.00 165.00 158.00 160.00 160.00 171.00 153.00 170.00 162.00 150.00
##  [91] 165.00 164.00 155.00 160.00 160.00 156.00 155.00 163.00 165.00 158.00
## [101] 160.00 167.00 160.00 167.00 156.00 172.00 163.00 167.00 150.00 165.00
## [111] 158.00 160.00 168.00 167.00 165.00 168.00 160.00 168.00 170.00 169.00
## [121] 150.00 170.00 157.00 150.00 160.00 160.00 165.00 150.00  58.00 159.00
## [131] 159.00 165.00 160.00 150.00 173.00 156.00 160.00 145.00 150.00 161.00
## [141] 160.00 163.00 165.00   1.68 160.00 167.00 167.00 156.00 167.00 157.00
## [151] 158.00 160.00 163.00 165.00 162.00 150.00 158.00 164.00 168.00 157.00
## [161] 170.00 165.00 160.00 170.00 165.00 160.00 158.00 160.00 163.00 160.00
## [171] 158.00 155.00 154.00 168.00 164.00 154.00 160.00 166.00 154.00 160.00
## [181] 174.00 165.00 165.00 170.00 168.00 155.00 155.00 161.00 172.00 165.00
## [191] 152.00 170.00 170.00 160.00 178.00 155.00 163.00 180.00 160.00 170.00
## [201] 163.00 150.00 173.00 159.00 162.00 168.00 160.00 165.00 158.00 164.00
## [211] 161.00 168.00 165.00 160.00 155.00 164.00 165.00 167.00 165.00 172.00
## [221] 166.00  50.00 165.00 159.00 160.00 156.00 163.00 158.00 158.00 158.00
## [231] 164.00 160.00 158.00 160.00 170.00 165.00 180.00 165.00 153.00 170.00
## [241] 165.00 173.00 155.00 165.00 154.00 168.00 170.00 164.00 168.00 159.00
## [251] 160.00 155.00 162.00 156.00 159.00 162.00 165.00 154.00 165.00 175.00
## [261] 165.00 154.00 168.00 170.00 174.00 163.00 169.00 163.00 162.00 164.00
## [271] 155.00 150.00 170.00 165.00 168.00 168.00 170.00 170.00 155.00 165.00
## [281] 156.00 160.00 168.00 155.00 160.00 156.00 157.00 160.00 168.00 153.00
## [291] 168.00 160.00 160.00 158.00 165.00 162.00 159.00 160.00 164.00 162.00
## [301] 165.00 157.00 160.00  62.00 160.00 155.00 160.00 162.00 170.00 160.00
## [311] 158.00 172.00 164.00 158.00 168.00 152.00 164.00 150.00 166.00 157.00
## [321] 159.00 167.00 165.00 160.00 157.00 157.00 157.00 168.00 166.00 165.00
## [331] 165.00     NA 160.00 165.00 168.00 169.00 169.00 159.00 167.00 172.00
## [341] 160.00 155.00 165.00 167.00 156.00 167.00 160.00 169.00 160.00 165.00
## [351] 164.00 161.00 150.00 160.00 163.00 160.00 162.00 160.00 160.00 160.00
## [361] 155.00 153.00 172.00 169.00 170.00 165.00 168.00 160.00 165.00 152.00
## [371] 172.00 158.00 160.00 157.00 169.00 158.00 167.00 166.00 160.00 158.00
## [381] 172.00 160.00 150.00 155.00 169.00 171.00 165.00 158.00 165.00 158.00
## [391] 150.00 160.00 165.00 173.00 167.00 150.00 165.00 175.00 160.00 164.00
## [401] 172.00 155.00 162.00 174.00 171.00     NA 168.00 163.00 160.00 160.00
## [411] 150.00 160.00 160.00 158.00 165.00 162.00 150.00 150.00 150.00 163.00
## [421] 158.00 159.00 177.00 171.00 173.00 158.00 168.00 167.00 170.00 147.00
## [431] 164.00 158.00 168.00 150.00 170.00 153.00 160.00 165.00 157.00 160.00
## [441] 168.00 160.00 158.00 157.00 160.00 157.00 166.00 158.00 167.00 165.00
## [451] 162.00 168.00 166.00 159.00 165.00 150.00 163.00   1.80 173.00 160.00
## [461] 160.00 157.00 160.00   1.65     NA   1.66 159.00 168.00 165.00 170.00
## [471] 170.00 166.00 172.00 166.00   1.78 158.00 155.00 170.00 157.00 160.00
## [481] 160.00 155.00 170.00 160.00 150.00 162.00 171.00 160.00 167.00 170.00
## [491] 168.00 163.00 158.00     NA     NA 168.00 172.00 160.00 158.00 170.00
## [501] 155.00 181.00 175.00 165.00 159.00 165.00 174.00 154.00 160.00 164.00
## [511] 167.00 163.00 155.00 157.00 160.00 155.00 160.00 173.00 164.00 153.00
## [521] 167.00 160.00 166.00 160.00 160.00 155.00 170.00 165.00 160.00 160.00
## [531] 152.00 163.00 171.00 187.00 170.00 169.00 155.00 160.00 172.00 167.00
## [541] 163.00 150.00 162.00 168.00 168.00 164.00 148.00 165.00 170.00 165.00
## [551] 163.00 165.00 160.00 160.00 160.00 175.00 161.00 153.00 162.00 162.00
## [561] 165.00 162.00 161.00 163.00 160.00

Correzione dati

In alcuni casi è possibile correggere i valori inseriti.

Nel caso dell’altezza si vede che l’altezza è stata inserita per la maggiori parte usando l’unità di misura cm mentre per alcune pazienti l’altezza sembra inserita in m

which(df$altezza < 100) # indici di riga

##  [1]   3  15  16  18  19  20  21  22  23  24 129 144 222 304 458 464 466
## [18] 475

df[which(df$altezza < 100), 'altezza'] # valori 

##  [1]  1.78  1.60  1.57  1.70  1.70  1.70  1.62  1.66  1.63  1.65 58.00
## [12]  1.68 50.00 62.00  1.80  1.65  1.66  1.78

alcune altezze (58,50,62) sono ancora strane. Vediamo le altre variabili corrispondenti alle medesime pazienti

df[which(df$altezza < 100), ] # valori 

##             codice2 peso altezza       BMI eta  CA eta_menarca estrogeni
## 3              <NA>   60    1.78 189370.03  47  NA          15        SI
## 15  01_02357_061820   48    1.60 187500.00  NA  NA          13        NO
## 16  01_18401_061820   62    1.57 251531.50   2  77          14        NO
## 18  01_65092_062220   58    1.70 200692.04  NA  64          13        SI
## 19  01_13605_062220   90    1.70 311418.69  30  74          13        NO
## 20  01_22471_062420   88    1.70 304498.27  71  74          12        NO
## 21  01_64368_062420   83    1.62 316262.76  71  74          12        NO
## 22  01_36480_062420   68    1.66 246770.21  71  68          12        NO
## 23  01_09521_062420   64    1.63 240882.23  72  70          11        NO
## 24  01_43451_062520   68    1.65 249770.43  16  68          14        NO
## 129 01_00001_072320   65   58.00    193.22  68 101          13        NO
## 144 01_27386_090229   62    1.68 219671.20  NA  98          15        NO
## 222 01_56772_093020   60   50.00    240.00  NA  98          13        NO
## 304 01_41804_020821  166   62.00    431.84  NA 102          12        NO
## 458 01_29523_110421  106    1.80 327160.49  NA 114          13        NO
## 464 01_52698_111221   66    1.65        NA  NA  98          14        NO
## 466 01_34228_111221   83    1.66 301204.82  NA  NA          13        NO
## 475 01_14220_112921   93    1.78 293523.55  NA 110          12        NO
##     contraccettivi nullipara gravidanza_30 allattamento menopausa
## 3               NO        NO            SI           SI        NO
## 15              NO        NO            SI           NO        SI
## 16              NO        NO            SI           NO        NO
## 18              NO        NO            NO           SI        SI
## 19              NO        SI          <NA>         <NA>        SI
## 20              NO        NO            SI           NO        SI
## 21              NO        NO            NO           NO        SI
## 22              NO        NO            SI           SI        SI
## 23              NO        NO            NO           NO        SI
## 24              NO        SI          <NA>         <NA>        SI
## 129             NO        SI          <NA>         <NA>        SI
## 144             NO        NO            NO           NO        NO
## 222             NO        NO            NO           SI        SI
## 304             SI        NO            SI           NO        NO
## 458             SI        NO            NO           NO        SI
## 464             NO        NO            NO           NO        SI
## 466             NO        NO            NO           SI        SI
## 475             NO        NO            NO           SI        NO
##     eta_menopausa alcol carne salumi verdura frutta densita
## 3              NA    NO    NO     SI      SI     SI       C
## 15             48    NO    NO     NO      SI     SI       B
## 16             NA    NO    NO     NO      NO     NO       C
## 18             47    NO    NO     SI      NO     NO       C
## 19             52    NO    NO     NO      NO     NO       A
## 20             48    NO    NO     NO      NO     NO       B
## 21             53    NO    NO     NO      NO     NO       B
## 22             52    NO    NO     NO      NO     NO       B
## 23             54    NO    NO     NO      NO     NO       B
## 24             52    NO    NO     NO      NO     NO       C
## 129            50    NO    NO     NO      NO     NO       B
## 144            NA    NO    SI     SI      NO     NO       D
## 222            58    NO    NO     NO      NO     NO       B
## 304            NA    NO    NO     NO      NO     NO       A
## 458            42    NO    NO     NO      NO     NO       B
## 464            53    NO    NO     NO      NO     NO       B
## 466            51    NO    NO     NO      SI     NO       C
## 475            NA    NO    NO     NO      NO     NO       C

Inputazione dati

Come si vede dalle slide precedenti è molto difficile

• assegnare un valore laddove sia mancante
• sostituire un valore ritenuto errato con uno ritenuto corretto:
  • richiede delle ipotesi sulla natura dell’errore commesso
• In definitiva la imputazione dei dati è una delle parti più importanti della raccolta dati e determina l’affidabilità di tutte le analisi successive per cui è necessario impegnarsi al massimo per la imputazione corretta e fare in modo che tutto il personale coinvolto sia consapevole delle problematiche legate agli errori di imputazione;
• Se possibile, è preferibile effettuare una serie di controlli software all’atto della imputazione per impedire che vengano inseriti dati palesemente scorretti:
  • controllo di tipo: numerico invece di carattere
  • di range di valori
  • se possibile, guidare l’input:
    • dati categorici: le possibili categorie sono selezionate dall’utente da un opportuno menu a tendina
    • date: selezione da un calendario
    • valori numerici: se si tratta di interi di cui si conosce il range di variazione si possono inserire usando uno slider o simili costrutti grafici.

Elimino dati errati

Facciamo alcune assunzioni che ci consentono di selezionare i dati che riteniamo buoni e di correggere laddove sia chiaro l’errore.

• il codice identificativo del paziente non può essere vuoto
• laddove le altezze sono minori di 2 assumiamo che ci sia stato un errore legato all’unità di misura e moltiplichiamo le altezze per 100 per ottenere il valore corretto
• nei casi in cui il valore sia intorno a 50 sommiamo 100 ( assumiamo un errore in cui l’utente ha inserito male la prima cifra)
• l’età del menarca deve essere al di sotto di 20 anni e al di sopra dei 2 anni.
• l’età della menopausa deve essere inferiore a 100 anni

df <- df[-which(is.na(df$codice2)),]
df <- df[-which(df$eta_menarca>20 | df$eta_menarca<5),]
df <- df[-which(df$eta_menopausa>100),]

df$altezza <-  ifelse(df$altezza < 2,df$altezza*100, df$altezza) 
df$altezza <-  ifelse(df$altezza > 40 & df$altezza < 70, df$altezza + 100, df$altezza) 


layout(matrix(1:8,4,2))
with(df,plot(peso))
with(df,plot(altezza))
with(df,plot(BMI))
with(df,plot(CA))
with(df,plot(eta))
with(df,plot(eta_menarca))
with(df,plot(eta_menopausa))

Istogramma

Come abbiamo visto nelle slide precedenti guardare la lista dei valori di una variabile non è semplice per cui si preferisce usare dei grafici. Finora abbiamo visto dei plot in cui l’asse x era la riga del paziente.

Per vedere la distribuzione statistica dei dati è necessario costruire un istogramma. Il numero di bin influenza il risultato finale.

Sovrapposta all’istogramma si può vedere - una distribuzione classica per esempio la gaussiana con pari media e varianza - una approssimazione non classica

Teniamo conto che le altezze presentano dati mancanti (NA) per cui nel calcolo della media si devono rimuovere i valori NA.

layout(matrix(1:4,2,2,byrow=T))
hist(df$altezza)
hist(df$altezza,30)

media <- mean(df$altezza,na.rm=TRUE)
deviazione <- sd(df$altezza,na.rm=TRUE)

hist(df$altezza,20,freq=FALSE)
x = pretty(df$altezza,n=100)
lines(x,dnorm(x,mean=media,sd=deviazione),col='red')
lines(density(na.omit(df$altezza)),col='blue')
legend(x='topright',legend=c('normale','approssimazione'),col=c('red','blue'),pch=20)

print(media)

## [1] 162.2165

print(deviazione)

## [1] 6.321985

Dal grafico sembra che l’istogramma sia normale i.e. gaussiano. Quanto possiamo avere fiducia che si tratti di una variabile gaussiana ?

Un test appropriato in questo caso è quello di Shapiro-Wilk.

shapiro.test(df$altezza)

## 
##  Shapiro-Wilk normality test
## 
## data:  df$altezza
## W = 0.9877, p-value = 0.0004133

Dato che il p-value < 0.05 convenzionalmente assumiamo che esiste sufficiente evidenza empirica che la distribuzione NON sia Normale o Gaussiana.

Per avere una confidenza ulteriore possiamo visualizzare un qq-plot che ci da un ulteriore prova di NON normalità della distribuzione. I punti nel grafico devono essere allineati lungo la linea retta se si tratta di normale.

qqnorm(df$altezza)
qqline(df$altezza)

Possiamo valutare la normalità anche del peso.

media <- mean(df$peso,na.rm=TRUE)
deviazione <- sd(df$peso,na.rm=TRUE)

layout(matrix(1:2,1,2))
hist(df$peso,20,freq=FALSE)
x = pretty(df$peso,n=100)
lines(x,dnorm(x,mean=media,sd=deviazione),col='red')
lines(density(na.omit(df$peso)),col='blue')
legend(x='topright',legend=c('normale','approssimazione'),col=c('red','blue'),pch=20)


qqnorm(df$peso)
qqline(df$peso)

Anche in questo caso si vede che c’è uno scostamento; il test di Shapiro-Wilk ci da:

shapiro.test(df$peso)

## 
##  Shapiro-Wilk normality test
## 
## data:  df$peso
## W = 0.90071, p-value < 2.2e-16

per cui non si tratta di distribuzione normale.

Coppie di variabili

Poichè ci aspettiamo che esista una associazione tra peso [kg] e altezza [cm],esploriamo queste due variabili simultaneamente.

E’ possibile aggiungere una retta di regressione al grafico. Il coefficiente di tale retta ci da una indicazione della pendenza della retta.

plot(peso ~ altezza,df)
abline(lm(peso ~ altezza,df),col='red')

coef(lm(peso ~ altezza,df))

## (Intercept)     altezza 
## -21.2417127   0.5529422

Una quantificazione della associazione tra le due variabili è fornita anche dal coefficiente di correlazione di Pearson o di Spearman. Quello di spearman è più appropriato in questo caso perchè le due variabili hanno unità di misura diverse e valori in range differenti.

cor(df$peso, df$altezza, use = 'complete', method='pearson')

## [1] 0.2541492

cor(df$peso, df$altezza, use = 'complete', method='spearman')

## [1] 0.2566899

Per quanto riguarda il BMI (Body Mass Index) ci aspettiamo che sia correlato a peso e altezza. In questo caso dobbiamo fare delle coppie di grafici.

pairs(df[,c('peso','altezza','BMI')])

Ci rendiamo conto però che i valori di BMI sono tutti schiacciati sullo 0 indipendentemente da altezza e peso, inoltre ci sono pochi valori che superano 15000 mentre da una rapida esplorazione dei siti web ci rendiamo conto che i valori di BMI non superano 40: ci deve essere un errore di calcolo.

Visto che la formula del BMI è nota sul web lo ricalcoliamo noi. E ripetiamo il grafico. Ora i valori di BMI sono compatibili con quelli della letteratura.

df <- within(df, BMIc <- peso / (altezza/100)^2 )
pairs(df[,c('peso','altezza','BMIc')])

Possiamo aggiungere la retta di regressione sul grafico peso-BMIc e calcolare la correlazione. Si vede che il valore di correlazione è molto alto.

plot(BMIc ~ scale(peso), df)
abline(modello<-lm(BMIc~scale(peso),df),col='red')

coef(modello)

## (Intercept) scale(peso) 
##   25.989484    4.746759

with(df, cor(peso,BMIc,use='complete',method='spearman') )

## [1] 0.8820422

Associazione tra variabili

Dato che ci sono pazienti in sovrappeso, è utile andare a capire se questa condizione è associata ad un regime alimentare particolare. Per esempio vogliamo vedere se le pazienti che mangiano carne sono in sovrappeso.

Facciamo il seguente grafico chiamato boxplot.

boxplot(peso ~ carne,df,notch=TRUE)
stripchart(peso ~ carne, df, add=T, vert=T, method='jitter',pch=20)

with(df, tapply(peso, list(carne), function(x){median(x,na.rm=T)}) ) 

##   NO   SI 
## 65.0 68.5

with(df, tapply(peso, list(carne), function(x){mad(x,na.rm=T)}) ) 

##      NO      SI 
## 10.3782 13.3434

Ci saremmo aspettati, nel caso in cui l’assunzione di carne sia legata al sovrappeso, che le pazienti ‘SI’ avessero una mediana di peso maggiore di quella delle pazienti che non ne fanno uso.

Tuttavia la differenza tra le mediane è circa 3 mentre la MAD è circa 10 per entrambe, cioè maggiore della differenza. Possiamo dire con un certo grado di confidenza che le due mediane sono differenti ?

A tale scopo si usa il test statistico di Wilcoxon:

wilcox.test(peso ~ carne, df)

## 
##  Wilcoxon rank sum test with continuity correction
## 
## data:  peso by carne
## W = 14492, p-value = 0.248
## alternative hypothesis: true location shift is not equal to 0

Dato che il p-value è maggiore di 0.05 convenzionalmente si assume che non abbiamo sufficienti informazioni per dire che le due mediane siano differenti.

Viceversa, se il p-value è minore di 0.05 convenzionalmente si assume che ci sia evidenza empirica della differenza delle mediane

Notiamo che, per così dire, il test conferma quanto avevamo già intuito dal grafico.

E’ sempre buona norma visualizzare i dati prima di effettuare i vari test statistici .

Analogamente possiamo visualizzare l’associazione con le altre variabili alimentari:

layout(matrix(1:4,2,2))
boxplot(peso ~ alcol,df)
boxplot(peso ~ salumi,df)
boxplot(peso ~ verdura,df)
boxplot(peso ~ frutta,df)

Oppure possiamo vedere varie combinazioni delle variabili:

layout(matrix(1:2,2,1))
boxplot(BMIc ~ alcol * carne, df)
stripchart(BMIc ~ alcol * carne, df, add=T, vert=T, method='jitter',pch=20)

boxplot(BMIc ~ salumi * verdura, df)
stripchart(BMIc ~ salumi * verdura, df, add=T, vert=T, method='jitter',pch=20)

Associazione con variabli con più livelli

La variabile densità prevede 5 livelli di densità della mammella A, AA, B, C, D.

boxplot(BMIc ~ densita, df)
stripchart(BMIc ~ densita, df, add=T, vert=T, method='jitter',pch=20)

Apparentemente il BMI è più alto per densità A e scende progressivamente per densità B, C, D.

Questa intuizione, che le mediane di BMIc dei vari gruppi sono diverse, può essere verificata usando il test Kruskal-Wallis.

Anche in questo caso il p-value < 0.05 indica che c’è una differenza tra le varie mediane.

kruskal.test(BMIc ~ densita, df)

## 
##  Kruskal-Wallis rank sum test
## 
## data:  BMIc by densita
## Kruskal-Wallis chi-squared = 79.662, df = 4, p-value < 2.2e-16

with(df,
     tapply(BMIc, list(densita), 
            function(x){median(x,na.rm=TRUE)}),
     )

##        A       AA        B        C        D 
## 28.28283 24.45606 26.03749 24.11487 21.82995

with(df,
          tapply(BMIc, list(densita),function(x){mad(x,na.rm=TRUE)})
     )

##        A       AA        B        C        D 
## 6.139905 5.477937 3.933546 3.551748 2.790094

Correlazioni tra variabili numeriche

Vediamo le possibili coppie di variabili numeriche.

pairs(df[,c('peso','altezza','BMIc','eta','CA','eta_menarca','eta_menopausa')])

summary(df[,c('peso','altezza','BMIc','eta','CA','eta_menarca','eta_menopausa')])

##       peso           altezza           BMIc            eta       
##  Min.   :  9.00   Min.   :145.0   Min.   :17.91   Min.   : 2.00  
##  1st Qu.: 60.00   1st Qu.:158.0   1st Qu.:22.57   1st Qu.:18.00  
##  Median : 66.00   Median :162.0   Median :25.10   Median :54.00  
##  Mean   : 68.34   Mean   :162.2   Mean   :26.03   Mean   :43.36  
##  3rd Qu.: 74.50   3rd Qu.:167.0   3rd Qu.:28.44   3rd Qu.:71.00  
##  Max.   :166.00   Max.   :181.0   Max.   :63.25   Max.   :72.00  
##  NA's   :6        NA's   :8       NA's   :8       NA's   :482    
##        CA         eta_menarca    eta_menopausa  
##  Min.   : 44.0   Min.   : 9.00   Min.   : 1.00  
##  1st Qu.: 97.0   1st Qu.:11.00   1st Qu.:45.00  
##  Median :101.0   Median :12.00   Median :50.00  
##  Mean   :101.2   Mean   :12.12   Mean   :45.78  
##  3rd Qu.:106.0   3rd Qu.:13.00   3rd Qu.:53.00  
##  Max.   :150.0   Max.   :17.00   Max.   :60.00  
##  NA's   :32      NA's   :8       NA's   :162

E’ possibile calcolare simultaneamente tutte le correlazioni. Si possono anche visualizzare i modo da individuare a colpo d’occhio le pù significative.

var_num <- df[,c('peso','altezza','BMIc','eta','CA','eta_menarca','eta_menopausa')]
M <- cor(var_num,method='spearman',use='pairwise.complete')
library(corrplot)

## corrplot 0.92 loaded

corrplot(M,method='ellipse')
corrplot(M,add=TRUE,method='number',type='upper',tl.pos='n')

Visualizzazione di dati categorici

Anche nel caso di variabili categoriche e possibile visualizzare i dati in modo da enfatizzare le relazioni.

Consideriamo la variabile densita. Possiamo vedere i dati in vari modi. Innanzitutto numericamente, cioè contando quanti elementi delle varie categorie sono presenti.

M_dens <- with(df, 
                table(densita)
              )
print(M_dens)

## densita
##   A  AA   B   C   D 
##  47   3 206 128  84

vediamo che la classe B è preponderante.

Possiamo visualizzare questa informazione per un maggiore impatto comunicativo.

barplot(M_dens,xlab='densita',ylab='# pazienti')

Abbiamo aggiunto etichette agli assi. E’ buona norma su qualunque grafico aggiungere etichette sugli assi.

Possiamo calcolare la percentuale di ciascuna classe. La funzione addmargins aggiunge la somma delle percentuali per farci capire che effettivamente è 1. La funzione format fa in modo che la rappresentazione dei numeri si fermi alle prime due cifre significative.

format(addmargins(prop.table(M_dens) * 100) , digits=2)

## densita
##        A       AA        B        C        D      Sum 
## " 10.04" "  0.64" " 44.02" " 27.35" " 17.95" "100.00"

Potremmo essere interessati ad incrociare la densita con la assunzione di carne.

M_dens_carne <- with(df,
                     table(densita,carne)
                     )
print(M_dens_carne)

##        carne
## densita  NO  SI
##      A   42   5
##      AA   3   0
##      B  171  33
##      C  109  19
##      D   66  16

Possiamo visualizzare questa informazione ancora con barplot in quattro modi.

layout(matrix(1:4,2,2))
barplot(M_dens_carne,legend=TRUE,beside=TRUE)
barplot(t(M_dens_carne),legend=TRUE,beside=TRUE)
barplot(M_dens_carne,legend=TRUE,beside=FALSE)
barplot(t(M_dens_carne),legend=TRUE,beside=FALSE)

Se invece dei valori assoluti ci interessano le percentuali:

M_dc_perc <- prop.table(M_dens_carne,margin=2) * 100
print(addmargins(M_dc_perc,margin=1))

##        carne
## densita          NO          SI
##     A    10.7416880   6.8493151
##     AA    0.7672634   0.0000000
##     B    43.7340153  45.2054795
##     C    27.8772379  26.0273973
##     D    16.8797954  21.9178082
##     Sum 100.0000000 100.0000000

layout(matrix(1:2,1,2))
barplot(M_dc_perc,legend=TRUE,beside=TRUE,ylab='%',xlab='consumo di carne')
barplot(M_dc_perc,legend=TRUE,beside=FALSE,ylab='%',xlab='consumo di carne')

Un altra modalità usata per graficare la tabella precedente è il mosaicplot:

mosaicplot(M_dens_carne, shade=TRUE)

Un test appropriato che da indicazioni sull’associazione tra densità e assunzione di carne è il \(\chi^2\) oppure il Fisher exact test:

chisq.test(M_dens_carne)

## Warning in chisq.test(M_dens_carne): Chi-squared approximation may be
## incorrect

## 
##  Pearson's Chi-squared test
## 
## data:  M_dens_carne
## X-squared = 2.4703, df = 4, p-value = 0.65

fisher.test(M_dens_carne)

## 
##  Fisher's Exact Test for Count Data
## 
## data:  M_dens_carne
## p-value = 0.7339
## alternative hypothesis: two.sided

il p-value > 0.05 ci dice che non c’è associazione.

Visualizzare più di 2 variabili

Spesso occorre considerare più di due variabili e la visualizzazione in tal caso non è possibile facilmente.

Consideriamo prima il caso di variabili numeriche.

Una prima modalità è di considerare tutte le coppie di variabili come si fa con pairs (le coppie possibili sono Nx(N-1) / 2, per esempio con N=5 abbiamo 5x4/2=10 coppie possibili).

Sono state sviluppate altre modalità.

stars(var_num[1:100,], col.stars = c(1:5)[df$gravidanza_30])

stars(iris[,1:4], col.stars = c(1:3)[iris$Species])

library(MASS)
parcoord(var_num[1:400,],col=c(1:6)[df$menopausa])

library(MASS)
parcoord(iris[,1:4],col=c(1:6)[iris$Species])

biplot(pcamod<-prcomp(scale(na.omit(var_num[,-c(4,6)]))))

plot(pcamod)

biplot(pcamod<-prcomp(scale(na.omit(iris[,1:4]))))

plot(pcamod)

library(caret)

## Loading required package: ggplot2

## Loading required package: lattice

featurePlot(x=iris[,c(1:4)],y=iris[,5], plot="pairs" , pch=20)

plot(cmdscale(dist(scale(iris[,-5]))),col=c(1:3)[iris[,5]],xlab='mds1',yla='mds2')

plot(hclust(dist(iris[,1:4])))

library(e1071)

M <- impute(var_num,what='median')

heatmap(scale(M))

heatmap(scale(as.matrix(iris[,1:4])))

library(lattice)
xyplot(CA ~ peso | densita, df, type=c('p','r'))

xyplot(peso ~ densita | verdura*frutta, df, type=c('p','r'), strip=strip.custom(strip.names=TRUE))

xyplot(BMIc ~ eta_menarca   | densita, df, type=c('p','r'))

library(lattice)
xyplot(CA ~ peso , df, type=c('p','r'), group=densita, auto.key = list(columns=5))

library(lattice)

xyplot(Sepal.Length ~ Sepal.Width , iris, type=c('p','r'), group=Species, auto.key = list(columns=3))

Pate II: Introduzione alle tecniche di Intelligenza artificiale

Misure di Performance

Ora che abbiamo imparato ad esplorare i dati ed abbiamo dimestichezza con essi vediamo qualche tecnica di machine learning elementare.

Nel seguito vedremo alcuni algoritmi di classificazione. Non è possibile entrare nel dettaglio di funzionamento di tali algoritmi nel breve spazio di queste lezioni. Ci limiteremo ad introdurre qualche breve cenno sui fondamenti e le istruzioni R per ottenere dei risultati.

Questi algoritmi basati su machine-learning cercano di assegnare un individuo ad una classe (classificazione). L’assegnazione è fatta dopo un apprendimento. L’apprendimento coinvolge un certo numero di individui di esempio (training-set).

Quando la macchina è addestrata, si cerca di misurare le prestazioni cioè la capacità di generalizzare il riconoscimento. Infatti una macchina può diventare così brava a riconoscere gli individui del training-set ma commettere errori quando cerca di classificare nuovi individui mai visti prima.

Per misurare le prestazioni si usa un set di individui che la macchina non ha mai visto (test-set). Si applica l’algoritmo di riconoscimento della macchina addestrata e si calcolano le prestazioni.

Allo scopo di calcolare le prestazioni si premettono le seguenti 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.

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 Sensitvity 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.

library(pROC)

# creo dataset binario
df$dens2 <- ifelse(df$densita %in% c('A','AA','B'), 'd1','d2')
  
plot(myroc <- roc(dens2 ~ peso,df))
text(.5,.1,paste('AUC',format(myroc$auc, digits=2)) ) 

# oppur enel caso di IRIS
copia <- iris
copia$myc <- ifelse(copia$Species=='setosa','c1','c2')

plot(myroc <- roc(myc ~ Sepal.Length,copia))
text(.5,.1,paste('AUC',format(myroc$auc, digits=2)) ) 

Pre-processing di base

In genere si scalano le features numeriche in modo da avere media 0 e varianza 1.

Inoltre si rimuovono le feature con alta correlazione tra loro e quelle che hanno varianza quasi zero.

Anche le PCA sono usate come tecnica di riduzione.

var_num_c <- var_num
var_num_c <- impute(var_num_c, what='median')

var_num_c = scale(var_num_c)
apply(var_num_c,2,mean)

##          peso       altezza          BMIc           eta            CA 
##  5.055416e-16 -5.231840e-16  2.131149e-16 -2.691955e-16 -2.577760e-16 
##   eta_menarca eta_menopausa 
## -5.677965e-17  2.025703e-16

apply(var_num_c,2,sd)

##          peso       altezza          BMIc           eta            CA 
##             1             1             1             1             1 
##   eta_menarca eta_menopausa 
##             1             1

nz <- nearZeroVar(var_num_c) # quali sono le featue a varianza zero ? 
var_num_c <- var_num_c[,-nz]

CC <- cor(var_num_c, use='complete')
corrplot(CC, method='ellipse')

hh <- findCorrelation(CC) # quali sono le cololle da elimianre ? 

var_num_c <- var_num_c[,-hh]
pca.pre <- prcomp(var_num_c) # tolgo la feature trovata prima

biplot(pca.pre)

var_num_c <- data.frame(var_num_c)
pairs(var_num_c)

Feature selection

Prima di effettuare una classificazione è utile individuare le feature che hanno un maggiore potere discriminativo. Ciò è utile principalmente per due motivi:

• feature non essenziali (non-informativi) possono ridurre le prestazioni del classificatore
• poche features sono più facilmente interpretabili dagli esseri umani (importante in medicina)

La rimozione di features non-informative o ridondanti può essere fatta essenzialmente con due approcci:

• metodi di tipo wrapper
  • il classificatore è condierato in combinzione con le feature: le feature si aggiungono e si tolgono provando una molteciplità di combinaizoni al fine di ottimizzare le prestazioni complessive
• metodi di tipo filter
  • le features sono valutate indipendentemente dal classificatore
  • la classe degli individui del data-set è tenuta in considerazione
  • in genere si usa la ROC, anche per multi-classe

library(caret)

filterVarImp(iris[,1:4],y=iris$Species)

##              setosa versicolor virginica
## Sepal.Length 0.9846     0.9326    0.9846
## Sepal.Width  0.9248     0.9248    0.8344
## Petal.Length 1.0000     1.0000    1.0000
## Petal.Width  1.0000     1.0000    1.0000

filterVarImp(var_num_c,df$densita)

##                       A        AA         B         C         D
## altezza       0.6241135 0.6241135 0.6241135 0.6241135 0.5408910
## BMIc          0.6312057 0.7303856 0.8544833 0.6312057 0.7303856
## CA            0.6118571 0.7562334 0.8234549 0.5992908 0.7562334
## eta_menarca   0.6489362 0.6489362 0.6489362 0.6941748 0.5866855
## eta_menopausa 0.5780142 0.5780142 0.5780142 0.6165049 0.5442985

Regressione logistica

Si usa per problemi a due classi. Stima la probabilità che una combinazione di parametri possa prevedere una classe. In particolare stima i coefficienti \(\beta_k\)

\(\log ( \frac{p}{1-p} ) = \sum_k \beta_k x_k\)

dove \(x_k\) sono le variabili predittori.

dum <- cbind(var_num_c, densita = df$dens2)
dum$densita <- ifelse(dum$densita == 'd1', 1,0)

modello.glm <- glm(densita ~ ., dum, family = binomial)


previsioni <- predict(modello.glm, type='response') 
previsioni <- ifelse(previsioni > 0.5, 1, 0)

confusionMatrix(as.factor(previsioni), as.factor(dum$densita))

## Confusion Matrix and Statistics
## 
##           Reference
## Prediction   0   1
##          0 156  90
##          1  81 166
##                                           
##                Accuracy : 0.6531          
##                  95% CI : (0.6093, 0.6951)
##     No Information Rate : 0.5193          
##     P-Value [Acc > NIR] : 1.264e-09       
##                                           
##                   Kappa : 0.3062          
##                                           
##  Mcnemar's Test P-Value : 0.5407          
##                                           
##             Sensitivity : 0.6582          
##             Specificity : 0.6484          
##          Pos Pred Value : 0.6341          
##          Neg Pred Value : 0.6721          
##              Prevalence : 0.4807          
##          Detection Rate : 0.3164          
##    Detection Prevalence : 0.4990          
##       Balanced Accuracy : 0.6533          
##                                           
##        'Positive' Class : 0               
## 

barplot(coef(modello.glm)) # importanza delle varie variabili

exp(coef(modello.glm))

##   (Intercept)       altezza          BMIc            CA   eta_menarca 
##     1.1195943     0.9203802     1.8144418     1.2569953     0.8990015 
## eta_menopausa 
##     1.0312485

plot(roc(response=dum$densita,predictor = previsioni))

## Setting levels: control = 0, case = 1

## Setting direction: controls < cases

Linear Discriminant Analisys (LDA)

Un algoritmo che vediamo è la LDA che è un analisi lineare (cioè in cui le features vengono solo moltiplicate per dei valori numerici e sommate tra loro) come la PCA.

Non entriamo nel dettaglio matematico ma ci limitiamo a vedere come si realizza in R. Immaginiamo che i nostri esemplari siano le pazienti con le varie densità e le features sono peso, altezza, eta, eta_menarca.

poichè la classe AA è fatta solo 3 esemplari li rimuoviamo dal data-set. Inoltre, una regola da tenere presente (che ha un presupposto matematico/teorico) è che per classificare N classi ci vogliono almeno N-1 features.

library(MASS)

features <- df[,c('peso','altezza','CA')]
classi <- df$densita

# togliamo i valori NA e di densita AA
indici_NA <- which(is.na(features$peso) | is.na(features$altezza) | is.na(features$CA) |  is.na(classi)  | classi=='AA')
features <- features[-indici_NA,]
classi <- classi[-indici_NA]
classi <- droplevels(classi)

modello.lda <- lda(x=features,grouping=classi)

modello.lda

## Call:
## lda(features, grouping = classi)
## 
## Prior probabilities of groups:
##         A         B         C         D 
## 0.1029748 0.4233410 0.2883295 0.1853547 
## 
## Group means:
##       peso  altezza        CA
## A 78.20000 161.7778 106.15556
## B 70.74595 161.7892 103.08108
## C 66.05159 162.4683  98.52381
## D 60.75309 163.1481  97.85185
## 
## Coefficients of linear discriminants:
##                 LD1         LD2          LD3
## peso    -0.07019014 -0.06488776 -0.010147164
## altezza  0.07541983  0.01674500 -0.148422700
## CA      -0.02202513  0.11988640 -0.007893631
## 
## Proportion of trace:
##    LD1    LD2    LD3 
## 0.9376 0.0575 0.0050

plot(modello.lda)

predizioni <- predict(modello.lda, type='class')

# matrice di confusione
table(predizioni$class, classi)

##    classi
##       A   B   C   D
##   A   4   2   1   0
##   B  38 155  90  37
##   C   3  22  32  29
##   D   0   6   3  15

Oppure con colori

pairs(predizioni$x, col=c(1:3)[classi])

Vediamo un altro data-set molto famoso il cosiddetto iris di Fisher. Questo data-set ha delle caratteristiche che consentono di illustrare meglio il funzionamento delle tecniche.

features <- iris[,c(1:4)]
classi <- iris[,5]

modello.lda <- lda(x=features,grouping=classi)

modello.lda

## Call:
## lda(features, grouping = classi)
## 
## Prior probabilities of groups:
##     setosa versicolor  virginica 
##  0.3333333  0.3333333  0.3333333 
## 
## Group means:
##            Sepal.Length Sepal.Width Petal.Length Petal.Width
## setosa            5.006       3.428        1.462       0.246
## versicolor        5.936       2.770        4.260       1.326
## virginica         6.588       2.974        5.552       2.026
## 
## Coefficients of linear discriminants:
##                     LD1         LD2
## Sepal.Length  0.8293776 -0.02410215
## Sepal.Width   1.5344731 -2.16452123
## Petal.Length -2.2012117  0.93192121
## Petal.Width  -2.8104603 -2.83918785
## 
## Proportion of trace:
##    LD1    LD2 
## 0.9912 0.0088

plot(modello.lda)

predizioni <- predict(modello.lda, type='class')

# matrice di confusione
table(predizioni$class, classi)

##             classi
##              setosa versicolor virginica
##   setosa         50          0         0
##   versicolor      0         48         1
##   virginica       0          2        49

Oppure possiamo visualizzare mediante colori

plot(predizioni$x, col=c(1:3)[classi])

pairs(features,col=c(1:3)[classi],pch=c(1:3)[predizioni$class])

library(caret)

confusionMatrix(predizioni$class, classi)

## Confusion Matrix and Statistics
## 
##             Reference
## Prediction   setosa versicolor virginica
##   setosa         50          0         0
##   versicolor      0         48         1
##   virginica       0          2        49
## 
## Overall Statistics
##                                           
##                Accuracy : 0.98            
##                  95% CI : (0.9427, 0.9959)
##     No Information Rate : 0.3333          
##     P-Value [Acc > NIR] : < 2.2e-16       
##                                           
##                   Kappa : 0.97            
##                                           
##  Mcnemar's Test P-Value : NA              
## 
## Statistics by Class:
## 
##                      Class: setosa Class: versicolor Class: virginica
## Sensitivity                 1.0000            0.9600           0.9800
## Specificity                 1.0000            0.9900           0.9800
## Pos Pred Value              1.0000            0.9796           0.9608
## Neg Pred Value              1.0000            0.9802           0.9899
## Prevalence                  0.3333            0.3333           0.3333
## Detection Rate              0.3333            0.3200           0.3267
## Detection Prevalence        0.3333            0.3267           0.3400
## Balanced Accuracy           1.0000            0.9750           0.9800

Il calcolo delle prestazioni con Cross-Validation può essere fatta con il codice seguente:

modello.lda <- lda(x=features,grouping=classi, CV =TRUE)

confusionMatrix(modello.lda$class,iris$Species)

## Confusion Matrix and Statistics
## 
##             Reference
## Prediction   setosa versicolor virginica
##   setosa         50          0         0
##   versicolor      0         48         1
##   virginica       0          2        49
## 
## Overall Statistics
##                                           
##                Accuracy : 0.98            
##                  95% CI : (0.9427, 0.9959)
##     No Information Rate : 0.3333          
##     P-Value [Acc > NIR] : < 2.2e-16       
##                                           
##                   Kappa : 0.97            
##                                           
##  Mcnemar's Test P-Value : NA              
## 
## Statistics by Class:
## 
##                      Class: setosa Class: versicolor Class: virginica
## Sensitivity                 1.0000            0.9600           0.9800
## Specificity                 1.0000            0.9900           0.9800
## Pos Pred Value              1.0000            0.9796           0.9608
## Neg Pred Value              1.0000            0.9802           0.9899
## Prevalence                  0.3333            0.3333           0.3333
## Detection Rate              0.3333            0.3200           0.3267
## Detection Prevalence        0.3333            0.3267           0.3400
## Balanced Accuracy           1.0000            0.9750           0.9800

In alternativa, con il paccehtto CARET si può eseguire la corss-validation nel seguente modo:

mioControl <- trainControl(method='cv',number=5)

lda.fit <- train(Species ~ . , data=iris,
                 method='lda',
                 trControl=mioControl)

lda.fit

## Linear Discriminant Analysis 
## 
## 150 samples
##   4 predictor
##   3 classes: 'setosa', 'versicolor', 'virginica' 
## 
## No pre-processing
## Resampling: Cross-Validated (5 fold) 
## Summary of sample sizes: 120, 120, 120, 120, 120 
## Resampling results:
## 
##   Accuracy  Kappa
##   0.98      0.97

confusionMatrix(lda.fit)

## Cross-Validated (5 fold) Confusion Matrix 
## 
## (entries are percentual average cell counts across resamples)
##  
##             Reference
## Prediction   setosa versicolor virginica
##   setosa       33.3        0.0       0.0
##   versicolor    0.0       32.0       0.7
##   virginica     0.0        1.3      32.7
##                           
##  Accuracy (average) : 0.98

K-nearest neighbours (k-NN)

Un altro algoritmo il k-NN richiede che vi siano una serie di individui prescelti che costituiscono un insieme di addestramento (train-test) e gli altri individui (insieme di test) sono classificati in base alla vicinanza agli individui del train-set.

Si prendono i k elementi più vicini all’individuo incognito e si vede la classe più frequente.

Il valore di k può influenzare le prestazioni della classificazione

La vicinanza è intesa come minima distanza Euclidea.

Come data-set usiamo prima un famoso data-set il Fisher-Iris

library(class)

set.seed(100)
indici <- sample(1:150, 50)# training set
train <- iris[indici, 1:4]
test  <- iris[-indici, 1:4]
cl    <- iris[indici,'Species']
cl_test <-   iris[-indici,'Species']

library(caret) # per matrice di confusione

modello.knn <- knn(train,test,cl,k=1)
confusionMatrix(modello.knn, cl_test)

## Confusion Matrix and Statistics
## 
##             Reference
## Prediction   setosa versicolor virginica
##   setosa         33          0         0
##   versicolor      0         29         1
##   virginica       0          4        33
## 
## Overall Statistics
##                                           
##                Accuracy : 0.95            
##                  95% CI : (0.8872, 0.9836)
##     No Information Rate : 0.34            
##     P-Value [Acc > NIR] : < 2.2e-16       
##                                           
##                   Kappa : 0.925           
##                                           
##  Mcnemar's Test P-Value : NA              
## 
## Statistics by Class:
## 
##                      Class: setosa Class: versicolor Class: virginica
## Sensitivity                   1.00            0.8788           0.9706
## Specificity                   1.00            0.9851           0.9394
## Pos Pred Value                1.00            0.9667           0.8919
## Neg Pred Value                1.00            0.9429           0.9841
## Prevalence                    0.33            0.3300           0.3400
## Detection Rate                0.33            0.2900           0.3300
## Detection Prevalence          0.33            0.3000           0.3700
## Balanced Accuracy             1.00            0.9319           0.9550

modello.knn <- knn(train,test,cl,k=3)
confusionMatrix(modello.knn, cl_test)

## Confusion Matrix and Statistics
## 
##             Reference
## Prediction   setosa versicolor virginica
##   setosa         33          0         0
##   versicolor      0         29         1
##   virginica       0          4        33
## 
## Overall Statistics
##                                           
##                Accuracy : 0.95            
##                  95% CI : (0.8872, 0.9836)
##     No Information Rate : 0.34            
##     P-Value [Acc > NIR] : < 2.2e-16       
##                                           
##                   Kappa : 0.925           
##                                           
##  Mcnemar's Test P-Value : NA              
## 
## Statistics by Class:
## 
##                      Class: setosa Class: versicolor Class: virginica
## Sensitivity                   1.00            0.8788           0.9706
## Specificity                   1.00            0.9851           0.9394
## Pos Pred Value                1.00            0.9667           0.8919
## Neg Pred Value                1.00            0.9429           0.9841
## Prevalence                    0.33            0.3300           0.3400
## Detection Rate                0.33            0.2900           0.3300
## Detection Prevalence          0.33            0.3000           0.3700
## Balanced Accuracy             1.00            0.9319           0.9550

modello.knn <- knn(train,test,cl,k=5)
confusionMatrix(modello.knn, cl_test)

## Confusion Matrix and Statistics
## 
##             Reference
## Prediction   setosa versicolor virginica
##   setosa         33          0         0
##   versicolor      0         30         1
##   virginica       0          3        33
## 
## Overall Statistics
##                                          
##                Accuracy : 0.96           
##                  95% CI : (0.9007, 0.989)
##     No Information Rate : 0.34           
##     P-Value [Acc > NIR] : < 2.2e-16      
##                                          
##                   Kappa : 0.94           
##                                          
##  Mcnemar's Test P-Value : NA             
## 
## Statistics by Class:
## 
##                      Class: setosa Class: versicolor Class: virginica
## Sensitivity                   1.00            0.9091           0.9706
## Specificity                   1.00            0.9851           0.9545
## Pos Pred Value                1.00            0.9677           0.9167
## Neg Pred Value                1.00            0.9565           0.9844
## Prevalence                    0.33            0.3300           0.3400
## Detection Rate                0.33            0.3000           0.3300
## Detection Prevalence          0.33            0.3100           0.3600
## Balanced Accuracy             1.00            0.9471           0.9626

pairs(train,col=c(1:3)[cl])

# nella figura seguente il colore è sempre la classe originale
# il pallino pieno indica i dati di train
plot(train[,c(3,4)],col=c(1:3)[cl], pch=20,xlim=range(c(train[,3],test[,3])), ylim=range(c(train[,4],test[,4])))
points(test[,c(3,4)],col=c(1:3)[cl_test],pch=c(1:3)[modello.knn])
legend('bottomright', legend=c('setosa tr','versic tr','virgin tr','setosa te','versic te','virgin te'), pch=c(20,20,20,1:3),col=c(1:3,1:3))

Classificationa and Regression Tree (CART)

A differenza del LDA che è una tecnica lineare in cui la linea di separazione tra le classi è una retta (o un piano nel caso di più dimensioni), i CART sono non lineari e la linea di separazione è piece-wise linear.

Gli alberi decisionali che si ottengono sono molto usati in ambito medicale.

features <- df[,c('peso','altezza','CA')]
classi <- df$densita

# togliamo i valori NA e di densita AA
indici_NA <- which(is.na(features$peso) | is.na(features$altezza) | is.na(features$CA) |  is.na(classi)  | classi=='AA')
features <- features[-indici_NA,]
classi <- classi[-indici_NA]
classi <- droplevels(classi)

d <- cbind(features,classi)

library(rpart)
modello.rpart <- rpart(classi ~ ., d)

par(mfrow = c(1,2), xpd = NA) # otherwise on some devices the text is clipped
plot(modello.rpart)
text(modello.rpart, use.n = TRUE)

library(rpart.plot)
rpart.plot(modello.rpart)

predizioni <- predict(modello.rpart,type='class')

confusionMatrix(predizioni, classi)

## Confusion Matrix and Statistics
## 
##           Reference
## Prediction   A   B   C   D
##          A   0   0   0   0
##          B  42 161  84  35
##          C   1  11  26  11
##          D   2  13  16  35
## 
## Overall Statistics
##                                           
##                Accuracy : 0.508           
##                  95% CI : (0.4601, 0.5558)
##     No Information Rate : 0.4233          
##     P-Value [Acc > NIR] : 0.0002212       
##                                           
##                   Kappa : 0.2163          
##                                           
##  Mcnemar's Test P-Value : < 2.2e-16       
## 
## Statistics by Class:
## 
##                      Class: A Class: B Class: C Class: D
## Sensitivity             0.000   0.8703   0.2063  0.43210
## Specificity             1.000   0.3611   0.9260  0.91292
## Pos Pred Value            NaN   0.5000   0.5306  0.53030
## Neg Pred Value          0.897   0.7913   0.7423  0.87601
## Prevalence              0.103   0.4233   0.2883  0.18535
## Detection Rate          0.000   0.3684   0.0595  0.08009
## Detection Prevalence    0.000   0.7368   0.1121  0.15103
## Balanced Accuracy       0.500   0.6157   0.5662  0.67251

features <- iris[,c(1:4)]
classi <- iris[,5]

d <- cbind(features,classi)

library(rpart)
modello.rpart <- rpart(classi ~ ., d)

par(mfrow = c(1,2), xpd = NA) # otherwise on some devices the text is clipped
plot(modello.rpart)
text(modello.rpart, use.n = TRUE)

rpart.plot(modello.rpart)

predizioni <- predict(modello.rpart, type='class')

confusionMatrix(predizioni, classi)

## Confusion Matrix and Statistics
## 
##             Reference
## Prediction   setosa versicolor virginica
##   setosa         50          0         0
##   versicolor      0         49         5
##   virginica       0          1        45
## 
## Overall Statistics
##                                          
##                Accuracy : 0.96           
##                  95% CI : (0.915, 0.9852)
##     No Information Rate : 0.3333         
##     P-Value [Acc > NIR] : < 2.2e-16      
##                                          
##                   Kappa : 0.94           
##                                          
##  Mcnemar's Test P-Value : NA             
## 
## Statistics by Class:
## 
##                      Class: setosa Class: versicolor Class: virginica
## Sensitivity                 1.0000            0.9800           0.9000
## Specificity                 1.0000            0.9500           0.9900
## Pos Pred Value              1.0000            0.9074           0.9783
## Neg Pred Value              1.0000            0.9896           0.9519
## Prevalence                  0.3333            0.3333           0.3333
## Detection Rate              0.3333            0.3267           0.3000
## Detection Prevalence        0.3333            0.3600           0.3067
## Balanced Accuracy           1.0000            0.9650           0.9450

Combinazione di alberi : Bagging e Random Forest

Per rendere più stabili e generalizzabili i classificatori ad albero ed evitare l’overfitting, si può adottare la tecnica del BAGGING (bootstrap aggregating) si può generare un certo numero di alberi adestrati su un training set leggermente diversi tra loro ed in ciascuno di questi fare l’addestramento. Ogni albero darà un suo risultato per la classificazione: la combinaizone di tutti questi risultati (per esempio con un meccanismo di voto) ci darà la classificazione finale.

La tecnica del Random Forest aggiunge al Bagging il fatto che i predittori (le feature ) non sono le stesse per gli alberi ma sono causalmente selezionate.

Il parametro mtry nel pacchetto CARET è il numero di predittori da selezionare casualmente.

mioControl <- trainControl(method='cv',number=5)


rf.fit <- train(Species ~ ., data = iris,
                method='rf',
                ntree=50,
                trControl=mioControl)
rf.fit

## Random Forest 
## 
## 150 samples
##   4 predictor
##   3 classes: 'setosa', 'versicolor', 'virginica' 
## 
## No pre-processing
## Resampling: Cross-Validated (5 fold) 
## Summary of sample sizes: 120, 120, 120, 120, 120 
## Resampling results across tuning parameters:
## 
##   mtry  Accuracy   Kappa
##   2     0.9266667  0.89 
##   3     0.9200000  0.88 
##   4     0.9200000  0.88 
## 
## Accuracy was used to select the optimal model using the largest value.
## The final value used for the model was mtry = 2.

Sui dati clinici i risultati sono i seguenti:

mioControl <- trainControl(method='cv',number=5)

dum <- cbind(var_num_c, densita = df$dens2)

rf.fit <- train(densita ~ ., data = dum,
                method='rf',
                ntree=50,
                trControl=mioControl)

rf.fit

## Random Forest 
## 
## 493 samples
##   5 predictor
##   2 classes: 'd1', 'd2' 
## 
## No pre-processing
## Resampling: Cross-Validated (5 fold) 
## Summary of sample sizes: 393, 395, 395, 395, 394 
## Resampling results across tuning parameters:
## 
##   mtry  Accuracy   Kappa    
##   2     0.6128044  0.2228975
##   3     0.5924152  0.1816424
##   5     0.5903335  0.1769111
## 
## Accuracy was used to select the optimal model using the largest value.
## The final value used for the model was mtry = 2.

confusionMatrix(rf.fit)

## Cross-Validated (5 fold) Confusion Matrix 
## 
## (entries are percentual average cell counts across resamples)
##  
##           Reference
## Prediction   d1   d2
##         d1 33.9 20.7
##         d2 18.1 27.4
##                             
##  Accuracy (average) : 0.6126

Neural Networks (NN)

Una rete neurale è una struttura non lineare che cerca di emulare il funzionamento dei neuroni.

Una rete neurale molto semplice è costituita da uno strato di neuroni di ingresso, uno strato nascosto ed uno strato di uscita.

copia <- iris

copia$Sepal.Length <- scale(copia$Sepal.Length)
copia$Sepal.Width <- scale(copia$Sepal.Width)
copia$Petal.Length <- scale(copia$Petal.Length)
copia$Petal.Width <- scale(copia$Petal.Width)

library(nnet)

modello.nnet <- nnet(Species ~ ., data=copia, size=5)

## # weights:  43
## initial  value 184.066113 
## iter  10 value 12.146570
## iter  20 value 3.161392
## iter  30 value 2.512449
## iter  40 value 2.502016
## final  value 2.502012 
## converged

library(NeuralNetTools)
plotnet(modello.nnet)

predizioni <- as.factor(predict(modello.nnet, type='class'))

confusionMatrix(predizioni, copia$Species)

## Confusion Matrix and Statistics
## 
##             Reference
## Prediction   setosa versicolor virginica
##   setosa         50          0         0
##   versicolor      0         49         0
##   virginica       0          1        50
## 
## Overall Statistics
##                                           
##                Accuracy : 0.9933          
##                  95% CI : (0.9634, 0.9998)
##     No Information Rate : 0.3333          
##     P-Value [Acc > NIR] : < 2.2e-16       
##                                           
##                   Kappa : 0.99            
##                                           
##  Mcnemar's Test P-Value : NA              
## 
## Statistics by Class:
## 
##                      Class: setosa Class: versicolor Class: virginica
## Sensitivity                 1.0000            0.9800           1.0000
## Specificity                 1.0000            1.0000           0.9900
## Pos Pred Value              1.0000            1.0000           0.9804
## Neg Pred Value              1.0000            0.9901           1.0000
## Prevalence                  0.3333            0.3333           0.3333
## Detection Rate              0.3333            0.3267           0.3333
## Detection Prevalence        0.3333            0.3267           0.3400
## Balanced Accuracy           1.0000            0.9900           0.9950

Il numero ottimale di neuroni dello strato interno non è determinabiel con una semplice formula. Si può procedere per tentativi come segue. La configurazione iniziale della rete è stabilita casualmente per cui ogni volta che si esegue la rete il risultato è diverso. Per ovviare a questo problema ho stabilito il ‘seme’ iniziale del ’generatore di casualità del computer in modo che il grafico visualizzato in ueste slide sia unico. Per esaminare più a fondo le prestazioni della rete bisogna svolgere ulteriori analisi che esulano daqeusto corso.

# divido in train e test
set.seed(1000)
indici <- sample(1:150,100)
train_data <- copia[indici,]
test_data  <- copia[-indici,]

# faccio variare il numero di neuroni da 1 a N
# ogni volta mi salvo l'accuratezzas
N <- 10
acc <- rep(NA,N) 

for (k in 2:N){
  modello.nnet <- nnet(Species ~ ., data=train_data, size=k, trace=FALSE)
  predizioni <- as.factor(predict(modello.nnet, newdata = test_data, type='class'))
  #print(length(levels(predizioni)))
  dum <- confusionMatrix(predizioni, test_data$Species)
  acc[k] <- dum$overall[3] # lower bound accuracy

}

plot(acc, ylim=c(0.1,1))
abline(h=0.95,lty=3)

Proviamo sul dataset clinico relativametne alla densita:

var_num_c2 <- as.data.frame(impute(var_num, what='median'))
clinici <- cbind(var_num_c2,densita=df$densita)
# scarto eta e eta_menopausa perchè troppi NA,
# elimino dove densita non nota e AA
# inserisco valore mediano laddove NA
clinici <- clinici[,-c(4,7)]
clinici <- clinici[-which(is.na(clinici$densita)|clinici$densita=='AA'),]
clinici$densita <- droplevels(clinici$densita)


# divido in train e test
set.seed(1000)
L <- dim(clinici)[1]
indici <- sample(1:L,.75 * L )
train_data <- clinici[indici,]
test_data  <- clinici[-indici,]

# faccio variare il numero di neuroni da 1 a N
# ogni volta mi salvo l'accuratezzas
N <- 10
acc <- rep(NA,N) 

for (k in 2:N){
  modello.nnet <- nnet(densita ~ ., data=train_data, size=k, trace=FALSE)
  predizioni <- as.factor(predict(modello.nnet, newdata = test_data, type='class'))
  dum <- confusionMatrix(predizioni, test_data$densita)
  acc[k] <- dum$overall[3] # lower bound accuracy

}



plot(acc, ylim=c(0.1,1))
abline(h=0.95,lty=3)

Support Vector Machine (SVM)

Le SVM è un algoritmo basato su una variante di quelli lineari: la curva (o superficie) di separazione è ancora una retta ( o piano) ma a differenza di LDA o NN le SVM producono un piano che in un certo senso è il migliore possibile, che massimizza la separazione tra le varie classi.

library(e1071)

set.seed(100)
indici <- sample(1:150, 70)
train <- copia[indici, 1:4]
test  <- copia[-indici, 1:4]
cl    <- copia[indici,'Species']
cl_test <-   copia[-indici,'Species']

library(caret) # per matrice di confusione

modello.svm <- svm(train,cl)
predizioni <- predict(modello.svm, newdata=test,type='class')
confusionMatrix(predizioni, cl_test)

## Confusion Matrix and Statistics
## 
##             Reference
## Prediction   setosa versicolor virginica
##   setosa         26          0         0
##   versicolor      0         24         2
##   virginica       0          1        27
## 
## Overall Statistics
##                                           
##                Accuracy : 0.9625          
##                  95% CI : (0.8943, 0.9922)
##     No Information Rate : 0.3625          
##     P-Value [Acc > NIR] : < 2.2e-16       
##                                           
##                   Kappa : 0.9437          
##                                           
##  Mcnemar's Test P-Value : NA              
## 
## Statistics by Class:
## 
##                      Class: setosa Class: versicolor Class: virginica
## Sensitivity                  1.000            0.9600           0.9310
## Specificity                  1.000            0.9636           0.9804
## Pos Pred Value               1.000            0.9231           0.9643
## Neg Pred Value               1.000            0.9815           0.9615
## Prevalence                   0.325            0.3125           0.3625
## Detection Rate               0.325            0.3000           0.3375
## Detection Prevalence         0.325            0.3250           0.3500
## Balanced Accuracy            1.000            0.9618           0.9557

plot(cmdscale(dist(train)),col=c(1:3)[cl], pch = c(1,3)[1:50 %in% modello.svm$index + 1])

Oppure con cross-validation mediante il pacchetto CARET:

mioControllo = trainControl(method='cv',number=5)

modello.svm2 <- train(Species ~ . , data=copia, method='svmLinear', trControl=mioControllo)

confusionMatrix(modello.svm2)

## Cross-Validated (5 fold) Confusion Matrix 
## 
## (entries are percentual average cell counts across resamples)
##  
##             Reference
## Prediction   setosa versicolor virginica
##   setosa       33.3        0.0       0.0
##   versicolor    0.0       30.7       1.3
##   virginica     0.0        2.7      32.0
##                           
##  Accuracy (average) : 0.96

Clustering o Apprendimento non-supervisionato

Sono algoritmi affinchè la macchina cerchi di imparare da sola i gruppi esistenti all’interno dei dati.

modello.kmeans <- kmeans(iris[,-5], 3)

pairs(iris[,-5], lower.panel= function(x,y,...){points(x,y,col = c(1:3)[modello.kmeans$cluster])} ,
      upper.panel = function(x,y,...){points(x,y,col = c(1:3)[iris[,5]])}) 

modello.kmeans <- kmeans(iris[,-5], 2)

pairs(iris[,-5], lower.panel= function(x,y,...){points(x,y,col = c(1:3)[modello.kmeans$cluster])} ,
      upper.panel = function(x,y,...){points(x,y,col = c(1:3)[iris[,5]])}) 

dati <- scale(na.omit(df[,c('peso','altezza', 'CA', 'eta_menarca')]))
nc <- 4
modello.kmeans.clinici <- kmeans(dati, nc)

pairs(dati,col=c(1:nc)[modello.kmeans.clinici$cluster])

plot(cmdscale(dist(dati)),col=c(1:nc)[modello.kmeans.clinici$cluster],xlab='mds1',ylab='mds2')

library(factoextra)
library(withr)

copia <- scale(iris[,c(1:4)])

fviz_nbclust(copia,kmeans)

copia.k <- kmeans(copia,2)
fviz_cluster(copia.k,copia)

dum <- scale(clinici[,-6])
fviz_nbclust(dum,kmeans)

dum.k <- kmeans(dum,2)
fviz_cluster(dum.k,dum)

Bibliografia

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
• 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