# Dijagnostika
# Uslovi koji se proveravaju:
# greske su normalne, nezavisne i sa konstantnom disperzijom
# korektnost modela: E(y)=X*beta
# opservacije koje odstupaju od ostalih, ili od modela

library(faraway)
data(savings)
g<-lm(sr~pop15+pop75+dpi+ddpi,savings)

# Provera uslova za gresku
#1. konstantna disperzija:
#potrebno je da disperzija reziduala ne zavisi ni od X-eva ni od Y
#najlakse se vidi preko grafika reziduala i ocena
plot(fitted(g),residuals(g),xlab="Fitted",ylab="Residuals") 
abline(h=0)
# apsolutne vrednosti reziduala i ocene (lakse se vide nepravilnosti)
plot(fitted(g),abs(residuals(g)),xlab="Fitted",ylab="|Residuals|")
# "numericki" nacin provere zavisnosti reziduala i ocena
summary(lm(abs(residuals(g)) ~ fitted(g)))

par(mfrow=c(3,3))
for(i in 1:9) plot(1:50,rnorm(50)) # konstantna disperzija
for(i in 1:9) plot(1:50,(1:50)*rnorm(50)) # izrazena nekonstantna disperzija (veæi indeksi imaju vecu disperziju) 
for(i in 1:9) plot(1:50,sqrt((1:50))*rnorm(50)) # slaba nekonstantna disperzija (manja razlika nego kod prethodnog)
for(i in 1:9) plot(1:50,cos((1:50)*pi/25)+rnorm(50))  # nije linearan model
par(mfrow=c(1,1))

# reziduali u odnosu na X
plot(savings$pop15,residuals(g),xlab="Population under 15",ylab="Residuals")  
plot(savings$pop75,residuals(g),xlab="Population over 75",ylab="Residuals")
#na prvom grafiku su grupisani reziduali u dve grupe (ispod i iznad 35)
#poredimo da li je znacajna razlika u disperziji izmedju grupa
#var.test je Fiserov test za jednakost disperzija
var.test(residuals(g)[savings$pop15>35],residuals(g)[savings$pop15<35])
# mala p vrednost => postoji razlika u disperziji (disperzija nije konstantna)

# resavanje nekonstantnosti disperzije:
# umesto y -> h(y), gde je h neka f-ja tako da je D(h(y))=const (najcesce log(y) ili sqrt(y))
data(gala)
gg<-lm(Species~Area+Elevation+Scruz+Nearest+Adjacent, gala)
#grafik reziduala i ocena
plot(fitted(gg),residuals(gg),xlab="Fitted",ylab="Residuals")

# probamo korenu transformaciju
gs<-lm(sqrt(Species)~Area+Elevation+Scruz+Nearest+Adjacent, gala)
# proveravamo graficki
plot(fitted(gs),residuals(gs),xlab="Fitted",ylab="Residuals",main="Square root Response")

#2. normalnost reziduala
#grafik sortiranih reziduala i Fi^{-1}(i/n+1), i=1,...,n 
qqnorm(residuals(g),ylab="Residuals")
qqline(residuals(g)) 
#histogram reziduala
hist(residuals(g))   

#Q-Q grafici poznatih raspodela
par(mfrow=c(3,3))
for(i in 1:9) qqnorm(rnorm(50)) #normalna
for(i in 1:9) qqnorm(exp(rnorm(50))) #lognormalna (asimetricna raspodela)
for(i in 1:9) qqnorm(rcauchy(50)) #kosijeva (raspodela sa teskim repom)
for(i in 1:9) qqnorm(runif(50)) #uniformna (raspodela sa lakim repom)
par(mfrow=c(1,1))
#ako dobijemo da ne vazi normalnost, pokusamo da izmenimo model
#mala odstupanja mogu da se tolerisu (pogotovo za velike uzorke)
#test normalnosti: Shapiro-Wilk test
shapiro.test(residuals(g))
#velika p-vrednost => prihvatamo H0
#najbolje je oslanjati se na QQ plot, a shapiro test koristiti samo kao pomocni metod za proveru normalnosti

#3. Nezavisnost reziduala
#proverava se ispitivanjem korelisanosti
#graficki: grafik reziduala u odnosu na vreme, kao i zavisnost susednih reziduala
#test: Durbin-Watson

data(airquality)
head(airquality)
#odnos svake dve promenljive
pairs(airquality,panel=panel.smooth)

#rad sa nedostajucim vrednostima:
# na.action=na.exclude => ne uzima u obzir nedostajuce vrednosti, ali pamti gde su bile
# na.action=na.omit => izbacujemo nedostajuce vrednosti i vise ih uopste ne vidimo

g<-lm(Ozone~Solar.R+Wind+Temp, airquality, na.action=na.exclude)
summary(g)
#proveravamo disperziju reziduala:
plot(fitted(g),residuals(g),xlab="Fitted",ylab="Residuals")

#vidimo da ima nelinearnosti i nekonstantnosti disperzije, pa transformisemo podatke
gl<-lm(log(Ozone)~Solar.R+Wind+Temp, airquality, na.action=na.exclude)
plot(fitted(gl),residuals(gl),xlab="Fitted",ylab="Residuals")

#Provera korelisanosti reziduala:
#u odnosu na vreme:
plot(residuals(gl),ylab="Residuals")
abline(h=0)
#kada bi bila velika korelisanost, imali bismo dugacke serije ispod ili iznad linije

#proveravamo korelisanost susednih:
plot(residuals(gl)[-153],residuals(gl)[-1], xlab=expression(hat(epsilon)[i]),ylab=expression(hat(epsilon)[i+1]))
# [-153] znaci svi osim 153-ceg
# residuals(gl)[-1] je shiftovano udesno od residuals(gl)[-153]
#ne vidi se da ima korelisanosti

#jos jedan nacin - linearni model reziduala u odnosu na same sebe
summary(lm(residuals(gl)[-1] ~ -1+residuals(gl)[-153]))  

install.packages("lmtest")
library(lmtest)
#Durbin Watsonov test korelisanosti
dwtest(Ozone~Solar.R+Wind+Temp, data=na.omit(airquality))
# dobija se da nisu korelisani