library(ISLR2)
attach(Wage)
?Wage

agelims <- range(age)
age.grid <- seq(from = agelims[1], to = agelims[2])
plot(age, wage, xlim = agelims, cex = .5, col = "darkgrey")
## -> nelinearna veza

### Polinomijalna regresija


fit <- lm(wage ~ poly(age, 4), data = Wage)   ## poly generira kovarijate h_0(age),...,h_4(age)
## gdje je {h_i} neka ortogonalna baza za prostor 4d polinoma.
coef(summary(fit)) ## lin. model kao i svaki drugi, samo nas ne zanimaju koeficijenti

## alternativno: 
fit2 <- lm(wage ~ age + I(age^2) + I(age^3) + I(age^4),
           data = Wage)
coef(fit2)
## ipak, predikcije su iste u oba modela

agelims <- range(age)
age.grid <- seq(from = agelims[1], to = agelims[2], by=0.1)
preds <- predict(fit, newdata = list(age = age.grid), se.fit = TRUE)
se.bands <- cbind(preds$fit + 2 * preds$se.fit, preds$fit - 2 * preds$se.fit)

# par(mfrow = c(1, 2), mar = c(4.5, 4.5, 1, 1),
#     oma = c(0, 0, 4, 0))
plot(age, wage, xlim = agelims, cex = .5, col = "darkgrey")
title("Degree-4 Polynomial", outer = T)
lines(age.grid, preds$fit, lwd = 2, col = "blue")
matlines(age.grid, se.bands, lwd = 1, col = "blue", lty = 3)
## -> siri intervali pri krajevima (manje podataka)

## izgleda kao da podaci dolaze iz dvije populacije (wage > ili < od 250) 
## -> kasnije cemo ukljuciti dodatne kovarijate, sada probajmo:

### Logisticka regresija -- model je P(wage>250 | age) = h(b0+...+b4*age^4) =: h(x^t b)
## gdje je h(z)=exp(z)/(1+exp(z))

fit <- glm(I(wage > 250) ~ poly(age, 4), data = Wage, family = binomial)


preds <- predict(fit, newdata = list(age = age.grid), se = T) ## type="link" po defaultu, daje x^t b_hat
pfit <- exp(preds$fit) / (1 + exp(preds$fit)) ## ili direktno ako stavimo type="response" u predict
se.bands.logit <- cbind(preds$fit + 2 * preds$se.fit,preds$fit - 2 * preds$se.fit)
## -> intervali pouzdanosti za x^t b_hat
se.bands <- exp(se.bands.logit) / (1 + exp(se.bands.logit))

plot(age, I((wage > 250)/5), xlim = agelims, type = "n",ylim = c(0, .2)) ## y-os skalirana na 0.2
points(jitter(age), I((wage > 250)/5), cex = .5, pch = "|", col = "darkgrey")
lines(age.grid, pfit, lwd = 2, col = "blue")
matlines(age.grid, se.bands, lwd = 1, col = "blue", lty = 3) 
### dosta siroki intervali (procjenjujemo jako malu vjerojatnost)

## splajnovi
library(splines)
fit <- lm(wage ~ bs(age, knots = c(25, 40, 60)), data = Wage)
?bs
## 3 cvora, d=3 -> 3+3+1=7 baznih funkcija (bs izostavlja slobodni clan-> to je bitno kod aditivnih modela)
dim(bs(age, knots = c(25, 40, 60)))

## -> to je opet obican linearan model sa kovarijatama izvedenim iz age
summary(fit)

pred <- predict(fit, newdata = list(age = age.grid), se = T)
plot(age, wage, cex = .5, col = "darkgrey")
lines(age.grid, pred$fit, lwd = 2)
lines(age.grid, pred$fit + 2 * pred$se, lty = "dashed")
lines(age.grid, pred$fit - 2 * pred$se, lty = "dashed")
abline(v= c(25, 40, 60), lty="dashed", col ="darkred") ## cvorovi

## umjesto bs(age, knots = c(25, 40, 60)) mogli smo zadati stupnjeve slobode
dim(bs(age, df = 6))
attr(bs(age, df = 6), "knots") 
## -> cvorovi postavljeni na kvantile od age (mozda prirodnije jer se prilagodava podacima)


## 1-splajnovi i step funkcije
plot(age, wage, cex = .5, col = "darkgrey")
lines(age.grid, pred$fit, lwd = 1)

# neprekidna po dijelovima linearna aproksimacija
fit1 <- lm(wage ~ bs(age, knots = c(25, 40, 60), degree=1),  data = Wage)
pred1 <- predict(fit1, newdata = list(age = age.grid), se = T)
lines(age.grid, pred1$fit, lwd = 1, col="darkblue")

# step funkcije
fit0 <- lm(wage ~ cut(age, c(agelims[1],25, 40, 60,agelims[2])), data = Wage)
pred0 <- predict(fit0, newdata = list(age = age.grid), se = T)
lines(age.grid, pred0$fit, lwd = 1, col="darkorange")


coef(summary(fit0))
## model je E[wage | age] = b0 + b1*1{age iz (25,40]} + ... + b3*1{age iz (60,80]}
## -> dakle b0 = prosjecna placa ljudi s age<=25 (bazna kategorija), b1 = prosjecna razlika u placi
## za ljude izmedu 25 i 40, u odnosu na one s <= 25 godina, itd. 

## prirodni splajn -- linearna aproksimacija na rubnim regijama (smanjenje varijance)
fit2 <- lm(wage ~ ns(age, df = 6), data = Wage) 
dim(ns(age, df = 6))
pred2 <- predict(fit2, newdata = list(age = age.grid),
                 se = T)
plot(age, wage, cex = .5, col = "darkgrey")
lines(age.grid, pred$fit, lwd = 2)
lines(age.grid, pred$fit + 2 * pred$se, lty = "dashed")
lines(age.grid, pred$fit - 2 * pred$se, lty = "dashed")
lines(age.grid, pred2$fit, col = "darkblue", lwd = 2)
lines(age.grid, pred2$fit + 2 * pred2$se, lty = "dashed", col = "darkblue")
lines(age.grid, pred2$fit - 2 * pred2$se, lty = "dashed", col = "darkblue")
### -> za isti broj stupnjeva slobode, veca fleksibilnost u sredisnjim regijama (tamo je vise podataka), 
### a opet manja varijanca na rubovima

### Smoothing spline (ponovno prirodan splajn, ali postoji samo jedan parametar koji kontrolira fleksibilnost 
### -> CV, i to LOOCV)
plot(age, wage, xlim = agelims, cex = .5, col = "darkgrey")
title("Smoothing Spline")
fit <- smooth.spline(age, wage, df = 16)
fit2 <- smooth.spline(age, wage, cv = TRUE) ## lambda, ili ekvivalentno df, odabrani unakrsnom validacijom

fit2$df # [1] 6.794596

lines(fit, col = "red", lwd = 2)
lines(fit2, col = "blue", lwd = 2)
legend("topright", legend = c("16 DF", "6.8 DF"), col = c("red", "blue"), lty = 1, lwd = 2, cex = .8)

## predict funkcija malo drugacija (nema newdata)
x=c(60,20, 40)
predict(fit2, x = x)

## Intervali pouzdanosti? Probajmo bootstrap!
## B puta cemo napraviti bootstrap uzorak duljine n iz originalnih (age_i,wage_i)-eva
## Aproksimirat cemo bootstrap uzorak sa smoothing splajnom te izracunati procjene u tockama iz age.grid
## Za svaku tocku iz age.grid na taj nacin cemo dobiti uzorak bootstrap procjenitelja duljine B (vidi poglavlje o bootstrapu)

n=length(age)
B=1000
boot.pred=matrix(NA, nrow=B, ncol=length(age.grid))
for (i in 1:B) {
  ind <- sample(1:n, size = n, replace = TRUE)
  x=age[ind]
  y=wage[ind]
  fit.boot=smooth.spline(x, y, cv = TRUE)
  boot.pred[i,]=predict(fit.boot, x= age.grid)$y 
}

boot.se=apply(boot.pred, 2, sd)

plot(age, wage, xlim = agelims, cex = .5, col = "darkgrey")
lines(age.grid, predict(fit2, x=age.grid)$y, lwd=2, col="darkred")
lines(age.grid, predict(fit2, x=age.grid)$y + 2*boot.se, lty="dashed", col="darkred")
lines(age.grid, predict(fit2, x=age.grid)$y - 2*boot.se, lty="dashed", col="darkred")

## drugi oblik bootstrap intervala pouzdanosti  koji se ne zasniva na pretpostavci normalnosti (*)
q1=apply(boot.pred, 2, quantile, probs=0.975)
q2=apply(boot.pred, 2, quantile, probs=0.025)
lines(age.grid, 2*predict(fit2, x=age.grid)$y  - q2, lty="dashed", col="darkblue")
lines(age.grid, 2*predict(fit2, x=age.grid)$y  - q1, lty="dashed", col="darkblue")
## -> slicni rezultati


### Lokalna regresija
## u praksi se velicina okoline tipicno zadaje preko tzv. span parametra
## to je postotak tocaka iz skupa za trening koje koristimo za prilagodbu lin. modela
## u nasoj notaciji to dobijemo ako stavimo
## h_lambda(x)=udaljenost od x do span*n najblizeg susjeda od x u {x_i:i=1,...,n}
## npr. u kNN metodi span je k/n -> span*n=k

plot(age, wage, xlim = agelims, cex = .5, col = "darkgrey")
title("Local Regression")
fit <- loess(wage ~ age, span = 0.05, data = Wage, degree = 1) 
## moze i kvadratna funkcija umjesto linearne, a za degree = 0 dobijemo NW procjenitelj
fit2 <- loess(wage ~ age, span = 0.2, data = Wage, degree = 1)
lines(age.grid, predict(fit, data.frame(age = age.grid)),
      col = "darkred", lwd = 2)
lines(age.grid, predict(fit2, data.frame(age = age.grid)),
      col = "darkblue", lwd = 2)
legend("topright", legend = c("Span = 0.05", "Span = 0.2"),
       col = c("red", "blue"), lty = 1, lwd = 2, cex = .8)

## gore se koristi tricube jezgra, tj. K_l(x0,t)=D((x0-t)/h_l(x0)) je 
tricubic=function(t){
  (1 - abs(t)^3)^3 * (abs(t)<=1)
}
t=seq(6,14, by=0.01)
plot(t, tricubic((10-t)/3), type="l") ## x0=10, h_l(x0)=3 
lines(t, tricubic((10-t)/1), lty="dashed") ## x0=10, h_l(x0)=1
## -> sirok "plato" oko t approx x0=10.


### usporedba s NW procjeniteljem (problem pristranosti na rubovima domene)

set.seed(1)
x<- runif(50, max = 3)
y <- 9 - x^2 + rnorm(50, sd = 0.3)
plot(x, y)
curve(9 - x^2, col = "grey", add = TRUE, lwd = 3) ## stvarna regresijska funkcija
grid.x <- seq(from = 0, to = 3, length.out = 300)
lok_lin <- loess(y ~ x, span=0.2, degree=1) ## 0.2*50=10
lines(grid.x, predict(lok_lin, newdata = grid.x), col="darkblue")
NW <- loess(y ~ x, span=0.2, degree=0) ## probaj span=0.1 -> manji problem na rubu, ali veci u unutrasnjosti
lines(grid.x, predict(NW, newdata = grid.x), col="darkorange")
## lok.lin. regresija tipicno je glada od NW procjenitelja 
## te ima manje problema s pristranosti na rubu (vidi desni rub!)

### GAMovi
## model je E[wage | year, age, education]=a + f1(year) + f2(age) + f3(education),
## pri cemu su f1 i f2 proizvoljne (potencijalno nelinearne) funkcije koje cemo procijeniti iz podataka 
## (ovo je jedna od glavnih prednosti GAMova), a 
## f3(education)=beta_i ako je education=i, za i=1,..,5

## ako cemo za procjenu funkcija f1 i f2 koristiti npr. prirodne splajnove, to je zapravo standardan linearan
## model sa tranformiranim kovarijatama (onda je automatski dostupna i cijela teorija):
gam1 <- lm(wage ~ ns(year, 4) + ns(age, 5) + education, data = Wage)
par(mfrow = c(1, 3))
plot.Gam(gam1, se = TRUE, col = "blue") ## crta fk(t) za k=1,2,3 (funkcije se centriraju
## tako da je prosjecna vrijednost na podacima=)
##Intepretacija?




## ako ne mozemo koristiti linearan model (npr. smoothing splajnovi ili lokalna regresija):
library(gam)
?gam

gam.m3 <- gam(wage ~ s(year, 4) + s(age, 5) + education, data = Wage) 
## s = smoothing splajn, 4 i 5 su stupnjevi slobode

par(mfrow = c(1, 3))
plot(gam.m3, se = TRUE, col = "blue")

## f1 ne izleda daleko od linearne -> to mozemo testirati

gam.m1 <- gam(wage ~ s(age, 5) + education, data = Wage)
gam.m2 <- gam(wage ~ year + s(age, 5) + education,data = Wage)
anova(gam.m1, gam.m2, gam.m3, test = "F") 
## H0 je uvijek da je manji model dovoljan:
## dodavanje linearne funkcije od year izgleda znacajno (tj. statisticki znacajno smanjuje RSS).
## s druge strane, cini se da nije potrebna nelinearna funkcija za year

summary(gam.m3) ## radi slicne testove

## predict radi kao i ranije

##umjesto smoothing splajnova, mozemo koristiti i lokalnu regresiju
gam.lo <- gam(
  wage ~ s(year, df = 4) + lo(age, span = 0.7, degree = 1) + education,
  data = Wage
)
## lo je kratica za loess
plot.Gam(gam.lo, se = TRUE, col = "darkblue")


## mozemo dodati i interakcije:
gam.lo.i <- gam(wage ~ lo(year, age, span = 0.5) + education, data = Wage)
library(akima)
par(mfrow=c(1,1))
plot(gam.lo.i)

### aditivna logisticka regresija
gam.lr <- gam(
  I(wage > 250) ~ year + s(age, df = 5) + education,
  family = binomial, data = Wage   ### kao kod obicne log. regresije (glm)
)
par(mfrow = c(1, 3))
plot(gam.lr, se = T, col = "darkblue")

## zapravo nema wage>250 u <HS kategoriji:
table(education, I(wage > 250))

## izbacimo tu kategoriju
gam.lr.s <- gam(
  I(wage > 250) ~ year + s(age, df = 5) + education,
  family = binomial, data = Wage,             
  subset = (education != "1. < HS Grad")
)
plot(gam.lr.s, se = T, col = "darkblue")