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=20) summary(regfit.full) # npr. CRBI se nalazi u prvih 6 modela, ali ne i u sedmom. kod FSS to nije slucaj. # ?plot.regsubsets # plot(regfit.full, scale="bic") ## forward 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) ## za p=1,..,6 sve metode daju iste modele, za p=7 to nije slucaj coef(regfit.full, 7) coef(regfit.fwd, 7) coef(regfit.bwd, 7) ## unakrsna validacija za bss (k=10) predict.regsubsets <- function(object, newdata, id, ...) { form <- as.formula(object$call[[2]]) mat <- model.matrix(form, newdata) coefi <- coef(object, id = id) 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 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 cv.errors <- matrix(NA, k, 19) 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, ], id = 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) library(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? ridge.pred <- predict(ridge.mod, s = 4, newx = x[test, ]) 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, ]) mean((ridge.pred - y.test)^2) # [1] 139856.6 # 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, ]) mean((lasso.pred - y.test)^2) # [1] 143673.6 ## slicno kao kod ridge regresije, 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 (vidi predavanja) ## obje metode daju puno bolje rezultate nego lin. regresija sa metodom najmanjih kvadrata (lambda=0) ## DZ - na trening skupu provesti best ss metodu, procijeniti testnu gresku odabranog modela na testnom skupu ## te usporediti sa ridge regresijom i lasso metodom ## PCR i PLS library(pls) set.seed(2) pcr.fit <- pcr(Salary ~ ., data = Hitters, scale = TRUE, validation="CV") summary(pcr.fit) validationplot(pcr.fit, val.type = "MSEP") set.seed(1) pcr.fit <- pcr(Salary ~ ., data = Hitters, subset = train, scale = TRUE, validation = "CV") validationplot(pcr.fit, val.type = "MSEP") which.min(pcr.fit$validation$PRESS) ## ali i manje se cini dovoljno dobro pcr.pred <- predict(pcr.fit, x[test, ], ncomp = 5) mean((pcr.pred - y.test)^2) ## [1] 142811.8 -> slicno kao ridge i lasso ## odabrani model je pcr.fit <- pcr(y ~ x, scale = TRUE, ncomp = 5) summary(pcr.fit) ## PLS vidjeti sami