# STATISTICA MULTIVARIATA
# I SEMESTRE 2019-2020
# DOCENTE: MINI V.
##############################################
####  analisi per componenti principali   ####
##############################################

##punti fondamentali dell'analisi da seguire##

1. motivazione della analisi (riduzione, ranking, multicollinearità ecc.)

2. struttura del database (ed eventuale esclusione di alcuni dati)
getwd()
voti=read.csv2("students.csv",header=T)
View(voti)
str(voti)
head(voti)
attach(voti)

esclusione della prima colonna
voti=voti[,2:4]


3. visualizzazione grafica
plot(voti)
pairs(x=voti,panel=panel.smooth)

4. matrice delle correlazioni
Creiamo una matrice di correlazione con solo due decimali, in modo da agevolare la lettura
res=cor(voti)
round(res,2)

5. analisi PER COMPONENTI PRINCIPALI: comando "prcomp"
La funzione prcomp comprende vari argomenti. 
L'analisi per componenti principali può essere eseguita su variabili originarie. 

pca=prcomp(voti)

Tuttavia, generalmente si richiede di 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.

pca1=prcomp(voti, scale=TRUE) --> standardizzazione delle variabili

RISULTATO OTTENUTO: 
si visualizza attraverso la dicitura del nome che si è attribuito all'analisi: 

pca1

si hanno tante componenti quante variabili originarie. Ogni componente è funzione lineare delle 
variabili originarie ciascuna moltiplicata per il proprio peso fattoriale. 

Le informazioni rilevanti si ricavano dalla funzione summary dell’output di prcomp. 

summary(pca1)

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.

 
6. quante componenti estrarre? (3 metodi)

(1) METODO: autovalori maggiori di 1 (o regola di Kaiser)
si chiede 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).

e_values=pca1$sdev^2
e_values

barplot(e_values)


(2) METODO: varianza totale spiegata superiore al 70%
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.

# 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? 

summary(pca1)

(3) METODO: visualizzazione grafica del "grafico a gomito"
Si 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.
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.


screeplot(pca1,type=c("lines"))

oppure

plot(pca1,type="lines")

grafico delle varianze: plot(pca1)


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.

7. DIAMO IL NOME ALLE COMPONENTI

pca1 # e guardiamo il factor loadings

8. correlazione tra le variabili originarie e le componenti

cor.pca1=cor(voti,pca1$scores)
cor.pca1

9. identificare i valori (scores) finali --> nuove variabili create per il database

yscores=pca1$x

se vogliamo tenere solo i primi 2:
final=yscores[,1:2]


10. rotazione: si migliora l'interpretazione?

ortogonale: 
varimax(pca1$x,normalize=TRUE,eps=1e-5)

obliqua:
promax(pca1$x,m=4)

si possono confrontare i valori 

11. Grafico delle due principali componenti  
biplot(pca1)