library(ISLR2) ?Hitters dim(Hitters) sum(is.na(Hitters$Salary)) Hitters <- na.omit(Hitters) dim(Hitters) ## best subset selection library(leaps) ?regsubsets regfit.full <- regsubsets(Salary ~ ., Hitters, nvmax=19, method="exhaustive") summary(regfit.full) # npr. CRBI se nalazi u prvih 6 modela, ali ne i u sedmom. ## npr. kod FSS to nikad nije slucaj. # ?plot.regsubsets # plot(regfit.full, scale="bic") ## forward/backward stepwise selection regfit.fwd <- regsubsets(Salary ~ ., data = Hitters, nvmax = 19, method = "forward") summary(regfit.fwd) regfit.bwd <- regsubsets(Salary ~ ., data = Hitters, nvmax = 19, method = "backward") summary(regfit.bwd) ### npr. predstanici modela sa kovarijata su coef(regfit.full, 7) coef(regfit.fwd, 7) coef(regfit.bwd, 7) ## unakrsna validacija za bss (k=10) predict.regsubsets <- function(object, newdata, broj_kovarijata, ...) { form <- as.formula(object$call[[2]]) mat <- model.matrix(form, newdata) coefi <- coef(object, id = broj_kovarijata) xvars <- names(coefi) mat[, xvars] %*% coefi } ## ne postoji predict funkcija za regsubsets(), ## buduci da su objekti koje koristimo klase regsubsets, ## potrebno je samo napisati predict n <- nrow(Hitters) k <- 10 ## broj grupa ("folds") set.seed(1) folds <- sample(rep(1:k, length = n)) ## replace=FALSE -> npr. folds[1]=10 znaci da ## prvi element skupa za ucenje stavljamo u 10-ti blok itd. cv.errors <- matrix(NA, k, 19) ### hiperparametar kompleksnosti je broj kovarijata u modelu for (j in 1:k) { best.fit <- regsubsets(Salary ~ ., data = Hitters[folds != j, ], nvmax = 19) ### kada izbacujemo j-ti blok, ### bitno je da odabir najboljeg modela ### za svaki broj kovarijata (sto je dio treninga) ### radimo na podacima bez j-tog bloka! ### Na taj nacin j-ti blok zaista igra ulogu nezavisnog testnog skupa. for (i in 1:19) { pred <- predict(best.fit, Hitters[folds == j, ], broj_kovarijata = i) cv.errors[j, i] <- mean((Hitters$Salary[folds == j] - pred)^2) ## sadrzi testnu gresku na j-tom bloku ## za najbolji model s i varijabli (na skupu bez j-tog bloka) ## kako bi dobili CV procjenu testne greske ## za najbolji model s i varijabli, uzimamo prosjek testnih ## gresaka na svih k blokova (vidi dolje) } } mean.cv.errors <- apply(cv.errors, 2, mean) plot(mean.cv.errors, type = "b") p0=which.min(mean.cv.errors) p0 ## 10 reg.best <- regsubsets(Salary ~ ., data = Hitters, nvmax = 19) coef(reg.best, p0) ## slicno bi bilo za forward i backward ss ############################# ridge i lasso ############# x <- model.matrix(Salary ~ ., Hitters)[, -1] ## generira matricu podataka X y <- Hitters$Salary dim(x) ######## ridge library(glmnet) ?glmnet grid <- 10^seq(10, -2, length = 100) ridge.mod <- glmnet(x, y, alpha = 0, lambda = grid) ## alpha=0 odgovara ridge regresiji, a alpha=1 je lasso (vidi help) ## ako ne zadamo vrijednosti za lambda, odredit ce se automatski ## (tako inace radimo) ## po deafultu se kovarijate standardiziraju (inace je kazna "nepostena") dim(coef(ridge.mod)) ## za svaki od 100 vrijednosti od lambda imamo 19+1=20 koeficijenata plot(ridge.mod, xvar="lambda", label=TRUE) ## veci lambda znaci vece smanjenje plot(ridge.mod, xvar="norm", label=TRUE) ## x-os je L1 norma koeficijenata (veca norma manji lambda) ## predict radi slicno kao i prije, ## npr. mozemo dobiti koeficijente za novu vrijednost za lambda, npr. 50, predict(ridge.mod, s = 50, type = "coefficients")[,1] ?predict.glmnet ## testna greska set.seed(1) train <- sample(1:nrow(x), nrow(x) / 2) test <- (-train) y.test <- y[test] # npr. za lambda=4 ridge.mod <- glmnet(x[train, ], y[train], alpha = 0, lambda = grid, thresh = 1e-12) ## zadnji parametar? (vidi help) ridge.pred <- predict(ridge.mod, s = 4, newx = x[test, ]) sqrt(mean((ridge.pred - y.test)^2)) # opcenito, koristimo unakrsnu validaciju set.seed(1) cv.out <- cv.glmnet(x[train, ], y[train], alpha = 0) plot(cv.out) ## na gornjoj x-osi su stupnjevi slobode, ## a dvije vertikalne crte oznacavaju ## lambda koji minimizira CV gresku, ## odnosno lambda dobiven uz one SE rule (veci lambda manja kompleksnost!) ?cv.glmnet ## vidi Value ?plot.cv.glmnet bestlam <- cv.out$lambda.min bestlam ## ili # bestlam <- cv.out$lambda.1se # bestlam ridge.pred <- predict(ridge.mod, s = bestlam, newx = x[test, ]) sqrt(mean((ridge.pred - y.test)^2)) # [1] 373.9741 summary(Hitters$Salary) ## mean = 535.9 sd(Hitters$Salary) ## 451.1187 # odabrani model je dakle out <- glmnet(x, y, alpha = 0) predict(out, type = "coefficients", s = bestlam)[,1] ########################## lasso lasso.mod <- glmnet(x[train, ], y[train], alpha = 1, lambda = grid) plot(lasso.mod, xvar="lambda", label=TRUE) ## za razliku od ridge regresije, lasso radi odabir kovarijata ridge.mod <- glmnet(x[train, ], y[train], alpha = 0, lambda = grid, thresh = 1e-12) plot(ridge.mod, xvar="lambda", label=TRUE) ## lambde nisu usporedive, ali stupnjevi slobode jesu # plot(lasso.mod, xvar="norm", label=TRUE) # plot(ridge.mod, xvar="norm", label=TRUE) set.seed(1) cv.out <- cv.glmnet(x[train, ], y[train], alpha = 1) plot(cv.out) bestlam <- cv.out$lambda.min lasso.pred <- predict(lasso.mod, s = bestlam, newx = x[test, ]) sqrt(mean((lasso.pred - y.test)^2)) # [1] 379.043 ## nesto losije nego ridge regresija, ## s tim da koristimo manje stupnjeva slobode (11 umjesto 19) ## te usput radimo i odabir kovarijata out <- glmnet(x, y, alpha = 1, lambda = grid) lasso.coef <- predict(out, type = "coefficients", s = bestlam)[,1] lasso.coef lasso.coef[lasso.coef != 0] ## tocno 11 nenul koeficijenata, a to je nepristran procjenitelj ## za stupnjeve slobode ## obje metode daju puno bolje rezultate nego lin. regresija ## sa metodom najmanjih kvadrata (lambda=0) ls=glmnet(x[train, ], y[train], alpha = 1, lambda = 0) ls_pred=predict(ls, s = 0, newx = x[test, ]) sqrt(mean((ls_pred - y.test)^2)) ## [1] 408.8121 ## DZ - na trening skupu provesti best ss metodu, ## procijeniti testnu gresku odabranog modela na testnom skupu, ## te usporediti sa ridge regresijom i lasso metodom