## Esercizio: analisi delle componenti principali

## "Ratios" example dal Biggeri: Caso 8.2

Tab8.15 <- read.table("Tab8.15.txt", sep="\t", h=T, dec=".")

Ratios <- Tab8.15[,-1]    # via l'indice della prima colonna
dimnames(Ratios)[[1]] <- c("Alfa", 2:14)

## ispezioniamo la struttura di correlazione tra le variabili
## (Tab. 8.16 p. 391)
cor(as.matrix(Ratios))

## analisi delle componenti principali in R base
## (NB: le variabili vengono scalate e centrate, ovvero
##  standardizzate; altrimenti "prevarrebbe" CR che ha più varianza)
ratios.pca <- prcomp(Ratios,
                     center = TRUE,
                     scale. = TRUE)

## varianza delle CP (Tab. 8.17, NB: nel sommario "standard" la
## prima riga è SD non varianza)
summary(ratios.pca)
## le varianze:
ratios.pca$sdev^2

## e i coefficienti delle 4 CP sono:
ratios.pca$rotation

## Riproduciamo la Fig. 8.11 (b): imprese nello spazio
## delle prime due CP

## scaliamo i dati, sottraendo la media da ogni variabile
Ratios.dem  <- apply(Ratios, 2, FUN=function(x) x-mean(x))

## moltiplichiamo (in senso matriciale, righe per colonne)
## per i coefficienti, ottenendo gli scores sulle 4 CP
ratios.scores <- as.matrix(Ratios.dem) %*% ratios.pca$rotation

## disegniamo la figura, evidenziando l'impresa Alfa
plot(ratios.scores[,"PC1"], ratios.scores[,"PC2"], pch=19, xlim=c(-0.3,0.3), ylim=c(-1.2,1.2))
abline(h=0, col=grey(0.8))
abline(v=0, col=grey(0.8))
points(ratios.scores["Alfa","PC1"], ratios.scores["Alfa","PC2"],
       pch=19, col="red")
text(ratios.scores["Alfa","PC1"] + 0.03, ratios.scores["Alfa","PC2"], "Alfa")


## Advanced PCA with FactoMineR:

## riproduciamo la Fig. 8.11 (a) sui dati di Tab. 8.15 in modo
## automatico. usando le funzionalità di questo package 

library(FactoMineR)
r.pca <- PCA(Ratios)

## E' piuttosto chiaro come -- nel piano definito dalle prime due CP 
## -- ROA e ROE siano estremamente simili, positivamente correlate 
## con CP1 e "poco" con CP2; ROD e CR opposte, molto correlate con 
## CP2 e poco con CP1; si veda il grafico "Variables Factor Map"  

## Per ispezionare i risultati
library(factoextra)

## Scores di ogni unità di osservazione (impresa) sulle 4 CP:
get_pca_ind(r.pca)$coord

## Coefficienti (a) delle 4 CP 
get_pca_var(r.pca)$coord