lunedì 8 maggio 2023

Analisi dei gruppi (clustering non gerarchico)

L'analisi dei gruppi si applica a dati multivariati ed è un metodo statistico di tassonomia numerica che riveste un ruolo importante nella analisi esplorativa dei dati.

L'obiettivo dell'analisi dei gruppi (cluster analysis o clustering) è concettualmente semplice: verificare la possibile esistenza, in un insieme di oggetti, di sottoinsiemi di oggetti particolarmente simili tra loro (gruppi/cluster).

Nonostante alla base dell'analisi dei gruppi vi sia un'idea semplice e logica, vi sono numerosi modi per realizzarla [1], qui ci occupiamo dell'implementazione del clustering non gerarchico con metodi esclusivi nelle due versioni:
→ clustering con il metodo di MacQueen (k-means);
→ clustering con il metodo di Rousseew (k-medoids).

Come dati impieghiamo i valori di BMI (indice di massa corporea) rilevati a livello europeo alcuni anni fa e pubblicati dall'Istat [2].

Per proseguire è necessario:
→ effettuare il download del file di dati bmi.csv
→ salvare il file nella cartella C:\Rdati\

Per il file di dati trovate link e modalità di download alla pagina Dati, ma potete anche semplicemente copiare i dati riportati qui sotto aggiungendo un ↵ Invio al termine dell'ultima riga e salvarli in C:\Rdati\ in un file di testo denominato bmi.csv (assicuratevi che il file sia effettivamente salvato con l'estensione .csv).

Nazione;sottopeso;normale;sovrappeso;obeso
Austria;2.4;49.6;33.3;14.7
Belgio;2.7;48.0;35.3;14.0
Bulgaria;2.2;43.8;39.2;14.8
Cipro;3.9;47.8;33.8;14.5
Croazia;1.9;40.7;38.7;18.7
Danimarca;2.2;50.0;32.9;14.9
Estonia;2.2;43.9;33.5;20.4
Finlandia;1.2;44.1;36.4;18.3
Francia;3.2;49.6;31.9;15.3
Germania;1.8;46.1;35.2;16.9
Grecia;1.9;41.3;39.4;17.3
Irlanda;1.9;42.3;37.0;18.7
Lettonia;1.7;41.8;35.2;21.3
Lituania;1.9;42.5;38.3;17.3
Lussemburgo;2.8;49.3;32.4;15.6
Malta;2.0;37.0;35.0;26.0
Olanda;1.6;49.0;36.0;13.3
Polonia;2.4;42.9;37.5;17.2
Portogallo;1.8;44.6;36.9;16.6
Regno Unito;2.1;42.2;35.6;20.1
Repubblica Ceca;1.1;42.1;37.6;19.3
Romania;1.3;42.9;46.4;9.4
Slovacchia;2.1;43.6;38.0;16.3
Slovenia;1.6;41.8;37.4;19.2
Spagna;2.2;45.4;35.7;16.7
Svezia;1.8;48.3;35.9;14.0
Ungheria;2.9;41.9;34.0;21.2

Inoltre è necessario scaricare dal CRAN il pacchetto aggiuntivo cluster [3], il pacchetto aggiuntivo factoextra [4] e il pacchetto aggiuntivo ggplot2 [5]. 

Iniziamo con il clustering con il metodo di MacQueen che impiega l'algoritmo k-means.

Copiate questo script, incollatelo nella Console di R e premete ↵ Invio:

# CLUSTERING (NON GERARCHICO) con il metodo di MacQueen (k-means)
#
library(cluster) # carica il pacchetto
mydata <- read.table("c:/Rdati/bmi.csv", header=TRUE, sep=";", row.names="Nazione") # importa i dati
z <- scale(mydata) # standardizza le variabili
windows() # apre e inizializza una nuova finestra grafica
#
myclust <- kmeans(z, 4, algorithm="MacQueen", nstart=50) # genera i 4 gruppi/cluster
clusplot(z, myclust$cluster, color=TRUE, labels=2, lines=0, main="Grafico dei cluster - metodo di MacQueen (k.means)", xlab="Componente principale 1", sub="", ylab="Componente principale 2", cex=0.6, col.txt="black", col.p="black") # traccia il grafico dei cluster 
#

Con prima cosa viene caricato il pacchetto cluster, quindi i dati sono importati in mydata

Con la funzione scale() della terza riga per ciascun dato viene calcolata la deviata normale standardizzata z. Questa funzione calcola per i dati di ciascuna colonna/variabile la media e la deviazione standard, poi calcola per ciascun dato x la corrispondente deviata normale standardizzata z come

z = (x – media) / deviazione standard

I valori di z sono salvati nel nuovo oggetto qui denominato per comodità mnemonica z.

Se nella Console di R digitate z 

> z 

sono mostrati i dati standardizzati, in coda ai quali sono riportate  la media 

attr(,"scaled:center")
 sottopeso    normale sovrappeso      obeso 
  2.103704  44.537037  36.240741  17.111111 

e la deviazione standard

attr(,"scaled:scale")
 sottopeso    normale sovrappeso      obeso 
 0.6111033  3.3717318  2.8989732  3.2322097 

impiegate per effettuare la standardizzazione dei dati.

Quindi viene aperta e inizializzata una nuova finestra grafica con windows().

La funzione kmeans() impiega i dati standardizzati z per generare 4 cluster, impiegando l'algoritmo "MacQueen" e 50 iterazioni (nstart=50) dell'algoritmo di scelta iniziale dei cluster. I risultati sono salvati in myclust.

I risultati del clustering salvati in myclust sono impiegatia per tracciare il grafico con la funzione clusplot() e i seguenti argomenti:
→ l'oggetto z contenente i dati standardizzati;
→ l'oggetto myclust$cluster contenente i cluster generati con kmeans();
→ color=TRUE per riportare i cluster in colore; 
→ labels=2 per riportare la numerazione assegnata ai cluster;
→ lines=0 per non riportare le linee che collegano i cluster;
→ sub="" che elimina il sottotitolo previsto di default dalla funzione;
→ cex=0.6 per rimpicciolire i caratteri del testo applicato;
→ col.txt che definisce il colore del testo che compare all'interno del grafico;
→ col.p che definisce il colore impiegato per rappresentare i punti all'interno del grafico.


Per semplicità nella funzione clusplot() sono stati lasciati i valori di default per numerosi altri argomenti, digitate help(clusplot) nella Console di R per un approfondimento del tema.

Se nella Console di R digitate mydata[order(-mydata$normale),] 

> mydata[order(-mydata$normale),]

sono mostrati i dati ordinati in ordine decrescente per la colonna che contiene la percentuale di soggetti con peso normale

                sottopeso normale sovrappeso obeso
Danimarca             2.2    50.0       32.9  14.9
Austria               2.4    49.6       33.3  14.7
Francia               3.2    49.6       31.9  15.3
Lussemburgo           2.8    49.3       32.4  15.6
Olanda                1.6    49.0       36.0  13.3
Svezia                1.8    48.3       35.9  14.0
Belgio                2.7    48.0       35.3  14.0
Cipro                 3.9    47.8       33.8  14.5
Germania              1.8    46.1       35.2  16.9
Spagna                2.2    45.4       35.7  16.7
Portogallo            1.8    44.6       36.9  16.6
Finlandia             1.2    44.1       36.4  18.3
Estonia               2.2    43.9       33.5  20.4
Bulgaria              2.2    43.8       39.2  14.8
Slovacchia            2.1    43.6       38.0  16.3
Polonia               2.4    42.9       37.5  17.2
Romania               1.3    42.9       46.4   9.4
Lituania              1.9    42.5       38.3  17.3
Irlanda               1.9    42.3       37.0  18.7
Regno Unito           2.1    42.2       35.6  20.1
Repubblica Ceca       1.1    42.1       37.6  19.3
Ungheria              2.9    41.9       34.0  21.2
Lettonia              1.7    41.8       35.2  21.3
Slovenia              1.6    41.8       37.4  19.2
Grecia                1.9    41.3       39.4  17.3
Croazia               1.9    40.7       38.7  18.7
Malta                 2.0    37.0       35.0  26.0
 
e potete vedere subito come - giusto per chiarire con un esempio - Danimarca, Austria, Francia, Lussemburgo, Olanda, Svezia, Belgio e Cipro, graficamente incluse nello stesso cluster, numericamente si distinguono dalle altre nazione per avere la maggior percentuale di soggetti con peso normale e contemporaneamente una ridotta percentuale di soggetti obesi.

Vediamo ora il clustering con il metodo di Rousseew che impiega l'algoritmo k-medoids.

# CLUSTERING NON GERARCHICO con il metodo di Rousseew (k-medoids)
#  
library(factoextra) # carica il pacchetto
library(cluster) # carica il pacchetto
mydata <- read.table("c:/Rdati/bmi.csv", header=TRUE, sep=";", row.names="Nazione") # importa i dati
windows() # apre e inizializza una nuova finestra grafica
#
myclust <- pam(mydata, 4, metric=c("euclidean"), stand=TRUE) # genera un oggetto pam che contiene i dati dei cluster
fviz_cluster(myclust, myclust$cluster, labelsize=7, main = "Grafico dei cluster - metodo di Rousseew (k-medoids)") # grafico dei cluster per le prime due componenti principali
#
windows() # apre e inizializza una nuova finestra grafica
distance <- get_dist(mydata, method="euclidean", stand=TRUE) # genera la matrice delle distanze 
fviz_dist(distance, order=TRUE, show_labels=TRUE, gradient = list(low="#00AFBB", mid="white", high="#FC4E07")) + theme(axis.text.x=element_text(angle=90, vjust=0.3))+ theme(axis.text.y=element_text(angle=0, vjust=0.3)) # grafico della matrice delle distanze
#
windows() # apre e inizializza una nuova finestra grafica
fviz_nbclust(mydata, FUNcluster=cluster::pam, method="wss") # calcola il numero ottimale di cluster 
#

I preliminari sono gli stessi dello script precedente a parte il fatto che in questo caso viene impiegato anche il pacchetto factoextra che a sua volta si basa sul pacchetto ggplot2 che deve anch'esso essere stato installato.

La prima funzione utilizzata è pam() che deriva il suo nome da "Partitioning Around Medoids", impiega l'algoritmo di clustering dei dati  k-medoids (un versione robusta dell'algoritmo k-means) e prevede come argomenti:
→ i dati mydata da impiegare;
→ il numero di cluster, in questo caso 4, in cui classificare i dati;
→ la metrica da impiegare che può essere "euclidean" o "manhattan";
stand=TRUE che impone la standardizzazione dei dati, che nel caso del clustering con il metodo di MacQueen che abbiamo visto sopra deve invece essere eseguita in via preliminare con la funzione scale().

Con la successiva funzione fviz_cluster() è generato il grafico dei cluster (myclust$cluster) ottenuti con l'algoritmo k-medoids e contenuti nell'oggetto myclust.


La matrice delle distanze distance generata con la funzione get_dist() sui dati (mydata), previa la loro standardizzazione (stand=TRUE), viene impiegata dalla funzione fviz_dist() per generare il grafico delle distanze tra gli oggetti clusterizzati, un grafico facoltativo e accessorio rispetto al precedente, ma che potrebbe interessare.

In diagonale compare in azzurro la distanza 0 della relazione di identità degli oggetti con se stessi. Il colore azzurro si va attenuando via via che inizia la dissimilarità tra gli oggetti, fino a trasformarsi in un rosso sempre più intenso al suo progressivo aumentare.

Il colore rosso, ad esempio, diventa la dominante nell'incrocio della Romania con tutte le altre Nazioni, a conferma della sua peculiarità: se guardate i dati, ha la percentuale in assoluto più elevata di sovrappeso (46.4%) e la percentuale in assoluto più bassa di obesi (9.4%). In termini di dissimilarità dalle restanti Nazioni, la Romania è immediatamente seguita da Malta, come risulta evidente anche nel grafico dei cluster.


Infine con la funzione fviz_nbclust() viene generato il grafico che in teoria consentirebbe di valutare il numero ottimale di cluster da imporre nella funzione pam().


Il numero ottimale di cluster dovrebbe corrispondere ad un punto in cui cambia la curvatura: qui vediamo il primo in corrispondenza di 2 cluster, e il secondo in corrispondenza di 4 cluster. Personalmente non ho mai trovato particolarmente utili gli automatismi di questo genere, preferisco ragionare sul quadro generale, ma lascio al lettore le valutazioni del caso.

Trovate il seguito e le strategie alternative di clustering e di analisi dei dati multivariati nei post:


----------

[1] Nonostante alla base dell'analisi dei gruppi vi sia un'idea semplice e logica, vi sono numerosi modi per realizzarla, infatti esistono:
metodi gerarchici, che danno luogo a una suddivisione ad albero (dendrogramma) in base alla distanza tra i singoli oggetti dell'insieme;
metodi non gerarchici, nei quali l'appartenenza di un oggetto (dell'insieme) ad uno specifico sottoinsieme/gruppo/cluster viene stabilita sulla base della sua distanza dal centro o dalla media dei dati o da un punto rappresentativo del cluster;
→ metodi bottom-up noti anche come metodi agglomerativi nei quali all'inizio del processo di classificazione ad ogni oggetto viene fatto corrispondere un cluster. In questo stadio gli oggetti sono considerati tutti dissimili tra di loro. Al passaggio successivo i due oggetti più simili sono raggruppati nello stesso cluster. Il numero dei cluster risulta quindi pari al numero di oggetti diminuito di uno. Il procedimento viene ripetuto ciclicamente, fino ad ottenere (all'ultimo passaggio) un unico cluster;
→ metodi top-down noti anche come metodi divisivi nei quali inizialmente tutti gli oggetti sono considerati come appartenenti ad un unico cluster, che viene via via suddiviso in cluster fino ad avere un numero di cluster uguale al numero degli oggetti;
metodi esclusivi, che prevedono che un oggetto possa appartenere esclusivamente a un cluster;
→ metodi non esclusivi (fuzzy) che prevedono che un oggetto possa appartenere, in modo quantitativamente diverso, a più di un cluster. 

[2] Vedere il post Indice di massa corporea (BMI).

[3] Trovate la documentazione nel manuale di riferimento del pacchetto: Package ‘cluster’. URL consultato il 06/01/2023.

[4] Trovate la documentazione nel manuale di riferimento del pacchetto: Package ‘factoextra’. URL consultato il 06/01/2023.

[5] Trovate la documentazione nel manuale di riferimento del pacchetto: Package ‘ggplot2’. URL consultato il 06/01/2023.

domenica 5 marzo 2023

Gestione dei dati mancanti

Per illustrare il tema, che è strettamente collegato all'ordinamento dei dati [1], vediamo tre cose:
→ come individuare i dati mancanti;
→ come i dati mancanti - identificati in R con la sigla NA (Not Available) - possono bloccare l'esecuzione di calcoli sulle variabili numeriche che li contengono e come in questo caso sia possibile omettere selettivamente dai calcoli i dati NA
→ come al bisogno sia possibile eliminare per intero da un vettore o matrice o tabella [2] i casi con dati NA.

Per proseguire dovete, seguendo le istruzioni fornite alla pagina Dati:
effettuare il download del file importa_csv.csv
salvare il file nella cartella C:\Rdati\

In alternativa copiate le otto righe riportate qui sotto, incollatele in un editor di file di testo aggiungendo un ↵ Invio al termine dell'ultima riga, e salvate il tutto in C:\Rdati\ in un file di testo denominato importa_csv.csv (attenzione: assicuratevi che il file sia effettivamente salvato con l'estensione .csv):

id;sesso;anni;peso_kg;altezza_m
MT;M;69;76;1,78
GF;F;56;63;
MC;F;53;71;1,60
SB;M;28;73;1,78
FE;F;61;54;1,54
AB;M;46;92;1,84
RF;F;31;81;1,56

Vediamo come individuare i dati mancanti. Copiate e incollate nella Console di R questo script e premete ↵ Invio:

# IDENTIFICA E CONTEGGIA I DATI MANCANTI 
#
# importa i dati, notare / invece di \ su windows
mydata <- read.table("C:/Rdati/importa_csv.csv", header=TRUE, sep=";", dec=",", row.names="id")
#
mydata # mostra i dati importati
#
is.na(mydata) # identifica i dati mancanti
#
colSums(is.na(mydata)) # conteggia i dati mancanti per colonna/variabile
#
which(is.na(mydata$altezza_m)) # identifica la posizione dei dati mancanti
#

Dopo avere importato i dati, richiamando il nome dell'oggetto mydata che li contiene questi sono mostrati con la sigla NA riportata automaticamente da R in corrispondenza di ciascuno dei dati mancanti.

> mydata # mostra i dati importati
   sesso anni peso_kg altezza_m
MT     M   69      76      1.78
GF     F   56      63        NA
MC     F   53      71      1.60
SB     M   28      73      1.78
FE     F   61      54      1.54
AB     M   46      92      1.84
RF     F   31      81      1.56

Nel caso di piccole tabelle come questa non serve altro: ma nel caso di estesi database ci servono funzioni che forniscano qualche soluzione più pratica e immediata. 

Impieghiamo quindi la funzione is.na() che identifica nella tabella mydata i dati mancanti riportando TRUE in corrispondenza di ciascuno di essi 

> is.na(mydata) # identifica i dati mancanti
   sesso  anni peso_kg altezza_m
MT FALSE FALSE   FALSE     FALSE
GF FALSE FALSE   FALSE      TRUE
MC FALSE FALSE   FALSE     FALSE
SB FALSE FALSE   FALSE     FALSE
FE FALSE FALSE   FALSE     FALSE
AB FALSE FALSE   FALSE     FALSE
RF FALSE FALSE   FALSE     FALSE

I dati della tabella mydata trasformati nei corrispettivi valori logici con la funzione is.na() possono allora essere impiegati come argomento della funzione colSums() che effettua il conteggio dei dati mancanti per ciascuna colonna/variabile

colSums(is.na(mydata)) # conteggia i dati mancanti per colonna/variabile
    sesso      anni   peso_kg altezza_m 
        0         0         0         1 

A questo punto abbiamo individuato in altezza_m la variabile con i dati mancanti e abbiamo conteggiato il loro numero

Non ci resta quindi che impiegare i dati della variabile mydata$altezza_m trasformati nei corrispettivi valori logici con is.na() come argomento della funzione which() per identificare la posizione dei  dati mancanti

> which(is.na(mydata$altezza_m)) # identifica la posizione dei dati mancanti
[1] 2

che nel nostro caso risultano essere il dato della riga numero 2 della variabile mydata$altezza_m.

Per l'effetto che i dati mancanti possono determinare sui calcoli copiate e incollate nella Console di R questo script e premete ↵ Invio:

# ESEMPIO DI MEDIA ARITMETICA NON COMPUTABILE A CAUSA DI DATI MANCANTI 
#
# importa i dati, notare / invece di \ su windows
mydata <- read.table("C:/Rdati/importa_csv.csv", header=TRUE, sep=";", dec=",", row.names="id")
#
mydata # mostra i dati importati
#
mean(mydata$peso_kg) # calcola la media aritmetica dei valori di peso
mean(mydata$altezza_m) # calcola la media aritmetica dei valori di altezza
#

Il risultato che compare nella Console di R è il seguente:

> mydata # mostra i dati importati 
   sesso anni peso_kg altezza_m
MT     M   69      76      1.78
GF     F   56      63        NA
MC     F   53      71      1.60
SB     M   28      73      1.78
FE     F   61      54      1.54
AB     M   46      92      1.84
RF     F   31      81      1.56
> #
> mean(mydata$peso_kg) # calcola la media aritmetica dei valori di peso 
[1] 72.85714
> mean(mydata$altezza_m) # calcola la media aritmetica dei valori di altezza 
[1] NA

Come potete vedere il calcolo della media effettuato con la funzione mean() ha avuto successo solamente per il peso. Per l'altezza in luogo del valore della media è stata riportata la sigla NA, evidente conseguenza del valore mancante dell'altezza nel caso GF.

Copiate e incollate nella Console di R queste due ulteriori righe di codice e premete ↵ Invio:

# la corretta gestione dei dati mancanti consente di calcolare la media dei valori di altezza 
#
mean(mydata$peso_kg) # ricalcola la media aritmetica dei valori di peso
mean(mydata$altezza_m, na.rm=TRUE) # ricalcola la media aritmetica dei valori di altezza
#

Il risultato dei calcoli della media che compare nella Console di R ora è questo:

> mean(mydata$peso_kg) # ricalcola la media aritmetica dei valori di peso
[1] 72.85714
> mean(mydata$altezza_m, na.rm=TRUE) # ricalcola la media aritmetica dei valori di altezza
[1] 1.683333

Quindi con l'aggiunta nella funzione mean() dell'argomento na.rm=TRUE, che rimuove dal calcolo i dati mancanti (in questo caso uno solo, ma potrebbero essere anche molti), è stato reso possibile il calcolo della media anche per l'altezza.

Con la funzione na.omit() è infine possibile eliminare definitivamente i casi con dati mancanti: anche un solo dato mancante/campo vuoto comporta la cancellazione dell'intera riga/caso. Copiate e incollate nella Console di R questo script e premete ↵ Invio:

# ELIMINA I CASI CON DATI MANCANTI 
#
# importa i dati, notare / invece di \ su windows
mydata <- read.table("C:/Rdati/importa_csv.csv", header=TRUE, sep=";", dec=",", row.names="id")
#
mydata # mostra i dati importati
#
newdata <- na.omit(mydata) # elimina da mydata i casi con dati mancanti
newdata # mostra i dati dopo eliminazione dei casi con dati mancanti
#
mean(newdata$peso_kg) # calcola la media aritmetica dei valori di peso
mean(newdata$altezza_m) # calcola la media aritmetica dei valori di altezza
#

Dai dati importati nella tabella mydata mediante la funzione na.omit() viene generata una nuova tabella denominata newdata dalla quale sono esclusi i casi con dati mancanti. Con l'eliminazione dalla tabella del caso GF i calcoli vanno subito a buon fine sia per il peso sia per l'altezza.

> mydata # mostra i dati importati 
   sesso anni peso_kg altezza_m
MT     M   69      76      1.78
GF     F   56      63        NA
MC     F   53      71      1.60
SB     M   28      73      1.78
FE     F   61      54      1.54
AB     M   46      92      1.84
RF     F   31      81      1.56
> #
> newdata <- na.omit(mydata) # elimina da mydata i casi con dati mancanti
> newdata # mostra i dati dopo eliminazione dei casi con dati mancanti
   sesso anni peso_kg altezza_m
MT     M   69      76      1.78
MC     F   53      71      1.60
SB     M   28      73      1.78
FE     F   61      54      1.54
AB     M   46      92      1.84
RF     F   31      81      1.56
> #  
> mean(newdata$peso_kg) # calcola la media aritmetica dei valori di peso 
[1] 74.5
> mean(newdata$altezza_m) # calcola la media aritmetica dei valori di altezza 
[1] 1.683333

Da notare che impiegando la funzione na.omit() la media dei valori di peso, che con lo script precedente era 72.85714 kg, ora è cambiata ed è diventata 74.5 kg in quanto in seguito all'eliminazione dell'intero caso GF è stato eliminato anche il suo valore di peso 63 kg che contribuiva a determinare la media.


----------

[1] Vedere il post Ordinamento dei dati.

[2] Parliamo di array o vettore nel caso di dati numerici monodimensionali, disposti su una sola riga, 

8
6
11
7

di matrice nel caso di dati numerici disposti su più righe e più colonne

8
9
15
14
6
7
18
12
11
8
17
13
7
4
19
17

e di tabella nei casi in cui il contenuto, disposto su più righe e più colonne, è rappresentato oltre che da dati numerici, anche da testo e/o operatori logici

M
7
9
VERO
F
3
12
VERO
F
5
10
FALSO