## Intervali poverenja

# Prvo cemo se baviti intervalima poverenja za koeficijente modela.
# Uzecemo primer proste regresije, jer ima samo 2 koeficijenta, sto se moze
# lako vizuelizovati.

# Napravicemo vestacke podatke
set.seed(888)
x <- runif(50, -3, 3)
y <- 1 + 3 * x + rnorm(50)
plot(x,y)
model <- lm(y ~ x)
abline(model)
summary(model)

# Zajednicki interval poverenja za koeficijente b_0 i b_1 je u obliku elipse
# Crtamo ga koristeci funkciju ellipse iz ellipse paketa
#install.packages("ellipse")
library(ellipse)
plot(ellipse(model), type = "l")
points(model$coef[1], model$coef[2])
points(1, 3, col="red")
# Stvarna vrednost je unutar intervala

# Nacrtajmo sada i pojedinacne intervale poverenja za b_0 i b_1
# Vrednosti krajeva intervala poverenja dobijamo funkcijom confint
confs <- confint(model)
confs
# Nacrtacemo funkcijom abline, linije za b_0 i za b_1 cime dobijamo pravougaonik
abline(v = confs[1, ]) # za intercept crtamo vertikalne linije
abline(h = confs[2, ]) # za x crtamo horizontalne linije

# Netacna ilustracija intervala poverenja koja daje neku intuiciju
# Ako uzmemo mnogo razlicitih uzoraka, svaki put dobijemo malo drugacije ocene
# beta i one ce se nalaziti u nekoj elipticnoj oblasti oko stvarnog beta, a
# tacnost intervala poverenja mozemo da vidimo kada uporedimo koliko vrednosti
# je izaslo iz intervala poverenja - za 200 ponavljanja ocekujemo da 10 (5%)
# ispadne iz 95% intervala poverenja.
# NAPOMENA. Ovo nije tacna interpretacija intervala poverenja, ali se dobija lep
# rezultat jer je ocena za seed 888 vrlo bliska teoretskom beta.

# Prava interpretacija je: Ako ponovimo uzorkovanje mnogo puta, u 95% slucajeva
# ce stvarna vrednost beta pripadati intervalu poverenja napravljenom na osnovu
# tog uzorka. Beta je broj, nije slucajan i nema raspodelu, vec je interval
# poverenja slucajan (jer je \hat{beta} slucajna) i za njega moze da se desi da
# obuhvati beta ili ne - za to je verovatnoca 95%.

# Kad smo to resili, da se vratimo intuitivnoj ilustraciji...
model_coefs <- replicate(200, {
  x <- runif(50, -3, 3)
  y <- 1 + 3 * x + rnorm(50)
  model <- lm(y ~ x)
  model$coef
})
points(t(model_coefs), col="darkgreen")

# I za ove intervale se moze raditi Bonferonijeva korekcija. Analogno testiranju
# hipoteza, nivo intervala poverenja povecamo sa 1-0.05 na 1-0.025,
# jer imamo dva prediktora
bonfi <- confint(model, level = 0.975)
abline(v = bonfi[1, ], col = "blue") # za intercept crtamo vertikalne linije
abline(h = bonfi[2, ], col = "blue") # za x crtamo horizontalne linije


# Intervali poverenja su jako povezani sa testiranjem hipoteza o koeficijentima.
# Naime, 95% interval poverenja za beta predstavlja skup razlicitih vrednosti
# beta_0 za koje se *ne bi* odbacila hipoteza da je beta = beta_0. Specijalno,
# ako 0 priprada intervalu poverenja za beta, onda bismo zakljucili da
# koeficijent beta nije znacajan. Ako je 0 van intervala poverenja, onda imamo
# statisticku znacajnost.

# Podseticemo se primera iz testiranja hipoteza na kojima smo videli razliku
# izmedju uzastopnih pojedinacnih testova i grupnih testova i videti kako se to
# odrazava na intervale poverenja.

# Prvo napravimo funkciju koja ce da crta sliku intervala za model.
draw_intervals <- function(model, ind = c(1, 2), bonfi = FALSE) {
  plot(ellipse(model, which = ind), type = "l")
  points(model$coef[ind[1]], model$coef[ind[2]])
  confs <- confint(model)
  abline(v = confs[ind[1], ])
  abline(h = confs[ind[2], ])
  
  # crtamo kooridinatni pocetak
  points(0, 0, col="darkgreen", cex=1.5, lwd = 1.5)
  abline(v = 0, col="darkgreen")
  abline(h = 0, col="darkgreen")
  
  if(bonfi) {
    p <- length(model$coef)
    bonfi <- confint(model, level = 1-(0.05/p))
    abline(v = bonfi[ind[1], ], col = "blue") # za intercept crtamo vertikalne linije
    abline(h = bonfi[ind[2], ], col = "blue") # za x crtamo horizontalne linije
  }
}

# testiramo na nasem modelu, nula nije na slici, jer je daleko, sve je znacajno
draw_intervals(model)
# Klasicni primeri kad je sve beznacajno, cisto radi osecaja, npr boston
boston <- MASS::Boston
boston_model <- lm(medv ~ . , boston)
summary(boston_model)
# indus je 4. a age 8. prediktor
draw_intervals(boston_model, c(4,8)) # (0,0) pripada svim intervalima
# uzmimo age i neki znacajan
draw_intervals(boston_model, c(8,9))
# age=0 pripada intervalu poverenja za age, dok dis=0 ne pripada intervalu
# poverenja za dis => age je beznacajan, a dis znacajan.

# Sada interesantniji primeri od ranijeg casa...

# Primer kada se skracuju dva prediktora, ali Bonferoni pomogne
set.seed(216)
x1 <- runif(50, -3, 3)
x2 <- runif(50, -3, 3)
x3 <- x2 + rnorm(50, sd = 1)
# Y pravimo da ne zavisi od x2 i x3, a posto su x2 i x3 korelirani, model ce ih
# oba uzeti u obzir, ali sa suprotnim znakovima
y <- 5 + rnorm(50)
model <- lm(y ~ x2 + x3)
summary(model) # jedva znacajni su x2 i x3
# Koordinatni pocetak je izvan pravougaonika, i izvan oba intervala za
# pojedinacne koeficijente sto znaci da su znacajna oba koeficijenta (sto
# odgovara p vrednosti 0.03), ali koord. pocetak pripada elipsi koja oslikava
# F-test, koji ipak tvrdi da su oba ova koeficijenta nul (p vrednost 0.093).

# Ali kad prosirimo intervale Bonferonijevom korekcijom...
draw_intervals(model, c(2,3), bonfi = TRUE)
# Koord. pocetak upada u plavi pravougaonik i tacno zakljucujemo da su
# koeficijenti uz x2 i x3 neznacajni, sto nam je F test sve vreme govorio.


# Primer kada je sve izgleda beznacajno ali F-test kaze da je znacajno
set.seed(216)
x1 <- runif(50, -3, 3)
x2 <- x1 + rnorm(50, sd = 0.1) # x2 pravimo da bude korelisan sa x1
# y uzmemo da zavisi samo od x2
y <- 1 + 3*x2 + rnorm(50)
# Kada napravimo model i sa x1 i sa x2...
model <- lm(y ~ x1 + x2)
summary(model) # nijedan koeficijent nije znacajan, ali F-test kaze da jesu
# Kako to izgleda na slici intervala poverenja...
draw_intervals(model, c(2,3))
# Pravougaonik sadrzi koord pocetak, sto ukazuje da su oba koeficijenta
# neznacajna, ako se gledaju pojedinacni testovi. Medjutim, (0,0) je daleko od
# elipse, koja ilustruje rezultat F-testa, pa zakljucujemo da neka veza izmedju
# prediktora i zavisne promenljive postoji, tj. verujemo F-testu i zajednickom
# intervalu poverenja. Zbog visoke koreliranosti, elipsa je vrlo izduzena.

# Kada izbacimo x1 iz price, dobijemo bolje rezultate...
model2 <- lm(y ~ x2)
summary(model2) # sve je znacajno
draw_intervals(model2, c(1,2)) # nijedan interval ne sadrzi nulu

# Vratimo ipak sliku prethodnog intervala poverenja
draw_intervals(model, c(2,3))
# Vrlo je interesantna informacija koju nosi ovaj interval u odnosu na prosto
# testiranje odredjene linearne hipoteze.

# Mozemo pomocu f-je linear_hypothesis sa 4. casa da testiramo hipotezu
# H0: b1 = 1.5, b2 = 1.5
linear_hypothesis(model, matrix(c(0,0, 1,0, 0,1), nrow=2), c(1.5, 1.5))
# Velika p-vrednost - ne odbacujemo H0, sto je ocekivano.
# Medjutim, mozemo da probamo i H0: b1 = -1, b2 = 4
linear_hypothesis(model, matrix(c(0,0, 1,0, 0,1), nrow=2), c(-1, 4))
# Opet je velika p-vr, pa i ovu hipotezu ne odbacujemo.
# Interval poverenja je vrlo izduzena elipsa, koja prati pravu x2 = 3 - x1, i on
# nam daje informaciju koje sve hipoteze H0 ce proci testiranje, tako da sa
# intervala vidimo da necemo odbaciti nijednu hipotezu o koeficijentima b1 i b2
# za koju je ispunjeno b1 + b2 = 3, ukoliko je b2 (okvirno) izmedju -1 i 5
# Testirajmo i za b2 = 5 i b1 = -2
linear_hypothesis(model, matrix(c(0,0, 1,0, 0,1), nrow=2), c(-2, 5))
# p-vr oko 0.05 - ne odbacujemo H0, ali vidimo da smo pri granici intervala pov.
# S druge strane, b1=5, b2 = -2, je ipak izvan elipse poverenja 
linear_hypothesis(model, matrix(c(0,0, 1,0, 0,1), nrow=2), c(-2, 5))

# Ocekivano, i hipoteza da su b1 + b2 = 3 prolazi test (bez navodjenja tacnih
# vrednosti za koeficijente)
linear_hypothesis(model, t(c(0, 1, 1)), 3)

# Bonus: Odredimo interval poverenja za b1 + b2
l <- c(0, 1, 1)
X <- model.matrix(model)
XtXi <- solve(t(X)%*%X)
n <- nrow(X)
p <- length(model$coefficients) - 1
estim <- t(l) %*% model$coefficients
stdev <- sqrt(sum(model$res^2)/(n-p-1) * t(l) %*% XtXi %*% l)
lwr <- estim - stdev * qt(0.975, n-p-1)
upr <- estim + stdev * qt(0.975, n-p-1)
c(lwr, upr)
# Sadrzi 3, sto smo i ocekivali

# Koristeci test statistiku za linearne hipoteze C*beta = gamma, moze se
# odrediti i zajednicki interval poverenja (elipsoid) za vise linearnih
# kombinacija koeficijenata.

## Intervali predvidjanja

# Ako dobijemo nov vektor prediktora x_0, nase predvidjanje za njega ce biti
# \hat{y_0} = x_0^T * \hat{beta} Kada pravimo interval poverenja bitno je da li
# pravimo interval poverenja za ocekivanje od y_0 ili za konkretno jedan
# podatak. Idejno, preciznije mozemo da ocenimo u koji interval ce upasti prosek
# vise y_0 dobijenih za isto x_0 (zbog greske epsilon, kojom unosimo slucajnost,
# ne mora svaki y_0 biti isti za isto x_0), nego interval u koji ce upasti
# konkretno jedan y_0 za neko x_0 (zbog varijacije u samim y_0, unete greskom
# epsilon)

# Prvo cemo pokazati intervale predvidjanja za ocekivanu vrednost

# Pocecemo sa vestackim uzorkom
set.seed(888)
x <- runif(50, -3, 3)
y <- 1 + 3 * x + rnorm(50)
model <- lm(y ~ x)
plot(x, y)
abline(model)
# Interval poverenja za ocekivanje se dobija funkcijom predict, uz argument
# interval = 'confidence' Mi cemo praviti traku poverenja, tj. uzeti intervale
# za razne vrednosti x i spojiti ih linijama da dobijemo krive.
x0 <- seq(-3, 3, length.out = 50)
intervals <- predict(model, data.frame(x=x0), interval = "confidence")
lines(x0, intervals[, 'upr'], type = "l", col="red")
lines(x0, intervals[, 'lwr'], type = "l", col="red")

# Jos jedna netacna ilustracija intervala poverenja koja daje neku intuiciju

# Ako ponovimo uzorkovanje mnogo puta, svaki put dobijemo malo drugacije ocene
# koeficijenata, sto menja pravu dobijenu modelom. Ta prava u svakoj tacki x
# daje ocenu ocekivanja y za x.
# Napravicemo mnogo modela za razlicite uzorke
models <- replicate(100, {
  x <- runif(50, -3, 3)
  y <- 1 + 3 * x + rnorm(50)
  lm(y ~ x)
}, simplify = FALSE)
# nacrtamo svaki od modela
lapply(models, abline, lwd = 0.5, col = "darkgrey")
# Ove razne prave stvaraju odredjenu traku koja se siri pri krajevima, i ta
# traka je vrlo slicna pravoj traci poverenja.


# Intervali predvidjanja za buducu vrednost
plot(x, y)
abline(model)
# Interval poverenja za buducu vrednost se dobija funkcijom predict, uz argument
# interval = 'prediction'
# Mi cemo praviti traku poverenja, tj uzeti intervale za
# razne vrednosti x i spojiti ih linijama da dobijemo krive.
x0 <- seq(-3, 3, length.out = 50)
p_intervals <- predict(model, data.frame(x=x0), interval = "prediction")
lines(x0, p_intervals[, 'upr'], type = "l", col="blue")
lines(x0, p_intervals[, 'lwr'], type = "l", col="blue")
# Vidimo da su oni mnogo siri od intervala za ocekivanje:
lines(x0, intervals[, 'upr'], type = "l", col="red")
lines(x0, intervals[, 'lwr'], type = "l", col="red")

# Mozemo nacrtati i intervale poverenja za medv ~ log(lstat) iz Boston podataka
boston <- MASS::Boston
log_model <- lm(medv ~ log(lstat), boston)
plot(medv ~ lstat, boston)

x0 <- seq(0, 40, length.out = 100)
intervals <- predict(log_model, data.frame(lstat=x0), interval = "confidence")
lines(x0, intervals[, 'fit'], type = "l", col="black")
lines(x0, intervals[, 'upr'], type = "l", col="red")
lines(x0, intervals[, 'lwr'], type = "l", col="red")
p_intervals <- predict(log_model, data.frame(lstat=x0), interval = "prediction")
lines(x0, p_intervals[, 'upr'], type = "l", col="blue")
lines(x0, p_intervals[, 'lwr'], type = "l", col="blue")