Nothing
autotune Lasso
In autotune: Faster and more Efficient Lasso (than 'glmnet' and 'scalreg') with Data-Driven Tuning
knitr::opts_chunk$set(
collapse = TRUE,
comment = "#>"
)
The R package autotune implements the Lasso with data-driven tuning for linear models.
Installation and Loading
Installation
The development version of the autotune package can be installed from GitHub using
# install.packages("devtools")
# Ensure that you have the Rcpp package installed with version >=1.0.13
# devtools::install_github("Tathagata-S/autotune")
When installing from GitHub, in order to build the package from source, you need to have the appropriate R development tools installed (Rtools on Windows, or these tools on Mac).
Load Package
After installation, the package can be loaded in the standard way:
library(autotune)
High Dimensional Regression
autotune performs lasso via the autotune_lasso() function.
We illustrate autotune lasso on simulated data using a linear model with $s=10$, $n=300$, $p=500$. Reader can also specify different configurations of data generating process in the following code chunk.
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))
Runtime of autotune Lasso
Given data (x, y), run the autotune lasso as follows with default $\alpha = 0.01$.
ptm <- proc.time()
fit.autotune <- autotune_lasso(x, y, alpha = 0.01)
proc.time() - ptm
The regression coefficients $\beta$s, intercept, final lambda and the sequence of estimated sigmas can be extracted from the fitted autotune lasso object as follows
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
Now, autotune also supports the S3 function coef which returns the estimated regression coefficients $\beta$s along with the intercept
head(coef(fit.autotune), n = 21)
Comparison with Lasso tuned via Cross Validation (CV Lasso)
We contrast our solution with Cross-Validation-tuned Lasso using cv.glmnet() in the glmnet package.
library(glmnet)
Runtime of CV Lasso
ptm2 <- proc.time()
fit.glmnet <- cv.glmnet(x, y)
proc.time()-ptm2
So, autotune Lasso shows faster runtimes as compared to CV Lasso.
Comparison of quality of tuning between autotune and CV Lasso
Now, we will visualize the quality of lambdas selected by autotune and CV with respect to the 10-fold CV MSE and true Relative MSE.
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)
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)
Across the lambda grid, we plot RMSE of solution path taken by different tuners.
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)))
Comparison on real data provided in Scaled Lasso's R package scalreg
The scalreg package was removed from CRAN in July 2026.
Comparing prediction errors of autotune and benchmarks: CV and Scaled Lasso
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)
Comparing via bootstrap
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)))
}
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)
#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)
)
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()
Try the autotune package in your browser
Any scripts or data that you put into this service are public.
autotune documentation built on Aug. 21, 2026, 5:18 p.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
knitr::opts_chunk$set( collapse = TRUE, comment = "#>" )
The R package autotune implements the Lasso with data-driven tuning for linear models.
Installation and Loading
Installation
The development version of the autotune package can be installed from GitHub using
# install.packages("devtools") # Ensure that you have the Rcpp package installed with version >=1.0.13 # devtools::install_github("Tathagata-S/autotune")
When installing from GitHub, in order to build the package from source, you need to have the appropriate R development tools installed (Rtools on Windows, or these tools on Mac).
Load Package
After installation, the package can be loaded in the standard way:
library(autotune)
High Dimensional Regression
autotune performs lasso via the autotune_lasso() function.
We illustrate autotune lasso on simulated data using a linear model with $s=10$, $n=300$, $p=500$. Reader can also specify different configurations of data generating process in the following code chunk.
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))
Runtime of autotune Lasso
Given data (x, y), run the autotune lasso as follows with default $\alpha = 0.01$.
ptm <- proc.time() fit.autotune <- autotune_lasso(x, y, alpha = 0.01) proc.time() - ptm
The regression coefficients $\beta$s, intercept, final lambda and the sequence of estimated sigmas can be extracted from the fitted autotune lasso object as follows
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
Now, autotune also supports the S3 function coef which returns the estimated regression coefficients $\beta$s along with the intercept
head(coef(fit.autotune), n = 21)
Comparison with Lasso tuned via Cross Validation (CV Lasso)
We contrast our solution with Cross-Validation-tuned Lasso using cv.glmnet() in the glmnet package.
library(glmnet)
Runtime of CV Lasso
ptm2 <- proc.time() fit.glmnet <- cv.glmnet(x, y) proc.time()-ptm2
So, autotune Lasso shows faster runtimes as compared to CV Lasso.
Comparison of quality of tuning between autotune and CV Lasso
Now, we will visualize the quality of lambdas selected by autotune and CV with respect to the 10-fold CV MSE and true Relative MSE.
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)
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)
Across the lambda grid, we plot RMSE of solution path taken by different tuners.
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)))
Comparison on real data provided in Scaled Lasso's R package scalreg
The scalreg package was removed from CRAN in July 2026.
Comparing prediction errors of autotune and benchmarks: CV and Scaled Lasso
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)
Comparing via bootstrap
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))) }
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)
#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) )
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()
Try the autotune package in your browser
Any scripts or data that you put into this service are public.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.