#Problemi sa prediktorima

#1. greske pri merenju
library(faraway)
data(cars)
plot(dist~speed, cars,ylab="distance")
g<-lm(dist~speed, cars)
summary(g)
abline(g)
#dodajemo "gresku merenja" prediktorima
ge1<-lm(dist~I(speed+rnorm(50)), cars)    
coef(ge1)
abline(ge1,lty=2)
#dodajemo jos vecu gresku
ge2<-lm(dist~I(speed+2*rnorm(50)), cars)  
coef(ge2)
abline(ge2,lty=3)
#5 puta veća greška (disperzija)
ge5<-lm(dist~I(speed+5*rnorm(50)), cars) 
coef(ge5)
abline(ge5,lty=4)
#=> nagib prave se smanjuje

#formiranje disperzija
vv<-rep(1:5/10,each=1000)
slopes<-numeric(5000)
#simulacija nagiba
for(i in 1:5000) 
  slopes[i]<-lm(dist~I(speed+sqrt(vv[i])*rnorm(50)), cars)$coef[2]    
#srednji nagibi za svaku disperziju 
betas<-c(coef(g)[2],colMeans(matrix(slopes,nrow=1000)))
variances<-c(0,1:5/10)+0.5
plot(variances,betas,xlim=c(0,1),ylim=c(3.86,4))    
#veza je poprilicno linearna; sto je veca disperzija, to se beta priblizava nuli (nagib)
gv<-lm(betas~variances)
coef(gv)
points(0,gv$coef[1],pch=3) #ekstrapolacija (krstic na grafiku)

#2. neodgovarajuci red velicine prediktora
#ako su X-evi razlicitog reda velicine, to moze da bude problem
#nije dobro ako su koeficijenti preblizu nuli ili prevelike vrednosti, najbolje bi bilo da u bete bliske jedinici
#linearna transformacija X-eva nece poremetiti model

data(savings)
g<-lm(sr~pop15+pop75+dpi+ddpi,savings)
summary(g)
#koeficijent uz dpi je vrlo mali, pa može da se skalira
#to znaci da dpi uzima velike vrednosti, pa cemo ga podeliti sa hiljadu
g<-lm(sr~pop15+pop75+I(dpi/1000)+ddpi,savings)
summary(g)
#nepromenjene vrednosti test statistika
#nema uticaja na strukturu i na kvalitet modela

#funkcija za skaliranje svih X-eve odjednom: scale; ovim se dobija EX=0, DX=1
scsav <- data.frame(scale(savings))
g<-lm(sr~ ., scsav)
summary(g)
#svi koeficijenti su obavezno po modulu manji od 1
#dobijamo intercept nula (jer je srednja vrednost svih x-eva nula, a regresiona "prava" prolazi kroz srednju vrednost)
#medjutim, interpretacija je komplikovanija jer se onda skalira i Y
#obicno se skalira samo jedno X

#3. Kolinearnost
#sta ako su X-evi korelisani (t(X)*X je singularna ili blizu singularna)?
#kako da proverimo ovo:
#1. korelaciona matrica za X
#2. za svako Xi: regresija pojedinacnog Xi u osnosu na sve ostale, ako se dobije veliko R vazi korelisanost
#3. sopstvene vrednosti matrice t(X)*X

#Rj^2 blisko 1 => variance inflation factor (VIF)=1/(1-Rj^2) velik

data(seatpos)
g<-lm(hipcenter ~ ., seatpos)
summary(g)
#=> nijedna promenljiva nije znacajna (p-vrednosti), a model je znacajan (F statistiku)
#postoji problem
round(cor(seatpos),3)   #na tri decimale
#ima dosta velikih korelacija

#drugi nacin:
x<-model.matrix(g)[,-1]  #matrica modela bez intercepta
e<-eigen(t(x) %*% x)     #sopstvene vrednosti, R daje sortirane vrednosti
e$val #sopstvene vrednosti - ako je neka mnogo mala, to je indikator korelisanosti
sqrt(e$val[1]/e$val) #condition numbers (skalirane sopstvene vrednosti)
#ovako transformisane, prevelike vrednosti(>=30) su problem, vidimo da ima vise ekstremnih slucajeva 
#provera regresijom, za prvi X (regresija prve kolone u odnosu na ostale) i prvo gledamo R^2:
Rsq<-summary(lm(x[,1] ~ x[,-1]))$r.squared
#odatle racunamo Variance inflation factor (koliko korelisanost utice na disperziju):
1/(1-Rsq) 
#ugradjena funkcija koja racuna VIF za sve prediktore
vif(x)
#lose je ako je VIF veliki

#ako se malo promeni X, doci ce do velikih promena u modelu
#ako malo promenimo y, a model matrica ostane ista
g<-lm(hipcenter+10*rnorm(38) ~ ., seatpos) #ima 38 podataka, pa zato dodajemo toliko brojeva   
summary(g)  #koeficijenti vise nisu onakvi kakvi su bili

#resenje je izbaciti neku od promenljivih pa proveriti ponovo
#posmatramo korelaciju izmedju promenljivih vezanih za duzine, ima ih ukupno 6, od trece do osme:
round(cor(x[,3:8]),2)  
#sve su medjusobno jako korelisane pa cemo izbaciti sve osim jedne, najosnovnije
g2<-lm(hipcenter~Age+Weight+Ht, seatpos)  
summary(g2)
#vidimo da su slicne R^2 za g i g2 => model je znacajan, ali sa mnogo manje promenljivih
#ali sada je visina postala znacajna promenljiva
#za problem kolinearnosti se koristi tzv. nazubljena regresija
