## Visestruka regresija

# Koristicemo podatke Boston iz MASS paketa, koji sadrze cene kuca u raznim
# delovima Bostona, uz promenljive koje mogu biti od uticaja, poput
# socioekonomske strukture stanovnistva, pristupacnosti putevima, blizini reke,
# broja soba, itd.
# Lep rad, iz kog potice ovaj skup podataka:
# Harrison, D. and Rubinfeld, D.L. (1978)
# "Hedonic prices and the demand for clean air"

?MASS::Boston
data <- MASS::Boston # ucitamo podatke u promenljivu data
summary(data) # neke sumarne statistike

# crtamo jedan po jedan grafik zavisnosti cene kuca (medv)
# od ostalih promenljivih
plot(medv ~ ., data) 

# nacrtacemo korelacionu matricu obelezja pomocu funkcije
# corrplot iz corrplot paketa
# install.packages("corrplot")
library(corrplot)
corrplot(cor(data))

cor(data)[, 'medv'] # gledamo korelacije sa medv promenljivom

# Prelazimo na modeliranje.
# Napravimo linearni model zavisnosti medv od svih ostalih obelezja 
# ('.' u formuli medv ~ . oznacava da koristimo sva obelezja)
full_model <- lm(medv ~ ., data = data)

# Pogledamo pregled modela, razmotrimo kako se interpretiraju koeficijenti,
# njihovi znakovi, norma, znacajnost itd.
summary(full_model)

# Uzimamo manji model radi demonstracije nedostataka modela
# pri razlicitim skalama podataka
model <- lm(medv ~ nox + dis + lstat, data)
summary(model) # nox ima mnogo veci koeficijent od otalih (osim intercepta)

# lstat je izrazen u procentima, uzecemo sada lstat/100 kao prediktor, da bismo
# dobili lstat izrazen kao proporcija, cime smanjujemo njegovu skalu.
# ( I(...) u "I(lstat/100)" se koristi kada transformisemo prediktore, da bi lm
# znao da treba da prvo primeni transformaciju, pa onda modelira. )
model2 <- lm(medv ~ nox + dis + I(lstat/100), data)
summary(model2) # sada je koeficijent uz lstat/100 ubedljivo najveci

# Znaci, za sustinski isti model dobili smo mnogo razlicite rezultate i moguce
# interpretacije, samo zbog skale na kojima su podaci. Zato se cesto svi podaci
# standardizuju (sa (X-m)/sigma, tj. (kolona-prosek_kolone)/stdev_kolone))
scaled <- data.frame(scale(data)) # scale daje matrix, pa menjamo u data.frame
# Obicno se ne standardizuje zavisna promenljiva y (kod nas 'medv').
scaled[, 'medv'] <- data[, 'medv']
# Prethodni model sa standardizovanim podacima
scaled_model <- lm(medv ~ nox + dis + lstat, data = scaled)
summary(scaled_model) # svi koeficijenti su istog reda velicine

# Nakon standardizacije, ((X^T)*X)/(n-1)  je jednaka korelacionoj matrici
scaled_full_model <- lm(medv ~ ., data = scaled)
X <- model.matrix(scaled_full_model)
(t(X)%*%X)[2:5, 2:5] / (506-1)
cor(data)[1:4, 1:4]

# Gledamo indikatorsku promenljivu, koja je interpretacija koeficijenta?
summary(full_model)
x1 <- data[5, ]
x2 <- x1
x2$chas <- TRUE # setujemo chas na 1 (TRUE)
# predvidjanje modela za x2 je 30.63026, a za x1 je 27.94352, dok je razlika
# 30.63026 - 27.94352 = 2.68674 ~ 2.687, sto je i koeficijent uz chas.
predict(full_model, rbind(x1,x2))


## Igranje sa R^2
# R^2 punog modela je 0.7406, dok je adjusted R^2 0.7338
summary(full_model) 

# Ako izbacimo neznacajne prediktore - age i indus, ostaje isti R^2
# (u formuli medv ~ . - indus - age, "oduzimanje" predstavlja izbacivanje
# prediktora, tj uzimamo sve (".") bez age ("-age") i bez indus ("-indus") )
model2 <- lm(medv ~ . - indus - age, data)

# R^2 je isti, ali adjusted R^2 se povecao, zato sto smo smanjili broj
# prediktora koje koristimo
summary(model2)
# R^2 ne opada sa brojem prediktora, ali popravljeni R^2 kaznjava koriscenje 
# velikog broja prediktora.

# Pogledacemo neke transformacije kojima mozemo da dobijemo bolji model
# Radi jednostavnosti, posmatramo samo zavisnost medv od lstat
plot(medv ~ lstat, data)
# provucemo linearni model, R^2 = 0.5441, adjR^2 = 0.5432
lin_model <- lm(medv ~ lstat, data)
summary(lin_model)
abline(lin_model, col="blue")

# Ovde zavisnost ocigledno nije linearna, pa dodajemo i kvadratni clan
# Ovde je model medv = a + b*lstat + c*lstat^2
# U formuli medv ~ ..., "+" oznacava dodavanje prediktora
quad_model <- lm(medv ~ lstat + I(lstat^2), data)
# R^2 = 0.6407, adjR^2 = 0.6393 (bolji od prave u svakom smislu)
summary(quad_model)
# Crtanje dobijene ocene modela radimo tako sto uzmemo niz xn tacaka za
# lstat-osu i uzmemo ocene dobijene modelom za te vrednosti lstat, pa crtamo 
# pomocu funkcije lines
xn <- seq(0, 40, 0.5)
# funkcija predict prima kao prvi argument model na osnovu kog predvidjamo, a
# kao drugi data.frame (obavezno data.frame) koji mora da sadrzi kolone sa 
# imenima prediktora (desna strana medv ~ ...)
yn <- predict(quad_model, data.frame(lstat = xn))
# Crtamo funkciju
lines(xn, yn, col="darkgreen") # lepse ocenjuje od prave

# Probamo i polinom 5-og stepena
# Formula moze da bude medv ~ lstat + I(lstat^2) + (stepeni 3,4) + I(lstat^5),
# ali krace je medv ~ poly(lstat, 5, raw=TRUE)
# U funkciji poly, 5 je stepen polinoma, a raw=TRUE oznacava da koristimo
# obicne polinome 1,x,x^2,... kao prediktore
poly_model <- lm(medv ~ poly(lstat, 5, raw=TRUE), data)
# R^2 = 0.6817, adjR^2 = 0.6785 (bolji i od kvadratne)
summary(poly_model)
yn <- predict(poly_model, data.frame(lstat = xn))
lines(xn, yn, col="red")

# Dobar model je i ako koristimo log(lstat) kao prediktor
log_model <- lm(medv ~ lstat + log(lstat), data)
summary(log_model)
yn <- predict(log_model, data.frame(lstat = xn))
lines(xn, yn, col="cyan", lwd=2)

# Vratimo se na pocetni veliki model, uzmemo log(lstat) pored lstat
new_model <- lm(medv ~ . - indus - age + log(lstat), data)
summary(new_model) # R^2 smo sa 0.74 povecali na 0.79

## Pretpostavke modela:
# 1. E(greske) = 0
# 2. cov(greske) = sigma*I
# 3. greske normalno raspodeljene

# Prvo proverimo normalnost reziduala
res <- new_model$residuals

# normalnost proveravamo qqplot-om
qqnorm(res) # qqplot
qqline(res) # teoretska linija koju treba da prate reziduali
# Reziduali nisu bas normalni, pa p-vrednosti testova mozda nisu precizne

# Uslove Gaus-Markova proveravamo grafikom zavisnosti reziduala od ocena y
plot(new_model$fitted.values, new_model$residuals)
abline(0, 0)
# Ocekivanje deluje da je blisko nuli, tj. da su reziduali centrirani oko nule
# Postoji blaga heteroskedasticnost, disperzija oko sredine grafika je nesto
# manja nego na krajevima
# Za uslove Gaus-Markova mozemo reci da su zadovoljeni, pa mozemo verovati
# ocenama standardnih gresaka i ocenama koeficijenata.


## ZASTO SU USLOVI GAUS-MARKOVA BITNI? Jer je pod njima ocena MNK BLUE.
# Pokazacemo neke efekte koji su moguci ako nisu ispunjeni uslovi G-M

# PRVI USLOV G-M: ocekivanje greske je 0
# Napravicemo vestacki model, y = 1 + 3x + eps, gde eps ima normalnu raspodelu
# sa ocekivanjem 5sin(x), a ne 0.
set.seed(216)
n <- 100
x <- seq(-2, 2, length.out = n)
y <- 1 + 3*x + rnorm(n, 5*sin(x), 1)
plot(x, y)
model <- lm(y ~ x)
summary(model)
abline(model)
# Ovaj model ne izgleda mnogo lose, ali pogledajmo efekat koji smo dobili
# na koeficijente modela, buduci da prvi G-M uslov nije ispunjen.
# Iznova cemo 1000 puta generisati uzorak (x,y) i provuci linearni model i uzeti
# koeficijente dobijenog modela.
model_coefs <- replicate(100, {
  x <- seq(-2, 2, length.out = n)
  y <- 1 + 3*x + rnorm(n, 5*sin(x), 1)
  model <- lm(y ~ x)
  model$coefficients
})
# Dobijamo 2*1000 matricu koeficijenata svakog od simuliranih modela.
model_coefs
# Ocekivanje (po zakonu v.b.) koeficijenata b_0 i b_1 je priblizno jednako
# proseku dobijenih koeficijenata iz simulacija
rowMeans(model_coefs)
# Ocekivanje ocena beta (\hat{beta}) dobijamo da je (1, 4.6) =/= (1, 3)

# odstupanje po formuli pristrasnosti koeficijenata beta je isto kao proseci
X <- model.matrix(model)
Eeps <- 5*sin(X[,2])
solve(t(X)%*%X)%*%t(X)%*%Eeps # = E(\hat{beta}) - beta, izvesti
# odstupanje proseka
rowMeans(model_coefs) - c(1, 3) # 1 i 3 su prave vrednosti koef

# Dakle, kada nije E(greske) = 0, i ocene beta su pristrasne

# Da smo nacrtali reziduale, videli bismo da postoji neko oscilovanje
plot(model$fit, model$res)
abline(0, 0)

# popravka, novi prediktor sin(x)
plot(x, y)
model <- lm(y~x+sin(x))
summary(model)
# crtamo predvidjanja
xn <- seq(-2, 2, 0.1)
yhat <- predict(model, data.frame(x = xn))
lines(xn, yhat)

# simuliramo model 1000 puta i uzimamo koeficijente
model_coefs <- replicate(100, {
  x <- seq(-2, 2, length.out = n)
  y <- 1 + 3*x + rnorm(n, 5*sin(x), 3)
  model <- lm(y~x+sin(x))
  model$coefficients
})
#model_coefs
rowMeans(model_coefs) # EB = (1,3,5) -> nepristrasna ocena B
# reziduali su sad lepsi
plot(model$fit, model$res)
abline(0, 0)


# TRECI USLOV G-M - Nekoreliranost gresaka
# Demonstriracemo efekat koji se dobija ako prosto ponovimo podatke 200 puta
# Model ostaje isti, jer nemamo nikakve nove podatke, ali smo zbog velicine
# podataka pouzdaniji u svoje ocene, pa svi koeficijenti postaju znacajni
# iako zapravo nisu (pored toga, i ocene su pogresne, ali zbog velicine podataka
# imamo laznu sigurnost u njih).
set.seed(216)
n <- 42
x <- seq(-3, 3, length.out = n)
y <- 1 + 0.1*x + rnorm(n)
plot(x, y)
model <- lm(y ~ x)
summary(model)
abline(model)

# isti model, ali svi koef sada znacajni
x2 <- rep(x, 200)
y2 <- rep(y, 200)
plot(x2, y2)
model2 <- lm(y2~x2)
summary(model2)
abline(model2)