# DIJAGNOSTIKA

# Bavicemo se dijagnostikom modela, tj. razmatranjem reziduala u cilju provere
# ispunjenosti pretpostavki modela.

# Izvrsicemo analizu nad podacima cars koji opisuju zavisnost duzine traga
# kocenja od brzine vozila.

cars <- SenSrivastava::E6.1
colnames(cars) <- c("dist", "speed")

plot(dist ~ speed, cars)

# Napravimo inicijalni model za nase podatke
m <- lm(dist ~ speed, cars)
summary(m)

abline(m)


# Simetricnost oko nule

# Za proveru prvog uslova G-M, moze se gledati grafik reziduala u odnosu na
# prediktore, ili u odnosu na predvidjene vrednosti modelom. Ako uocimo neke
# pravilnosti, to znaci da nisu uracunali svu zavisnost i trebalo bi razmisliti
# o dodavanju prediktora.
plot(m$res ~ cars$speed)
plot(m$res ~ m$fit)
# Ovde postoji neka dodatna zavisnost
# Laksi nacin je koristeci funkciju residualPlot(s) iz paketa "car"
library(car)
residualPlots(m)
# Jasno se vidi da postoji neka kvadratna zavisnost reziduala od prediktora,
# pa cemo da dodamo kvadratni clan u model
m <- update(m, . ~ . + I(speed^2))
summary(m)
# Pogledajmo opet grafike reziduala
residualPlots(m)
# Sada su ravni i simetricni oko nule, sto znaci da je prvi uslov G-M ispunjen.


# Testiranje normalnosti reziduala

# Najlaksi nacin za testiranje normalnosti je gledati qqplot. Ranije smo
# koristili funkcije qqnorm i qqline, a sada cemo koristiti funkciju qqPlot iz
# "car" paketa koja daje lepsi prikaz, uz intervale poverenja.
library(car)
set.seed(216)
x <- rnorm(50)
qqPlot(x)
# Neka odstupanja od normalnosti postoje, jer uzorak ne moze biti savrseno
# normalan, ali vidimo da je u granicama tolerancije.

# Neki uzorak koji je nenormalan...
set.seed(216)
x <- rexp(50)
qqPlot(x)
# Uocava se jasno odstupanje od prave linije

# Zanima nas normalnost reziduala naseg modela. Ispravnije je crtati
# standardizovane reziduale kada se gleda normalnost, nego obicne reziduale, da
# bi se uzela u obzir disperzija reziduala (koji su *ocene* gresaka, pa nemaju
# istu disperziju kao greske).
qqPlot(rstandard(m))
# Reziduali izgledaju normalno, sumnjiva je samo tacka 23.

# Testiranje normalnosti se moze izvrsiti i pomocu testova normalnosti. Jedan od
# njih je Shapiro-Wilk test.
shapiro.test(rstandard(m))
# U paketu "nortest" su implementirani i drugi testovi normalnosti, poput
# Anderson-Darlingovom i Kramer-fon Mizesovog.
library(nortest)
ad.test(rstandard(m))
cvm.test(rstandard(m))
# Shapiro-Wilk test ima visoku p vrednost, pa ne bismo odbacili normalnost, ali
# ostala dva testa odbacuju hipotezu o normalnosti. Za sada cemo nastaviti sa
# daljom analizom reziduala, pa cemo se vratiti na normalnost kada dobijemo
# bolji model.

# Problem Heteroskedasticnosti

# Druga pretpostavka Gaus-markova je da greske imaju jednaku disperziju, sto
# nazivamo homoskedasticnosu. Ova pretpostavka je cesto narusena i bitno je to
# detektovati i korigovati ovaj problem u sto vecoj meri.
# Heteroskedasticnost proveravamo graficki sa grafika zavisnosti reziduala od
# predvidjenih vrednosti. Umesto reziduala korisno je gledati i kvadrate
# reziduala ili njihove norme. Takodje, korektnije je gledati standardizovane
# reziduale nego obicne.
plot(rstandard(m) ~ fitted(m))
# uocavamo da se rasprsenost povecava sa rastom ocene y
plot(rstandard(m)^2 ~ fitted(m))
# sa grafika kvadrata reziduala se jasnije vidi porast disperzije
plot(abs(rstandard(m)) ~ fitted(m))
# takodje i apsolutne vrednosti vise govore od obicnih reziduala.

# Jedan od formalnih testova homoskedasticnosti je Brojs-paganov test
# On je implementiran u "lmtest" paketu
library(lmtest)
bptest(m)
# p vrednost je jako mala, pa odbacujemo hipotzu o homoskedasticnoti.
# Treba pomenuti da, iako ovaj test nekad prihvati hipotezu o
# homoskedasticnosti, ukoliko grafici reziduala indikuju da postoji problem,
# treba pristupiti daljoj analizi i ispraviti heteroskedasticnost, iako nije
# formalno znacajna po ovom testu.

# Resavanje heteroskedasticnosti

# Resavanju heteroskedasticnosti se moze pristupiti na vise nacina. Prvi koji
# cemo isprobati je transformisanje zavisne promenljive. Neka intuicija bi bila
# da, ako disperzija raste sa porastom vrednosti prediktora, uzimanje logaritma
# ili korena od zavisne promenljive bi moglo, zbog konkavnosti, da ublazi
# odstupanja of proseka pri velikim vrednostima prediktora.
# U nasem slucaju, ocigledan je porast disperzije reziduala sa rastom ocena y,
# pa ima smisla probati logaritamsku ili korenu transformaciju
plot(abs(rstandard(m)) ~ fitted(m))
# Probacemo prvo logaritam
plot(log(dist) ~ speed, cars)
# izgleda mnogo bolje, napravimo model
mlog <- lm(log(dist) ~ speed + I(speed^2), cars)
summary(mlog)
plot(abs(rstandard(mlog)) ~ fitted(mlog))
# Sada jeste blaza heteroskedasticnost, ali sada opada, previse smo korigovali.
bptest(mlog) # test potvrdjuje heteroskedasticnost
# Probacemo ipak i korenu transformaciju
plot(sqrt(dist) ~ speed, cars)
# ovo izgleda jos lepse. Pravimo model
msq <- lm(sqrt(dist) ~ speed + I(speed^2), cars)
summary(msq) # R^2 = 0.9251, za sad najbolji model
plot(abs(rstandard(msq)) ~ fitted(msq))
# Bolje izgledaju reziduali, mada postoje blage naznake heteroskedasticnosti
bptest(msq) # bptest ne odbacuje homoskedasticnost

# Probacemo jos jednu transformaciju, a to je Box-Cox transformacija, koja ima
# za cilj da podatke nacini normalnijim i homoskedasticnijim. Ona zavisi od
# parametra lambda, koji se odredjuje metodom maksimalne verodostojnosti.
# Parametar lambda dobijamo pozivanjem boxcox funkcije iz paketa MASS i
# nalazenjem onog lambda za koje je najveca vrednost funkcije verodostojnosti
library(MASS)
bc <- boxcox(m)
bc$x[which.max(bc$y)]
# U ovom slucaju je lambda oko 0.35
BC <- function(y) (y^0.35 - 1)/0.35
# Vidimo kako izgleda grafik, pa napravimo model...
plot(BC(dist) ~ speed, cars)
mbc <- lm(BC(dist) ~ speed + I(speed^2), cars)
summary(mbc)
plot(abs(rstandard(mbc)) ~ fitted(mbc))
# Sada grafik reziduala izgleda jako dobro, prisutnost heteroskedasticnosti je
# prakticno zanemarljiva.
bptest(mbc) # bptest daje vrlo visoku p vrednost, sto potvrdjuje nase zapazanje