Archivio tag: acp

Quante componenti estrarre? La parallel analysis

L’analisi in componenti principali (ACP) è una celebre tecnica di analisi multivariata che, dato un gruppo di variabili quantitative misurate sulle medesime unità statistiche, consente di identificare un insieme ridotto di nuove variabili in grado di rappresentare i dati in maniera più sintetica. Date n variabili osservate, fra loro più o meno correlate, la ACP le trasforma in altrettante variabili fra loro linearmente indipendenti. Queste nuove variabili sono dette “componenti” e vengono estratte in ordine d’importanza: da quella che spiega la quota maggiore di varianza dei dati osservati a quella che ne spiega la quota minore.

Lo scopo della ACP è ridurre i dati trasformandoli, cercando di veicolare la medesima quantità d’informazione utilizzando un numero inferiore di variabili. Ne abbiamo già parlato in un precedente post, al quale si rimanda per ulteriori approfondimenti.

Quello che vorrei affrontare in questo articolo è un vecchissimo problema che ancora sembra affliggere molte ricerche, ovvero la determinazione del numero di componenti da estrarre. Una volta individuate le componenti, bisogna scegliere quante tenerne: se le tenessimo tutte non avremmo certo ridotto i dati, ma se ne tenessimo troppo poche perderemmo preziose informazioni. Si tratta di un problema non da poco, anche perché spesso è attraverso la ACP che viene definita la dimensionalità di strumenti di misura come i test psicologici utilizzati nella prassi psicodiagnostica.

Prima di addentrarci nel problema, vediamo di procurarci dei dati da utilizzare per svolgere qualche analisi. Il pacchetto Flury contiene tutti i dataset utilizzati in un volume del 1997 di B. Flury dedicato alla statistica multivariata; la libreria va installata perché non è distribuita con la versione standard di R. Noi utilizzeremo il dataset microtus, che contiene la registrazione di otto variabili biometriche su 288 roditori appartenenti appunto al genere “Microtus”.

install.pakcages("Flury") # Installa la libreria Flury
library(Flury)            # Carica la libreria Flury
data(microtus)            # Attiva il dataset

head(microtus)            # Visualizza l'anteprima del dataset
      Group M1Left M2Left M3Left Foramen Pbone Length Height Rostrum
1 multiplex   2078   1649   1708    3868  5463   2355    805     475
2 multiplex   1929   1551   1550    3825  4741   2305    760     450
3 multiplex   1888   1613   1674    4440  4807   2388    775     460
4 multiplex   2020   1670   1829    3800  4974   2370    766     460
5 multiplex   2223   1814   1933    4222  5460   2470    815     475
6 multiplex   2190   1800   2066    4662  4860   2535    838     521

La variabile Group (prima colonna) classifica i roditori in tre gruppi: “multiplex”, “subterraneus” e “unknown”. Noi ci interesseremo del gruppo “subterraneus”, per cui isoliamo i dati di questi esemplari in un nuovo dataset:


subterraneus <- microtus[microtus$Group=="subterraneus", -1] [/code]

R offre diversi comandi per realizzare una ACP. Qui utilizzeremo la funzione principal del pacchetto psych, anche questo da installare manualmente. Utilizzare la funzione principal è molto semplice: oltre alla matrice di dati osservati, dobbiamo specificare anche il numero di componenti da estrarre. Dato che ancora non abbiamo deciso quante componenti utilizzare (non abbiamo nessun aggancio teorico), ne estraiamo tante quante sono le variabili osservate.


install.packages(“psych”) # Installa la libreria psych
library(psych) # Carica la libreria psych
k <- ncol(subterraneus) # Numero di componenti (= num. variabili) PCA <- principal(subterraneus, nfactors=k, rotate="none") [/code]

Il problema che sorge a questo punto è capire quante sono effettivamente le componenti che andrebbero estratte. Come già visto, una cosa che possiamo fare è visionare gli autovalori:

> PCA$values
[1] 4.02 1.21 0.85 0.68 0.46 0.38 0.27 0.12

Gli autovalori rappresentano la varianza delle singole componenti, e il rapporto tra ogni autovalore e la somma degli autovalori di tutte le componenti fornisce la varianza spiegata da ogni componente. Del 100% della variabilità originaria, quanto è spiegato da ogni componente? Lo si può capire così:

> PCA$values / sum(PCA$values)
[1] 0.50 0.15 0.11 0.08 0.06 0.05 0.03 0.02

Una regola molto utilizzata per stabilire il numero di componenti è quella di estrarne tante quanti sono gli autovalori maggiori di 1. Nel nostro caso dovremmo tenere solo le prime due componenti, che complessivamente spiegano il 65% della varianza (0.50 + 0.15). Siamo però davvero sicuri che tenere le prime due componenti sia la soluzione migliore? L’autovalore della seconda componente vale 1.21, un valore molto vicino alla fatidica soglia 1: non è che magari vi è così vicino da dover essere escluso?

Ragionando in ottica di analisi in componenti principali potremmo anche tenere la seconda componente: avrà pure un autovalore troppo vicino a 1, ma spiega pur sempre un 15% di varianza. Molto spesso, però, l’analisi degli autovalori calcolati dalla ACP è utilizzata per stabilire il numero di componenti in un’analisi fattoriale, la quale oltre che alla percentuale di varianza spiegata deve guardare soprattutto all’economia della soluzione trovata, alla generalizzabilità del risultato e alla congruenza con assunzioni teoriche più o meno fondate.

C’è un tipo di analisi molto raffinata che aiuta nella scelta del numero di componenti da estrarre: la parallel analysis. Non si tratta di una tecnica innovativa, anzi; essa risale agli anni ’60 ed è basata su un’osservazione: talvolta può capitare che alcuni campioni di dati producano degli autovalori maggiori di 1 per puro effetto del caso. La differenza tra un autovalore di poco superiore a 1 e un altro di poco inferiore può essere molto sottile, così sottile che può essere stato il caso a determinarla.

La parallel analysis consiste nel confrontare gli autovalori calcolati sui dati osservati non solo con il valore soglia 1, ma anche con degli autovalori calcolati su matrici di dati casuali. L’idea è che una componente, perché sia considerata significativa, non solo debba presentare un autovalore maggiore di 1, ma il suo autovalore deve anche essere superiore al corrispettivo autovalore prodotto da dei dati casuali. In sostanza, la componente deve dimostrare di essere più forte del caso.

La tecnica prevede di realizzare, parallelamente all’analisi in componenti principali, delle altre analisi su matrici di dati generati casualmente. Definito un numero n di iterazioni, si producono n matrici di dati casuali dalle caratteristiche simili a quelle dei dati osservati: stesso numero di unità campionarie e stesso numero di variabili. Su ogni matrice si esegue un’analisi in componenti principali e si calcolano gli autovalori; come risultato si avrà, per ogni componente, una serie di autovalori di numerosità n.

Per esempio, nel dataset subterraneus abbiamo 8 variabili misurate su 46 unità. Possiamo definire un n = 1000 e quindi costruire 1000 matrici di dati casuali con 46 righe e 8 colonne. Su ogni matrice si esegue un’analisi in componenti principali e si calcolano gli autovalori, ottenendo così 1000 autovalori per ognuna delle 8 componenti possibili.

Al termine dell’operazione, per ognuna delle 8 componenti, sui 1000 autovalori si calcola un indice riassuntivo che può essere la media o la mediana. Avremo come risultato 8 medie (o mediane), che rappresentano gli autovalori generati dai dati casuali; questi indici rispondono alla domanda: «Quali autovalori avremmo ottenuto se i dati osservati fossero stati determinati completamente dal caso?».

A questo punto si confrontano autovalori “reali” e autovalori “casuali”: il numero di componenti da estrarre corrisponde al numero di autovalori reali che risultano più elevati degli autovalori casuali.

Realizzare una parallel analysis in R è abbastanza semplice. La libreria psych contiene l’ottima funzione fa.parallel, ma per diversi motivi (che magari approfondiremo in un altro post) io preferisco utilizzare uno script realizzato da me artigianalmente. Nel mio account GitHub ho caricato uno script R che consente di realizzare in maniera completamente automatica la procedura descritta sopra. Per scaricare lo script è sufficiente andare in basso a destra nella pagina e cliccare su “Donwload ZIP” oppure cliccare direttamente su questo link. Lo script su chiama parallelPCA.R e va estratto dall’archivio zip.

Vediamo ora di applicare la parallel analysis sul dataset microtus. Lo script deve preventivamente caricato in R; il caricamento può essere eseguito direttamente da menù: File > Sorgente codice R… per renderlo disponibile nell’area di lavoro. Quindi, utilizzando il comando parallelPCA si esegue la parallel analysis.


parpca <- parallelPCA(subterraneus) [/code]

L’output è prevalentemente di tipo grafico, e mostra uno scree-graph al quale sono sovrapposti, in grigio, i quantili delle distribuzioni degli autovalori dei dati casuali: i puntini rappresentano le mediane mentre le barre rappresentano il 2.5° e 97.5° percentile.

Gli stessi dati sono riportati nell’output testuale: autovalori reali, medie, intervalli di confidenza e quantili degli autovalori calcolati sui dati casuali.

> parpca
Parallel Analysis by Random Data 
Correlation matrix: pearson 

Observed values:
        c1   c2   c3   c4   c5   c6   c7   c8
eigen 4.02 1.21 0.85 0.68 0.46 0.38 0.27 0.12

Parallel Confidence Intervals:
         c1   c2   c3   c4   c5   c6   c7   c8
CI Sup 1.93 1.61 1.38 1.19 1.02 0.85 0.69 0.52
Mean   1.68 1.40 1.20 1.03 0.89 0.74 0.60 0.45
CI Inf 1.43 1.19 1.02 0.88 0.76 0.63 0.51 0.39

Parallel Quantiles:
        c1   c2   c3   c4   c5   c6   c7   c8
97.5% 2.00 1.59 1.34 1.16 1.01 0.87 0.73 0.59
50%   1.67 1.40 1.20 1.03 0.89 0.74 0.60 0.45
2.5%  1.45 1.23 1.06 0.91 0.76 0.61 0.47 0.32

A questo punto abbiamo un’arma in più per decidere il numero di componenti da estrarre. Dei due autovalori maggiori di 1, solo uno si dimostra effettivamente “più forte” del caso: già il secondo autovalore prodotto dai dati reali (1.21) è inferiore alla mediana degli autovalori prodotti dai dati casuali (1.40), e gli altri di conseguenza a seguire.

Quindi, quante componenti estrarre? Premesso che certe decisioni dovrebbe essere la teoria a guidarle, la parallel analysis ci indica come potenzialmente rischiosa l’assunzione di due componenti. Se sotto analisi ci fosse un questionario psicometrico e questa analisi fosse il preludio di un’analisi fattoriale, con le componenti da estrarre che rappresentano dei costrutti latenti, io sceglierei un unico fattore.

Bibliografia essenziale

Horn, J. L. (1965). A rationale and test for the number of factors in factor analysis. Psychometrika, 30(2), 179-185.

Zwick, W. R., & Velicer, W. F. (1986). Comparison of five rules for determining the number of components to retain. Psychological Bulletin, 99, 432-442.

Franklin, S. B., Gibson, D. J., Robertson, P. A., Pohlmann, J. T., & Fralish, J. S. (1995). Parallel analysis: a method for determining significant principal components. Journal of Vegetation Science, 6(1), 99-106.

Crawford, A. V., Green, S. B., Levy, R., Lo, W. J., Scott, L., Svetina, D., & Thompson, M. S. (2010). Evaluation of parallel analysis methods for determining the number of factors. Educational and Psychological Measurement, 70(6), 885–901.

La scelta del numero di componenti nell’Analisi in Componenti Principali

Nei differenti ambiti di ricerca, accade spesso di dover analizzare contemporaneamente un numero molto elevato di variabili. In questi casi, per comprendere meglio le strutture latenti all’interno dei dati, potrebbe risultare molto utile sostituire le p variabili originarie con un numero ridotto di variabili artificiali (o componenti principali) che garantiscono la sintesi e la facilità di lettura, garantendo la minor perdita possibile di informazione. Una metodologia statistica multivariata che permette di fare ciò è l’Analisi delle Componenti Principali o ACP. In breve, si suppone che le variabili artificiali siano una combinazione lineare delle variabili osservate.

Partendo da una matrice dei dati, l’applicazione dell’ACP consente di sostituire alle p variabili (tra loro correlate) un nuovo insieme di variabili artificiali dette COMPONENTI PRINCIPALI che:

1.sono tra loro incorrelate (ortogonali);
2.sono elencate in ordine decrescente rispetto alla loro varianza (nel senso che la prima CP sarà la combinazione lineare di massima varianza, la seconda CP sarà la seconda combinazione lineare in termini di varianza e cosi via).

Alla fine si otterranno tante combinazioni lineari quante sono le variabili originarie e spetterà al ricercatore scegliere il numero di Componenti principali ritenuto idoneo e interpretarle. Si pongono pertanto diversi problemi che impongono allo studioso delle decisioni che influenzeranno l’esito della ricerca. Tuttavia, esistono delle “regole” che aiutano il ricercatore a minimizzare la componente soggettiva della ricerca, in termini di scelta del numero di CP ed interpretazione delle stesse, ed essere il più obiettivo possibile.

Per quanto riguarda la scelta del numero di componenti da utilizzare, si fa riferimento a tre criteri:
1. Quota di varianza totale spiegata
2. Scree-graph
3. Eigenvalue one o Regola di Kaiser

Solitamente, per individuare il numero di componenti da utilizzare, si cerca di seguire congiuntamente i tre criteri.

Secondo il primo criterio, si deve considerare un numero di CP tale che esse tengano conto di una percentuale sufficientemente elevata di varianza totale (ad esempio, almeno il 70%). Nel definire la percentuale minima di varianza accettabile, occorre tener conto del numero di variabili originarie; pertanto al crescere del numero di variabili potrà essere accettata una percentuale minore di varianza spiegata.
Il secondo criterio, fa uso di un grafico chiamato scree-graph degli autovalori in funzione del numero di CP. Poiché gli autovalori sono decrescenti, il grafico assume la forma di una spezzata con pendenza sempre negativa (vedi sotto).

Analizzando il grafico, si potrà individuare un punto nel quale si manifesta una brusca variazione di pendenza, in corrispondenza della quale si individua il numero k di CP da considerare. Tuttavia, può accadere che la diminuzione degli autovalori sia graduale e il grafico non evidenzi salti evidenti. Inoltre, in letteratura si distinguono differenti posizioni in relazione all’inclusione (Cattel, 1966) o all’esclusione (Harman, 1976) della CP in corrispondenza del gomito.
Il terzo criterio, suggerisce di considerare tutte le CP il cui autovalore è maggiore di 1. La “ratio” di questo criterio deriva dal fatto che l’autovalore di una CP è uguale alla sua varianza e che operando su variabili standardizzate queste hanno varianza unitaria.
 Pertanto, si decide di mantenere una CP solo se essa spiega una quota di varianza totale maggiore di quella di una singola variabile (D’Urso, 2008).
Per chiarire meglio il tutto, facciamo un esempio utilizzando parte del dataset “microtus” presente all’interno della libreria “Flury” e prendendo in considerazione solamente il gruppo “subterraneus”.

library(Flury)
data(microtus)
subterraneus <- microtus[microtus$Group=="subterraneus", -1]

Il dataset è composto da 45 casi e 8 variabili, espresse su unità di misura differenti. L’ACP si applica direttamente sul nostro dataset attraverso la funzione “prcomp” presente nel pacchetto “stats”.

acp <- prcomp(subterraneus, scale=TRUE)

La funzione prcomp comprende vari argomenti, per un approfondimento dei quali si rimanda all’help. Nel nostro caso sarà sufficiente specificare due argomenti: il primo è rappresentato dalla matrice dei dati; il secondo argomento, scale=TRUE, ha la funzione di normalizzare le variabili; pertanto la nostra analisi delle componenti principali sarà condotta sulla matrice di correlazione.
L’output della funzione è il seguente:

Standard deviations:
[1] 2.0054992 1.0994572 0.9225566 0.8246091 0.6753426 0.6199444 0.5229512 0.3523907

Rotation:
               PC1         PC2          PC3         PC4         PC5          PC6         PC7         PC8
M1Left  -0.3088583 -0.10891920  0.545170129 -0.67358259 -0.04817153  0.109379942 -0.35688311  0.01300575
M2Left  -0.3521551 -0.12119956 -0.604983110  0.03425009 -0.32677472  0.388151140 -0.48524914 -0.03478596
M3Left  -0.3740423 -0.01511945  0.322459880  0.39303305 -0.65913792 -0.300524226  0.08555331  0.26327065
Foramen -0.3411396 -0.55597610 -0.016570304  0.04551210  0.31006428  0.319768896  0.42584195  0.43867421
Pbone   -0.2600426  0.61212751 -0.297386290 -0.44021781 -0.11785081  0.007275735  0.44230634  0.25665302
Length  -0.4664496 -0.08624148  0.009645834  0.02441142  0.03028856 -0.044189223  0.35170048 -0.80481979
Height  -0.2748288  0.53074157  0.316756946  0.43889866  0.34728294  0.443174798 -0.17858683  0.03021476
Rostrum -0.4044517  0.01688066 -0.208957907  0.03933720  0.47421882 -0.668130405 -0.31190640  0.149518

Le informazioni rilevanti si ricavano dalla funzione summary dell’output di prcomp. La funzione restituisce una matrice che riporta: nella prima riga le radici quadrate degli autovalori, nella seconda le proporzioni di varianza spiegata da ciascuna componente e nella terza riga le percentuali di varianza cumulata.

summary(acp)
Importance of components:
                          PC1    PC2    PC3    PC4     PC5     PC6     PC7     PC8
Standard deviation     2.0055 1.0995 0.9226 0.8246 0.67534 0.61994 0.52295 0.35239
Proportion of Variance 0.5028 0.1511 0.1064 0.0850 0.05701 0.04804 0.03418 0.01552
Cumulative Proportion  0.5028 0.6539 0.7602 0.8452 0.90225 0.95029 0.98448 1.00000

Per scegliere il numero di componenti da utilizzare ci possiamo rifare alle tre regole precedentemente descritte.
Secondo la prima regola, sembrerebbe sensato prendere in considerazione le prime due CP, poiché spiegano il 65% della varianza totale, quota che può essere valutata congrua, considerando la riduzione di variabili da 8 a 2.
Anche secondo la regola di Kaiser, dovremo prendere in considerazione le prime due CP, poiché risultano le uniche con autovalore maggiore di uno.
Infine, per utilizzare come metodo di giudizio anche lo scree-graph, non dobbiamo fare altro che digitare il seguente comando:

plot(acp, type="lines")

Analizzando il nostro grafico, possiamo notare una brusca variazione di pendenza in corrispondenza della seconda CP. Pertanto anche l’analisi dello scree-graph conferma quanto visto con i due precedenti criteri e pertanto scegliamo le prime due componenti principali.

Come avrete notato, nonostante l’utilizzo di tre criteri differenti, rimane comunque un certo margine di soggettività nella scelta, che talvolta può risultare determinante nel prendere in considerazione una CP in più o in meno, con evidenti distorsioni nei risultati della ricerca.