Ricerca Sociale con R

Una wiki per l'analisi dei dati con R

Strumenti Utente

Strumenti Sito


Il Laboratorio di Analisi dei Dati con R, dell'Università di Teramo su piattaforma Meet, inizia il 9 aprile 2021 - Iscrizione - email
r:analisi_esplorativa:analisi_in_componenti_principali

L'analisi in componenti principali

L'analisi in componenti principali è una tecnica fattoriale di analisi multivariata dei dati.

Una presentazione dettagliata della tecnica è disponibile su Wikipedia.

Testi in italiano, con esempi in campo sociologico sono (fra gli altri):

La funzione prcomp

prcomp(x, retx = TRUE, center = TRUE, scale. = FALSE,
       tol = NULL, ...)

Per mostrare i comandi e gli output, utilizziamo il Datasets di esempio “USArrests”.

Come prima cosa, andrebbe sempre controllata la linearità delle relazioni fra le variabili contenute nel dataset, ad esempio con uno Scatterplot (grafico a dispersione).

 plot(USArrests)
&nolink |

Utilizziamo la funzione prcomp, con l'opzione scale=TRUE, ovvero normalizzando le variabili (varianza = 1).

> PCA <- prcomp(USArrests, scale = TRUE)
> print(PCA)
Standard deviations:
[1] 1.5748783 0.9948694 0.5971291 0.4164494

Rotation:
                PC1        PC2        PC3         PC4
Murder   -0.5358995  0.4181809 -0.3412327  0.64922780
Assault  -0.5831836  0.1879856 -0.2681484 -0.74340748
UrbanPop -0.2781909 -0.8728062 -0.3780158  0.13387773
Rape     -0.5434321 -0.1673186  0.8177779  0.08902432

L'output prodotto con print(PCA) consiste nelle deviazioni standard delle componenti e nella matrice di rotazione (con i variable loadings)

> summary(PCA)
Importance of components:
                          PC1    PC2     PC3     PC4
Standard deviation     1.5749 0.9949 0.59713 0.41645
Proportion of Variance 0.6201 0.2474 0.08914 0.04336
Cumulative Proportion  0.6201 0.8675 0.95664 1.00000

Il comando summary() in questo caso restituisce una serie di valori che consentono di valutare la rilevanza relativa di ciascuna componente, ovvero: la deviazione standard di ciascuna componente; la proporzione di varianza spiegata, e le proporzioni cumulate.

Per scegliere di utilizzare solo alcune variabili e non altre:

> prcomp(~Murder + Assault + Rape, data = USArrests, 
+     scale = TRUE)
Standard deviations:
[1] 1.5357670 0.6767949 0.4282154

Rotation:
               PC1        PC2        PC3
Murder  -0.5826006  0.5339532 -0.6127565
Assault -0.6079818  0.2140236  0.7645600
Rape    -0.5393836 -0.8179779 -0.1999436

Output

L'oggetto (output) PCA creato dall'analisi in componenti principali contiene i seguenti valori:

  • sdev = deviazioni standard delle componenti principali
  • rotation = la matrice di rotazione con i variable loadings (consente di interpretare i fattori, in base al peso che hanno le variabili nel definirli)
  • x = le coordinate fattoriali dei casi (se retx=TRUE, come di default)
  • center = il valore corrispondente all'attuale origine dell'asse, per ciascuna variabile (se center=TRUE, come di default)
  • scale = come center, ma per la varianza di ciascuna variabile (se scale=TRUE)

Ciascuno dei valori contenuti nell'output può infatti essere visualizzato chiamando, ad es., per ottenere le coordinate fattoriali:

> PCA$x
                       PC1         PC2         PC3          PC4
Alabama        -0.97566045  1.12200121 -0.43980366  0.154696581
Alaska         -1.93053788  1.06242692  2.01950027 -0.434175454
Arizona        -1.74544285 -0.73845954  0.05423025 -0.826264240
Arkansas        0.13999894  1.10854226  0.11342217 -0.180973554
....

La matrice delle correlazioni fra i nuovi valori e i vecchi si ottiene con il comando cor() applicato alle due matrici.

> cor.PCA=cor(USArrests, PCA$x)
> cor.PCA
                PC1        PC2        PC3         PC4
Murder   -0.8439764  0.4160354 -0.2037600  0.27037052
Assault  -0.9184432  0.1870211 -0.1601192 -0.30959159
UrbanPop -0.4381168 -0.8683282 -0.2257242  0.05575330
Rape     -0.8558394 -0.1664602  0.4883190  0.03707412

Grafici

La funzione plot() restituisce un grafico come questo:

> plot(PCA)
&nolink |

mentre la funzione “screeplot()” rappresenta graficamente i fattori in questo modo:

> screeplot(PCA,type=c("lines")) 
&nolink |

aiutandoci a stabilire in numero di componenti principali da tenere in considerazione per rappresentare il fenomeno.

La funzione biplot() fornisce invece una rappresentazione sintetica dei casi sul piano fattoriale, utilizzando la matrice PCA$x, che può eventualmente essere esportata per includere le coordinate fattoriali come variabili nella matrice originaria:

> biplot(PCA)
&nolink |

Script

prcomp-ex.R
plot(USArrests)
PCA<-prcomp(USArrests, scale = TRUE)
print(PCA)
summary(PCA)
plot(PCA)
biplot(PCA)
PCA$x
cor.PCA=cor(USArrests, PCA$x)
cor.PCA

La funzione princomp e RCommander

La funzione princomp, calcola l'ACP con la risoluzione matriciale, ovvero calcolando gli autovalori delle matrice delle correlazioni o delle covarianze (e non della matrice originaria).

princomp(x, cor = FALSE, scores = TRUE, covmat = NULL,
         subset = rep(TRUE, nrow(as.matrix(x))), ...)

Nella console di R

> PCA1 <- princomp(USArrests, cor = TRUE)
> print(PCA1)
Call:
princomp(x = USArrests, cor = TRUE)

Standard deviations:
   Comp.1    Comp.2    Comp.3    Comp.4 
1.5748783 0.9948694 0.5971291 0.4164494 

4  variables and  50 observations.

L'output prodotto con print(PCA1) non contiene la matrice delle componenti (variable loadings), che può essere però chiamata con il comando:

> loadings(PCA1)

Loadings:
         Comp.1 Comp.2 Comp.3 Comp.4
Murder   -0.536  0.418 -0.341  0.649
Assault  -0.583  0.188 -0.268 -0.743
UrbanPop -0.278 -0.873 -0.378  0.134
Rape     -0.543 -0.167  0.818       

               Comp.1 Comp.2 Comp.3 Comp.4
SS loadings      1.00   1.00   1.00   1.00
Proportion Var   0.25   0.25   0.25   0.25
Cumulative Var   0.25   0.50   0.75   1.00

oppure anche con

> PCA1$loadings

In questo caso, i pesi fattoriali (variable loadings) si riferiscono agli autovalori della matrice delle correlazioni.

Il comando summary():

> summary(PCA1)
Importance of components:
                          Comp.1    Comp.2    Comp.3     Comp.4
Standard deviation     1.5748783 0.9948694 0.5971291 0.41644938
Proportion of Variance 0.6200604 0.2474413 0.0891408 0.04335752
Cumulative Proportion  0.6200604 0.8675017 0.9566425 1.00000000

Anche i comandi plot() screeplot() e biplot() restituiscono grafici corrispondenti.

L'oggetto (output) PCA1 è però diverso, in quanto diverso è il metodo di calcolo delle componenti. Esso quindi contiene:

  • sdev = deviazioni standard delle componenti principali
  • loadings (anziché rotation) = la matrice con i variable loadings
  • scores (anziché x) = le coordinate fattoriali dei casi (può essere omessa, con retx=FALSE)
  • center = il valore corrispondente all'attuale origine dell'asse, per ciascuna variabile (la media)
  • scale = il fattore di scala delle variabili

Quindi, per ottenere le coordinate fattoriali, il comando sarà:

> PCA1$scores
                    Comp.1      Comp.2      Comp.3       Comp.4
Alabama        -0.98556588  1.13339238 -0.44426879  0.156267145
Alaska         -1.95013775  1.07321326  2.04000333 -0.438583440
Arizona        -1.76316354 -0.74595678  0.05478082 -0.834652924
Arkansas        0.14142029  1.11979678  0.11457369 -0.182810896
....

Ed è con la matrice PCA1$scores che si potrà anche calcolare la matrice delle correlazioni variabili-componenti:

> cor.PCA1=cor(USArrests, PCA1$scores)
> cor.PCA1
             Comp.1     Comp.2     Comp.3      Comp.4
Murder   -0.8439764  0.4160354 -0.2037600  0.27037052
Assault  -0.9184432  0.1870211 -0.1601192 -0.30959159
UrbanPop -0.4381168 -0.8683282 -0.2257242  0.05575330
Rape     -0.8558394 -0.1664602  0.4883190  0.03707412

Infine, volendo ruotare gli autovalori [o alcuni di essi] con il metodo varimax o promax:

> varimax(PCA1$loadings[,1:2])

oppure

> promax(PCA1$loadings[,1:2])

restituisce direttamente i valori degli autovalori ruotati.

Per scegliere di utilizzare alcune variabili e non altre:

> princomp(~Murder + Assault + UrbanPop, data = USArrests, na.action = na.exclude, 
+     cor = TRUE)
Call:
princomp(formula = ~Murder + Assault + UrbanPop, data = USArrests, 
    na.action = na.exclude, cor = TRUE)

Standard deviations:
   Comp.1    Comp.2    Comp.3 
1.3656547 0.9795415 0.4189100 

3  variables and  50 observations.

- Per la rotazione dei fattori e il calcolo dei punteggi individuali, vedi anche ACP con rotazione delle componenti

Script

princomp-ex.R
PCA1<-princomp(USArrests, cor = TRUE)
print(PCA1)
loadings(PCA1)
summary(PCA1)
PCA1$scores
plot(PCA1)
biplot(PCA1)
PCA1$scores
cor.PCA1=cor(USArrests, PCA1$scores)
cor.PCA1
varimax(PCA1$loadings[,1:2])
promax(PCA1$loadings[,1:2])

In RCommander

Utilizzando lo stesso esempio di prima, possiamo vedere come viene utilizzata in RCommander.

Come prima cosa, carichiamo il dataset “USArrests” (Dati / Dati presenti nei pacchetti / Leggi i dati da un pacchetto caricato → scrivere il nome del dataset).

Poi selezioniamo la funzione “Analisi delle componenti principali”, nel Menu “Statistiche / Analisi dimensionale”.

&nolink |


Infine selezioniamo le opzioni dalla finestra di dialogo:

&nolink |


Rcmdr>  .PC <- princomp(~Assault+Murder+Rape+UrbanPop, cor=TRUE, data=USArrests)

L'output automatico di RCommander è composto da:

- La matrice dei component loadings

Rcmdr>  unclass(loadings(.PC))  # component loadings
Comp.1     Comp.2     Comp.3      Comp.4
Assault  -0.5831836  0.1879856  0.2681484  0.74340748
Murder   -0.5358995  0.4181809  0.3412327 -0.64922780
Rape     -0.5434321 -0.1673186 -0.8177779 -0.08902432
UrbanPop -0.2781909 -0.8728062  0.3780158 -0.13387773

- La varianza delle componenti, ovvero la deviazione standard al quadrato (sd^2)

Rcmdr>  .PC$sd^2  # component variances
Comp.1    Comp.2    Comp.3    Comp.4 
2.4802416 0.9897652 0.3565632 0.1734301

- La funzione summary(), che corrisponde a quella vista sopra

Rcmdr>  summary(.PC) # proportions of variance
Importance of components:
  Comp.1    Comp.2    Comp.3     Comp.4
Standard deviation     1.5748783 0.9948694 0.5971291 0.41644938
Proportion of Variance 0.6200604 0.2474413 0.0891408 0.04335752
Cumulative Proportion  0.6200604 0.8675017 0.9566425 1.00000000

- E infine, se lo si è selezionato, il grafico delle componenti

Rcmdr>  screeplot(.PC)
&nolink |


Attenzione: RCommander elimina l'oggetto contenente gli output alla fine della procedura:

Rcmdr>  remove(.PC)

La funzione PCA in FactoMineR

L'analisi in componenti principali in FactoMineR adotta un approccio ancora diverso. Come in princomp, le componenti vengono estratte dalla matrice delle correlazioni; diversamente che negli altri casi, però, l'analisi delle variabili e dei casi, nonché l'interpretazione dei fattori, viene condotta in base all'associazione variabili-componenti o individui-componenti.

In questo modo, è possibile rappresentare sui piani fattoriali le variabili eventualmente escluse dall'analisi (supplementari, o illustrative), anche categoriali, e non solo quantitative. Le misure di associazione e di qualità della rappresentazione utilizzate sono infatti — oltre alle correlazioni fra variabili attive e componenti — i coseni quadri e i valori test, che possono essere calcolati anche per le modalità delle variabili categoriali.

FactoMineR ha un plugin per RCommander, che fornisce una interfaccia a menu per l'esecuzione dei comandi.

Lo script utilizzato in questo esempio è disponibile all'indirizzo http://factominer.free.fr/docs/code_pca.r

In primo luogo, carichiamo la libreria necessaria, e il dataset di esempio:

> library(FactoMineR)
> data(decathlon)

Eseguiamo poi il comando, indicando le variabili supplementari quantitative (quanti.sup) e quelle qualitative (quali.sup). Le altre variabili contenute nel dataset saranno automaticamente utilizzate come attive — ovvero per il calcolo delle componenti. Questo significa che il file deve contenere tutte e solo le variabili necessarie all'analisi.

> res <- PCA(decathlon, quanti.sup = 11:12, quali.sup = 13)

L'output produce anche il grafico dei casi

&nolink |Fig. 1

che corrisponde al comamdo plot(res, choix=“ind”, axes=1:2)

e quello delle variabili

&nolink |

che corrisponde al comamdo plot(res, choix=“var”, axes=1:2)

Possiamo omettere i grafici con l'opzione graph = FALSE.

Gli oggetti contenuti in res sono:

> res
**Results for the Principal Component Analysis (PCA)**
The analysis was performed on 41 individuals, described by 13 variables
*The results are available in the following objects:

   name                description                                              
1  "$eig"              "eigenvalues"                                            
2  "$var"              "results for the variables"                              
3  "$var$coord"        "coord. for the variables"                               
4  "$var$cor"          "correlations variables - dimensions"                    
5  "$var$cos2"         "cos2 for the variables"                                 
6  "$var$contrib"      "contributions of the variables"                         
7  "$ind"              "results for the individuals"                            
8  "$ind$coord"        "coord. for the individuals"                             
9  "$ind$cos2"         "cos2 for the individuals"                               
10 "$ind$contrib"      "contributions of the individuals"                       
11 "$quanti.sup"       "results for the supplementary quantitative variables"   
12 "$quanti.sup$coord" "coord. for the supplementary quantitative variables"    
13 "$quanti.sup$cor"   "correlations suppl. quantitative variables - dimensions"
14 "$quali.sup"        "results for the supplementary categorical variables"    
15 "$quali.sup$coord"  "coord. for the supplementary categories"                
16 "$quali.sup$v.test" "v-test of the supplementary categories"                 
17 "$call"             "summary statistics"                                     
18 "$call$centre"      "mean of the variables"                                  
19 "$call$ecart.type"  "standard error of the variables"                        
20 "$call$row.w"       "weights for the individuals"                            
21 "$call$col.w"       "weights for the variables"  

Anche in questo caso, possiamo utilizzare i singoli oggetti contenuti in res per approfondire l'analisi.

Autovalori

Autovalori e varianza spiegata (% e %cum):

> res$eig
        eigenvalue percentage of variance
comp 1   3.2719055              32.719055
comp 2   1.7371310              17.371310
comp 3   1.4049167              14.049167
comp 4   1.0568504              10.568504
comp 5   0.6847735               6.847735
comp 6   0.5992687               5.992687
comp 7   0.4512353               4.512353
comp 8   0.3968766               3.968766
comp 9   0.2148149               2.148149
comp 10  0.1822275               1.822275
        cumulative percentage of variance
comp 1                           32.71906
comp 2                           50.09037
comp 3                           64.13953
comp 4                           74.70804
comp 5                           81.55577
comp 6                           87.54846
comp 7                           92.06081
comp 8                           96.02958
comp 9                           98.17773
comp 10                         100.00000

Grafico degli autovalori:

> barplot(res$eig[, 1])
&nolink |

Variabili

Medie e deviazioni standard (scarto tipo) delle variabili attive:

> round(cbind(res$call$centre, res$call$ecart.type), 
+     2)
              [,1]  [,2]
100m         11.00  0.26
Long.jump     7.26  0.31
Shot.put     14.48  0.81
High.jump     1.98  0.09
400m         49.62  1.14
110m.hurdle  14.61  0.47
Discus       44.33  3.34
Pole.vault    4.76  0.27
Javeline     58.32  4.77
1500m       279.02 11.53

(cbind serve ad unire i dataframe, round serve ad arrotondare i decimali)

Coordinate fattoriali delle prime quattro componenti:

> res$var$coord[, 1:4]
                  Dim.1      Dim.2       Dim.3       Dim.4
100m        -0.77471983  0.1871420 -0.18440714 -0.03781826
Long.jump    0.74189974 -0.3454213  0.18221105  0.10178564
Shot.put     0.62250255  0.5983033 -0.02337844  0.19059161
High.jump    0.57194530  0.3502936 -0.25951193 -0.13559420
400m        -0.67960994  0.5694378  0.13146970  0.02930198
110m.hurdle -0.74624532  0.2287933 -0.09263738  0.29083103
Discus       0.55246652  0.6063134  0.04295225 -0.25967143
Pole.vault   0.05034151 -0.1803569  0.69175665  0.55153397
Javeline     0.27711085  0.3169891 -0.38965541  0.71227728
1500m       -0.05807706  0.4742238  0.78214280 -0.16108904

Correlazioni delle variabili con le prime quattro componenti:

> res$var$cor[, 1:4]
                  Dim.1      Dim.2       Dim.3       Dim.4
100m        -0.77471983  0.1871420 -0.18440714 -0.03781826
Long.jump    0.74189974 -0.3454213  0.18221105  0.10178564
Shot.put     0.62250255  0.5983033 -0.02337844  0.19059161
High.jump    0.57194530  0.3502936 -0.25951193 -0.13559420
400m        -0.67960994  0.5694378  0.13146970  0.02930198
110m.hurdle -0.74624532  0.2287933 -0.09263738  0.29083103
Discus       0.55246652  0.6063134  0.04295225 -0.25967143
Pole.vault   0.05034151 -0.1803569  0.69175665  0.55153397
Javeline     0.27711085  0.3169891 -0.38965541  0.71227728
1500m       -0.05807706  0.4742238  0.78214280 -0.16108904

Casi

Coordinate, contributi e coseni quadri (vedi: Analisi delle corrispondenze):

> round(cbind(res$ind$coord[, 1:4], res$ind$contrib[, 
+     1:4], res$ind$cos2[, 1:4]), 2)
            Dim.1 Dim.2 Dim.3 Dim.4 Dim.1 Dim.2 Dim.3 Dim.4
SEBRLE       0.79  0.77  0.83  1.17  0.47  0.84  1.19  3.18
CLAY         1.23  0.57  2.14 -0.35  1.14  0.46  7.96  0.29
KARPOV       1.36  0.48  1.96 -1.86  1.38  0.33  6.64  7.95
BERNARD     -0.61 -0.87  0.89  2.22  0.28  1.07  1.37 11.38
YURKOV      -0.59  2.13 -1.23  0.87  0.26  6.38  2.61  1.76
....
            Dim.1 Dim.2 Dim.3 Dim.4
SEBRLE       0.11  0.11  0.12  0.25
CLAY         0.12  0.03  0.37  0.01
KARPOV       0.16  0.02  0.33  0.30
BERNARD      0.05  0.10  0.10  0.65
YURKOV       0.04  0.50  0.16  0.08
....

Caratterizzazione degli assi (prime 3 componenti):

> dimdesc(res)
$Dim.1
$Dim.1$quanti
            correlation      p.value
Points        0.9561543 0.000000e+00
Long.jump     0.7418997 2.849886e-08
Shot.put      0.6225026 1.388321e-05
High.jump     0.5719453 9.362285e-05
Discus        0.5524665 1.802220e-04
Rank         -0.6705104 1.616348e-06
400m         -0.6796099 1.028175e-06
110m.hurdle  -0.7462453 2.136962e-08
100m         -0.7747198 2.778467e-09


$Dim.2
$Dim.2$quanti
          correlation      p.value
Discus      0.6063134 2.650745e-05
Shot.put    0.5983033 3.603567e-05
400m        0.5694378 1.020941e-04
1500m       0.4742238 1.734405e-03
High.jump   0.3502936 2.475025e-02
Javeline    0.3169891 4.344974e-02
Long.jump  -0.3454213 2.696969e-02


$Dim.3
$Dim.3$quanti
           correlation      p.value
1500m        0.7821428 1.554450e-09
Pole.vault   0.6917567 5.480172e-07
Javeline    -0.3896554 1.179331e-02

Infine, FactoMineR include una funzione che consente di esportare l'intero output in un file csv:

> write.infile(res,file="my_FactoMineR_results.csv")

Script

FM-PCA-ex.R
data(decathlon)
res <- PCA(decathlon,quanti.sup=11:12,quali.sup=13)
res
res$eig
barplot(res$eig[,1])
round(cbind(res$call$centre,res$call$ecart.type),2)
res$var$coord[,1:4]
res$var$cor[,1:4]
round(cbind(res$ind$coord[,1:4],res$ind$contrib[,1:4],res$ind$cos2[,1:4]),2)
dimdesc(res)
write.infile(res,file="my_FactoMineR_results.csv")
r/analisi_esplorativa/analisi_in_componenti_principali.txt · Ultima modifica: 04/02/2021 10:16 da admin