PCA e EFA (notebook 3)

Statistica II - 750AA

Introduzione

In questo notebook, esploreremo sia l’Analisi delle Componenti Principali (PCA) che l’Analisi Fattoriale Esplorativa (EFA) su diversi dataset: alcuni generati, altri pre-caricati di R (come iris e mtcars), e un dataset reale (ANSES-CIQUAL food composition table 2020).

# Alcuni pacchetti che useremo: se non presenti installatele con il comando install.packages()
library(readxl)
library(tidyr, stringr)
library(ggplot2, ggcorrplot)
library(corrplot)
library(ggfortify)

Esempi con dataset generati

Come nel notebook precedente, generiamo dei dati partendo da campioni gaussiani di medie e deviazioni standard specificate.

for( i in 2:d ){
  data_gen <- cbind( data_gen,  rnorm(n, mean=0, sd=1))
  names_data_gen <- c(names_data_gen, print(paste0("feature_", as.character(i))))
}

[1] "feature_2"
[1] "feature_3"
[1] "feature_4"
[1] "feature_5"

In questo dataset la varianza campionaria dovrebbe essere all’incirca diagonale con degli \(1\) sulla diagonale – perché le colonne sono state generate in modo indipendente e abbiamo specificato la deviazione standard teorica sd=1.

I comandi per la varianza, covarianza e correlazione sono rispettivamente var(), cov() e cor().

var_data_gen

            feature_1   feature_2   feature_3  feature_4   feature_5
feature_1  0.90319778 -0.04827194  0.10746503 0.06946613 -0.21129440
feature_2 -0.04827194  1.03936852 -0.06219997 0.14722352 -0.13245614
feature_3  0.10746503 -0.06219997  1.01835020 0.03320205 -0.25230152
feature_4  0.06946613  0.14722352  0.03320205 0.83073817  0.05959285
feature_5 -0.21129440 -0.13245614 -0.25230152 0.05959285  1.15805656

cor_data_gen

            feature_1   feature_2   feature_3  feature_4   feature_5
feature_1  1.00000000 -0.04982170  0.11205394 0.08019543 -0.20660064
feature_2 -0.04982170  1.00000000 -0.06045843 0.15843836 -0.12073201
feature_3  0.11205394 -0.06045843  1.00000000 0.03609810 -0.23233060
feature_4  0.08019543  0.15843836  0.03609810 1.00000000  0.06075712
feature_5 -0.20660064 -0.12073201 -0.23233060 0.06075712  1.00000000

Con cov() e cor() è possibile anche calcolare covarianze tra coppie di data frames.

cor(data_gen[-3], data_gen[3])

            feature_3
feature_1  0.11205394
feature_2 -0.06045843
feature_4  0.03609810
feature_5 -0.23233060

Per visualizzare le correlazioni mediante una mappa di calore possiamo usare il pacchetto corrplot oppure ggcorrplot (che estende ggplot2).

# installare corrplot se non è già presente
library(corrplot)

# plot di base
corrplot(cor_data_gen)

Ci sono moltissime opzioni per migliorare il grafico (leggere l’help per informazioni).

corrplot(cor_data_gen, method="pie", type="lower", diag=FALSE, number.cex=1, number.digits = 2, title="Dataset generato (features indipendenti)")

Similmente, usando ggcorrplot si producono grafici con ggplot.

Dataset con correlazioni

Le matrici trovate di correlazione sono poco interessanti, perché il dataset non presenta associazioni interessanti tra le varie features. Per introdurre artificialmente queste correlazioni, trasformiamo le features con una mappa opportuna.

In questo caso il correlogramma è più interessante

Dalla teoria sappiamo che \(\operatorname{var}(xA) = A^T \operatorname{var}(x) A \approx A^T A\) essendo \(\operatorname{var}(x)\approx I_{d\times d}\). Verifichiamolo numericamente. Calcoliamo la varianza campionaria.

var(data_gen_A)

           feature_1  feature_2  feature_3  feature_4 feature_5
feature_1  2.3673428  0.7354934  1.6560901 -1.0357747 -1.973163
feature_2  0.7354934  5.3050634 -0.3901245 -1.5062156  1.309246
feature_3  1.6560901 -0.3901245  2.7275217 -0.8481627 -1.888022
feature_4 -1.0357747 -1.5062156 -0.8481627  1.9154925  1.637760
feature_5 -1.9731634  1.3092457 -1.8880223  1.6377603  3.822365

Per calcolare la trasposta di una matrice usiamo la funzione t().

cov_A_teorica <- t(A) %*% A
cov_A_teorica

     [,1] [,2] [,3] [,4] [,5]
[1,]    3    1    2   -1   -2
[2,]    1    4    0   -1    1
[3,]    2    0    3   -1   -2
[4,]   -1   -1   -1    2    2
[5,]   -2    1   -2    2    4

PCA sul Dataset Generato

Per eseguire l’analisi delle componenti principali, possiamo usare il comando prcomp() (consigliato) oppure il comando princomp(). Si veda https://cran.r-project.org/web/packages/LearnPCA/vignettes/Vig_07_Functions_PCA.pdf per una discussione delle differenze.

È fortemente consigliato per la PCA di standardizzare i dati. Possiamo farlo prima con il comando scale()

# usiamo direttamente il dataset generato con correlazioni indotte da A
data_gen <- data_gen_A

# standardizziamo
data_gen_standard <- scale(data_gen, center=TRUE, scale=TRUE)

# Notiamo che la media è nulla per ciascuna feature
summary(data_gen_standard)

   feature_1          feature_2         feature_3      
 Min.   :-2.19265   Min.   :-2.4748   Min.   :-2.6179  
 1st Qu.:-0.74504   1st Qu.:-0.6652   1st Qu.:-0.8311  
 Median : 0.02771   Median : 0.0973   Median : 0.1511  
 Mean   : 0.00000   Mean   : 0.0000   Mean   : 0.0000  
 3rd Qu.: 0.71495   3rd Qu.: 0.6608   3rd Qu.: 0.7401  
 Max.   : 2.48814   Max.   : 2.1326   Max.   : 2.2966  
   feature_4         feature_5       
 Min.   :-2.1960   Min.   :-2.24490  
 1st Qu.:-0.7081   1st Qu.:-0.69485  
 Median :-0.1085   Median :-0.05175  
 Mean   : 0.0000   Mean   : 0.00000  
 3rd Qu.: 0.6783   3rd Qu.: 0.64878  
 Max.   : 2.5069   Max.   : 2.19352  

# Notiamo che la varianza coincide con la correlazione
var(data_gen_standard)

           feature_1  feature_2  feature_3  feature_4  feature_5
feature_1  1.0000000  0.2075405  0.6517324 -0.4864011 -0.6559434
feature_2  0.2075405  1.0000000 -0.1025590 -0.4724998  0.2907433
feature_3  0.6517324 -0.1025590  1.0000000 -0.3710691 -0.5847321
feature_4 -0.4864011 -0.4724998 -0.3710691  1.0000000  0.6052628
feature_5 -0.6559434  0.2907433 -0.5847321  0.6052628  1.0000000

cor(data_gen_standard)

           feature_1  feature_2  feature_3  feature_4  feature_5
feature_1  1.0000000  0.2075405  0.6517324 -0.4864011 -0.6559434
feature_2  0.2075405  1.0000000 -0.1025590 -0.4724998  0.2907433
feature_3  0.6517324 -0.1025590  1.0000000 -0.3710691 -0.5847321
feature_4 -0.4864011 -0.4724998 -0.3710691  1.0000000  0.6052628
feature_5 -0.6559434  0.2907433 -0.5847321  0.6052628  1.0000000

A questo punto possiamo usare il comando prcomp().

pca_gen <- prcomp(data_gen_standard)
summary(pca_gen)

Importance of components:
                          PC1    PC2     PC3     PC4
Standard deviation     1.6071 1.0881 0.35344 0.32929
Proportion of Variance 0.6457 0.2960 0.03123 0.02711
Cumulative Proportion  0.6457 0.9417 0.97289 1.00000

Oppure è equivalente usare l’opzione scale.TRUE sul data frame originale.

summary(pca_gen)

Importance of components:
                          PC1    PC2    PC3     PC4     PC5
Standard deviation     1.6421 1.1700 0.7642 0.57752 0.13055
Proportion of Variance 0.5393 0.2738 0.1168 0.06671 0.00341
Cumulative Proportion  0.5393 0.8131 0.9299 0.99659 1.00000

Leggiamo dall’output di summary() la deviazione standard (\(\sqrt{\lambda_j}\)) per ciascuna componente principale, oltre alla proporzione di varianza spiegata e la proporzione cumulativa. Possiamo visualizzarli tramite dei plot. Cominciamo con lo scree plot.

variances_pca_gen <- (pca_gen$sdev)**2

# uno scree plot con il comando plot
plot(variances_pca_gen, type='l')
points(variances_pca_gen)

# uno scree plot con ggplot2

scree_data <- data.frame(PC = 1:d, Variance = variances_pca_gen)
ggplot(scree_data, aes(x = PC, y = Variance)) +
  geom_line() +
  geom_point() +
  labs(title = "Scree Plot: dataset generato", x = "Componenti Principali", y = "Autovalori")+
  scale_x_continuous(label = ~ scales::comma(.x, accuracy = 1))

Presentiamo il plot della varianza cumulata.

# calcoliamo la percentuale di varianza cumulativa

cum_var <- data.frame(PC = 1:d, Variance = cumsum(variances_pca_gen)/sum(variances_pca_gen))

# plot di base
plot(cum_var, type='l')
points(cum_var)

abline(h=0.95, col="red", lty=3)

# stesso plot ma con ggplot2

ggplot(cum_var, aes(x = PC, y = Variance)) +
  geom_line() +
  geom_point() +
  labs(title = "Varianza cumulata: dataset generato (indipendenti)", x = "Componenti Principali", y = "percentuale di varianza")+
  scale_x_continuous(label = ~ scales::comma(.x, accuracy = 1)) +
  geom_hline(yintercept = 0.95, linetype="dashed",  color = "red", size=1)

Presentiamo infine il biplot in cui sia gli scores che i loadings vengono rappresentati sullo stesso piano. Il comando di base è biplot(). Per grafici più personalizzabili e chiari consigliamo di usare autoplot() dal pacchetto ggfortify.

Osserviamo la matrice delle rotazioni (detta anche dei loadings).

pca_gen$rotation

                  PC1         PC2        PC3          PC4         PC5
feature_1 -0.52764223 -0.04329013 -0.4143550 -0.658173590  0.33887623
feature_2 -0.07517918 -0.83412884 -0.2190672 -0.004680497 -0.50056451
feature_3 -0.47915887  0.23391487 -0.5286863  0.653921203 -0.09256595
feature_4  0.46305473  0.39627489 -0.5409181 -0.309175053 -0.49027027
feature_5  0.52146798 -0.30100678 -0.4563090  0.208766530  0.62101906

EFA sul Dataset Generato

Passiamo ora all’Analisi Fattoriale Esplorativa (EFA). Anche qui abbiamo diversi comandi: factanal() del pacchetto stats (di base) e fa() del pacchetto psych, più avanzato.

# cerchiamo un solo fattore (che potrebbe essere vicino alla PC1)

efa_gen <- factanal(data_gen, factors=2)

# library(psych)
# efa_gen<- fa(data_gen, fators=3)

summary(efa_gen)

             Length Class    Mode     
converged     1     -none-   logical  
loadings     10     loadings numeric  
uniquenesses  5     -none-   numeric  
correlation  25     -none-   numeric  
criteria      3     -none-   numeric  
factors       1     -none-   numeric  
dof           1     -none-   numeric  
method        1     -none-   character
rotmat        4     -none-   numeric  
STATISTIC     1     -none-   numeric  
PVAL          1     -none-   numeric  
n.obs         1     -none-   numeric  
call          3     -none-   call     

Controlliamo in che misura l’equazione matriciale della EFA: \(\operatorname{var}(x) = L^T L + \Psi\) è numericamente valida.

cor(data_gen)

           feature_1  feature_2  feature_3  feature_4  feature_5
feature_1  1.0000000  0.2075405  0.6517324 -0.4864011 -0.6559434
feature_2  0.2075405  1.0000000 -0.1025590 -0.4724998  0.2907433
feature_3  0.6517324 -0.1025590  1.0000000 -0.3710691 -0.5847321
feature_4 -0.4864011 -0.4724998 -0.3710691  1.0000000  0.6052628
feature_5 -0.6559434  0.2907433 -0.5847321  0.6052628  1.0000000

 efa_gen$loadings[1:d] %*% t(efa_gen$loadings[1:d]) +diag(efa_gen$uniquenesses, nrow=d, ncol=d)

            [,1]        [,2]        [,3]        [,4]        [,5]
[1,]  1.00000022  0.01177468  0.59398615 -0.51862035 -0.67765244
[2,]  0.01177468  0.99999347  0.01046037 -0.00913314 -0.01193377
[3,]  0.59398615  0.01046037  0.99999970 -0.46073068 -0.60201122
[4,] -0.51862035 -0.00913314 -0.46073068  0.99999751  0.52562719
[5,] -0.67765244 -0.01193377 -0.60201122  0.52562719  0.99999987

Confrontiamo il vettore dei loadings con la prima componente principale riscalata per la sua deviazione standard.

loading_PCA <- pca_gen$sdev[1]* pca_gen$rotation[,1]

print("Loadings PCA: ")

[1] "Loadings PCA: "

loading_PCA

 feature_1  feature_2  feature_3  feature_4  feature_5 
-0.8664501 -0.1234530 -0.7868348  0.7603899  0.8563113 

print("Loadings EFA: ")

[1] "Loadings EFA: "

efa_gen$loadings[1:d]

[1]  0.81769128  0.01439991  0.72641859 -0.63424957 -0.82873873

Esercizio: effettuare una PCA su un dataset non standardizzato, verificare che le PC tendono a catturare le componenti che hanno varianza maggiore (invece di rilevare le correlazioni).

Esempi con dataset precaricati

Cominciamo con il dataset iris.

head(iris)

summary(iris)

  Sepal.Length    Sepal.Width     Petal.Length    Petal.Width   
 Min.   :4.300   Min.   :2.000   Min.   :1.000   Min.   :0.100  
 1st Qu.:5.100   1st Qu.:2.800   1st Qu.:1.600   1st Qu.:0.300  
 Median :5.800   Median :3.000   Median :4.350   Median :1.300  
 Mean   :5.843   Mean   :3.057   Mean   :3.758   Mean   :1.199  
 3rd Qu.:6.400   3rd Qu.:3.300   3rd Qu.:5.100   3rd Qu.:1.800  
 Max.   :7.900   Max.   :4.400   Max.   :6.900   Max.   :2.500  
       Species  
 setosa    :50  
 versicolor:50  
 virginica :50  
                
                
                

pairs(iris[-5], col=iris[,5])

Studiamo la correlazione.

ggcorrplot(cor(iris[-5]), type="lower")

Vediamo una forte correlazione positiva tra le lunghezze, larghezze dei petali e lunghezza del sepalo, mentre una correlazione debole negativa tra queste e la larghezza del sepalo.

Suggerimento: alcuni grafici comuni (scatterplot, heatmap della correlazione) si possono condensare in uno solo, ed eseguire con comandi diretti del pacchetto GGally.

library(GGally)

ggpairs(iris[-5])

PCA sul dataset Iris

pca_iris <- prcomp(iris[-5], scale. = TRUE)
summary(pca_iris)

Importance of components:
                          PC1    PC2     PC3     PC4
Standard deviation     1.7084 0.9560 0.38309 0.14393
Proportion of Variance 0.7296 0.2285 0.03669 0.00518
Cumulative Proportion  0.7296 0.9581 0.99482 1.00000

Notiamo che le prime due PC spiegano il \(95\%\) della varianza.

scree_data <- data.frame(PC = 1:4, Variance =pca_iris$sdev**2)
ggplot(scree_data, aes(x = PC, y = Variance)) +
  geom_line() +
  geom_point() +
  labs(title = "Scree Plot: dataset iris", x = "Componenti Principali", y = "Autovalori")+
  scale_x_continuous(label = ~ scales::comma(.x, accuracy = 1))

Mostriamo il biplot ed intepretiamo il risultato della PCA.

autoplot(pca_iris, data = iris, colour="Species", x=1, y=2, # selezionare eventualmente altre componenti principali ad esempio x=2, y=3
         label=FALSE, loadings = TRUE, loadings.label = TRUE, loadings.label.colour="blue", loadings.label.repel=TRUE,
         loadings.label.size = 3, 
         main = "Biplot PCA dataset generato")

Il biplot mostra che la PC1 è associata alle tre variabili petali e lunghezza dei sepali, mentre la PC2 è associata alla larghezza dei sepali. Vediamolo attraverso la matrice dei loadings.

round(pca_iris$rotation[, 1:2], 2)

               PC1   PC2
Sepal.Length  0.52 -0.38
Sepal.Width  -0.27 -0.92
Petal.Length  0.58 -0.02
Petal.Width   0.56 -0.07

EFA sul Dataset Iris

Possiamo eseguire la EFA sullo stesso dataset.

library(psych)
efa_iris <- fa(iris[-5], nfactors = 2)

# Riassunto dei risultati: loadings
efa_iris$loadings

Loadings:
             MR1    MR2   
Sepal.Length  0.940  0.184
Sepal.Width          0.978
Petal.Length  0.962 -0.121
Petal.Width   0.925       

                 MR1   MR2
SS loadings    2.664 1.010
Proportion Var 0.666 0.253
Cumulative Var 0.666 0.919

efa_iris$communalities

Sepal.Length  Sepal.Width Petal.Length  Petal.Width 
   0.8099449    0.9682729    0.9950000    0.9022485 

L’analisi rivela la stessa situazione: tre variabili possono essere espresse in termini del primo fattore, la rimanente tramite il secondo fattore.

PCA sul dataset mtcars

Riprendiamo il dataset mtcars e studiamolo mediante PCA e poi EFA.

head(mtcars)

# togliamo le features 8 e 9 perché non numeriche (sono indicatrici)

mtcars_num <- mtcars[-(8:9)]

cor_mtcars <- cor(mtcars_num)

ggcorrplot(cor_mtcars, method="square", type="upper", title="Dataset generato (features dipendenti)", hc.order = TRUE, lab=FALSE)

Riconosciamo due blocchi di variabili correlate in senso opposto: il primo è carb, wt, hp, cyl, disp, qsec, mpg, drat, gear.

Vediamo cosa evidenza la PCA.

pca_mtcars <- prcomp(mtcars_num, scale. = TRUE)

# Riassunto dei risultati
summary(pca_mtcars)

Importance of components:
                          PC1    PC2     PC3     PC4     PC5     PC6
Standard deviation     2.3782 1.4429 0.71008 0.51481 0.42797 0.35184
Proportion of Variance 0.6284 0.2313 0.05602 0.02945 0.02035 0.01375
Cumulative Proportion  0.6284 0.8598 0.91581 0.94525 0.96560 0.97936
                           PC7    PC8     PC9
Standard deviation     0.32413 0.2419 0.14896
Proportion of Variance 0.01167 0.0065 0.00247
Cumulative Proportion  0.99103 0.9975 1.00000

# Scree plot

scree_mtcars <- data.frame(PC = 1:ncol(mtcars_num), Variance =pca_mtcars$sdev**2)

ggplot(scree_mtcars, aes(x = PC, y = Variance)) +
  geom_line() +
  geom_point() +
  labs(title = "Scree Plot: dataset mtcars", x = "Componenti Principali", y = "Autovalori")+
  scale_x_continuous(label = ~ scales::comma(.x, accuracy = 1))

NA
NA

Riconosciamo un gomito a \(k=3\), mentre il criterio di Kaiser darebbe \(k=2\). Il grafico della varianza spiegata cumulata darebbe \(k=4\) (quasi \(k=3\)) per arrivare al \(95\%\).

cum_var <- data.frame(PC = 1:length(pca_mtcars$sdev), Variance = cumsum(pca_mtcars$sdev^2)/sum(pca_mtcars$sdev^2))

ggplot(cum_var, aes(x = PC, y = Variance)) +
  geom_line() +
  geom_point() +
  labs(title = "Varianza cumulata: dataset generato (indipendenti)", x = "Componenti Principali", y = "percentuale di varianza")+
  scale_x_continuous(label = ~ scales::comma(.x, accuracy = 1)) +
  geom_hline(yintercept = 0.95, linetype="dashed",  color = "red", size=1)

Visualizziamo i biplot per le PC1, PC2 e pure PC3.

autoplot(pca_mtcars, data = mtcars, colour="cyl", x=1, y=2, # selezionare eventualmente altre componenti principali ad esempio x=2, y=3
         label=FALSE, loadings = TRUE, loadings.label = TRUE, label.size=3, shape=20,  loadings.label.colour="red", loadings.label.repel=TRUE, loadings.label.size = 4, 
         main = "Biplot PCA dataset mtcars")

NA
NA

Abbiamo quindi una PC 1 in cui i due gruppi di variabili sono ben separati. Visto l’allineamento di mpg e cyl sembrerebbe una feature legata al consumo. La PC 2 invece è più allineata con gear e qsec, quindi più legata alle prestazioni (accelerazione da ferma ecc.).

Possiamo confermare la nostra ipotesi visualizzando il plot degli scores relativi a PC1 e PC2 con i nomi dei modelli.

autoplot(pca_mtcars, data = mtcars, colour="mpg", size="qsec", x=1, y=2,
         label=TRUE, label.repel=TRUE, loadings = FALSE, loadings.label = FALSE, label.size=3, shape=20,  loadings.label.colour="red", loadings.label.repel=TRUE, loadings.label.size = 4, 
         main = "Biplot PCA dataset mtcars")

NA
NA

Consideriamo anche la terza componente principale (plottiamo solo PC2 e PC3). Questa sembra essere anche legata al peso wt e al numero di carburatori carb.

autoplot(pca_mtcars, data = mtcars, colour="cyl", size="wt", x=2, y=3,
         label=FALSE, loadings = TRUE, loadings.label = TRUE, label.size=3, shape=20,  loadings.label.colour="red", loadings.label.repel=TRUE, loadings.label.size = 4, 
         main = "Biplot PCA dataset mtcars")

NA
NA

EFA sul Dataset mtcars

Eseguiamo ora una EFA sullo stesso dataset.

# Eseguiamo l'EFA
efa_mtcars <- factanal(mtcars_num, factors = 2)

# Riassunto dei risultati
print(efa_mtcars)

Call:
factanal(x = mtcars_num, factors = 2)

Uniquenesses:
  mpg   cyl  disp    hp  drat    wt  qsec  gear  carb 
0.169 0.097 0.086 0.079 0.294 0.189 0.338 0.144 0.242 

Loadings:
     Factor1 Factor2
mpg  -0.757  -0.509 
cyl   0.764   0.566 
disp  0.839   0.458 
hp    0.474   0.835 
drat -0.838         
wt    0.837   0.332 
qsec         -0.811 
gear -0.861   0.339 
carb          0.869 

               Factor1 Factor2
SS loadings      4.235   3.127
Proportion Var   0.471   0.347
Cumulative Var   0.471   0.818

Test of the hypothesis that 2 factors are sufficient.
The chi square statistic is 52.32 on 19 degrees of freedom.
The p-value is 5.92e-05 

Un dataset reale

La tabella di composizione alimentare ANSES-CIQUAL, disponibile in formato Excel, è gestito dal CIQUAL all’interno dell’Osservatorio dell’Alimentazione, un’unità di ANSES (l’agenzia francese per la sicurezza alimentare, ambientale e della salute occupazionale). La sua missione è raccogliere, valutare e rendere disponibili i dati sulla composizione nutrizionale degli alimenti consumati in Francia. Questo file contiene la composizione di 3185 alimenti per 67 componenti, come carboidrati, zuccheri, proteine, grassi, vitamine e minerali, fornendo informazioni dettagliate per ogni 100 g della parte commestibile degli alimenti.

È importante notare che, in caso di valori mancanti, viene utilizzato un trattino al posto del numero, e tali valori non dovrebbero essere considerati come zero (cosa che invece faremo per semplicità). Inoltre, il termine traccia viene utilizzato quando un componente è rilevato ma non può essere quantificato con precisione, oppure quando si stima che il contenuto di un alimento sia molto basso, ma non nullo. In questi casi, il contenuto medio viene pubblicato come traccia.

Con queste premesse, carichiamo il file dal formato Excel.

library(readxl)

ciqual <- read_excel("datasets/Table Ciqual 2020_ENG_2020 07 07.xls", 
     col_types = c("skip", "skip", "skip", 
         "text", "text", "skip", "skip", "text", 
         "skip", "text", "text", "text", "text", 
         "text", "text", "text", "text", "text", 
         "text", "text", "text", "text", "text", 
         "text", "text", "text", "text", "text", 
         "text", "text", "text", "text", "text", 
         "text", "skip", "skip", "skip", "skip", 
         "skip", "skip", "skip", "skip", "skip", 
         "skip", "skip", "skip", "skip", "skip",
         "text", "text", "text", "text", "text", 
         "text", "text", "text", "text", "text", 
         "text", "text", "text", "text", "text", 
         "text", "text", "text", "text", "text", 
         "text", "text", "text", "text", "text", 
         "text", "text", "text"))

head(ciqual)

NA

Ripuliamo il dataset rimuovendo caratteri non standard come righe, e i dati mancanti ecc.

ciqual_tidy <- ciqual

ciqual_tidy[ciqual=="-"]<- "0"
ciqual_tidy[] <- lapply(ciqual_tidy, function(x) gsub("< ", "", x))
ciqual_tidy[] <- lapply(ciqual_tidy, function(x) gsub(",", ".", x))
ciqual_tidy[] <- lapply(ciqual_tidy, function(x) gsub("traces", "0", x))
ciqual_tidy$alim_nom_eng <- ciqual$alim_nom_eng

ciqual_tidy[-(1:3)] <- lapply(ciqual_tidy[-(1:3)], function(x) as.numeric(x))
ciqual_tidy <- drop_na(ciqual_tidy)
names_ciqual_short <- substr(colnames(ciqual_tidy), 1, 10)

head(ciqual_tidy)

NA

Selezioniamo alcune colonne per effettuare PCA e EFA.

c_select <- c(1, 7:9, 11:14, 17, 30, 40, 41)

# oppure selezioniamo tutte le colonne eccetto le prime 3

#c_select <- c(-(2:3))
 
ciqual_reduced <- ciqual_tidy[c_select]


colnames(ciqual_reduced)<- names_ciqual_short[c_select]

cor_ciqual <- cor(ciqual_reduced[-1])

ggcorrplot(cor_ciqual, type="upper", hc.order=TRUE, lab=FALSE)

Notiamo una forte correlazione positiva tra carboidrati e zuccheri, tra grassi ed energia, tra sale e sodio, e una correlazione negativa (più o meno forte) tra l’acqua e tutto il resto.

PCA su dataset reale

Passiamo alla PCA.

pca_result <- prcomp(ciqual_reduced[-1], center=TRUE, scale.=TRUE)


# Riassunto dei risultati
summary(pca_result)

Importance of components:
                          PC1    PC2    PC3     PC4     PC5     PC6
Standard deviation     1.5497 1.4509 1.3704 1.04099 1.00218 0.92773
Proportion of Variance 0.2183 0.1914 0.1707 0.09852 0.09131 0.07824
Cumulative Proportion  0.2183 0.4097 0.5804 0.67895 0.77026 0.84850
                           PC7     PC8     PC9    PC10    PC11
Standard deviation     0.85564 0.63403 0.58519 0.38596 0.20239
Proportion of Variance 0.06656 0.03654 0.03113 0.01354 0.00372
Cumulative Proportion  0.91506 0.95160 0.98273 0.99628 1.00000

# Scree plot

scree_ciqual <- data.frame(PC = 1:ncol(ciqual_reduced[-1]), Variance =pca_result$sdev**2)

ggplot(scree_ciqual, aes(x = PC, y = Variance)) +
  geom_line() +
  geom_point() +
  labs(title = "Scree Plot: dataset CIQUAL", x = "Componenti Principali", y = "Autovalori")+
  scale_x_continuous(label = ~ scales::comma(.x, accuracy = 1))

NA
NA
NA
NA

Possiamo selezionare le prime tre componenti principali basandoci sul criterio di Kaiser e scree plot.

Osserviamo ora il biplot.

autoplot(pca_result, data = ciqual_reduced, col="alim_grp_n",
         loadings = TRUE, loadings.label = TRUE, loadings.label.colour="black", loadings.label.repel=TRUE,
         loadings.label.size = 3, 
         main = "Biplot del Dataset ciqual reduced PC1-PC2")

La PC1 sembra essere legata alla presenza o meno di acqua e ai grassi e carboidrati. La PC2 invece distingue piuttosto se l’alimento è proteico/grasso/salato oppure ricco di zuccheri (dolce?). Visualizziamo solamente il plot degli scores senza i loadings.

autoplot(pca_result, data = ciqual_reduced, col="alim_grp_n",
         loadings = FALSE, loadings.label = FALSE, 
         main = "Biplot del Dataset ciqual reduced PC1-PC2")

Quali sono gli alimenti classificati come miscellaneous?

head(ciqual_tidy[ciqual_tidy$alim_grp_nom_eng =="miscellaneous", ])

NA

Sono le salse: quindi PC2 potrebbe indicare un’asse dolce vs salato. Vediamo ora la PC3.

autoplot(pca_result, data = ciqual_reduced, col="alim_grp_n", x=2, y=3,
         loadings = TRUE, loadings.label = TRUE, loadings.label.colour="black", loadings.label.repel=TRUE,
         loadings.label.size = 3, 
         main = "Biplot del Dataset ciqual reduced PC2-PC3")

La PC3 sembra distinguere tra contenuto di sale e grassi/proteine.

Esercizo: In effetti le salse sono dei veri e propri outliers, anche perché stiamo confrontando gli indicatori nutrizionali per 100 g di prodotto. Ripetere tutta l’analisi delle componenti principali rimuovendo le righe corrispondenti al gruppo miscellaneous. I risultati si prestano a interepretazoni diverse?

EFA su dataset reale

Eseguiamo infine l’analisi fattoriale esplorativa.

# Eseguiamo l'EFA
efa_ciqual <- factanal(ciqual_reduced[-1], factors = 3)

# Riassunto dei risultati
print(efa_ciqual)

Call:
factanal(x = ciqual_reduced[-1], factors = 3)

Uniquenesses:
Energy, N  Water (g/1 Protein (g Carbohydra Fat (g/100 Sugars (g/ 
     0.324      0.427      0.880      0.005      0.336      0.509 
fructose ( lactose (g Salt (g/10 Selenium ( Sodium (mg 
     0.921      0.984      0.076      0.996      0.005 

Loadings:
           Factor1 Factor2 Factor3
Energy, N           0.139   0.810 
Water (g/1 -0.147  -0.531  -0.519 
Protein (g         -0.193   0.287 
Carbohydra          0.997         
Fat (g/100                  0.812 
Sugars (g/          0.697         
fructose (          0.250  -0.127 
lactose (g          0.123         
Salt (g/10  0.960                 
Selenium (                        
Sodium (mg  0.997                 

               Factor1 Factor2 Factor3
SS loadings      1.939   1.906   1.690
Proportion Var   0.176   0.173   0.154
Cumulative Var   0.176   0.350   0.503

Test of the hypothesis that 3 factors are sufficient.
The chi square statistic is 1392.81 on 25 degrees of freedom.
The p-value is 1.31e-278 

Troviamo 3 fattori: uno legato al sale/sodio, uno legato ai carboidrati/acqua, e un terzo legato all’energia/grassi/proteine. Ovviamente la varianza cumulata è ancora insufficiente (il test rifiuta l’ipotesi che 3 fattori siano sufficienti). Cosa accade se aumentiamo il numero di fattori latenti?

Esercizi

1. Generare una matrice \(A\) simmetrica e semidefinita positiva con 3 righe e 3 colonne. Implementare il calcolo della radice quarta di \(A\), ovvero di una matrice \(Q\) tale che \(Q^4=A\).

2.Generare 100 coppie di dati Gaussiani che abbiano la bisettrice del primo quadrante come asse principale e varianza 9 lungo tale direzione. Dire quale deve essere l’altro asse principale.

3. Si generi una tabella di 3 colonne e 25 righe. Le prime due colonne siano Gaussiane (media e varianza a piacere), la terza colonna sia la somma delle prime due colonne più un piccolo rumore. Ottenere le varianze spiegate delle componenti principali della tabella, rappresentarle in un grafico illustrativo e visualizzare il grafico delle prime due componenti principali.

4. Generare un vettore Gaussiano di 5000 punti con componenti principali di varianze 9, 4 e 1 ed in modo che i suoi assi principali siano le bisettrici dei quadranti del piano xy e l’asse delle quote. Ottenere il grafico della proiezione dei punti sul piano xz.

5. Generare un campione Gaussiano quadri-dimensionale di numerosità 100, in modo che questo risulti approssimativamente bidimensionale e che la prima componente principale sia la somma delle prime due componenti.

6. Creare una tabella con 5 colonne e 350 righe, e popolarla di valori in modo che il biplot dell’analisi delle componenti principali per la tabella creata presenti alcune frecce molto corte.

7. Generare una tabella di dati in modo che, analizzata attraverso l’analisi delle componenti principali, riveli un accumulo di dati sul piano principale che però si rivela ingannevole se visto attraverso gli altri piani.

8. Considerare la tabella di dati iris. Studiare la corrispondente tabella ottenuta attraverso i logaritmi dei dati originari (ovvero, indagare rispetto a una dipendenza di tipo potenza tra i fattori).

• Valutare l’opportunità di standardizzare la tabella ottenuta.

• Ricavare l’analisi delle componenti principali per la nuova tabella di dati.

• Verificare, colorando i punti rappresentati sui piani principali in accordo con la corrispondente specie, la collocazione dei dati in dipendenza dalla specie.