inst/doc/autotune-lasso-vignette.R

In autotune: Faster and more Efficient Lasso (than 'glmnet' and 'scalreg') with Data-Driven Tuning

## ----include = FALSE----------------------------------------------------------
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)

## ----install-github, eval = FALSE---------------------------------------------
# # install.packages("devtools")
# # Ensure that you have the Rcpp package installed with version >=1.0.13
# # devtools::install_github("Tathagata-S/autotune")

## ----setup--------------------------------------------------------------------
library(autotune)

## -----------------------------------------------------------------------------
set.seed(10)
n = 300
p = 500
s = 10
beta = c(rep(1, s), rep(0, p - s))
x = matrix(rnorm(n * p), ncol = p)
# Maunal sigma allocation
# y = x %*% beta + rnorm(n, sd = 1)
# Dynamic sigma allocation with snr specified
snr = 2
y = x %*% beta + rnorm(n, sd = sqrt(var(x%*%beta)/snr))

## -----------------------------------------------------------------------------
ptm <- proc.time()
fit.autotune <- autotune_lasso(x, y, alpha = 0.01)
proc.time() - ptm

## -----------------------------------------------------------------------------
b.autotune <- fit.autotune$beta
intercept.autotune <- fit.autotune$a0
lambda.autotune <- fit.autotune$lambda
sigma.seq.autotune <- fit.autotune$CD.path.details$sigma_sq_seq
sigma.estimate.autotune <- fit.autotune$sigma_sq

## -----------------------------------------------------------------------------
head(coef(fit.autotune), n = 21)

## -----------------------------------------------------------------------------
library(glmnet)

## -----------------------------------------------------------------------------
ptm2 <- proc.time()
fit.glmnet <- cv.glmnet(x, y)
proc.time()-ptm2

## ----fig.width=7.5, fig.height=7.5--------------------------------------------
plot(fit.glmnet, cex.lab = 1.5, ylab = "10-fold CV MSE")

# plotting log of lambdas selected by CV
abline(v = log(c(fit.glmnet$lambda.min, fit.glmnet$lambda.1se)), lty = "dashed",col = c(rgb(0,1,0), rgb(0,0.75,0)), lwd = 3)

# plotting log of lambda selected by autotune
abline(v = log(lambda.autotune), col = "blue", lty = "dashed", lwd = 3)
legend("topleft", inset = 0.05,             legend = c("Autotune LASSO",
             "CV_LASSO(min)",
             "CV_LASSO(1se)"),
  col = c("blue", rgb(0,1,0), rgb(0,0.75,0)),
  lty = "dashed", lwd = 2, cex = 1)

## ----fig.width=7.5, fig.height=7.5--------------------------------------------
mse_glmnet <- apply(fit.glmnet$glmnet.fit$beta, 2, function(x) mean((x- beta)^2)/mean(beta^2))
ymin = 0.95 * min(mse_glmnet)
ymax = max(1, max(mse_glmnet))
oldpar <- par(mgp = c(3.5, 1, 0), mar = c(5, 4, 2.5, 0.5) + 1.4)
plot(log(fit.glmnet$lambda), mse_glmnet, 
     type = 'b', col = "red", cex.lab = 2, cex=2, cex.axis = 2,
     ylim = c(ymin, ymax),
     ylab = "Relative MSE w.r.t. true coefficients", xlab = expression(paste("log(",lambda,")")))

# plotting log of lambdas selected by CV
abline(v = log(c(fit.glmnet$lambda.min, fit.glmnet$lambda.1se)), lty = "dashed",col = c(rgb(0,1,0), rgb(0,0.75,0)), lwd = 3)

# plotting log of lambdas selected by autotune
abline(v = log(lambda.autotune), col = "blue", lty = "dashed", lwd = 3)
legend(
  "topleft", inset = 0.05,                  
  legend = c("Autotune LASSO",
             "CV_LASSO(min)",
             "CV_LASSO(1se)"),
  col = c("blue", rgb(0,1,0), rgb(0,0.75,0)),
  lty = "dashed", lwd = 2, cex = 1)
par(oldpar)

## ----fig.width=7.5, fig.height=7.5--------------------------------------------
temp <- fit.autotune$CD.path.details$lambda0
seq.lambdas.autotune <- temp * c(var(y), sigma.seq.autotune)

final.rmse.autotune <- (mean(b.autotune - beta)^2) / (mean(beta)^2)
intermediate.rmse.autotune <- approx(x = fit.glmnet$lambda, y = mse_glmnet, xout = seq.lambdas.autotune)
intermediate.rmse.autotune$y[1] <- 1
rmse.path.autotune <- c(intermediate.rmse.autotune$y, final.rmse.autotune)

oldpar <- par(mgp = c(3.5, 1, 0), mar = c(5, 4, 2.5, 0.5) + 1.4)
plot(log(fit.glmnet$lambda), mse_glmnet, 
     type = 'b', col = "red", cex.lab = 2, cex=2, cex.axis = 2,
     ylim = c(min(ymin, 0.95 * final.rmse.autotune), ymax),
     ylab = "Relative MSE w.r.t. the true coefficients", xlab = expression(paste("log(",lambda,")")))

lines(c(log(seq.lambdas.autotune), log(fit.autotune$lambda)), 
      rmse.path.autotune, 
      col = "blue", 
      lwd = 2)                  
points(log(seq.lambdas.autotune), 
       intermediate.rmse.autotune$y, 
       col = "blue", 
       pch = 1,               
       lwd = 2,
       cex = 2) 
points(log(fit.autotune$lambda),
      final.rmse.autotune,
      col = "blue",
      pch = 11,
      lwd = 3,
      cex = 3
    )
par(oldpar)

## -----------------------------------------------------------------------------
library(AUC)

## -----------------------------------------------------------------------------
b.glmnet = coef(fit.glmnet, s = "lambda.min")[-1]

## -----------------------------------------------------------------------------
auc(roc(abs(b.autotune), as.factor(beta != 0)))
auc(roc(abs(b.glmnet), as.factor(beta != 0)))

## -----------------------------------------------------------------------------
data("sp500")
attach(sp500)

scaled_lasso_benchmark <- read.csv(
  system.file(
    "extdata",
    "sp500-scalreg-1.0.1.csv",
    package = "autotune"
  )
)
scaled_lasso_observed <- scaled_lasso_benchmark[
  which(!scaled_lasso_benchmark$bootstrap),
  ,
  drop = FALSE
]

n <- 190
rang <- 252 - n

X_train = sp500.percent[1:n, 3:(dim(sp500.percent)[2])]
Y_train = sp500.percent[1:n, 1]

x_train = scale(X_train)
y_train = Y_train - mean(Y_train)

X_test = sp500.percent[(1:rang) + n, 3: (dim(sp500.percent)[2])]
Y_test = sp500.percent[(1:rang) + n, 1]
x_test <- scale(X_test)
y_test <- Y_test - mean(Y_test)
# Default value of beta_iter_max is 40
ans_autotune <- autotune_lasso(x_train, y_train, beta_iter_max = 40, trace_it = TRUE)
pred_err_autotune <- mean( (y_test - x_test %*% ans_autotune$beta)^2 )

pred_err_scallas <- scaled_lasso_observed$test_mse

cv_fit <- cv.glmnet(x_train, y_train, alpha = 1, intercept = F)
pred_err_cv_min <- mean( (y_test - predict(cv_fit, newx = x_test, s = "lambda.min"))^2 )
pred_err_cv_1se <- mean( (y_test - predict(cv_fit, newx = x_test, s = "lambda.1se"))^2 )


pred_err_autotune
pred_err_scallas
pred_err_cv_min
pred_err_cv_1se

scaled_lasso_observed$nonzero_coefs
sum(ans_autotune$beta != 0)
sum(coef(cv_fit, s = "lambda.min") != 0)
sum(coef(cv_fit, s = "lambda.1se") != 0)
detach(sp500)

## -----------------------------------------------------------------------------
real_life_prediction_comparison <- function(x, y, x_test, y_test, scaled_lasso_result, alpha = 0.01, beta_iter_max = 40, plot_regu = TRUE, trace_it = TRUE){

ans_autotune <- autotune_lasso(x, y, alpha = alpha, beta_iter_max = beta_iter_max, trace_it = trace_it)
# ans_glmnet <- glmnet(x, y, alpha = 1, intercept = FALSE)
cv_fit <- cv.glmnet(x, y, alpha = 1, intercept = F)
vary <- var(y_test)
mse_glmnet <- apply(cv_fit$glmnet.fit$beta, 2, function(b) mean((y_test - x_test%*%b)^2)/vary)

# log_glmnet_lambdas <- log(fit.glmnet$lambda)
# log_autotune_lambdas <- log(intermediate_lambdas_used)



if(plot_regu) {
  intermediate_lambdas_used <- c( c(var(y)/2, ans_autotune$CD.path.details$sigma_sq_seq) * ans_autotune$CD.path.details$lambda0)
  
  autotune_lambdas <- c(intermediate_lambdas_used, ans_autotune$lambda)
  
  final.rmse.autotune <- mean((y_test - x_test%*%ans_autotune$beta)^2) / vary
  intermediate.rmse.autotune <- approx(x = log(cv_fit$lambda), y = mse_glmnet, xout = log(intermediate_lambdas_used), rule = 2)
  intermediate.rmse.autotune$y[1] <- 1
  rmse.path.autotune <- c(intermediate.rmse.autotune$y, final.rmse.autotune)
  
  
  ysquaredmean <- mean(y^2)
  cv_fit$cvm <- cv_fit$cvm/ysquaredmean
  cv_fit$cvsd <- cv_fit$cvsd/ysquaredmean
  cv_fit$cvup <- cv_fit$cvm + cv_fit$cvsd
  cv_fit$cvlo <- cv_fit$cvm - cv_fit$cvsd
  
  cv_fit_lowerlim <- min(cv_fit$cvm - cv_fit$cvsd, mse_glmnet)
  cv_fit_upperlim <- max(cv_fit$cvm + cv_fit$cvsd, mse_glmnet)
  
  ymin = 0.95 * min(cv_fit_lowerlim , min(rmse.path.autotune))
  ymax = max(cv_fit_upperlim, rmse.path.autotune)
  ymax_trunc = ymin + 0.8 * (ymax - ymin)  # Truncate vertical lines here
  
  ydiff = ymax - ymin
  oldpar <- par(mgp = c(3.5, 1, 0), mar = c(6, 5, 4, 2) + 0.1)
  on.exit(par(oldpar), add = TRUE)
  
  par(mar = c(5, 7, 4, 2))
  plot(cv_fit,
       ylim = c(ymin, max(cv_fit_upperlim, rmse.path.autotune)),
       cex.lab = 2,
       cex.axis = 2,
       ylab = expression(
    "||" * (Y[test] - X[test] * hat(beta)) * "||"[2]^2
    ~ "/" ~
    "||" * Y[test]  * "||"[2]^2
  ))
  # expression(paste(sqrt("RMSE")))
  
  points(log(cv_fit$lambda), mse_glmnet,
         type = 'b', col = "orange", lwd = 3)
  
  
  lines(rep(log(ans_autotune$lambda), 2),
        c(ymin, ymin + 0.92 * ydiff), col = "blue", lwd = 0.5, lty = "dashed")
  
  lines(log(autotune_lambdas), rmse.path.autotune,
        col = "blue", lwd = 1.5)
  
  points(log(autotune_lambdas), rmse.path.autotune,
         col = "blue", pch = 1, lwd = 1.5, cex = 1)
  
  points(log(ans_autotune$lambda), final.rmse.autotune,
         col = "blue", pch = 11, lwd = 2, cex = 2)
  
  
  
  vlines <- c(
    log(cv_fit$lambda.min),
    log(cv_fit$lambda.1se)
  )
  vline_labels <- c("CV(min)", "CV(1se)")
  vline_colors <- c(rgb(0,1,0), rgb(0,0.75,0))
  vline_ltys <- c("dashed", "dashed")
  vline_lwds <- c(3, 3)
  
  for (i in seq_along(vlines)) {
    lines(rep(vlines[i], 2),
          c(ymin, ymin + (0.55 + i * 0.1) * ydiff),
          col = vline_colors[i], lty = vline_ltys[i],
          lwd = vline_lwds[i])
  
    text(x = vlines[i], y = ymin + (0.57 + i * 0.1) * ydiff,
         labels = vline_labels[i], srt = -30, adj = 0,
         xpd = TRUE, cex = 1.4, col = vline_colors[i])
  }
  
  text(x = log(ans_autotune$lambda),
       y =  ymin + 0.95 * ydiff,
       labels = "Autotune", srt = -30, adj = 0,
       xpd = TRUE, cex = 1.4, col = "blue")
  
  if (!is.null(cv_fit$cvsd)) {
    arrows(x0 = log(cv_fit$lambda), 
           y0 = cv_fit$cvm - cv_fit$cvsd,
           x1 = log(cv_fit$lambda),
           y1 = cv_fit$cvm + cv_fit$cvsd,
           angle = 90, code = 3, length = 0.05, col = "red", lwd = 1.5)
  }
  
  
  legend("topright",
         inset = c(0.05, 0.02),
         legend = c("CV Error", "Test Error"),
         col = c("red", "orange"),
         pt.cex = c(1.5, 2), pch = c(19, 19),
         lwd = c(1, 3),
         text.col = c("red", "orange"),
         horiz = FALSE,
         bty = "n")

}

# mse_autotune <- mean( (y_test - x_test %*% ans_autotune$beta)^2 )
# mse_cv_min <- mean( (y_test - predict(cv_fit, newx = x_test, s = "lambda.min"))^2 )
# mse_cv_1se <- mean( (y_test - predict(cv_fit, newx = x_test, s = "lambda.1se"))^2 )


oos_scallas <- scaled_lasso_result$out_of_sample_r2
oos_autotune <- 1 -  mean( (y_test - x_test %*% ans_autotune$beta)^2 ) / vary
oos_cv_min <- 1 - mean( (y_test - predict(cv_fit, newx = x_test, s = "lambda.min"))^2 ) / vary
oos_cv_1se <- 1 - mean( (y_test - predict(cv_fit, newx = x_test, s = "lambda.1se"))^2 ) / vary

return(data.frame(out_of_sample_autotune = oos_autotune,
out_of_sample_scaled_lasso = oos_scallas,
out_of_sample_cv_min = oos_cv_min,
out_of_sample_cv_1se = oos_cv_1se,
nonzero_coefs_autotune = sum(ans_autotune$beta != 0),
nonzero_coefs_scaled_lasso = scaled_lasso_result$nonzero_coefs,
nonzero_coefs_cvmin = sum(coef(cv_fit, s = "lambda.min") != 0),
nonzero_coefs_cv1se = sum(coef(cv_fit, s = "lambda.1se") != 0)))
}

## ----fig.width=7.5, fig.height=7.5--------------------------------------------
ans <- real_life_prediction_comparison(
  x_train,
  y_train,
  x_test,
  y_test,
  scaled_lasso_result = scaled_lasso_observed,
  plot_regu = TRUE
)
names(ans)

## -----------------------------------------------------------------------------
library(dplyr)

## ----fig.width=7.5, fig.height=7.5--------------------------------------------
#Number of bootstraps
B <- 20


result_list <- list()

result_list[[1]] <- data.frame(
  n = n,
  bootstrap = FALSE,
  real_life_prediction_comparison(
    x_train,
    y_train,
    x_test,
    y_test,
    scaled_lasso_result = scaled_lasso_observed,
    plot_regu = FALSE,
    trace_it = FALSE
  )[1:4]
)

for( j in (1:B)+1) {
  set.seed(j)
  boot.sample <- sample(n, size = n, replace = TRUE)
  boot.x <- x_train[boot.sample, ]
  boot.y <- y_train[boot.sample]
  scaled_lasso_result <- scaled_lasso_benchmark[
    which(scaled_lasso_benchmark$seed == j),
    ,
    drop = FALSE
  ]
  result_list[[j]] <- c(
    n = n,
    bootstrap = TRUE,
    real_life_prediction_comparison(
      boot.x,
      boot.y,
      x_test,
      y_test,
      scaled_lasso_result = scaled_lasso_result,
      plot_regu = FALSE,
      trace_it = FALSE
    )[1:4]
  )
}

mat <- do.call(rbind, result_list)     
df <- as.data.frame(mat, stringsAsFactors = FALSE)

colnames(df) <- c("n", "bootstrap", "Autotune",
                  "Scaled", "CVmin",
                  "CV1se")
df <- df %>%
  mutate(
    n = as.integer(n),
    bootstrap = as.logical(bootstrap), 
    across(!c(n, bootstrap), as.numeric)
  )

## ----fig.width=7.5, fig.height=7.5--------------------------------------------
alg_cols <- colnames(df)[3:6]
boot_only <- df %>% filter(bootstrap == TRUE)
boxlist <- lapply(alg_cols, function(cn) boot_only[[cn]])
names(boxlist) <- alg_cols

# Uncomment the following line and line 448 for saving the boxplot in your working directory

# pdf(file.path(paste0("OOS_boxplot_for_n=",df[1,1],"_B=", B,".pdf")), width = 8, height = 9)

oldpar <- par(mar = c(4, 6, 4, 2))
boxplot(boxlist, main = paste0("Bootstrapped Out-of-Sample R^2 for n = ", df[1,1]), ylab = expression(paste("Out-of-Sample ",R^2)), ylim = c(max(0.98 * min(df[alg_cols]), 0), 1),
        cex.lab = 2,     
  cex.axis = 1.8,    
  cex.main = 1.5)

nonboot <- df %>% filter(bootstrap == FALSE)
if(nrow(nonboot) > 0){
  points(1:4, as.numeric(nonboot[1, alg_cols]), pch = 18, col = "red", cex = 2)
}

legend(
  "topright",
  legend = expression(paste("Observed O.O.S. ",R^2)),
  col = "red",
  pch = 18,
  pt.cex = 2,
  cex = 2,          
  bty = "n"           
)

par(oldpar)
# dev.off()