## Analiza glavnih komponenti

# Kod višestruke regresije u opštem slučaju prediktori su korelisani. Analiza
# glavnih komponenti (Principal Component Analysis) koristi se kao alat u
# ispitivanju podataka i pri pravljenju modela; u pretprocesiranju, ispitivanju
# skupa opservacija i skupa prediktora. Ovom statističkom procedurom se,
# korišćenjem linearnih transformacija, početni skup prediktora transformiše u
# novi skup nekorelisanih prediktora. Pritom se identifikuju oni prediktori koji
# su najbitniji za podatke, koji opisuju najviše varijabilnosti, i oni koji nisu
# toliko bitni. Te nebitne prediktore onda možemo izbaciti i sprovesti linearnu
# regresiju sa transformisanim prediktorima i sa manje njih. Na taj način se PCA
# koristi za smanjivanje dimenzionalnosti i povećanje interpretabilnosti.

# Generisemo neke podatke
set.seed(216)
x1 <- rnorm(100)
x2 <- 1 + x1 + rnorm(100, sd = 0.5)

# Vidimo da je najveca varijabilnost podataka duz prave y = 1 + x
par(pty="s") # hocemo kvadratni grafik
plot(x1, x2, xlim = c(-3, 3), ylim = c(-2, 4))

arrows(c(mean(x1) - 1, mean(x1) + 1), c(mean(x2) - 1, mean(x2) - 1),
       c(mean(x1) + 1, mean(x1) - 1), c(mean(x2) + 1, mean(x2) + 1),
       col="blue", lwd = 2)

# Osnovna ideja je da se rotira skup prediktora tako da se pravci najvece
# varijabilnosti nalaze na osama i da prediktori postanu ortogonalni.

# To dobijamo uzimanjem Z = XU, gde je U matrica sacinjena od sopstvenih vektora
# matrice X^T * X. Tako ce Z^T * Z biti dijagonalna matrica sa sopstvenim
# vrednostima od X^T * X na dijagonali.
# Pravimo matriu X
X <- cbind(x1 = x1, x2 = x2)
# Obavezno je centriranje podataka da bi se dobili ortogonalni prediktori
X <- scale(X)
# Racunamo X^T * X
XtX <- t(X) %*% X
# Trazimo sopstvene vrednosti i sopstvene vektore
e <- eigen(XtX)
colnames(e$vectors) <- c("EV1", "EV2")
rownames(e$vectors) <- c("x1", "x2")

# Transformisemo X tako da dobijemo matricu Z sacinjenu od novih prediktora,
# koji su ortogonalni.
Z <- X %*% e$vectors
plot(Z)
t(Z) %*% Z
cor(Z)

# Jedna od najbitnijih upotreba za PCA je smanjivanje dimenzije. Naime, ako
# imamo mnogo prediktora, korisno je smanjiti broj prediktora sa kojima se radi
# analiza, radi interpretabilnosti. medjutim, nekad se ne mogu prosto izbaciti
# prediktori iz modela, vec neke linearne kombinacije istih uticu na model.
# Analiza glavnih komponenti tezi da nadje ove linearne kombinacije prediktora
# iz X koje najvise doprinose varijabilnosti podataka.

# Uradicemo analizu glavnih komponenti na bazi longley
X <- longley[, -7] # kolona 7 je Employed, koja je zavisna promenljiva

# Pre analize glavnih komponenti obavezno je standardizovati prediktore! 
Xs <- scale(X)
# Ako se ne standardizuju, problem koji se moze javiti je da prediktor koji je u
# intervalu [-100, 100] izgleda da doprinosi varijabilnosti vise nego neki koji
# je na [-0.1, 0.1], iako su ravnopravni kad se skaliraju.

# Jedan nacin za odredjivanje glavnih komponenti je koriscenje prcomp funkcije
pca <- prcomp(X, scale = TRUE)
par(pty="m") # necemo vise kvadratni grafik

# Veoma je koristan grafik "screeplot" kod PCA, koji pokazuje udeo objasnjene
# varijabilnosti svake od glavnih komponenti. Okvirno bi trebalo uzeti onoliko
# prvih glavnih komponenti koje objasnjavaju 80% varijabilnosti. U ovom slucaju
# ce to biti prve dve.
# Nekoliko nacina da se nacrta ovaj grafik su sledeci...
plot(pca, type="l")
screeplot(pca)
# Dok je najlepsi iz paketa factoextra
library(factoextra)
fviz_eig(pca)
# Pored screeplota, numericki mozemo da vidimo udeo disperzije svake od
# komponenata koristeci funkciju iz factoextra paketa...
get_eigenvalue(pca)
# ili prosto summary
summary(pca)

# Kada smo odredili glavne komponente, treba transformisati prediktore u novu
# bazu odredjenu glavnim komponentama. Rucni metod je mnozenje matricom rotacije
pca_data <- Xs %*% pca$rotation
plot(pca_data)
cor(pca_data)
# Jednostavniji metod je koriscenjem predict funkcije koja radi slicno kao i za
# linearne modele.
pca_data <- predict(pca, longley)
plot(pca_data)
cor(pca_data)

# Uporedicemo modele dobijene sa svim prediktorima i one sa glavnim komponentama
model <- lm(Employed ~ . , longley)
summary(model) # R^2 = 0.9955

# Kada uzmemo svih 6 glavnih komponenti, dobijamo sustinski isti model, samo su
# prediktori u drugoj bazi vektorskog prostora.
model_pca <- lm(longley$Employed ~ pca_data)
summary(model_pca) # R^2 = 0.9955

# ocekivano je da prve glavne komponente budu znacajnije od kasnijih, ali to
# nije uvek slucaj! Ovde vidimo da je 5. glavna komponenta znacajna dok 4. nije.
# napravimo model bez 4..
model_pca2 <- lm(longley$Employed ~ pca_data[, c(1:3,5)])
summary(model_pca2) # R^2 = 0.9938
# dok sa 4. umesto 5. komponente izgleda ovako
model_pca22 <- lm(longley$Employed ~ pca_data[, 1:4])
summary(model_pca22) # R^2 = 0.98618
# Bitna je razlika u modelu, dok ako uzmemo samo prve 3 komponente imamo slican
# model, tj. 4. komponenta zaista ne doprinosi modelu.
model_pca3 <- lm(longley$Employed ~ pca_data[, 1:3])
summary(model_pca3) # R^2 = 0.986

# Ako bismo pratili pravilo o 80% varijabilnosti, uzeli bismo samo prve dve
# glavne komponente, ali R^2 pada sa 0.98 na 0.92, sto nije beznacajno, pa je
# bolje uzeti i trecu komponentu ako je bitnija prediktivna moc od
# interpretacije.
model_pca4 <- lm(longley$Employed ~ pca_data[, 1:2])
summary(model_pca4) # R^2 = 0.9289