#Netipicne opservacije
#tacke koje ne odgovaraju modelu ili uticu na promenu modela
#1. Leverages (tezinske tacke) - "skrecu" model ka sebi vise od ostalih
library(faraway)
data("savings")
g<-lm(sr~pop15+pop75+dpi+ddpi, savings)
ginf<-influence(g)
ginf$hat #hi-ovi    
sum(ginf$hat) #daje broj koeficijenata modela
#graficko ispitivanje uticajnih opservacija, prikaz usecene normalne raspodele
countries<-row.names(savings)
halfnorm(lm.influence(g)$hat,labs=countries,ylab="Leverages") #porede se reziduali i kvantili polunormalne raspodele
#labs=countries znaci da se dodaju imena za dve najvece tacke (za vise od 2 tacke dodaje se argument nlab)

#studentizovani reziduali
gs<-summary(g)
gs$sig
stud<-residuals(g)/(gs$sig*sqrt(1-ginf$hat))  #standardizujemo reziduale
#ocekujemo da prate normalnu raspodelu
qqnorm(stud)
abline(0,1)
#reziduali su normalni

#2. Autlejeri - tacke koje odudaraju od modela vise od ostalih:
#postoje autlejeri koji menjaju oblik prave i oni koji ne menjaju
#sa grafika vidimo potencijalni autlejere, napravimo model sa i bez njega, i uporedimo
#ako je razlika velika, jeste autlejer

#Dzeknajf
jack<-rstudent(g) #dzeknajf reziduali
jack[which.max(abs(jack))] #daje najvecu vrednost reziduala, pa je to kandidat za autlejer
qt(.05/(50*2),44)   #delimo alfa sa n i jos delimo za 2 zbog dvostranog testa, 44=n-p-1
#nema autlejera

#Autlejeri mogu medusobno da se cuvaju
#Autlejeri u jednom modelu ne moraju biti autlejeri u drugom modelu
#Veliki reziduali mogu biti posledica nenormalnosti raspodele
#Uticaj pojedinacne tacke zavisi i od velicine baze

#baza sa vise autlejera
data(star)
plot(star$temp,star$light,xlab="log(Temperature)",ylab="log(Light Intensity)")
#vidimo da 4 tacke odudaraju od ostatka, inace izgleda da postoji pozitivna korelacija izmedju x i y
#pravimo model
ga<-lm(light~temp, star)
abline(ga) #regresiona prava
#cetiri tacke levo su tezinske, pa su odvukle liniju na tu stranu (linija nema smisao)
#u kom rangu se krecu standardizovani (jackknife) reziduali:
range(rstudent(ga))
qt(1-0.05/(2*47),44)#kriticna oblast:(ovde za n=47, p=2)
#po testu, nema reziduala, iako mi vidimo da ih ima (ova 4 autlajera se dobro stite medjusobno)

#model kada iskljucimo ove tacke
ga<-lm(light~temp, data=star, subset=(temp>3.6))
abline(ga,lty=2)#dobijamo isprekidanu liniju koja mnogo bolje opisuje linearnu vezu

#3. Uticajne tacke: njihovim uklanjanjem iz modela, on se drasticno menja
#mogu i ne moraju biti tezinske i mogu a ne moraju biti autlejeri (ali obicno su jedno od ta dva)
#metoda Kukove statistike meri Kukovo rastojanje od modela
#mera uticaja je razlika u koeficijentima modela sa i bez te tacke
#utvrdjuje i efekat reziduala (odstupanje Y-ona od svoje ocene) i tezinski efekat
cook<-cooks.distance(g)
#tacke za koje je Kukovo rastojanje veliko su uticajne
halfnorm(cook,3,labs=countries,ylab="Cook's distances")
# Libija je vrlo uticajna tacka, pa je izbacujemo i pravimo novi model
gl<-lm(sr~pop15+pop75+dpi+ddpi,savings,subset=(cook < max(cook))) #izbacen je max(cook)
#poredimo stari i novi model:
summary(gl)
summary(g)
#koeficijent uz ddpi, bez Libije je 0.6102790, a sa libijom 0.4096949 (promena za 50%)
#nije dobro kada su koeficijenti toliko osetljivi na pojedinacne opservacije

# hocemo da vidimo sta se desava kad izbacimo neke druge promenljive
plot(ginf$coef[,2],ylab="Change in pop15 coef") #coef[,2] znaci ocene za drugu promenljivu pop15

identify(1:50,ginf$coef[,2],countries)
#dodaju se imena na grafik, tamo gde kliknemo; kad se klikne finish, onda ce se ispisati imena
#=>Japan najvise odskace; sta se desava ako ga izostavimo?
gj<-lm(sr~pop15+pop75+dpi+ddpi,savings,subset=(countries != "Japan"))
summary(gj)
#koeficijenti su se znacajno promenili, znaci da je Japan uticajan
#ddpi vise nije znacajno, R^2 se dosta smanjio

#Provera strukture modela
d<-residuals(lm(sr~pop75+dpi+ddpi,savings))
m<-residuals(lm(pop15~pop75+dpi+ddpi,savings))
plot(m,d,xlab="pop15 residuals",ylab="Savings residuals")
#poredimo nagibe na grafiku i orginalni
coef(lm(d ~ m))
coef(g)
abline(0,coef(g)['pop15'])
#nagib na grafiku i nagib za pop15 u regresiji su isti
#grafik parcijalnih reziduala
plot(savings$pop15,residuals(g)+coef(g)['pop15']*savings$pop15,xlab="pop'n under 15", ylab="Savings(Adjusted)")
abline(0,coef(g)['pop15'])
prplot(g,1) #ugradjena funkcija
#imamo dve grupe na grafiku, pa to treba ispitati
g1<-lm(sr~pop15+pop75+dpi+ddpi,savings,subset=(pop15 > 35))
g2<-lm(sr~pop15+pop75+dpi+ddpi,savings,subset=(pop15 < 35))
summary(g1) #ne postoji veza
summary(g2) #postoji veza
