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 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")