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\
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:
----------
[1] Vedere Using R for Multivariate Analysis. URL consultato il 02/01/2023.
[2] Vedere Principal component analysis. URL consultato il 02/01/2023.
[3] Vedere Principal Component Analysis in R. URL consultato il 02/01/2023.
[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().
Nessun commento:
Posta un commento