#TRANSFORMACIJE
#pokusavamo ako su ugrozene pretpostavke modela

#1.transformacie za Y (obicno logY ili Y^a, gde je obicno a<1)
#transformacija treba da se dobro slaze sa podacima i da dobijena veza ima smisla, da moze da se objasni
#Box-Cox metoda - najpoznatija metoda za transformacija Y, samo za pozitivne podatke
#umesto y nalazimo funkciju g(y), metodom maksimalne verodostojnosti

library(faraway)
library(MASS) #u njemu se nalazi funkcija boxcox
data(savings)
g<-lm(sr~pop15+pop75+dpi+ddpi,savings)
boxcox(g,plotit=T) #vidimo da je max oko 1, ali ne moze precizno da se odredi interval povenja
boxcox(g,plotit=T,lambda=seq(0.5,1.5,by=0.1)) #lakše se čita 95%-ni interval
#isprekidane linije su 95% interval poverenja
#maksimum je lambda=1, sto znaci da nema nikakve potrebe za transformacijom

#primer sa transformacijom
data(gala)
g<-lm(Species~Area+Elevation+Nearest+Scruz+Adjacent,gala)
boxcox(g,plotit=T)
boxcox(g,lambda=seq(0.0,1.0,by=0.05),plotit=T)   
#lambda bi trebalo da bude blizu 0.3 => treci koren iz Y

#Problemi kod Box-Cox metoda:
#autlejeri mogu da poremete model, smatra se da lambda nije vece od 1 jer se u suprotnom javljaju autlejeri
#ako ima negativnih Y-ona na pocetku, prvo treba dodati neku konstantu tako da svi Y-oni postanu pozitivni, a onda izvrsiti metod
#ako je opseg vrednosti velik, Box-Cox nece imati stvarni efekat

#2.Transformacije za X
#2.1.Regresija slomljenog stapa 
#ako sa grafika vidimo da razlicitim delovima podataka odgovaraju razliciti modeli
#primer toga smo imali u bazi savings za pop15>35 i pop15<35
#uticaj se vidi samo ako posmatrano jedno x i jedno y
plot(sr ~ pop15,savings,xlab="Pop'n under 15",ylab="Savings Rate")
abline(v=35,lty=5) 
g1<-lm(sr~pop15, savings, subset=(pop15 < 35))
g2<-lm(sr~pop15, savings, subset=(pop15 > 35))
#crtamo odvojeno modele na istom grafiku
segments(20,g1$coef[1]+g1$coef[2]*20,35,g1$coef[1]+g1$coef[2]*35)
segments(48,g2$coef[1]+g2$coef[2]*48,35,g2$coef[1]+g2$coef[2]*35)
#segments je za crtanje duzi, zadaje se koordinatama pocetne tacke i koordinatama krajnje tacke (4 agrumenta)

#prekid je nezgodan jer u zavisnosti od izbora prekidne tacke, dobijace se razlicite prave
#preporucuje se da "slomljeni stap" bude neprekidna funkcija
#jedino ako imamo dobar razlog da uzmemo neko drugo c tako da ne bude neprekidno, onda cemo ga uzeti i ostaviti prekid
#npr. temperatura vode, pa nam je nula stepeni neka prelazna tacka 
#definisemo dve funkcije
lhs<-function(x) ifelse(x < 35,35-x,0)
rhs<-function(x) ifelse(x < 35,0,x-35)
gb<-lm(sr~lhs(pop15)+rhs(pop15), savings) #model sa prethodnim funkcijama
x<-seq(20,48,by=1) 
py<-gb$coef[1]+gb$coef[2]*lhs(x)+gb$coef[3]*rhs(x)
lines(x,py,lty=2) #ucrtavanje prave modela
#dobija se neprekidan stap, a pojedinacne funkcije se spajaju u tacki 35

#2.2.Polinovi
#umesto samo X u modelu figurisu i X^2, X^3,...
#kako se bira stepen polinoma?
#dodajemo sve vece stepene (probamo d=1 (linearna regresija), vidimo da nije dobro pa uzmemo d=2 itd), do prvog koji ne bude statisticki znacajan
#Ili dodamo odmah polinom velikog stepena pa oduzimamo (od najveceg!) dok ne oduzmemo sve koji nisu znacajni
# d=1    
summary(lm(sr~ddpi,savings)) #=> koeficijent je znacajan 
# d=2
summary(lm(sr~ddpi+I(ddpi^2),savings)) #=> oba koeficijenta su znacajna
# d=3
summary(lm(sr~ddpi+I(ddpi^2)+I(ddpi^3),savings))  #=>nijedan koef nije znacajan
#treci stepen je narusio model, pa nam nije potreban

#prikaz znacajnosti manjeg reda polinoma
#pravimo linearnu transformaciju oduzimanjem 10
savings<-data.frame(savings,mddpi=savings$ddpi-10)   
summary(lm(sr~mddpi+I(mddpi^2),savings))     
#kvadratni clan jeste znacajan ali linearni nije, ali mi znamo da je znacajan (nije mogla znacajnost prediktora da se promeni linearnom transformacijom)

#zbog numericke stabilnosti dobro je korisititi ortogonalne polinome
g<-lm(sr~poly(ddpi,4),savings) #poly daje ortogonalne polinome
summary(g)
#pokazuje da je kvadratni model najbolji

#polinom vise promenljivih drugog stepena
g<-lm(sr~polym(pop15,ddpi,degree=2),savings)
summary(g)

#2.3.Regresioni splajnovi
#posmatramo funkciju koja ima prekide u oblasti deifnisanosti, a na intervalima na kojima je definisana je glatka
#potrebno je resenje koje se ponasa lokalno kao polinom, a globalno je kao slomljeni stap
#tj serija izlomljenih polinoma, nesto sastavljeno od razlicitih polinoma
#glatka je na cvorovima

#pretpostavimo da je pravi model y=sin(2*pi*x^3)^3+e i generisemo podatke
funky<-function(x) sin(2*pi*x^3)^3
x<-seq(0,1,by=0.01)
y<-funky(x)+0.1*rnorm(101)
matplot(x,cbind(y,funky(x)),type="pl",ylab="y",pch=18,lty=1)
#gledamo kako polinomi reda 4 i reda 12 odgovaraju podacima
g4<-lm(y~poly(x,4))
g12<-lm(y~poly(x,12))
matplot(x,cbind(y,g4$fit,g12$fit),type="pll",ylab="y",pch=18,lty=c(1,2))
#polinom 4. stepena je crvena isprekidana linija => ne odgovara podacima
#polinom 12.stepena je zelena linija => odgovara podacima

#kreiramo B-splajn bazu
#potrebno je imati 3 dodatna cvora na pocetku i na kraju kako bi baza bila dobra
library(splines)
knots <- c(0,0,0,0,0.2,0.4,0.5,0.6,0.7,0.8,0.85,0.9,1,1,1,1)
#zbugsnuti cvorovi oko 0.85 zbog minimuma koji zelimo da istaknemo
#tamo gde ima vise promena, se stavlja veci broj cvorova
#pravi se funkcija od k-tog do (k+3)-ceg cvora zbog argumenta ord koji je po defaultu 4 (help(splineDesign))

bx<-splineDesign(knots,x)
gs<-lm(y~bx)
matplot(x,bx,type="l") #osnovna funkcija
matplot(x,cbind(y,gs$fit),type="pl",ylab="y",pch=18,lty=1) #prilagodjena osnovna funkcija
#vidimo da je ovo jos bolje od polinoma 12. stepena
#za splajnove je vrlo tesko interpretirati dobijeni rezultat, pa su dobri za predvidjanje ali ne i za objasnjavanje fenomena

#IZBOR PROMENLJIVIH
#biramo promenljive koje daju najbolji model
#zelimo da model sto bolje odgovara podacima ali i da sto bolje i lakse moze da se interpretira
#medjutim, ta dva zahteva su u suprotnosti jedan s drugim
#postoje dva osnovna nacina za izbor promenljivih: testiranje i kriterijumi

#procedure na osnovu testiranja
#najjednostavnija: eliminacija unazad (backward elimination)
data(state)
statedata<-data.frame(state.x77,row.names=state.abb) #pravimo bazu
g<-lm(Life.Exp ~ ., data=statedata) #y=life expectancy, X-sve ostale promenljive
summary(g)
#za neke promenljive nije ocekivano da nisu znacajne (na primer income)

#povrsina drzave je najmanje znacajna, pa nju prvu izbacujemo
g<-update(g, . ~ . - Area) #menjamo model, umesto da pravimo novi
summary(g)

#sledece izbacujemo nepismenost (HS.Grad(procenat sa zavrsenom srednjom skolom) i pismenost su korelisani)
g<-update(g, . ~ . - Illiteracy)
summary(g)

#izbacujemo prihod
g<-update(g, . ~ . - Income)
summary(g)

#izbacujemo populaciju (posto je na granici ne mora da se izbaci u slucaju da je njena interpretacija jednostavna)
g<-update(g, . ~ . - Population) 
summary(g)
#ostale su samo znacajne promenljive i tu stajemo sa izbacivanjem
#prilikom ovakvih izbacivanja skoro uopste se nije promenio Multiple R-squared, sto opravdava izbacivanje
#ako iz poslednjeg modela izbacimo jos neku promenljivu, onda se Multiple R^2 dosta smanjuje

#promenljive izbacene iz modela mogu biti u vezi sa y
summary(lm(Life.Exp ~ Illiteracy+Murder+Frost, statedata))
#nepismenost je u ovom modelu znacajna, ali HS Grad u modelu daje bolje rezultate

#Procedure na osnovu kriterijuma 
#daje sto manji model sa sto boljim slaganjem
#1.AIC
#Akaike Information Criterion: AIC = -2 max lnL + 2p
#Bayes Information Criterion (BIC) - bolji za manje modele: BIC = -2max lnL + p ln(n)
#kod regresionih modela razlika u konstanti moze da se zanemari

g<-lm(Life.Exp ~ ., data=statedata)
step(g) #AIC kriterijum
#dobija se isti model kao eliminacijom unazad (sa ukljucenom populacijom)

#2.Poprsavljeni (adjusted) R^2
install.packages("leaps")
library(leaps)
#pretpostavimo da znamo koliko parametara hocemo (p),onda je za svako p najbolji model:
b<-regsubsets(Life.Exp~.,data=statedata)     
(rs<-summary(b))
#najbolji model sa jednim prediktorom je y~Murder
#najbolji model sa dva prediktora je y~Murder+HS.Grad ...

plot(2:8,rs$adjr2,xlab="No. of Parameters",ylab="Adjusted R-square")
#najveci Adjusted R-square za 5 parametara (intercept+4 prediktora), pa je on najbolji

#3.Cp kriterijum
#uz pomoc Cp-a odredimo najbolji model, fiksiramo broj prediktora i dobijemo najbolji model od mogucih
plot(2:8,rs$cp,xlab="No. of Parameters",ylab="Cp Statistic")
abline(0,1) #Cp=p linija
#ako bi izabrali 2 ili 3 promenljive, Cp bi bilo dosta veliko 
#dobijemo da je potrebno da imamo 4 promenljive: y i 3 prediktora (moze se uzeti i 5 jer je 4 na granici) 
#iz prethodne tabele za cetvrti model citamo promenljive koje imaju *:Population, Murder, HS.Grad, Frost
#iste rezultate smo dobijali i ranije

#Metodi za izbor promenljivih su dosta osetljivi na autlajere i uticajne tacke
#proveravamo koja tacka ima najvecu tezinu
h<-lm.influence(g)$hat
names(h)<-state.abb
rev(sort(h))
b<-regsubsets(Life.Exp~.,data=statedata,subset=(state.abb!="AK"))   #izbacimo aljasku
#Aljaska ima najvecu tezinu, izbacujemo je da vidimo da li ima promena
rs<-summary(b)
rs$which[which.max(rs$adjr),] #koje treba zadrzati po adjusted R^2 kriterijumu?
#Area je sada ukljucena u model
rs$adjr2

stripchart(data.frame(scale(statedata)),vertical=TRUE,method="jitter")
#jitter - stvara "buku" (da se tacke pomere malo levo ili desno, da bi grafik bio pregledniji, da bude manje poklapanja))
#scale() - standardizuje podatke 

#transformacija promenljivih Population i Area
b<-regsubsets(Life.Exp~log(Population)+Income+Illiteracy+Murder+HS.Grad+Frost+log(Area),statedata)
rs<-summary(b)
rs$which[which.max(rs$adjr),]
rs$adjr2
#dobili smo jos bolji model (veci R^2), koji se jos bolje slaze sa podacima i ima 4 prediktora

#NEDOSTAJUCI PODACI 
#1.ako je baza dovoljno velika, izbrisati opservacije sa nedostajucim vrednostima
#2.dopuniti bazu predvidjenom vrednoscu na osnovu ostalih podataka ili srednjom vrednoscu
#3.EM metod

data(chmiss)
chmiss
g<-lm(involact ~ .,chmiss)
summary(g)
#R pravi model zanemarujuci slucajeve sa nedostajucim vrednostima, ali DF je 21 (skoro polovina podataka nedostaje)
#zamena nedostajucih vrednosti srednjom vrednoscu
cmeans<-apply(chmiss,2,mean,na.rm=T)
cmeans
mchm<-chmiss
for(i in c(1,2,3,4,6)) mchm[is.na(chmiss[,i]),i]<-cmeans[i]
g<-lm(involact ~ ., mchm)
summary(g)
#predvidjanje nedostajucih vrednosti
gr<-lm(race~fire+theft+age+income,chmiss)
chmiss[is.na(chmiss$race),]
predict(gr,chmiss[is.na(chmiss$race),])
#predvidjena vrednost je negativna a ne bi trebala da bude, pa vrsimo transformaciju
logit<-function(x) log(x/(1-x))
ilogit<-function(x) exp(x)/(1+exp(x)) #inverzna transformacija
gr<-lm(logit(race/100) ~ fire+theft+age+income,chmiss)
ilogit(predict(gr,chmiss[is.na(chmiss$race),]))*100
#posto je koriscena baza u kojoj su namerno izbacene vrednosti mogu se vredvidjene vrednosti uporediti sa stvarnim
data(chredlin)
chredlin$race[is.na(chmiss$race)]
#prve dve predvidjene vrednosti su dobre, dok ostale i nisu bas
#uspesnost predvidjanja zavisi od kolinearnosti izmedju prediktora