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

Regressione lineare multipla

La regressione lineare semplice (regressione bivariata) costruita sul piano cartesiano

y = a + b·x

può essere estesa a uno spazio n-dimensionale: abbiamo così la regressione lineare multipla nella quale la variabile dipendente y viene fatta dipendere non più da una, bensì da k variabili indipendenti, indicate come x₁, x₂, x₃ … xₖ secondo l'equazione

y = a + b₁·x₁ + b₂·x₂ + b₃·x₃ + ... + bₖ·xₖ

Il risultato è rappresentato quindi da una intercetta a e da tanti coefficienti angolari b (indicati con b₁, b₂, b₃ ... bₖ) quante sono le variabili indipendenti (x₁, x₂, x₃ ... xₖ) che contribuiscono a determinare il valore della variabile dipendente y.

Come esempio impieghiamo il set di dati ais nel quale, tra i dati ematologici rilevati in un gruppo di atleti australiani, abbiamo tre grandezze:

→ la concentrazione degli eritrociti (globuli rossi) espressa in milioni per microlitro di sangue (10⁶/µL) e contenuta nella variabile rcc;

→ il valore ematòcrito cioè la frazione del volume del sangue occupata dagli eritrociti espressa in percentuale (%) e contenuta nella variabile hc;

→ la concentrazione dell'emoglobina espressa in grammi per decilitro di sangue (g/dL) e contenuta nella variabile hg.

Le tre grandezze sono strettamente legate tra loro dal punto di vista biologico, in quanto gli eritrociti hanno ciascuno un certo volume e contengono ciascuno una data quantità di emoglobina: pertanto all'aumentare del numero di eritrociti aumenta la frazione percentuale del volume del sangue occupata dagli eritrociti (il valore ematòcrito) e aumenta la quantità di emoglobina. Mentre la diminuzione dell'ematòcrito e la diminuzione dell'emoglobina sono indice della diminuzione del numero degli eritrociti.

L'esempio è stato scomposto in una serie di script indipendenti - onde facilitare il riutilizzo del codice per analizzare i propri dati - che impiegano i pacchetti aggiuntivi gvlma, car e rgl, che devono essere scaricati dal CRAN. Accertatevi di avere installato inoltre il pacchetto DAAG che contiene i dati impiegati nello script, o in alternativa  procedete come indicato nel post Il set di dati ais. 

Iniziamo con una analisi esplorativa preliminare dei dati, copiate questo primo script, incollatelo nella Console di R e premete ↵ Invio:

# REGRESSIONE LINEARE MULTIPLA - analisi esplorativa preliminare dei dati

#

library(DAAG) # carica il pacchetto che include il set di dati ais

library(car) # carica il pacchetto per la grafica

mydata <- ais[,c(1,3,4)] # le colonne 1, 3, 4 contengono le variabili rcc, hc, hg

windows() # apre e inizializza una nuova finestra grafica

#

scatterplotMatrix(~rcc+hc+hg, col="black", pch=20, regLine=list(method=lm, lty=1, lwd=2, col="chartreuse3"), smooth=FALSE, diagonal=list(method="histogram", breaks="FD"), main="Istogrammi e grafici xy con rette di regressione", data=mydata) # istogrammi e grafici xy con rette di regressione

#

Dopo avere caricato i pacchetti necessari, copiato nell'oggetto mydata le tre variabili che ci interessano, aperto e inizializzato una nuova finestra grafica, l'analisi preliminare dei dati viene effettuata con la funzione scatterplotMatrix() nella quale:

→ ~rcc+hc+hg specifica le variabili da rappresentare, che sono la concentrazione degli eritrociti rcc, il valore ematòcrito hc e la concentrazone di emoglobia hg;

→ col="..." specifica il colore dei grafici; 

→ pch=... specifica il tipo di simbolo da impiegare (un piccolo cerchio pieno);

→ per la retta di regressione regLine sono impiegati il modello lineare (lm), una linea continua (lty=1) con uno spessore doppio (lwd=2) e un colore (col="chartreuse3"), mentre l'argomento smooth=FALSE fa si che non venga aggiunta ai grafici la rappresentazione della funzione non lineare prevista di default nel pacchetto car;

→ l'argomento diagonal=list() indica la rappresentazione della variabile riportata nella diagonale, che qui è l'istogramma (method="histogram") con il numero di classi (breaks="FD") calcolato mediante la regola di Freedman-Diaconis. 

Da notare che le espressioni che è possibile impiegare per l'argomento  diagonal=... sono:

→ list(method="boxplot");

→ list(method="density", bw="nrd0", adjust=1, kernel="gaussian", na.rm=TRUE);

→ list(method="histogram");

→ list(method="oned");

→ list(method="qqplot").

Vista la relazione lineare molto buona che intercorre tra le variabili sembra ragionevole continuare ed effettuare il calcolo della regressione lineare multipla, copiate questo secondo script, incollatelo nella Console di R e premete ↵ Invio:

# REGRESSIONE LINEARE MULTIPLA - calcolo della regressione lineare multipla

#  

library(DAAG) # carica il pacchetto che include il set di dati ais 

mydata <- ais[,c(1,3,4)] # le colonne 1, 3, 4 contengono le variabili rcc, hc, hg

reglin <- lm(rcc~hg+hc, data=mydata) # calcolo della regressione 

#

coefficients(reglin) # intercetta e coefficienti angolari 

confint(reglin, level=0.95) # intervalli di confidenza al 95% 

#  

La funzione lm(...) viene impiegata per calcolare la regressione lineare multipla specificando con l'espressione rcc~hg+hc la variabile dipendente (rcc), le due variabili indipendenti (hg e hc) e l'equazione 

rcc = a + b₁·hg + b₂·hc

quindi possiamo ricavare con la funzione coefficients() i valori dell'intercetta (a), del coefficiente angolare di hg (b₁) e del coefficiente angolare di hc (b₂):

(Intercept)          hg          hc 

-0.24269563  0.03283747  0.10403395 

e riportare l'equazione della regressione lineare multipla come

rcc = -0.24269563 + 0.03283747·hg + 0.10403395·hc

Successivamente con la funzione confint() sono riportati gli intervalli di confidenza al 95% dell'intercetta e dei due coefficienti angolari:

                  2.5 %     97.5 %

(Intercept) -0.53163317 0.04624192

hg          -0.02459817 0.09027311

hc           0.08267073 0.12539718

Da notare che gli intervalli di confidenza forniscono un test di significatività per intercetta e coefficienti angolari:

→ dato un intervallo di confidenza al 95% che va da -0.53163317 a 0.04624192 (che quindi include il valore 0) l'intercetta non è significativamente diversa da 0, cosa molto logica in quanto al diminuire fino ad azzerarsi dei valori di emoglobina ed ematòcrito deve necessariamente diminuire e azzerarsi il numero degli eritrociti;

→ con un intervallo di confidenza al 95% che va da -0.02459817 a 0.09027311 (che quindi include il valore 0) il coefficiente angolare della concentrazione di emoglobina hg non è significativamente diversa da 0, cosa meno logica in quanto non solo dal punto di vista biologico deve esistere una proporzionalità tra numero degli eritrociti e concentrazione dell'emoglobina, ma anche nell'analisi esplorativa preliminare dei dati riportata con lo script precedente risulta evidente la proporzionalità tra numero degli eritrociti e concentrazione dell'emoglobina che è una dato di fatto dal punto di vista biologico;

→ con un intervallo di confidenza al 95% che va da 0.08267073 a 0.12539718 (che quindi non include il valore 0) il coefficiente angolare del valore ematòcrito hc risulta significativamente diverso da 0.

Come detto all'inizio, oltre ad effettuare una analisi esplorativa preliminare dei dati è opportuno sottoporre i risultati ad una valutazione critica a posteriori, che effettuiamo ora sotto forma di una diagnostica numerica della regressione lineare multipla, copiate il terzo script, incollatelo nella Console di R e premete ↵ Invio:

# REGRESSIONE LINEARE MULTIPLA - diagnostica numerica della regressione

#

library(gvlma) # carica il pacchetto per la diagnostica numerica 

library(car) # carica il pacchetto per la grafica 

library(DAAG) # carica il pacchetto che include il set di dati ais

mydata <- ais[,c(1,3,4)] # le colonne 1, 3, 4 contengono le variabili rcc, hc, hg

reglin <- lm(rcc~hg+hc, data=mydata) # calcolo della retta di regressione lineare multipla

#

summary(gvlma(reglin)) #  test globale per l'assunto di linearità 

outlierTest(reglin) # valore p di Bonferroni per la presenza di dati aberranti (outliers) 

#

Ecco la prima parte del riepilogo dei risultati fornito dalla funzione summary():

Call:

lm(formula = rcc ~ hg + hc, data = mydata)

Residuals:

     Min       1Q   Median       3Q      Max 

-0.55791 -0.11284 -0.00699  0.09950  0.53595 

Coefficients:

            Estimate Std. Error t value Pr(>|t|)    

(Intercept) -0.24270    0.14652  -1.656   0.0992 .  

hg           0.03284    0.02913   1.127   0.2609    

hc           0.10403    0.01083   9.603   <2e-16 ***

---

Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.1744 on 199 degrees of freedom

Multiple R-squared:  0.8565,    Adjusted R-squared:  0.855 

F-statistic: 593.8 on 2 and 199 DF,  p-value: < 2.2e-16

Il valore di probabilità Pr(>|t|) ottenuto con il test t di Student conferma (ovviamente) i risultati della significatività ricavati dagli intervalli di confidenza riportati sopra, per i quali valgono le considerazioni già fatte. E se la significatività statistica non è coerente con una realtà biologica assodata, dovrà ovviamente essere quest'ultima a prevalere: per cui assumiamo che l'equazione calcolata sia, fino a prova contraria, una descrizione sufficientemente adeguata della relazione di funzione che intercorre tra le tre variabili esaminate.

Che questo sia vero ce lo conferma la funzione gvlma(), acronimo di global validation of linear models assumptions, che fornisce questi risultati:

→ i dati sono distribuiti in modo lineare (Global Stat), assunto accettabile;

→ i dati non presentano asimmetria (Skewness), assunto accettabile;

→ i dati non presentano curtosi (Skewness), assunto accettabile;

→ nell'ambito dei valori osservati la varianza è costante (Heteroscedasticity), assunto accettabile.

ASSESSMENT OF THE LINEAR MODEL ASSUMPTIONS

USING THE GLOBAL TEST ON 4 DEGREES-OF-FREEDOM:

Level of Significance =  0.05 

Call:

 gvlma(x = reglin) 

                    Value p-value                   Decision

Global Stat        6.5801 0.15981    Assumptions acceptable.

Skewness           0.5149 0.47304    Assumptions acceptable.

Kurtosis           0.7559 0.38461    Assumptions acceptable.

Link Function      3.9844 0.04592 Assumptions NOT satisfied!

Heteroscedasticity 1.3249 0.24972    Assumptions acceptable.

L'assenza nei dati di asimmetria e di curtosi, e la loro omoschedasticità, sono condizioni sufficienti per garantire la loro adeguata descrizione mediante una regressione lineare. Il primo test (Global Stat) di valutazione globale degli assunti di linearità, che sostanzialmente comprende tutti e tre i test precedenti in forma compatta, li conferma. In questo contesto il test rimanente e l'unico discordante (Link Function) può essere ignorato.

Con la funzione outlierTest() viene  effettuato il test di Bonferroni, che non evidenzia dati anomali (outliers) ma segnala il dato 78 come quello che si discosta maggiormente dagli altri:

No Studentized residuals with Bonferroni p < 0.05

Largest |rstudent|:

   rstudent unadjusted p-value Bonferroni p

78 -3.28855          0.0011917      0.24073

In tema di valutazione critica a posteriori, aggiungiamo ora alla diagnostica numerica anche una diagnostica grafica della regressione lineare multipla, per effettuarla copiate questo quarto script, incollatelo nella Console di R e premete ↵ Invio:

# REGRESSIONE LINEARE MULTIPLA - diagnostica grafica della regressione

#

library(car) # carica il pacchetto per la grafica

library(DAAG) # carica il pacchetto che include il set di dati ais

mydata <- ais[,c(1,3,4)] # le colonne 1, 3, 4 contengono le variabili rcc, hc, hg

reglin <- lm(rcc~hg+hc, data=mydata) # calcolo della retta di regressione lineare multipla

#

windows() # apre e inizializza una nuova finestra grafica

qqPlot(reglin, main="Grafico dei quantili per i residui", xlab="t-quantili", ylab="Residui studentizzati") # grafico globale per linearità e dati aberranti

#

windows() # apre e inizializza una nuova finestra grafica

ceresPlots(reglin, ask=FALSE, main="Grafico di Ceres", ylab="Residui di Ceres (rcc)") # grafico di Ceres per la valutazione separata della linearità di hg e hc vs. rcc

#

windows() # apre e inizializza una nuova finestra grafica

plot(mydata$rcc, reglin$fitted.values, xlim=c(4,7), ylim=c(4,7)) # grafico valori originari vs valori stimati

abline(0, 1, lty=2) # linea di indentità y = x

#

La funzione qqPlot() riporta una distribuzione dei residui lineare, evidenziando nel contempo i dati 78 e 161 come quelli che dei quali sarebbe utile controllare la validità, ovvero che sono situati in intervalli di valori per i quali potrebbe essere opportuno acquisire più dati. 

La funzione ceresPlot() genera due grafici di Ceres separati per emoglobina ed ematòcrito, e consente di  visualizzare una relazione lineare tra ematòcrito ed eritrociti migliore di quella che intercorre tra emoglobina ed eritrociti, confermando in un certo senso i valori di significatività riportati sopra.

La diagnostica grafica della regressione lineare multipla si conclude con questo grafico minimalista, ma che entra nell'essenza dei risultati.

Il grafico, generato con la funzione plot(), riporta in ascisse i valori originari della concentrazione degli eritrociti (mydata$rcc) e in ordinate i corrispondenti valori reglin$fitted.values ricalcolati da hg e hc impiegando l'equazione della regressione lineare multipla ove

reglin$fitted.values = -0.24269563 + 0.03283747·hg + 0.10403395·hc

La linea di identità y = x che viene riportata evidenzia, senza la necessità di ulteriori statistiche, che l'informazione acquisita misurando la concentrazione dell'emoglobina e l'ematòcrito può essere compressa in un unico valore, che fornisce una stima mediamente adeguata della concentrazione degli eritrociti.

Nota: con "mediamente adeguata" si intende che, data la misura diretta della concentrazione degli eritrociti riportata in ascisse, la sua misura indiretta - ottenuta misurando emoglobina ed ematòcrito ed applicando la regressione lineare multipla - comporta una perdita di informazione [aumento dell'incertezza della misura] che si riflette nella "dispersione media" dei valori sull'asse delle ordinate. Se non ci fosse perdita di informazione i punti sarebbero perfettamente allineati sulla retta y = x.

Ora, visto che con tre variabili la cosa è ancora fattibile, non resta che aggiungere un grafico tridimensionale (3D).

Copiate quest'ultimo script, incollatelo nella Console di R e premete ↵ Invio:

# GRAFICO 3D ANIMATO con il pacchetto rgl

#  

library(rgl) # carica il pacchetto per la grafica 3D

library(DAAG) # carica il pacchetto che include il set di dati ais 

mydata <- ais[,c(1,3,4)] # le colonne 1, 3, 4 contengono le variabili rcc, hc, hg

windows() # apre e inizializza una nuova finestra grafica

#  

# grafico 3D animato

#  

plot3d(mydata$rcc, y=mydata$hg, z=mydata$hc, type="p", col="red", size=3) # realizza il grafico

play3d(spin3d(axis=c(0,0,1), rpm=10), duration=10) # anima il grafico

#  

Dopo avere caricato il pacchetto necessario con il comando library(rgl) il grafico viene realizzato con la funzione plot3d().

Quindi con la funzione play3d() viene realizzata una animazione della durata di 10 secondi (duration=10) con i parametri specificati dalla funzione spin3d() e cioè l'asse sul quale fare ruotare l'animazione, nel nostro caso  l'asse delle z (0,0,1), e il numero di rotazioni al minuto (rpm=10).

Quando l'animazione che compare nella finestra "RGL device" aperta nella Console di R dallo script si ferma, se “afferrate” il grafico 3D facendo click con il tasto sinistro del mouse e tenendolo premuto senza rilasciarlo potete ruotarlo a vostro piacimento muovendo il mouse. 

Anche se non aggiunge nulla di nuovo, il grafico [1] ci conforta documentando la linearità della relazione tra le tre variabili.

Conclusione: mettere in relazione tra di loro troppe variabili contemporaneamente può portare a risultati di difficile comprensione. Tuttavia una analisi esplorativa preliminare dei dati e una attenta valutazione critica a posteriori dei risultati possono aiutare a identificare i casi nei quali l'impiego della regressione lineare multipla può essere appropriato.

----------

[1] Trovate altre modalità di rappresentazione sotto forma di grafici tridimensionali nel post Grafici 3D.