## Testiranje hipoteza o parametrima modela

# Prvo, u prostoj regresiji y = b_0 + b_1x testiramo da li je koeficijent b_1
# statisticki znacajan, tj. testiramo hipotezu H0: b_1 = 0 protiv alternative
# H1: b_1 =/= 0. U sustini, testiramo da li postoji nekakva linearna veza
# izmedju x i y ili je dovoljno da model bude samo y = b_0, tj. da li je model
# podjednako deskriptivan kao i prosek od y.

# Generisimo podatke za koje su x i y linearno povezani
set.seed(216)
x <- runif(100, -3, 3)
y <- 2.16 + 10*x + rnorm(100, sd = 5)
plot(x, y)
# Napravimo linearan model za y ~ x
model <- lm(y ~ x)
summary(model) # 3 zvezdiice, p vrednost < 2*10^{-16} => koef. uz x je znacajan
abline(model)
abline(mean(y), 0, col="red") # mnogo smo bolji od proseka, sto kaze i R^2


# Generisimo podatke tako da x i y nemaju veze jedan s drugim
set.seed(216)
x <- runif(100, -3, 3)
y <- 2.16 + rnorm(100)
plot(x, y)
# Napravimo linearan model za y ~ x
model <- lm(y ~ x)
summary(model) # beznacajan koeficijent uz x, p vrednost 0.96, nema zvezdica
abline(model)
abline(mean(y), 0, col="red") # Prakticno isti model kao i prosek, R^2 skoro 0

# Ali, testovi mogu biti i nesigurni
set.seed(216)
x <- runif(100, -3, 3)
y <- 2.16 + 0.105*x + rnorm(100)
plot(x, y)
# Napravimo linearan model za y ~ x
model <- lm(y ~ x)
summary(model)
abline(model)
# p vr. oko 5%, oko praga znacajnosti, nije ocigledno sta raditi, ali ako smo
# pred testiranje odabrali prag znacajnosti 0.05, ne odbacujemo H0: b_1 = 0, jer
# nemamo dokaza da tvrdimo da je b_1 razlicit od 0 (p vr > alfa).

# Ako povecamo uzorak, mozda dobijemo sigurnije rezultate...
set.seed(216)
x <- runif(200, -3, 3)
y <- 2.16 + 0.105*x + rnorm(200)
plot(x, y)
# Napravimo linearan model za y ~ x
model <- lm(y ~ x)
summary(model)
# Sad vec imamo statisticku znacajnost za prag znacajnosti 0.05, jer je p=0.0016
# sto je manje od 0.05. sada mozemo da odbacimo hipotezu da je b_1 = 0.
abline(model)


## Testiranje u visestrukoj regresiji

# Prvo sto testiramo je da li postoji uopste linearna veza izmedju nasih
# prediktora i zavisne promenljive. Da li je nas model sustinski isti kao
# y = b_0 + eps, tj. objasnjava y isto koliko i prosek od y?

# Testiramo hipotezu H0 da su svi koeficijenti b_1,b_2,...,b_p jednaki 0
# (ne gledamo b_0), protiv H1 da je bar jedan od njih razlicit od 0.

# Ovaj test se zove F test i rezultat ovog testiranja je u poslednjem redu u
# summary, F-statistic i p-value za nju su bitne vrednosti.

# Radi demonstracije, modeliramo prvo sum da vidimo kad je F-stat malo i 
# statisticki neznacajno, tj p-vrednost velika i ne odbacujemo hipotezu H0 da
# su svi koeficijenti (osim b_0) nule
set.seed(216)
x1 <- runif(50, -3, 3)
x2 <- runif(50, -3, 3)
x3 <- runif(50, -3, 3)
y <- 5 + rnorm(50)
model <- lm(y ~ x1 + x2 + x3)
summary(model)

# Ako je model zaista linearan, dobicemo statisticki znacajne rezultate
set.seed(216)
x1 <- runif(50, -3, 3)
x2 <- runif(50, -3, 3)
x3 <- runif(50, -3, 3)
y <- 5 + x1 + x2 + x3 + rnorm(50)
model <- lm(y ~ x1 + x2 + x3)
summary(model)

## Testiranje hipoteze za jedan koeficijent

# Vratimo se Boston skupu podataka
boston <- MASS::Boston
# napravimo prvo model sa svim prediktorima
full_model <- lm(medv ~ ., boston)
# vecina je znacajna, osim indus i age, koji imaju visoke p-vrednosti
summary(full_model)

# videli smo prosli puta da mozemo da izbacimo indus i age i da se nista bitno
# ne promeni u modelu
model2 <- lm(medv ~ . - indus - age, boston)
summary(model2)

# Medjutim, nije pametno izbacivati vise koeficijanata odjednom na osnovu
# pojedinacnih testova znacajnosti, jer se izbacivanjem jednog koeficijenta
# mogu dobiti razni efekti. Nekoliko primera sledi...

# Ukoliko bismo iz modela full_model izbacili prediktor lstat, koji je znacajan,
model_no_lstat <- lm(medv ~ . - lstat, boston)
# odjednom prediktor age, koji je malopre bio potpuno neznacajan, postaje
# statisticki znacajan, stavise sa jako malom p-vrednoscu (3 zvezdice ***)
summary(model_no_lstat)
# Neka interpretacija ovog efekta bi bila to sto, kako lstat oznacava procenat
# populacije "nize klase", izbacivanjem njega iz modela prediktor age (procenat
# kuca stare gradnje) postaje bitan, jer je ljudi sa manjim primanjima cesce
# zive u starijim zgradama, pa su lstat i age povezane i objasnjavaju isti
# efekat, pa kada se izbaci lstat, age preuzima ulogu indikatora klase.

# Dalje, interesantan je i efekat izbacivanja znacajnog prediktora rad...
model_no_rad <- lm(medv ~ . - rad, boston)
# ... nakon cega prediktor tax postaje beznacajan, iako je imao 2 zvezdice u
# punom modelu.
summary(model_no_rad)

# Tako da vidimo da se znacajnost prediktora obavezno mora gledati u okviru
# celog modela, uzimajuci u obzir njihove povezanosti, jer se izbacivanjem 
# jednog mogu bitno promeniti rezultati modela. Testovi pojedinacnih prediktora
# se ne smeju gledati kao testiranje postojanja ikakvog uticaja prediktora na
# zavisnu promenljivu, vec testiranje uticaja tog prediktora na y, pod uslovom
# da su ostali prediktori prisutni.

# Potencijalno objasnjenje zasto se gubi znacajnost tax nakon izbacivanja rad je
# opisano sledecim vestackim primerom.
set.seed(216)
# Uzimamo tri prediktora, pri cemu je x3 visoko korelisan sa x2
x1 <- runif(50, -3, 3)
x2 <- runif(50, -3, 3)
x3 <- x2 + rnorm(50, sd = 1)
# Y pravimo tako da x2 i x3 imaju suprotan uticaj na Y, tako da se njihovi
# uticaji ponistavaju
y <- 5 + x1 + x2 - x3 + rnorm(50)
# Kada modeliramo y po svim prediktorima...
model <- lm(y ~ x1 + x2 + x3)
# ... dobijamo da su sva tri znacajna, sa vrlo malim p-vrednostima
summary(model)
# Medjutim, ako izbacimo x3, tj. modelujemo samo preko x1 i x2...
model2 <- lm(y ~ x1 + x2)
# ... odjednom i x2 postaje beznacajan za model i mogli smo da uzmemo samo x1
summary(model2)
# Kada smo imali i x2 i x3 u modelu, bilo je neophodno ukljuciti i jedan i drugi
# da bi uticaj x2 mogao da bude ponisten uticajem x3, medjutim kada ne 
# ne razmatramo x3, x2 vise nema sta da ponistava, i postaje nebitan, jer samo
# donosi svoj doprinos koji se ne ponistava postojanjem x3.


# Sada cemo proci primer kada su u pocetku oba koeficijenta beznacajna, ali
# ne smemo ih zajedno izbaciti, vec izbacivanjem jednog, drugi postaje znacajan
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)
# ... oba imaju povece p-vrednosti i nisu znacajni. Interesantna je vrednost
# koeficijenata - uz x1 je 1.14 a uz x2 1.84, koji u zbiru daju oko 3, sto je
# pravi koeficijent uz x2, ali posto su visoko korelisani x1 i x2, oni 'podele'
# taj uticaj.
summary(model)

# Izbacimo li x2 iz modela...
model2 <- lm(y ~ x1)
# ...x1 pokupi sav uticaj i postaje znacajan. Dodatno, koeficijent postaje 3
summary(model2)
# Ista stvar vazi i ako izbacimo x1
model3 <- lm(y ~ x2)
summary(model3)
# Ovaj efekat se dobija jer i x1 i x2 objasnjavaju isti uticaj na y, pa iako je
# taj uticaj znacajan (sto vidimo kad uzmemo samo jedan od prediktora), kad se
# oba prediktora ukljuce u model, mogu oba postati beznacajna, ili je mogao da
# jedan od njih pokupi sav uticaj.

# Pogledamo li malo bolje rezultate modelasa x1 i x2
summary(model)
# Vidimo da, iako su oba pojedinacno beznacajna, F-test na dnu nam ipak daje
# p-vrednost skoro nula - tj. ukazuje da se ne mogu svi koeficijenti izbaciti, 
# vec da postoji jaka veza izmedju y i prediktora.
# Uz to, i R^2 nam je dosta veliki, sto ukazuje na istu stvar.

# Ovo ukazuje na korisnost F-testa i potvrdjuje pravilo da je bolje gledati
# rezultate testiranja za grupu koeficijenata nego vise pojedinacnih testova,
# jer grupni testovi bolje regulisu broj koeficijenata, a i smanjuju verovatnocu
# pravljenja greske prve vrste nego pri vise testiranja pojedinacnih prediktora.

# Znacaj F-testa naspram pojedinacnih testova se ogleda i na sledecem primeru.
# Ako napravimo 100 prediktora koji su svi sum, i y generisemo potpuno nezavisno
# od njih, pukom srecom ocekujemo da cemo dobiti oko 5 statisticki znacajnih
# koeficijenata (sa nivoom znacajnosti 5%), iako nikakva veza ne postoji.
set.seed(216)
xs <- sapply(1:100, function(i) runif(1000, -3, 3))
xs <- data.frame(xs)
y <- rnorm(1000)
data <- data.frame(cbind(xs, y))
# Kada napravimo model...
model <- lm(y ~ ., data)
# ... dobijemo upravo ocekivanih 5 koeficijenata sa bar 1 zvezdicom (sto 
# odgovara stat. znacajnom rezultatu za nivo znacajnosti 5%). Tako da bismo
# gledajuci pojedinacne testove pomislili da postoji veza izmedju y i tih 5
# prediktora, ali znamo da je sve u podacima cist sum. Jedan od prediktora je
# znacajan i na nivou 1% (2 zvezdice), iako ne postoji nikakva veza.
# Medjutim, F-test opet daje tacnu ocenu - p vrednost mu je dosta velika, sto
# nam kaze da ne mozemo da odbacimo hipotezu da su svi koeficijenti nule. Stoga
# bismo (u praksi) u ovom slucaju poverovali F-testu i napravili neku dublju
# analizu da utvrdimo postoji li veza izmedju y i nekih prediktora.
summary(model)

# Iz ovog razloga, sto ako imamo veci broj prediktora, verovatnoca da slucajno
# dobijemo znacajan rezultat postaje veca, uvodi se Bonferonijeva korekcija pri
# testiranju vise prediktora. Naime, ako testiramo p prediktora, testiranje
# pojedinacnih koeficijenata vrsimo sa nivoom znacajnosti alfa/p, umesto alfa,
# tj smanjujemo prihvatljivu gresku prve vrste na svakom testu. Prakticno, ovo
# mozemo i da postignemo ako p vrednost testa pomnozimo sa p (br. pred) i 
# posmatramo odnos sa nivoom znacajnosti alfa.

# Na primer, ako ponovimo nas model gde se skracuju dva prediktora
set.seed(216)
x1 <- runif(50, -3, 3)
x2 <- runif(50, -3, 3)
x3 <- x2 + rnorm(50, sd = 1)
# Y pravimo da zavisi samo od x1, a posto su x2 i x3 korelirani, model ce ih
# oba uzeti u obzir, ali sa suprotnim znakovima
y <- 5 + x1 + rnorm(50)
# Kada modeliramo y po svim prediktorima...
model <- lm(y ~ x1 + x2 + x3)
# ... dobijamo da su sva tri znacajna. Ali x2 i x3 imaju p-vrednosti 0.03, sto
# kada pomnozimo sa 3, sto je br. prediktora, dobijamo korigovane p-vrednosti
# od 0.09, sto vise nije znacajno sa nivoom zncajnosti 0.05.
summary(model)

## Poredjenje 2 modela - testiranje znacajnosti grupe koeficijenata

# Nekad zelimo da proverimo da li neka grupa k koeficijenata b_i1,...,b_ik moze
# da se anulira ili je statisticki znacajna. To bi na primer hteli da uradimo
# pre nego sto izbacimo age i indus iz modela za Boston.

# Testiramo H0: b_age=0, b_indus=0 protiv hipoteze da je bar jedan koef. nenula.
# Postupak je sledeci
# Definisemo ceo model
full_model <- lm(medv ~ ., boston)
# definisemo model bez indus i age
sub_model <- lm(medv ~ . - indus - age, boston)
# testiranje vrsimo funkcijom anova. Konvencija je da se prvo stavi manji model
anova(sub_model, full_model)
# Ona nam daje rezultat testiranja da li je razlika u sumi kvadrata reziduala
# izmedju punog modela i modela bez indus i age znacajna, sto se svodi na nase
# testiranje da li su i age i indus nule ili ne. Zanima nas kolona F i Pr(>F).
# Prva nam daje vrednost F statistike, a druga nam daje p-vrednost testa.
# U ovom slucaju, p vrednost je jako velika, pa ne odbacujemo H0 da su i age i
# indus neznacajni, pa je opravdan nas potez izbacivanja njih iz modela.

# Ako bismo dodali i npr. rad na listu za izbacivanje, rezultati su drukciji
sub_model2 <- lm(medv ~ . - indus - age - rad, boston)
anova(sub_model2, full_model)
# Sada je p vrednost jako mala i odbacujemo H0 da su sva tri prediktora 
# statisticki neznacajna