library(ggplot2) library(dplyr) library(glmnet) # glmnet regularization (after 1., there is no order, but numbers are easier to reference): # 1. generalizes * regression, where * is a distribution; see argument family of ?glmnet # 2. can (which does not mean will) reduce increased std. error of beta due to multicollinearity # 3. can (which does not mean will) fit a model with more predictors than data # 4. can (which does not mean will) see through useless predictors # 5. can (which does not mean will) reduce MSE(\hat{y}) # 6. can (which does not mean will) reduce MSE(\hat{\beta}) # fake data ---- ## re 4. N <- 51 x <- seq(0, 1, length.out = N) y <- -2 + 5*x + rnorm(N) / 5 df <- data.frame(x = x, y = y) ggplot(df, aes(x, y)) + geom_point() # don't do this! there are better ways to do this sort of idea # see Reference 2. Chapter 5 fit <- lm(y ~ x + I(x^2) + I(x^3) + I(x^4) + I(x^5) + I(x^6) + I(x^7) + I(x^8) + I(x^9) + I(x^10) + I(x^11) + I(x^12) + I(x^13) + I(x^14) + I(x^15) + I(x^16) + I(x^17) + I(x^18) + I(x^19) + I(x^20), data = df) df$yhat <- predict(fit) beta <- coef(fit) ggplot(df) + geom_point(aes(x, y), alpha = 0.1) + geom_line(aes(x, yhat), linewidth = 3, color = "blue") X <- model.matrix(fit)[, -1] y <- df$y fitr <- cv.glmnet(X, y, standardize = FALSE) betar <- coef(fitr, s = "lambda.1se") df$yhatr <- predict(fitr, newx = X, s = "lambda.1se") ggplot(df) + geom_point(aes(x, y), alpha = 0.1) + geom_line(aes(x, yhat), color = "blue") + geom_line(aes(x, yhatr), size = 2, color = "orange") # abalone ---- ## re 1. ---- df <- read.csv("https://roualdes.sfo3.digitaloceanspaces.com/data/abalone.csv") fit <- lm(age ~ length + diameter + height + whole_weight + shucked_weight + viscera_weight + shell_weight, data = df) coef(fit) X <- model.matrix(fit)[,-1] y <- df$age fitr <- glmnet(X, y, lambda=0, standardize = FALSE, thresh = 1e-12) coef(fitr) ## re 2. ---- summary(lm(age ~ length, data = df)) summary(lm(age ~ diameter, data = df)) summary(fit2 <- lm(age ~ length + diameter, data = df)) cor(X) ### bootstrapping a glmnet model ---- Xb <- model.matrix(fit2)[,-1] y <- df$age fitr <- cv.glmnet(Xb, y) coef(fitr) lambda <- fitr$lambda.1se boot_reg <- function(data, idx) { X <- data$X[idx,] y <- data$y[idx] fit <- glmnet(X, y, lambda = data$lambda) coef(fit) } data <- list(X = Xb, y = y, lambda = lambda) bootstrap <- function(data, stat, R = 1001, ...) { if(is.vector(data)) { N <- length(data) } if(is.data.frame(data)) { N <- nrow(data) } if(is.list(data)) { N <- nrow(data$X) } m0 <- stat(data, 1:N, ...) l <- length(m0) ms <- matrix(rep(NA, R * l), ncol = l) for (r in 1:R) { idx <- sample(1:N, N, replace = TRUE) ms[r, ] <- as.vector(stat(data, idx, ...)) } return(list(ms = ms, m0 = m0)) } b <- bootstrap(data, boot_reg, R = 1000) apply(b$ms, 2, sd) # compare third element to Std. Error for diameter from age ~ diameter model above ## re 1. & 4. ---- fit <- lm(age ~ length + diameter + height + whole_weight + shucked_weight + viscera_weight + shell_weight, data = df) X <- model.matrix(fit)[,-1] y <- df$age plot(glmnet(X, y)) fitr <- cv.glmnet(X, y) plot(fitr) ## order of co coef(glmnet(X, y, lambda = fitr$lambda[2])) ## re 5. ---- # ¡You should not blindly use regularization! fitr <- cv.glmnet(X, y) mean((df$age - predict(fit))^2) # lm() mean((df$age - predict(fitr, newx = X))^2) # cv.glmnet() ### classical cross-validation ---- M <- nrow(X) trdx <- sample(1:M, round(0.7 * M)) trainX <- X[trdx,] trainY <- y[trdx] trainDF <- as.data.frame(trainX) trainDF$age <- trainY testX <- X[-trdx,] testY <- y[-trdx] fit <- lm(age ~ length + diameter + height + whole_weight + shucked_weight + viscera_weight + shell_weight, data = trainDF) fitr <- cv.glmnet(trainX, trainY) mean((testY - predict(fit, newdata = as.data.frame(testX)))^2) mean((testY - predict(fitr, newx = testX))^2) # references # OER: Open Educational Resources are both free and publicly accessible # OER 1. https://uc-r.github.io/regularized_regression # OER 2. https://hastie.su.domains/ElemStatLearn/ starting in Section 3.4 and then scattered throughout much of the rest of the book