La funzione principal() del pacchetto psych facilita l'analisi delle componenti principali con dati sociali e psicometrici, anche provenienti da surveys.

Come si è visto in L'analisi in componenti principali, per ruotare le componenti estratte, è possibile usare funzioni come varimax() alla matrice dei pesi componenziali ottenuta con princomp() e prcomp(), ma sarebbe poi necessario ricalcolare i punteggi individuali.

Questa funzione, invece:

• consente di scegliere il metodo di rotazione (rotate, di default 'varimax') e ricalcola i punteggi (scores = TRUE);
• consente di pesare i casi (argomento weight);
• può operare l'imputazione dei casi mancanti (missing = TRUE; per default è FALSE), coma la mediana o con la media (argomento `impute`)

In comune con le funzioni per l'analisi fattoriale classica (esplorativa):

• richiede di indicare il numero di componenti (qui vengono chiamati fattori) da estrarre;
• fornisce misure di valutazione del modello (residui e fitness), dove per 'modello' si deve intendere il numero di componenti scelte.

Le funzioni di factoextra (viste in L'analisi in componenti principali) non sono però utilizzabili con i risultati prodotti.

I grafici sono costruiti sulla libreria base graphics (vedi i grafici con R)

library(tidyverse)
library(psych)

Useremo il dataset decathlon di FactoMineR:

data(decathlon, package = "FactoMineR")
# prime righe del dataset
head(decathlon)

##          100m Long.jump Shot.put High.jump  400m 110m.hurdle Discus Pole.vault
## SEBRLE  11.04      7.58    14.83      2.07 49.81       14.69  43.75       5.02
## CLAY    10.76      7.40    14.26      1.86 49.37       14.05  50.72       4.92
## KARPOV  11.02      7.30    14.77      2.04 48.37       14.09  48.95       4.92
## BERNARD 11.02      7.23    14.25      1.92 48.93       14.99  40.87       5.32
## YURKOV  11.34      7.09    15.19      2.10 50.42       15.31  46.26       4.72
## WARNERS 11.11      7.60    14.31      1.98 48.68       14.23  41.10       4.92
...

Statistiche descrittive e matrice delle correlazioni

Iniziamo sempre dallo studio delle variabili e delle correlazioni.

La funzione describe() fa parte del pacchetto psych e produce una tabella come quella che segue

describe(decathlon[1:10],
         # non inserisco le misure di asimmetria
         skew = FALSE, 
         # aggiungo la differenza interquartile
         IQR = TRUE) 

##             vars  n   mean    sd    min    max range   se   IQR
## 100m           1 41  11.00  0.26  10.44  11.64  1.20 0.04  0.29
## Long.jump      2 41   7.26  0.32   6.61   7.96  1.35 0.05  0.45
## Shot.put       3 41  14.48  0.82  12.68  16.36  3.68 0.13  1.09
## High.jump      4 41   1.98  0.09   1.85   2.15  0.30 0.01  0.12
## 400m           5 41  49.62  1.15  46.81  53.20  6.39 0.18  1.37
## 110m.hurdle    6 41  14.61  0.47  13.97  15.67  1.70 0.07  0.77
## Discus         7 41  44.33  3.38  37.92  51.65 13.73 0.53  4.17
## Pole.vault     8 41   4.76  0.28   4.20   5.40  1.20 0.04  0.42
## Javeline       9 41  58.32  4.83  50.31  70.52 20.21 0.75  5.62
## 1500m         10 41 279.02 11.67 262.10 317.00 54.90 1.82 14.08
...

Passiamo a guardare la matrice delle correlazioni:

cor(decathlon[1:10]) %>% 
  round(2)

##              100m Long.jump Shot.put High.jump  400m 110m.hurdle Discus
## 100m         1.00     -0.60    -0.36     -0.25  0.52        0.58  -0.22
## Long.jump   -0.60      1.00     0.18      0.29 -0.60       -0.51   0.19
## Shot.put    -0.36      0.18     1.00      0.49 -0.14       -0.25   0.62
## High.jump   -0.25      0.29     0.49      1.00 -0.19       -0.28   0.37
## 400m         0.52     -0.60    -0.14     -0.19  1.00        0.55  -0.12
## 110m.hurdle  0.58     -0.51    -0.25     -0.28  0.55        1.00  -0.33
## Discus      -0.22      0.19     0.62      0.37 -0.12       -0.33   1.00
## Pole.vault  -0.08      0.20     0.06     -0.16 -0.08        0.00  -0.15
## Javeline    -0.16      0.12     0.37      0.17  0.00        0.01   0.16
## 1500m       -0.06     -0.03     0.12     -0.04  0.41        0.04   0.26
...

Quando le variabili sono numerose, può essere comodo usare una heatmap per avere una visione d'insieme del dataset.

# colori
col <- colorRampPalette(c("blue", "white", "red"))(20)
 
# heatmap
heatmap(cor(decathlon[1:10]),
        col = col, 
        symm = TRUE)

\

La rotazione Varimax massimizza la somma della varianza dei pesi componenziali al quadrato. Si ottengono dunque, per ciascuna componente, pesi più elevati per un numero inferiore di variabili.

# rotazione varimax
PCA.rot <- principal(decathlon[1:10], 
                     rotate="varimax", 
                     nfactors=4, 
                     scores=TRUE)
PCA.rot

## Principal Components Analysis
## Call: principal(r = decathlon[1:10], nfactors = 4, rotate = "varimax", 
##     scores = TRUE)
## Standardized loadings (pattern matrix) based upon correlation matrix
##               RC1   RC2   RC3   RC4   h2   u2 com
## 100m         0.75 -0.27 -0.16 -0.07 0.67 0.33 1.4
## Long.jump   -0.81  0.11  0.17  0.09 0.71 0.29 1.2
## Shot.put    -0.16  0.77  0.13  0.37 0.78 0.22 1.6
## High.jump   -0.25  0.62 -0.25  0.16 0.54 0.46 1.8
## 400m         0.87  0.08  0.19 -0.08 0.80 0.20 1.1
## 110m.hurdle  0.76 -0.30  0.07  0.19 0.70 0.30 1.5
## Discus      -0.12  0.85 -0.01 -0.06 0.74 0.26 1.0
## Pole.vault  -0.19 -0.24  0.84  0.15 0.82 0.18 1.3
## Javeline     0.02  0.23  0.01  0.89 0.84 0.16 1.1
## 1500m        0.24  0.40  0.68 -0.43 0.87 0.13 2.7
## 
##                        RC1  RC2  RC3  RC4
## SS loadings           2.75 2.16 1.34 1.22
## Proportion Var        0.28 0.22 0.13 0.12
## Cumulative Var        0.28 0.49 0.63 0.75
## Proportion Explained  0.37 0.29 0.18 0.16
## Cumulative Proportion 0.37 0.66 0.84 1.00
## 
## Mean item complexity =  1.5
## Test of the hypothesis that 4 components are sufficient.
## 
## The root mean square of the residuals (RMSR) is  0.08 
##  with the empirical chi square  22.28  with prob <  0.022 
## 
## Fit based upon off diagonal values = 0.93

print.psych

Il pacchetto prevede una funzione che consente di "stampare" gli output selezionando e riorganizzandone gli elementi: in questo caso, possiamo filtrare i factor loadings al di sopra di una certa soglia (cut, qui 0,3) e riordinare gli items (sort):

# output ordinato
print.psych(PCA.rot, sort = T, cut = 0.3)

## Principal Components Analysis
## Call: principal(r = decathlon[1:10], nfactors = 4, rotate = "varimax", 
##     scores = TRUE)
## Standardized loadings (pattern matrix) based upon correlation matrix
##             item   RC1   RC2   RC3   RC4   h2   u2 com
## 400m           5  0.87                   0.80 0.20 1.1
## Long.jump      2 -0.81                   0.71 0.29 1.2
## 110m.hurdle    6  0.76                   0.70 0.30 1.5
## 100m           1  0.75                   0.67 0.33 1.4
## Discus         7        0.85             0.74 0.26 1.0
## Shot.put       3        0.77        0.37 0.78 0.22 1.6
## High.jump      4        0.62             0.54 0.46 1.8
## Pole.vault     8              0.84       0.82 0.18 1.3
## 1500m         10        0.40  0.68 -0.43 0.87 0.13 2.7
## Javeline       9                    0.89 0.84 0.16 1.1
## 
##                        RC1  RC2  RC3  RC4
## SS loadings           2.75 2.16 1.34 1.22
## Proportion Var        0.28 0.22 0.13 0.12
## Cumulative Var        0.28 0.49 0.63 0.75
## Proportion Explained  0.37 0.29 0.18 0.16
## Cumulative Proportion 0.37 0.66 0.84 1.00
## 
## Mean item complexity =  1.5
## Test of the hypothesis that 4 components are sufficient.
## 
## The root mean square of the residuals (RMSR) is  0.08 
##  with the empirical chi square  22.28  with prob <  0.022 
## 
## Fit based upon off diagonal values = 0.93

Elementi dell'output

• h2 = comunalità (varianza spiegata dalla componente; sarebbe uguale a 1 se avessimo considerato tutte le componenti)
• u2 = unicità di ciascuna variabile (varianza unica della variabile: $1-h2$)
• com = complessità (complexity): «funzione del numero dei fattori che descrivono una variabile in una determinata soluzione fattoriale» (Richard J. Hofmann 1973, p. 1, traduzione mia); varia da 1 a $nf$ (numero dei fattori).
• root mean square of the residuals: deve approssimarsi a 0;
• test del chi quadrato, dove i gradi di libertà del modello sono calcolati a partire dal numero di variabili e di componenti indicate, ovvero: n corr items ($p (p-1)/2$) -- n corr items-fattori -- n corr fattori ($p nf + nf*(nf-1)/2$)

10*9 / 2 - (10*4 - 4*3/2)

## [1] 11

• complexity: media dei valori della colonna com; qualità della riproduzione delle variabili da parte del modello;
• fit: riduzione/riproduzione della matrice delle correlazioni da parte del modello.

# autovalori
PCA.rot$values

##  [1] 3.2719055 1.7371310 1.4049167 1.0568504 0.6847735 0.5992687 0.4512353
##  [8] 0.3968766 0.2148149 0.1822275

plot(PCA.rot$values, type = "b")

Qui sono evidenziati gli autovalori in prima ipotesi significativi, ovvero quelli che:

• precedono il “gomito” della curva;
• sono maggiori di 1, ovvero spiegano più varianza di una singola variabile (criterio di Kaiser)

Confrontando i risultati con e senza rotazione, osserveremo che i pesi componenziali si distribuiscono diversamente:

# senza rotazione
PCA <- principal(decathlon[1:10], 
                 rotate="none", 
                 nfactors=4, 
                 scores=TRUE)
PCA$loadings

## Loadings:
##             PC1    PC2    PC3    PC4   
## 100m        -0.775  0.187 -0.184       
## Long.jump    0.742 -0.345  0.182  0.102
## Shot.put     0.623  0.598         0.191
## High.jump    0.572  0.350 -0.260 -0.136
## 400m        -0.680  0.569  0.131       
## 110m.hurdle -0.746  0.229         0.291
## Discus       0.552  0.606        -0.260
## Pole.vault         -0.180  0.692  0.552
## Javeline     0.277  0.317 -0.390  0.712
## 1500m               0.474  0.782 -0.161
## 
##                  PC1   PC2   PC3   PC4
## SS loadings    3.272 1.737 1.405 1.057
## Proportion Var 0.327 0.174 0.140 0.106
## Cumulative Var 0.327 0.501 0.641 0.747

PCA.rot$loadings

## Loadings:
##             RC1    RC2    RC3    RC4   
## 100m         0.753 -0.270 -0.160       
## Long.jump   -0.815  0.110  0.169       
## Shot.put    -0.156  0.774  0.134  0.375
## High.jump   -0.247  0.623 -0.250  0.155
## 400m         0.869         0.194       
## 110m.hurdle  0.757 -0.299         0.186
## Discus      -0.116  0.852              
## Pole.vault  -0.187 -0.238  0.839  0.149
## Javeline            0.228         0.886
## 1500m        0.243  0.398  0.678 -0.434
## 
##                  RC1   RC2   RC3   RC4
## SS loadings    2.752 2.162 1.341 1.217
## Proportion Var 0.275 0.216 0.134 0.122
## Cumulative Var 0.275 0.491 0.625 0.747

La varianza spiegata resta la stessa (così come le misure di adeguatezza del modello a 4 componenti), ma le componenti dovrebbero essere più facili da interpretare, almeno dal punto di vista della ripartizione dei pesi. Resta da valutare l'adeguatezza delle componenti stessi.

Da questo output si comprende intuitivamente il senso dell'indice di complessità (com), tanto più alto quando più numerose sono le componenti su cui la varianza dell'item stesso per così dire si distribuisce.

Vedi anche Analisi delle corrispondenze

Per calcolare i contributi possiamo procedere elevando i pesi componenziali al quadrato e calcolato le proporzioni o percentuali di colonna (la somma dei pesi al quadrato per colonna è infatti pari alla varianza spiegata da ciascun fattore):

# contributi
as.table(PCA.rot$loadings)^2 %>% 
  proportions(1) %>% round(2)

##              RC1  RC2  RC3  RC4
## 100m        0.85 0.11 0.04 0.01
## Long.jump   0.93 0.02 0.04 0.01
## Shot.put    0.03 0.77 0.02 0.18
## High.jump   0.11 0.73 0.12 0.04
## 400m        0.94 0.01 0.05 0.01
## 110m.hurdle 0.82 0.13 0.01 0.05
## Discus      0.02 0.98 0.00 0.00
## Pole.vault  0.04 0.07 0.86 0.03
## Javeline    0.00 0.06 0.00 0.94
## 1500m       0.07 0.18 0.53 0.22

\

I coseni quadri sono i pesi componenziali al quadrato, la cui somma di riga è uguale a 1 (uniqueness).

# senza rotazione
PCA$loadings^2 %>% as.table() %>% round(2)

##              PC1  PC2  PC3  PC4
## 100m        0.60 0.04 0.03 0.00
## Long.jump   0.55 0.12 0.03 0.01
## Shot.put    0.39 0.36 0.00 0.04
## High.jump   0.33 0.12 0.07 0.02
## 400m        0.46 0.32 0.02 0.00
## 110m.hurdle 0.56 0.05 0.01 0.08
## Discus      0.31 0.37 0.00 0.07
## Pole.vault  0.00 0.03 0.48 0.30
## Javeline    0.08 0.10 0.15 0.51
## 1500m       0.00 0.22 0.61 0.03

# con rotazione varimax
PCA.rot$loadings^2 %>% as.table() %>% round(2)

##              RC1  RC2  RC3  RC4
## 100m        0.57 0.07 0.03 0.00
## Long.jump   0.66 0.01 0.03 0.01
## Shot.put    0.02 0.60 0.02 0.14
## High.jump   0.06 0.39 0.06 0.02
## 400m        0.75 0.01 0.04 0.01
## 110m.hurdle 0.57 0.09 0.01 0.03
## Discus      0.01 0.73 0.00 0.00
## Pole.vault  0.04 0.06 0.70 0.02
## Javeline    0.00 0.05 0.00 0.78
## 1500m       0.06 0.16 0.46 0.19

La variabile '100m' è rappresentata leggermente peggio, ma 'Long.jump' è rappresentata decisamente meglio sulla prima componente, passando da 0.55 a 0.66. 400m passa da 0.46 a 0.75

I grafici

Di default, il comando plot() produce lo scatterplot delle variabili per tutte le componenti:

plot(PCA.rot)

Per scegliere le componenti da rappresentare, usiamo l'argomento choose; per le etichette, l'argomento labels (che richiede un vettore carattere).

plot(PCA.rot,
     # componenti
     choose = c(1,2),
     # titolo
     main = "Varimax",
     # etichette con i nomi delle righe
     labels = names(decathlon), 
     # limiti degli assi
     xlim = c(-1, +1), ylim = c(-1, +1), 
     # dimensione del carattere
     cex = 0.80)

Il comando biplot() produce il grafico PCA classico, sempre per tutte le componenti:

# grafico delle componenti
biplot(PCA.rot)

\ E può essere modificato (ad esempio) come segue:

# margini
par(mar=c(4, 0.1,4,0.1))
# grafico delle prime due componenti
biplot(PCA.rot, choose = c(1,2),
       main = "Varimax")
abline(h = 0, lty = 2)
abline(v = 0, lty = 2)

Possiamo infine prendere in considerazioni i residui, ovvero le differenze fra la matrice originaria delle correlazioni e la matrice dei pesi componenziali.

# matrice dei residui (uniqueness)
factor.residuals(r = cor(decathlon[,1:10]), 
                 f = PCA.rot2$loading)

La funzione resid(), stampa i residui in forma di matrice triangolare:

resid(PCA.rot)

##             100m  Lng.j Sht.p Hgh.j 400m  110m. Discs Pl.vl Javln 1500m
## 100m         0.33                                                      
## Long.jump    0.08  0.29                                                
## Shot.put     0.02 -0.09  0.22                                          
## High.jump    0.08  0.05 -0.06  0.46                                    
## 400m        -0.09  0.07 -0.06  0.04  0.20                              
## 110m.hurdle -0.05  0.11  0.02  0.08 -0.09  0.30                        
## Discus       0.09  0.01 -0.04 -0.18 -0.09  0.03  0.26                  
## Pole.vault   0.14 -0.08  0.05  0.13 -0.05 -0.02  0.04  0.18            
## Javeline    -0.05  0.02 -0.13 -0.10  0.04 -0.10  0.01 -0.11  0.16      
## 1500m       -0.06  0.05 -0.08  0.00  0.00  0.01 -0.07 -0.12  0.10  0.13

resid(PCA.rot) %>% 
  # senza ordinamento e dendrogramma
  heatmap(Rowv =NA, Colv =NA)

\

Si tratta di un metodo di rotazione non ortogonale, che in questo caso potrebbe essere più adatto, dato che le diverse capacità atletiche potrebbero non essere indipendenti.

PCA.rot2 <- principal(decathlon[1:10], 
                     rotate="promax", 
                     nfactors=4, 
                     scores=TRUE)

PCA.rot2$loadings

## Loadings:
##             RC1    RC2    RC3    RC4   
## 100m         0.737 -0.195 -0.183       
## Long.jump   -0.818         0.213       
## Shot.put            0.785         0.400
## High.jump   -0.159  0.619 -0.288  0.135
## 400m         0.903  0.170  0.132       
## 110m.hurdle  0.767 -0.222         0.267
## Discus              0.862              
## Pole.vault  -0.202 -0.271  0.893  0.183
## Javeline     0.145  0.254         0.918
## 1500m        0.261  0.420  0.594 -0.363
## 
##                  RC1   RC2   RC3   RC4
## SS loadings    2.771 2.173 1.356 1.260
## Proportion Var 0.277 0.217 0.136 0.126
## Cumulative Var 0.277 0.494 0.630 0.756

ACP_psych.R

# pacchetti
library(tidyverse)
library(psych)
 
# dati 
data(decathlon, package = "FactoMineR")
 
describe(decathlon[1:10],
         skew = FALSE, IQR = TRUE) 
 
# matrice delle correlazioni
cor(decathlon[1:10]) %>% round(2)
 
# heatmap
 
# colori dal blu al rosso
col <- colorRampPalette(c("blue", "white", "red"))(20)
 
heatmap(cor(decathlon[1:10]),
        col = col, 
        symm = TRUE)
 
# rotazione varimax
PCA.rot <- principal(decathlon[1:10], 
                     rotate="varimax", 
                     nfactors=4, 
                     scores=TRUE)
PCA.rot
 
# autovalori
plot(PCA.rot$values, type = "b")
 
# senza rotazione
PCA <- principal(decathlon[1:10], 
                 rotate="none", 
                 nfactors=4, 
                 scores=TRUE)
# confronto
PCA$loadings
PCA.rot$loadings
 
# output ordinato
print.psych(PCA.rot, sort = T, cut = 0.3)
 
# contributi
as.table(PCA.rot$loadings)^2 %>% 
  proportions(1) %>% round(2)
 
# cos2 senza rotazione
PCA$loadings^2 %>% as.table() %>% round(2)
 
# cos2 con rotazione varimax
PCA.rot$loadings^2 %>% as.table() %>% round(2)
 
# grafico
plot(PCA.rot)
 
plot(PCA.rot,
     # componenti
     choose = c(1,2),
     # titolo
     main = "Varimax",
     # etichette con i nomi delle righe
     labels = names(decathlon), 
     # limiti degli assi
     xlim = c(-1, +1), ylim = c(-1, +1), 
     # dimensione del carattere
     cex = 0.80)
 
# grafico delle componenti
biplot(PCA.rot)
 
# grafico delle prime due componenti
biplot(PCA.rot, choose = c(1,2),
       main = "Varimax")
abline(h = 0, lty = 2)
abline(v = 0, lty = 2)
 
# matrice dei residui (uniqueness)
factor.residuals(r = cor(decathlon[,1:10]),
                 f = PCA.rot2$loading)
 
resid(PCA.rot)
 
resid(PCA.rot) %>% 
  # senza ordinamento e dendrogramma
  heatmap(Rowv =NA, Colv =NA)
 
 
# promax
PCA.rot2 <- principal(decathlon[1:10], 
                     rotate="promax", 
                     nfactors=4, 
                     scores=TRUE)
 
PCA.rot2$loadings