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().

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 quindi impiegati 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:

→ Analisi delle componenti principali 

→ Analisi dei gruppi (clustering gerarchico) 

→ Analisi dei gruppi (clustering non esclusivo) 

----------

[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’.

https://cran.r-project.org/web/packages/cluster/cluster.pdf

[4] Trovate la documentazione nel manuale di riferimento del pacchetto: Package ‘factoextra’.

https://cran.r-project.org/web/packages/factoextra/factoextra.pdf

[5] Trovate la documentazione nel manuale di riferimento del pacchetto: Package ‘ggplot2’.

https://cran.r-project.org/web/packages/ggplot2/ggplot2.pdf