Analisi delle distribuzioni dei dati nei gruppi/cluster

Se trovate utile l'impiego delle tecniche di clustering (o analisi dei gruppi) [1] come metodo per l'analisi esplorativa dei dati – ma anche se le affrontate per la prima volta – potrebbe interessarvi questo esempio.

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

# ANALIZZARE LE DISTRIBUZIONI DEI DATI NEI CLUSTER con ggpairs  

#  

library(cluster) # carica il pacchetto per il clustering

library(factoextra) # carica un pacchetto addizionale per il clustering

library(GGally) # carica il pacchetto che estende le funzioni di ggplot2

#

# preparazione dei dati

#

mydata <- read.table("c:/Rdati/ema.csv", header=TRUE, sep=";", dec=".") # importa i dati dei 643 casi

newdata <- subset(mydata, select=c("gr", "hb", "mcv", "hb_a2", "ferro")) # seleziona le variabili per il clustering

#

# clustering non gerarchico con il metodo di Rousseew (k-medoids)

#

myclust <- pam(newdata, 3, metric=c("euclidean"), stand=TRUE) # genera un oggetto 'pam' che contiene i risultati del clustering

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

#

Vediamo di cosa si tratta. Abbiamo nel file ema.csv una raccolta di 643 casi dei quali abbiamo raccolto i risultati di tre analisi del sangue richieste a scopo diagnostico:

→ esame emocromocitometrico – del quale ci interessano in particolare i valori di: (i) gr la concentrazione degli eritrociti nel sangue (globuli rossi) espressa in migliaia di miliardi per litro di sangue 10¹²/L o in milioni per microlitro di sangue 10⁶/µL (le due espressioni sono numericamente equivalenti); (ii) hb la concentrazione dell'emoglobina espressa in grammi per decilitro di sangue g/dL; (iii) mcv il volume medio dei globuli rossi espresso in milionesimi di miliardesimo di litro (10-¹⁵) o femtolitri (fL);

→ emoglobina A2 – hb_A2 la variante fisiologica dell'emoglobina presente nell'adulto espressa come percentuale dell'emoglobina totale;

→ ferro – ferro la concentrazione del ferro espressa in microgrammi per decilitro di siero µg/dL. 

Sappiamo che i 643 casi comprendono:

→ soggetti sani, impiegati come controlli, con i risultati delle tre analisi all'interno dei valori fisiologici (valori "normali");

→ soggetti con alterazioni dell'esame emocromocitometrico e con carenza di ferro (confermata dalla diminuita concentrazione del ferro nel siero);

→  soggetti con alterazioni dell'esame emocromocitometrico e portatori del gene della beta-talassemia (confermato dall'aumento della percentuale di emoglobina A2);

→  soggetti con alterazioni dell'esame emocromocitometrico ma non classificabili ulteriormente con i dati disponibili.

Siamo interessati a una analisi esplorativa dei dati e, seguendo lo script, procediamo come segue.

1) Preparazione dei dati. 

Installiamo i tre pacchetti che ci servono, cluster, factoextra e GGally, quindi ci procuriamo i dati impiegati nell'esempio (per il file di dati trovate link e istruzioni alla pagina di Download):

→ effettuando il download del file di dati ema.csv

→ salvando il file nella cartella C:\Rdati\

Con la funzione read.table() importiamo i dati nell'oggetto mydata. Quindi con la funzione subset() riportiamo nell'oggetto newdata le sole variabili che ci interessano, cioè i risultati delle analisi indicate sopra (643 record), in cinque colonne:

> newdata

      gr   hb  mcv hb_a2 ferro

1   4.90 13.3 82.8   1.8   106

2   4.66 10.8 73.6   2.6   148

3   5.43 11.5 66.5   4.8   104

4   5.41 10.8 73.4   2.5    74

..........

642 5.28 11.0 65.4   1.9    11

643 4.03  7.2 55.0   1.9    21

2) Clustering non gerarchico (analisi dei gruppi). 

A questo punto eseguiamo il clustering non gerarchico, lo facciamo impiegando il metodo di Rousseew. 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 newdata da impiegare;

→ il numero di cluster, in questo caso 3, in cui classificare i dati;

→ la metrica da impiegare che è "euclidean" (in alternativa può essere "manhattan");

→ stand=TRUE che impone la standardizzazione dei dati, che viene eseguita automaticamente dalla funzione.

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

Date le premesse riportate all'inizio abbiamo impostato 3 come numero dei gruppi (cluster) in cui classificare i dati allo scopo di separare (i) sani, (ii) carenza di ferro, (iii) beta-talassemia eterozigote. Consapevoli del fatto che data la casistica, la quantità e il tipo delle informazioni di cui disponiamo, le altre possibili e meno frequenti cause di alterazione dell'esame emocromocitometrico (e ve ne sono molte) rimarranno inestricabilmente sovrapposte tra loro e/o con qualcuno di questi tre gruppi (a causa dei vincoli imposti dalla informazione disponibile che è, come sempre, limitata).

3) Analisi delle distribuzioni dei dati nei cluster.

Eccoci arrivati al punto. Non ci basta avere, per ciascun cluster, l'indicazione dei soggetti (numerati in ordine progressivo in mydata) appartenenti a ciascun cluster. Vogliamo avere anche l'indicazione della corrispondenza tra i cluster e le distribuzioni dei risultati delle grandezze misurate. Lo facciamo in una nuova finestra grafica aperta con windows() associando a ciascun caso (record) della tabella newdata con la funzione cbind() il cluster di appartenenza myclust$clustering estratto dall'oggetto 'pam' myclust che contiene i risultati del clustering e che viene riportato aggiungendo una sesta colonna alla tabella.

#

# analisi delle distribuzioni dei dati nei cluster

#

windows() # apre e inizializza una nuova finestra grafica

newclust <- cbind(newdata, myclust$clustering) # riporta per ciascun caso il cluster di appartenenza

ggpairs(newclust, columns=1:5, aes(color=as.factor(myclust$clustering), alpha=0.5), title="Scatterplot e kernel density plot dei dati dei cluster") # scatterplot e kernel density plot dei dati nei cluster

# 

Ora abbiamo la tabella che possiamo impiegare per generare il grafico della distribuzione dei dati dei cluster con la funzione ggpairs() che impiega gli argomenti:

→ newclust per indicare la tabella che contiene i dati;

→ 1:5 per indicare che i dati da rappresentare sono contenuti nelle colonne da 1 a 5 (la sesta contiene il cluster di appartenenza);

→ as.factor() in quanto le variabili numeriche 1, 2 e 3 della colonna 6 che indicano il cluster di appartenenza devono essere trasformate in variabili fattore;

→ alpha=0.5 per rappresentare i colori con una trasparenza del 50%.

Da notare che al bisogno potete anche esportare la tabella newclust, che nelle prime cinque colonne contiene i dati e nella sesta il cluster di appartenenza come qui indicato 

> newclust

      gr   hb  mcv hb_a2 ferro myclust$clustering

1   4.90 13.3 82.8   1.8   106                  1

2   4.66 10.8 73.6   2.6   148                  1

3   5.43 11.5 66.5   4.8   104                  2

4   5.41 10.8 73.4   2.5    74                  3

..........

642 5.28 11.0 65.4   1.9    11                  3

643 4.03  7.2 55.0   1.9    21                  3

e salvarla sotto forma di file [2].

Ed ecco il grafico risultante.

La conclusione? 

1) l'analisi dei gruppi fornisce una sintesi nella quale le grandezze soggiacenti non sono più riconoscibili singolarmente in quanto le componenti principali sono una combinazione lineare di tutte le variabili originali (nel primo dei grafici la prima componente principale in ascisse è la variabile formata dalla combinazione lineare di tutte le variabili originali che spiega la maggior quantità di varianza, la seconda componente principale in ordinate è la variabile che spiega la maggior quantità di varianza in ciò che rimane una volta rimosso l'effetto della prima componente principale);

2) la distribuzione dei dati all'interno dei cluster fornisce una rappresentazione di dettaglio che consente di ricollegare questi ultimi con la distribuzione originaria delle variabili la cui combinazione ha dato loro origine, e consente di riconoscere come:

→ il cluster 1 (rosso) contiene i soggetti con emoglobina A2 e ferro normali;

→ il cluster 2 (verde) contiene i soggetti con emoglobina A2 aumentata; 

→ il cluster 3 (azzurro) contiene i soggetti con ferro diminuito.

Pur con distribuzioni che in parte si sovrappongono – un fatto determinato dalla incompletezza della informazione disponibile rispetto alla complessità dei casi – l'esempio illustra come:

→ il clustering fornisce una vista sintetica dell'informazione globale contenuta nei dati;

→ possiamo aggiungere all'analisi dei cluster una vista di dettaglio delle distribuzioni dei singoli dati all'interno dei cluster;

→ le due viste degli stessi dati si integrano sinergicamente l'una con l'altra. 

----------

[1] Fate click su clustering nelle Parole chiave.

[2] Vedere il post Esportazione dei dati in un file .csv.

Analisi delle componenti principali

L'analisi delle componenti principali (Principal Component Analysis da cui l'acronimo PCA), insieme all'analisi dei gruppi (cluster analysis o clustering), fa parte delle tecniche di analisi esplorativa dei dati e viene impiegata nell'analisi di dati multivariati [1]. 

Dato un campione descritto da n variabili, le sue componenti principali  [2, 3]:

→ sono ciascuna una combinazione lineare delle n variabili;

→ sono nello stesso numero n delle variabili, essendo la prima componente principale la variabile formata dalla combinazione lineare delle variabili originali che spiega la maggior quantità di varianza, la seconda componente principale la variabile che spiega la maggior quantità di varianza in ciò che rimane una volta rimosso l'effetto della prima componente principale, e così via, fino a spiegare tutta la varianza osservata;

→ sono non correlate, come dimostrato dal fatto che i coefficienti di correlazione tra le componenti principali sono uguali a 0 (zero);

→ concentrano di solito la maggior parte dell'informazione del campione in 2 componenti principali, con una perdita di informazione contenuta;

→ consentono una migliore introspezione nei dati in quanto 2 dimensioni sono facilmente rappresentabili e analizzabili.

Per un corretto utilizzo della PCA è opportuno ricordare che:

→ la tecnica prevede l'assunto che i dati abbiano una distribuzione gaussiana;

→ tra le variabili deve esistere una correlazione lineare;

→ le componenti principali sono una combinazione lineare delle variabili analizzate;

→ le componenti principali non sono invarianti rispetto alla scala, pertanto è necessario standardizzare le variabili che descrivono il campione sul quale sono calcolate le componenti principali;

→ a causa della standardizzazione di cui al punto precedente le nuove variabili calcolate non possono essere interpretate come i dati originari, dei quali sono una "sintesi" lineare  [2, 3].

Come dati impieghiamo i valori di BMI (indice di massa corporea) rilevati a livello europeo alcuni anni fa e pubblicati dall'Istat. Si tratta di una tabella che, per ciascuna delle Nazioni elencate, riporta la percentuale di soggetti sottopeso, con peso normale, sovrappeso e obesi [4].

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

Infine  è necessario scaricare dal CRAN il pacchetto aggiuntivo car e il pacchetto aggiuntivo cluster.

Copiate e incollate questo primo script nella Console di R e premete ↵ Invio:

# ANALISI DELLE COMPONENTI PRINCIPALI - ispezione dei dati

#

library(car) # carica il pacchetto necessario per la grafica

mydata <- read.table("c:/Rdati/bmi.csv", header=TRUE, sep=";", row.names="Nazione") # importa i dati

windows() # apre e inizializza una nuova finestra grafica

#

scatterplotMatrix(~sottopeso+normale+sovrappeso+obeso, regLine=list(method=lm, lty=1, lwd=2, col="red"), smooth=FALSE, diagonal=list(method="density", kernel="gaussian", adjust=1), col = "black", main="Matrice dei grafici xy", data=mydata) # traccia il grafico di dispersione xy che incrocia tutte le variabili 

#  

Le prime tre righe servono solamente a caricare il pacchetto necessario per la grafica e a importare i dati e si commentano da sole.

Per effettuare una prima ispezione dei dati impieghiamo la funzione scatterplotMatrix che genera i grafici xy che incrociano le quattro variabili numeriche della tabella, aggiunge con regLine=-... la retta di regressione e riporta sulla diagonale (diagonal=...) i kernel density plot (method="density") delle distribuzioni dei dati.

Il grafico documenta una proporzionalità inversa tra la percentuale di soggetti con peso normale e le percentuali di soggetti sovrappeso e obesi (all'aumentare della percentuale di soggetti con peso normale le ultime due diminuiscono), e una proporzionalità diretta tra soggetti con peso normale e soggetti sottopeso.

L'idea ora è di semplificare il quadro impiegando la PCA. Per effettuarne l'analisi numerica copiate e incollate questo secondo script nella Console di R e premete ↵ Invio:

# ANALISI DELLE COMPONENTI PRINCIPALI - analisi numerica

#  

mydata <- read.table("c:/Rdati/bmi.csv", header=TRUE, sep=";", row.names="Nazione") # importa i dati

#

options(scipen=999) # esprime i numeri in formato fisso

#

z <- scale(mydata) # calcola la deviata normale standardizzata dei dati

cbind(mydata, z) # mostra i dati originari e la loro deviata normale standardizzata z

#

mypca <- princomp(z, cor=TRUE) # calcola le componenti principali sulle variabili standardizzate

cbind(mydata, mypca$scores[,1:4]) # mostra i dati originari e le componenti principali 

#

round(cor(cbind(mydata, mypca$scores[,1:4])), digits=3) # calcola i coefficienti di correlazione delle  componenti principali e dei dati originari

summary(mypca) # sintetizza l'importanza delle componenti principali 

mypca$sdev^2 # si valuta il criterio di Kraiser che prevede varianza > 1

#

options(scipen=0) # ripristina la notazione scientifica

#

Dopo avere importato i dati, con la funzione options() e l'argomento scipen=999 imponiamo di rappresentare i risultati numerici in formato fisso (di default R prevede la notazione scientifica) per facilitarne la lettura [5].

Le componenti principali non sono invarianti rispetto alla scala: pertanto è necessario standardizzare le variabili che descrivono il campione sul quale sono calcolate le componenti principali. Questo viene fatto calcolando per i dati di ciascuna variabile la media e la deviazione standard, poi calcolando per ciascuno dato x la corrispondente deviata normale standardizzata z come

z = (x – media) / deviazione standard

Con la funzione scale() [6] sono calcolate automaticamente le deviate normali standardizzate z dei dati e con la funzione cbind() accanto ai dati originari (mydata) sono riportati in una tabella le rispettive deviate normali standardizzate (z).

> cbind(mydata, z) # mostra i dati originari e la loro deviata normale standardizzata z

                sottopeso normale sovrappeso obeso    sottopeso     normale  sovrappeso       obeso

Austria               2.4    49.6       33.3  14.7  0.484854650  1.50159124 -1.01440770 -0.74596370

Belgio                2.7    48.0       35.3  14.0  0.975769982  1.02705765 -0.32450826 -0.96253381

Bulgaria              2.2    43.8       39.2  14.8  0.157577761 -0.21859302  1.02079566 -0.71502512

Cipro                 3.9    47.8       33.8  14.5  2.939431313  0.96774095 -0.84193284 -0.80784088

Croazia               1.9    40.7       38.7  18.7 -0.333337572 -1.13800184  0.84832080  0.49157977

Danimarca             2.2    50.0       32.9  14.9  0.157577761  1.62022463 -1.15238759 -0.68408653

Estonia               2.2    43.9       33.5  20.4  0.157577761 -0.18893467 -0.94541776  1.01753574

Finlandia             1.2    44.1       36.4  18.3 -1.478806681 -0.12961797  0.05493644  0.36782542

Francia               3.2    49.6       31.9  15.3  1.793962204  1.50159124 -1.49733732 -0.56033218

Germania              1.8    46.1       35.2  16.9 -0.496976016  0.46354901 -0.35900323 -0.06531479

Grecia                1.9    41.3       39.4  17.3 -0.333337572 -0.96005175  1.08978561  0.05843955

Irlanda               1.9    42.3       37.0  18.7 -0.333337572 -0.66346826  0.26190627  0.49157977

Lettonia              1.7    41.8       35.2  21.3 -0.660614460 -0.81176000 -0.35900323  1.29598302

Lituania              1.9    42.5       38.3  17.3 -0.333337572 -0.60415156  0.71034091  0.05843955

Lussemburgo           2.8    49.3       32.4  15.6  1.139408427  1.41261619 -1.32486245 -0.46751642

Malta                 2.0    37.0       35.0  26.0 -0.169699127 -2.23536077 -0.42799317  2.75009660

Olanda                1.6    49.0       36.0  13.3 -0.824252904  1.32364114 -0.08304345 -1.17910392

Polonia               2.4    42.9       37.5  17.2  0.484854650 -0.48551816  0.43438113  0.02750097

Portogallo            1.8    44.6       36.9  16.6 -0.496976016  0.01867378  0.22741130 -0.15813055

Regno Unito           2.1    42.2       35.6  20.1 -0.006060683 -0.69312661 -0.22102334  0.92471998

Repubblica Ceca       1.1    42.1       37.6  19.3 -1.642445126 -0.72278496  0.46887610  0.67721129

Romania               1.3    42.9       46.4   9.4 -1.315168237 -0.48551816  3.50443367 -2.38570880

Slovacchia            2.1    43.6       38.0  16.3 -0.006060683 -0.27790972  0.60685599 -0.25094631

Slovenia              1.6    41.8       37.4  19.2 -0.824252904 -0.81176000  0.39988616  0.64627270

Spagna                2.2    45.4       35.7  16.7  0.157577761  0.25594057 -0.18652837 -0.12719197

Svezia                1.8    48.3       35.9  14.0 -0.496976016  1.11603270 -0.11753842 -0.96253381

Ungheria              2.9    41.9       34.0  21.2  1.303046871 -0.78210165 -0.77294290  1.26504444

Se invece nella Console di R digitate z 

> z 

sono mostrati i soli 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.

La funzione princomp() calcola le componenti principali sulle variabili standardizzate (z) impiegando la matrice di correlazione (cor=TRUE).

Quindi mediante la funzione cbind()  i dati originari (mydata) sono combinati con le quattro colonne contenenti i valori della prima componente principale (mypca$scores[,1]), della seconda componente principale (mypca$scores[,2]), della terza componente principale (mypca$scores[,3]) e della quarta componente principale (mypca$scores[,4]) e sono così visualizzati:

> cbind(mydata, mypca$scores[,1:4]) # mostra i dati originari e le componenti principali 

                sottopeso normale sovrappeso obeso     Comp.1      Comp.2      Comp.3         Comp.4

Austria               2.4    49.6       33.3  14.7  1.9568592 -0.25836896 -0.57968660 -0.00008584140

Belgio                2.7    48.0       35.3  14.0  1.6573003 -0.53681742  0.34757492  0.00090365951

Bulgaria              2.2    43.8       39.2  14.8 -0.3065287 -1.00537387  0.76399423  0.00304859104

Cipro                 3.9    47.8       33.8  14.5  2.8676299  0.45690930  1.71080991 -0.00033835804

Croazia               1.9    40.7       38.7  18.7 -1.4782942  0.10560338  0.51186402  0.00259366777

Danimarca             2.2    50.0       32.9  14.9  1.9030622 -0.26177070 -0.96268268 -0.00023905961

Estonia               2.2    43.9       33.5  20.4  0.0772148  1.39278775 -0.34702817 -0.00035012011

Finlandia             1.2    44.1       36.4  18.3 -1.0242159 -0.14902653 -1.16654614  0.00167352686

Francia               3.2    49.6       31.9  15.3  2.8277112  0.53285016  0.24082546 -0.00114263428

Germania              1.8    46.1       35.2  16.9  0.2246564 -0.11800412 -0.74397641  0.00094526687

Grecia                1.9    41.3       39.4  17.3 -1.3355279 -0.40371977  0.60210826 -0.01534678459

Irlanda               1.9    42.3       37.0  18.7 -0.8930939  0.30605551  0.05586142 -0.01670407396

Lettonia              1.7    41.8       35.2  21.3 -1.1340925  1.20428268 -0.52824564  0.00064173257

Lituania              1.9    42.5       38.3  17.3 -0.9256566 -0.28759995  0.28705908  0.00249320079

Lussemburgo           2.8    49.3       32.4  15.6  2.3063902  0.33950282 -0.18750396  0.01764811897

Malta                 2.0    37.0       35.0  26.0 -2.2351068  2.86564834  0.23542898  0.00007011548

Olanda                1.6    49.0       36.0  13.3  0.8396220 -1.43885040 -1.09559132 -0.01670331763

Polonia               2.4    42.9       37.5  17.2 -0.2687562  0.04968007  0.78163159  0.00190832184

Portogallo            1.8    44.6       36.9  16.6 -0.3104836 -0.39590962 -0.28827071 -0.01655924560

Regno Unito           2.1    42.2       35.6  20.1 -0.6497827  1.00578052  0.06207928  0.00080150891

Repubblica Ceca       1.1    42.1       37.6  19.3 -1.7935828 -0.03541460 -0.91094768  0.02066918164

Romania               1.3    42.9       46.4   9.4 -1.8998797 -4.01115995  1.00294656  0.00752978219

Slovacchia            2.1    43.6       38.0  16.3 -0.3885290 -0.45558544  0.41177524  0.00234357074

Slovenia              1.6    41.8       37.4  19.2 -1.3688625  0.24227630 -0.23901380  0.00198400366

Spagna                2.2    45.4       35.7  16.7  0.3802247 -0.01464965 -0.04421772  0.00108011741

Svezia                1.8    48.3       35.9  14.0  0.8261779 -1.10531918 -0.78871246  0.00153297266

Ungheria              2.9    41.9       34.0  21.2  0.1455442  1.97619331  0.86846432 -0.00039790369

Alla riga di codice successiva con la funzione cor() sono calcolati i coefficienti di correlazione tra tutte le variabili, arrotondando i risultati a tre decimali con la funzione round( .... , digits=3):

> round(cor(cbind(mydata, mypca$scores[,1:4])), digits=3) # calcola i coefficienti di correlazione delle  componenti principali e dei dati originari

           sottopeso normale sovrappeso  obeso Comp.1 Comp.2 Comp.3 Comp.4

sottopeso      1.000   0.412     -0.567 -0.108  0.755  0.340  0.560  0.001

normale        0.412   1.000     -0.481 -0.688  0.902 -0.327 -0.283  0.005

sovrappeso    -0.567  -0.481      1.000 -0.290 -0.679 -0.666  0.309  0.004

obeso         -0.108  -0.688     -0.290  1.000 -0.473  0.877 -0.088  0.005

Comp.1         0.755   0.902     -0.679 -0.473  1.000  0.000  0.000  0.000

Comp.2         0.340  -0.327     -0.666  0.877  0.000  1.000  0.000  0.000

Comp.3         0.560  -0.283      0.309 -0.088  0.000  0.000  1.000  0.000

Comp.4         0.001   0.005      0.004  0.005  0.000  0.000  0.000  1.000

L'analisi numerica si conclude con la sintesi fornita dalla funzione summary():

> summary(mypca) # sintetizza dell'importanza delle componenti principali 

Importance of components:

                          Comp.1    Comp.2    Comp.3        Comp.4

Standard deviation     1.4380589 1.1980885 0.7046275 0.00840918288

Proportion of Variance 0.5170033 0.3588540 0.1241250 0.00001767859

Cumulative Proportion  0.5170033 0.8758574 0.9999823 1.00000000000

che fornisce una indicazione importante. 

Con la quota di varianza (Proportion of Variance) fornita da ciascuna delle quattro componenti principali vediamo che la prima componente principale (Comp.1) spiega il 51.7% della variabilità del campione, la seconda componente principale (Comp.2) spiega il 35.9% della variabilità del campione, la terza componente principale (Comp.3) spiega il 12.4% della variabilità del campione, che sommati danno il 100%, essendo trascurabile la variabilità spiegata dalla quarta componente principale (Comp.4).

La quota cumulativa di varianza (Cumulative Proportion) ci conferma che nella prima e nella seconda componente principale è contenuta la maggior parte della variabilità del campione (87.6%). Se siamo disposti a sacrificare una piccola quota dell'informazione (varianza) contenuta nel campione (il 12.4%), una analisi del campione può essere effettuata impiegando due sole variabili, la prima componente principale (Comp.1) e la seconda componente principale (Comp.2), invece delle quattro variabili originali (sottopeso, normale, sovrappeso, obeso), consentendo quindi una migliore introspezione nei dati.

La penultima riga di codice (mypca$sdev^2) riporta la varianza spiegata dalle componenti principali, per consentire di applicare il criterio di Kraiser, che prevede che siano significative le componenti principali con una varianza maggiore di 1. Che si confermano essere le prime due:

> mypca$sdev ^ 2 # si valuta il criterio di Kraiser che prevede varianza > 1

       Comp.1        Comp.2        Comp.3        Comp.4 

2.06801325718 1.43541616920 0.49649985926 0.00007071436 

Infine con options(scipen=0) viene ripristinata la notazione scientifica.

Dato che quindi possiamo considerare adeguata un'analisi del campione che impiega solamente la prima e la seconda componente principale, passiamo all'analisi grafica rappresentando le due componenti in un grafico cartesiano xy con la funzione biplot().

Copiate e incollate questo terzo script nella Console di R e premete ↵ Invio:

# ANALISI DELLE COMPONENTI PRINCIPALI - analisi grafica

#  

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

#

mypca <- princomp(z, cor=TRUE) # calcola le componenti principali sulle variabili standardizzate

#

biplot(mypca, xlim=c(-0.4,0.4), xlab="Componente principale 1", ylab="Componente principale 2", cex=0.6) # traccia il grafico con i vettori delle componenti principali

# 

Da notare che nel biplot:

→ l'asse inferiore riporta la scala nella quale sono espressi i valori della componente principale 1

→ l'asse di sinistra riporta la scala nella quale sono espressi i valori della componente principale 2

→ l'asse superiore riporta la scala nella quale sono espressi i pesi dei vettori sulla componente principale 1

→ l'asse di destra riporta la scala nella quale sono espressi i pesi dei vettori sulla componente principale 2

→ la proiezione di un vettore sull'asse di una specifica componente principale fornisce il peso che la variabile rappresentata nel vettore ha su quella componente principale;

→ vettori separati tra loro da angoli piccoli indicano variabili con correlazione positiva;

→ vettori che divergono molto tra loro (fino a 180 gradi) indicano variabili con correlazione negativa.

Dall'ispezione del biplot si ricava che:

→ alcune Nazioni (Svezia, Olanda, Belgio, Danimarca, Austria, Lussemburgo, Francia, Cipro) si proiettano prevalentemente nella dimensione "normale" e in alcuni casi (Francia e Cipro) con una possibile rilevanza della dimensione "sottopeso";

→ altre Nazioni (Regno Unito, Lettonia, Estonia, Ungheria, Malta) si proiettano prevalentemente nella dimensione "obeso", con un dato particolarmente rilevante che riguarda Malta;

→ le Nazioni rimanenti si proiettano nella dimensione del "sovrappeso", con un dato particolarmente rilevante che riguarda la Romania.

Sulla base dei rilievi forniti dall'analisi delle componenti principali è possibile ricondurre le singole osservazioni ai dati originari, mediante i quali, per esempio, si può confermare che Malta, con il 26% di obesi, e la Romania, con il 46.4% di sovrappeso e il 9.4% di obesi, si distaccano in modo importante dalle altre Nazioni.

La cosa molto interessante è che l'analisi delle componenti principali fornisce il collegamento con le tecniche complementari, e precisamente con l'analisi dei gruppi (clustering). Lo vediamo applicando ai dati di BMI il clustering non gerarchico con il metodo di MacQueen e l'algoritmo k-means.

Copiate e incollate questo quarto e ultimo script 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 

#

Dopo i soliti preliminari, con la funzione kmeans() sono generati 4 cluster impiegando l'algoritmo "MacQueen".

Quindi viene impiegata per tracciare il grafico la funzione clusplot() con i seguenti argomenti:

→ l'oggetto z contenente i dati (standardizzati);

→ l'oggetto myclust$cluster contenente i cluster;

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

La distribuzione delle Nazioni che si vede in questo grafico si sovrappone perfettamente a quella realizzata con l'analisi delle componenti principali e in aggiunta identifica quattro cluster ben separati l'uno dall'altro e ben individuati.

In un caso ancora abbastanza semplice come questo possiamo anche pensare di validare i risultati fin qui ottenuti rianalizzando i dati originali: per questo li ordiniamo in senso decrescente prima per la prevalenza di peso normale, poi per la prevalenza di obesi, infine per la prevalenza di sovrappeso, con questo risultato

che, considerata la Romania con i suoi valori estremi come un caso a sé stante, conferma l'esistenza nei dati di sottoinsiemi numericamente ben caratterizzati, corrispondenti a quelli rilevati con la cluster analysis.

In conclusione nel caso di dati multivariati l'analisi delle componenti principali (PCA), pur con gli aspetti un po' delicati ricordati all'inizio:

→ consente di ridurre il numero delle variabili;

→ consente di semplificare la rappresentazione grafica dei dati; 

→ comporta una perdita di informazione tuttavia controllata e misurabile;

→ fornisce le basi per realizzare il clustering.

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

→ Analisi dei gruppi (clustering gerarchico) 

→ Analisi dei gruppi (clustering non gerarchico) 

→ Analisi dei gruppi (clustering non esclusivo) 

----------

[1] Vedere: Using R for Multivariate Analysis.

https://little-book-of-r-for-multivariate-analysis.readthedocs.io/en/latest/src/multivariateanalysis.html

 
[2] Vedere: Principal component analysis.

https://en.wikipedia.org/wiki/Principal_component_analysis

[3] Vedere: Principal Component Analysis in R.

https://www.datacamp.com/tutorial/pca-analysis-r

[4] I dati sono illustrati nel post Indice di massa corporea (BMI). 

[5] Vedere il post Esprimere in numeri in formato fisso.

[6] Digitate help(scale) nella Console di R per la documentazione della funzione scale().