## Multikolinearnost

# Ako su neki prediktori linearno zavisni, matrica X^T * X ce biti singularna,
# tj. inverz te matrice ne postoji, pa stoga ne mozemo odrediti ocenu metodom
# najmanjih kvadrata na uobicajen nacin. Tacnije, ocena metodom najmanjih
# kvadrata nece biti jedinstvena. Resenje za ovaj problem je da se odstrani
# jedan od prediktora iz pomenute linearne zavistnosti da bismo ostali sa
# linearno nezavisnim prediktorima. Ovaj problem se lako otkriva prosto zbog
# singularnosti matrice X^T * X.

# Podmukliji problem je prisustvo jake linearne korelisanosti nekih prediktora,
# iako nijedan prediktor nije tacno linearna kombinacija ostalih. Tada matrica
# X^T * X moze da bude "skoro singularna", sto dovodi do numerickih
# nestabilnosti i problema sa odredjivanjem tacnih koeficijenata modela metodom
# najmanjih kvadrata.

# Efekte multikolinearnosti smo vidajli tokom kursa vise puta, posebno kada smo
# pravili vestacke primere kod kojih smo namerno uzimali da je jedan
# prediktorjednak drugom uz dodat sum. Dobili smo da su ocene koeficijenata
# pogresne, kao i da testovi ne daju tacne rezultate, pa nekad oznace sve
# prediktore kao beznacajne.

# Sada cemo se baviti nacinima detekcije ovog problema.
# Uzecemo za primer podatke longley i cilj je modelovati promenljivu Employed na
# osnovu ostalih prediktora.

# Napravimo model
model <- lm(Employed ~ ., longley)
summary(model)
# Nekoliko prediktora su oznaceni kao neznacajni, iako bi imalo smisla da uticu
# na zavisnu promenljivu. Provericemo da li postoji kolinearnost medju
# prediktorima koja bi objasnila ovakve rezultate.

# Prvi metod, kojim se mogu otkriti medjusobne korelacije parova prediktora je
# gledanje korelacione matrice.
library(corrplot)
corrplot(cor(longley))

# U ovom skupu podataka vidimo mnogo veoma jakih linearnih korelacija, sto se
# vidi i na grafiku parova promenljivih.
plot(longley)

# Ovo ukazuje da neki prediktori "rade isti posao", tj. buduci da su visoko
# korelirani, objasnjavaju isti uticaj na zavisnu promenljivu, pa bi mogao da se
# izabere jedan, a ostali izbace. Pogledacemo i druge metode pre izbacivanja.

# Drugi metod kojim se moze otkriti multikolinearnost je posmatranje
# uslovljenosti (condition number) matrice X^T * X. Ako su lambda_1, ...,
# lambda_p sopstvene vrednosti X^T * X, sortirane opadajuce, condition number za
# i-tu sopstvenu vrednost se dobija kao k_i = sqrt(lambda_1 / lambda_i). Velike
# vrednosti (>30) za k_i se smatraju indikatorima prisustva multikolinearnosti.
X <- model.matrix(model)[, -1] # izbacujemo kolonu za intercept
eig <- eigen(t(X) %*% X)
sqrt(eig$values[1] / eig$values)
# Vidimo da ima nekoliko vrednosti vecih od 30, sto znaci da problem postoji u
# vise prediktora, a ne da postoji samo jedna linearna kombinacija kojom se
# treba zabaviti.

# Treci metod nalazenja multikolinearnosti je merenje faktora inflacije
# disperzije (variance inflation factor - vif). Ideja je da se izracuna R_j^2 -
# koeficijent determinacije modela zavisnosti j-og prediktora u odnosu na
# ostale, pa se definise vif = 1 / (1 - R_j^2). Ako je vrednost R_j^2 blizu 1,
# tj prediktor j se moze dobro modelovati preko ostalih, onda ce vif biti jako
# veliko. Vrednost za koju treba obratiti paznju je vif > 5 ili 10.
vif(model)
# Vidimo da postoji jako velika prisutnost inflacije disperzije, stpo nam
# ukazuje da su ocene standardnih gresaka koeficijenata veoma preuvelicane, sto
# dovodi do "beznacajnih" rezultata testova, odnosno velikih p vrednosti i
# sirokih intervala poverenja.

# Razmotrimo sada koje promenljive da izbacimo
# Pogledajmo opet korelacije
corrplot(cor(longley))
# Promenljive Unemployed i Armed.Forces nisu jako korelisane sa ostalim, pa cemo
# ih sacuvati, dok ako pogledamo korelacije ostalih prediktora vidimo izuzetno
# visoke korelacije.
cor(X[, -c(3,4)])
# Dakle ovi prediktori objasnjavaju slicnu stvar, pa je mozda dovoljno uzeti
# samo jednu od njih, recimo Year.
model2 <- lm(Employed ~ Armed.Forces + Unemployed + Year, longley)
summary(model2)
summary(model)
# R^2 se jako malo smanjio, pa imamo slican kvalitet modela, a mnogo manje
# prediktora, pri cemu su sada svi znacajni.
vif(model2)
# Vidimo da su se sada i vrednosti faktora inflacije disperzije smirile.

# Interesantan primer multikolinearnosti je polinomijalna regresija. Vidimo to
# na primeru i kako R resava problem.
set.seed(216)
x <- runif(50, -3,3)
y <- 1 + x +x^2 +x^3 + rnorm(50, sd = 3)
plot(y ~ x)

m <- lm(y ~ x + I(x^2) + I(x^3))
summary(m)
vif(m)
cor(model.matrix(m))
# Korelacija izmedju x i x^3 je jako velika, pa je x postao beznacajan. Vidimo
# da postoji i znacajna inflacija disperzije. Ipak, nije korektno izbacivati
# odredjene clanove polinoma, buduci da nas je zanima da li je polinom treceg
# stepena dobar model ili je dovoljno ici do drugog.

# Zato funkcija poly u R-u kao rezultat vraca ortogonalne polinome, koji su
# konstruisani tako da prvi bude stepena 1, drugi stepena 2, a treci stepena 3,
# pri cemu ne postoji korelisanost izmedju njih, pa se iz njih moze videti da li
# postoji uticaj odredjenog reda polinoma na model. Ukoliko se prosledi
# parametar raw=TRUE funkciji poly onda dobijamo "obicne" polinome kao ranije.
m <- lm(y ~ poly(x, 4))
summary(m)
cor(model.matrix(m))