}
#- Concatenate results:
results_from_simulation <- list()
results_from_simulation$betas <- rbind(oracle_beta, #1
lasso_beta, #2
lasso_beta2,
quantile_beta[which.min(quantile_cvm), ],
quantile_beta)
results_from_simulation$opt   <- c(tau.seq[which.min(quantile_cvm)])
results_from_simulation$mse   <- c(oracle_mse,
lasso_mse,
lasso_mse2,
quantile_mse[which.min(quantile_cvm)],
quantile_mse)
results_from_simulation$meanY <- mean(y)
out[[b]] <- results_from_simulation
} # End simulation
# Save simulation results
save(simParams, out, beta0, file = paste0("results/results_s2_OLS_ar_", simParams$model, ".RData"))
}
rm(list = ls())
library(xrnet)
library(readxl)
library(glmnet)
library(ggplot2)
library(dplyr)
library(xtable)
#- Source important files
simSheet <- read_xlsx("simulationList.xlsx", sheet = 1)
# Make sure to comment and un-comment the correct model you are interested in.
#- Section 3.2 (Low Dimensional)
# idx = 1: n = 400, p = 50
# idx = 2: n = 800, p = 50
# idx = 3: n = 400, p = 200
# idx = 4: n = 800, p = 200
# idx = 5: Model 2
# idx = 7: Model 3
# idx = 8: Model 4
# OLS_ar_: idx = 1 - 4 (Table 2); idx = 5, 7, 8 (Table 3)
# OLS_cs_: idx = 1, 3 (Table S4)
# OLS_ind_: idx =  1, 3 (Table S4)
idx <- 1
load(file = paste0("results/results_s2_OLS_ar_", idx, ".RData"))
#load(file = paste0("results/results_s2_OLS_cs_", idx, ".RData"))
#load(file = paste0("results/results_s2_OLS_ind_", idx, ".RData"))
#- Section 3.2 (High Dimensional)
# idx = 1: n = 400, p = 2000; model 1
# idx = 2: n = 400, p = 2000; model 2
# OLS_HD_ar_: idx = 1 - 2 (Table 4)
#load(file = paste0("results/results_s2_OLS_HD_ind_", idx, ".RData"))
#load(file = paste0("results/results_s2_OLS_HD_ar_", idx, ".RData"))
#- Section 3.4:
# idx = 1, 5 (Table 5)
#load(file = paste0("results/results_s4_logistic_", idx, ".RData"))
beta0
B <- length(out)
source(paste0("sourceFiles/evaluationMeasures.R"))
# Store results in these matrices
bias_matrix <- matrix(NA, nrow = B, ncol = 23) # Calculate bias
mse_matrix  <- matrix(NA, nrow = B, ncol = 23) # Calculate prediction measure (PMSE or AUC)
fp_matrix   <- matrix(NA, nrow = B, ncol = 23) # Calculate # of FP
fn_matrix   <- matrix(NA, nrow = B, ncol = 23) # Calculate # of FN
beta_res    <- list()
opt_quant   <- matrix(NA, nrow = B, ncol = 1) # Calculate optimal tau
mean_y      <- numeric() # Calculate mean value of outcome (useful for binary outcome)
#- Get simulation results
for(b in 1:B) {
bias_matrix[b, ] <- apply(out[[b]]$betas, 1, function(x) sqrt(sum((x - beta0)^2)))
beta_res[[b]]    <- apply(out[[b]]$betas, 1, function(x) getBetaCounts(x, beta0))
fn_matrix[b, ]   <- apply(out[[b]]$betas, 1, function(x) getModelSelection(x, beta0))[1, ]
fp_matrix[b, ]   <- apply(out[[b]]$betas, 1, function(x) getModelSelection(x, beta0))[2, ]
mse_matrix[b, ]  <- out[[b]]$mse
opt_quant[b, ]   <- out[[b]]$opt
mean_y[b]        <- out[[b]]$meanY
}
beta_matrix <- matrix(NA, nrow = 23, ncol = 5)
for(i in 1:5) {
tmp <- matrix(NA, nrow = B, ncol = 23)
for(b in 1:B) {
tmp[b, ] <- beta_res[[b]][i, ]
}
beta_matrix[, i] <- apply(tmp, 2,  mean)
}
rm(i, b)
#####- Output results
res_table <- cbind(apply(bias_matrix, 2, mean),
apply(fp_matrix, 2, mean),
apply(fn_matrix, 2, mean),
beta_matrix,
apply(mse_matrix, 2, mean))
res_table <- res_table[c(1, 2, 3, 4, 15), ];
colnames(res_table) <- c("MSB", "FP", "FN", paste0("beta_", which(beta0!=0)), "PMSE")
rownames(res_table) <- c("Oracle", "LASSO", "nLASSO", "asymLASSO", "asymLASSO (tau = 0.5)")
xtable(res_table)
apply(opt_quant, 2, mean)
#for(i in 1:B) {
#  print(out[[i]]$betas[4, c(1, 4)])
#}
res_table[, 1]
rm(list = ls())
library(xrnet)
library(readxl)
library(glmnet)
library(ggplot2)
library(dplyr)
library(xtable)
#- Source important files
simSheet <- read_xlsx("simulationList.xlsx", sheet = 1)
# Make sure to comment and un-comment the correct model you are interested in.
#- Section 3.2 (Low Dimensional)
# idx = 1: n = 400, p = 50
# idx = 2: n = 800, p = 50
# idx = 3: n = 400, p = 200
# idx = 4: n = 800, p = 200
# idx = 5: Model 2
# idx = 7: Model 3
# idx = 8: Model 4
# OLS_ar_: idx = 1 - 4 (Table 2); idx = 5, 7, 8 (Table 3)
# OLS_cs_: idx = 1, 3 (Table S4)
# OLS_ind_: idx =  1, 3 (Table S4)
idx <- 1
load(file = paste0("results/results_s2_OLS_ar_", idx, ".RData"))
#load(file = paste0("results/results_s2_OLS_cs_", idx, ".RData"))
#load(file = paste0("results/results_s2_OLS_ind_", idx, ".RData"))
#- Section 3.2 (High Dimensional)
# idx = 1: n = 400, p = 2000; model 1
# idx = 2: n = 400, p = 2000; model 2
# OLS_HD_ar_: idx = 1 - 2 (Table 4)
#load(file = paste0("results/results_s2_OLS_HD_ind_", idx, ".RData"))
#load(file = paste0("results/results_s2_OLS_HD_ar_", idx, ".RData"))
#- Section 3.4:
# idx = 1, 5 (Table 5)
#load(file = paste0("results/results_s4_logistic_", idx, ".RData"))
beta0
B <- length(out)
source(paste0("sourceFiles/evaluationMeasures.R"))
# Store results in these matrices
bias_matrix <- matrix(NA, nrow = B, ncol = 23) # Calculate bias
mse_matrix  <- matrix(NA, nrow = B, ncol = 23) # Calculate prediction measure (PMSE or AUC)
fp_matrix   <- matrix(NA, nrow = B, ncol = 23) # Calculate # of FP
fn_matrix   <- matrix(NA, nrow = B, ncol = 23) # Calculate # of FN
beta_res    <- list()
opt_quant   <- matrix(NA, nrow = B, ncol = 1) # Calculate optimal tau
mean_y      <- numeric() # Calculate mean value of outcome (useful for binary outcome)
#- Get simulation results
for(b in 1:B) {
bias_matrix[b, ] <- apply(out[[b]]$betas, 1, function(x) sqrt(sum((x - beta0)^2)))
beta_res[[b]]    <- apply(out[[b]]$betas, 1, function(x) getBetaCounts(x, beta0))
fn_matrix[b, ]   <- apply(out[[b]]$betas, 1, function(x) getModelSelection(x, beta0))[1, ]
fp_matrix[b, ]   <- apply(out[[b]]$betas, 1, function(x) getModelSelection(x, beta0))[2, ]
mse_matrix[b, ]  <- out[[b]]$mse
opt_quant[b, ]   <- out[[b]]$opt
mean_y[b]        <- out[[b]]$meanY
}
beta_matrix <- matrix(NA, nrow = 23, ncol = 5)
for(i in 1:5) {
tmp <- matrix(NA, nrow = B, ncol = 23)
for(b in 1:B) {
tmp[b, ] <- beta_res[[b]][i, ]
}
beta_matrix[, i] <- apply(tmp, 2,  mean)
}
rm(i, b)
#####- Output results
res_table <- cbind(apply(bias_matrix, 2, mean),
apply(fp_matrix, 2, mean),
apply(fn_matrix, 2, mean),
beta_matrix,
apply(mse_matrix, 2, mean))
res_table <- res_table[c(1, 2, 3, 4, 15), ];
colnames(res_table) <- c("MSB", "FP", "FN", paste0("beta_", which(beta0!=0)), "PMSE")
rownames(res_table) <- c("Oracle", "LASSO", "nLASSO", "asymLASSO", "asymLASSO (tau = 0.5)")
xtable(res_table)
apply(opt_quant, 2, mean)
#for(i in 1:B) {
#  print(out[[i]]$betas[4, c(1, 4)])
#}
rm(list = ls())
library(xrnet)
library(readxl)
library(glmnet)
library(ggplot2)
library(dplyr)
library(xtable)
#- Source important files
simSheet <- read_xlsx("simulationList.xlsx", sheet = 1)
# Make sure to comment and un-comment the correct model you are interested in.
#- Section 3.2 (Low Dimensional)
# idx = 1: n = 400, p = 50
# idx = 2: n = 800, p = 50
# idx = 3: n = 400, p = 200
# idx = 4: n = 800, p = 200
# idx = 5: Model 2
# idx = 7: Model 3
# idx = 8: Model 4
# OLS_ar_: idx = 1 - 4 (Table 2); idx = 5, 7, 8 (Table 3)
# OLS_cs_: idx = 1, 3 (Table S4)
# OLS_ind_: idx =  1, 3 (Table S4)
idx <- 2
load(file = paste0("results/results_s2_OLS_ar_", idx, ".RData"))
#load(file = paste0("results/results_s2_OLS_cs_", idx, ".RData"))
#load(file = paste0("results/results_s2_OLS_ind_", idx, ".RData"))
#- Section 3.2 (High Dimensional)
# idx = 1: n = 400, p = 2000; model 1
# idx = 2: n = 400, p = 2000; model 2
# OLS_HD_ar_: idx = 1 - 2 (Table 4)
#load(file = paste0("results/results_s2_OLS_HD_ind_", idx, ".RData"))
#load(file = paste0("results/results_s2_OLS_HD_ar_", idx, ".RData"))
#- Section 3.4:
# idx = 1, 5 (Table 5)
#load(file = paste0("results/results_s4_logistic_", idx, ".RData"))
beta0
B <- length(out)
source(paste0("sourceFiles/evaluationMeasures.R"))
# Store results in these matrices
bias_matrix <- matrix(NA, nrow = B, ncol = 23) # Calculate bias
mse_matrix  <- matrix(NA, nrow = B, ncol = 23) # Calculate prediction measure (PMSE or AUC)
fp_matrix   <- matrix(NA, nrow = B, ncol = 23) # Calculate # of FP
fn_matrix   <- matrix(NA, nrow = B, ncol = 23) # Calculate # of FN
beta_res    <- list()
opt_quant   <- matrix(NA, nrow = B, ncol = 1) # Calculate optimal tau
mean_y      <- numeric() # Calculate mean value of outcome (useful for binary outcome)
#- Get simulation results
for(b in 1:B) {
bias_matrix[b, ] <- apply(out[[b]]$betas, 1, function(x) sqrt(sum((x - beta0)^2)))
beta_res[[b]]    <- apply(out[[b]]$betas, 1, function(x) getBetaCounts(x, beta0))
fn_matrix[b, ]   <- apply(out[[b]]$betas, 1, function(x) getModelSelection(x, beta0))[1, ]
fp_matrix[b, ]   <- apply(out[[b]]$betas, 1, function(x) getModelSelection(x, beta0))[2, ]
mse_matrix[b, ]  <- out[[b]]$mse
opt_quant[b, ]   <- out[[b]]$opt
mean_y[b]        <- out[[b]]$meanY
}
beta_matrix <- matrix(NA, nrow = 23, ncol = 5)
for(i in 1:5) {
tmp <- matrix(NA, nrow = B, ncol = 23)
for(b in 1:B) {
tmp[b, ] <- beta_res[[b]][i, ]
}
beta_matrix[, i] <- apply(tmp, 2,  mean)
}
rm(i, b)
#####- Output results
res_table <- cbind(apply(bias_matrix, 2, mean),
apply(fp_matrix, 2, mean),
apply(fn_matrix, 2, mean),
beta_matrix,
apply(mse_matrix, 2, mean))
res_table <- res_table[c(1, 2, 3, 4, 15), ];
colnames(res_table) <- c("MSB", "FP", "FN", paste0("beta_", which(beta0!=0)), "PMSE")
rownames(res_table) <- c("Oracle", "LASSO", "nLASSO", "asymLASSO", "asymLASSO (tau = 0.5)")
xtable(res_table)
apply(opt_quant, 2, mean)
rm(list = ls())
library(xrnet)
library(readxl)
library(glmnet)
library(ggplot2)
library(dplyr)
library(xtable)
#- Source important files
simSheet <- read_xlsx("simulationList.xlsx", sheet = 1)
# Make sure to comment and un-comment the correct model you are interested in.
#- Section 3.2 (Low Dimensional)
# idx = 1: n = 400, p = 50
# idx = 2: n = 800, p = 50
# idx = 3: n = 400, p = 200
# idx = 4: n = 800, p = 200
# idx = 5: Model 2
# idx = 7: Model 3
# idx = 8: Model 4
# OLS_ar_: idx = 1 - 4 (Table 2); idx = 5, 7, 8 (Table 3)
# OLS_cs_: idx = 1, 3 (Table S4)
# OLS_ind_: idx =  1, 3 (Table S4)
idx <- 2
load(file = paste0("results/results_s2_OLS_ar_", idx, ".RData"))
#load(file = paste0("results/results_s2_OLS_cs_", idx, ".RData"))
#load(file = paste0("results/results_s2_OLS_ind_", idx, ".RData"))
#- Section 3.2 (High Dimensional)
# idx = 1: n = 400, p = 2000; model 1
# idx = 2: n = 400, p = 2000; model 2
# OLS_HD_ar_: idx = 1 - 2 (Table 4)
#load(file = paste0("results/results_s2_OLS_HD_ind_", idx, ".RData"))
#load(file = paste0("results/results_s2_OLS_HD_ar_", idx, ".RData"))
#- Section 3.4:
# idx = 1, 5 (Table 5)
#load(file = paste0("results/results_s4_logistic_", idx, ".RData"))
beta0
B <- length(out)
source(paste0("sourceFiles/evaluationMeasures.R"))
# Store results in these matrices
bias_matrix <- matrix(NA, nrow = B, ncol = 23) # Calculate bias
mse_matrix  <- matrix(NA, nrow = B, ncol = 23) # Calculate prediction measure (PMSE or AUC)
fp_matrix   <- matrix(NA, nrow = B, ncol = 23) # Calculate # of FP
fn_matrix   <- matrix(NA, nrow = B, ncol = 23) # Calculate # of FN
beta_res    <- list()
opt_quant   <- matrix(NA, nrow = B, ncol = 1) # Calculate optimal tau
mean_y      <- numeric() # Calculate mean value of outcome (useful for binary outcome)
#- Get simulation results
for(b in 1:B) {
bias_matrix[b, ] <- apply(out[[b]]$betas, 1, function(x) sqrt(sum((x - beta0)^2)))
beta_res[[b]]    <- apply(out[[b]]$betas, 1, function(x) getBetaCounts(x, beta0))
fn_matrix[b, ]   <- apply(out[[b]]$betas, 1, function(x) getModelSelection(x, beta0))[1, ]
fp_matrix[b, ]   <- apply(out[[b]]$betas, 1, function(x) getModelSelection(x, beta0))[2, ]
mse_matrix[b, ]  <- out[[b]]$mse
opt_quant[b, ]   <- out[[b]]$opt
mean_y[b]        <- out[[b]]$meanY
}
beta_matrix <- matrix(NA, nrow = 23, ncol = 5)
for(i in 1:5) {
tmp <- matrix(NA, nrow = B, ncol = 23)
for(b in 1:B) {
tmp[b, ] <- beta_res[[b]][i, ]
}
beta_matrix[, i] <- apply(tmp, 2,  mean)
}
rm(i, b)
#####- Output results
res_table <- cbind(apply(bias_matrix, 2, mean),
apply(fp_matrix, 2, mean),
apply(fn_matrix, 2, mean),
beta_matrix,
apply(mse_matrix, 2, mean))
res_table <- res_table[c(1, 2, 3, 4), ];
colnames(res_table) <- c("MSB", "FP", "FN", paste0("beta_", which(beta0!=0)), "PMSE")
rownames(res_table) <- c("Oracle", "LASSO", "nLASSO", "asymLASSO")
xtable(res_table)
apply(opt_quant, 2, mean)
rm(list = ls())
library(data.table)
library(dplyr)
library(hdrm)
library(glmnet)
library(xrnet)
library(doParallel)
library(pROC)
library(VennDiagram)
dat <- readData(bcTCGA)
X   <- cbind(dat$y, dat$X)
X   <- scale(X, scale = TRUE)
colnames(X)[1] <- "BRCA1"
variantId <- colnames(X)
colnames(X)[1] <- "BRCA1"
variantId <- colnames(X)
tau.seq <- seq(0.05, 0.95, by = 0.05)
nCores <- detectCores() - 2; nCores
y.idx <- which(variantId == "BRCA1")
set.seed(2021)
foldId <- sample(1:10, dim(X)[1], replace = TRUE)
# Fit LASSO
registerDoParallel(nCores)
fit.lasso  <- cv.glmnet(X[, -y.idx], X[, y.idx],
family = "gaussian",
alpha = 1,
type.measure = "deviance",
nfolds = 10,
foldid = foldId,
nlambda = 20,
parallel = TRUE)
stopImplicitCluster()
l1_beta <- coef(fit.lasso, s = "lambda.min")[-1]
registerDoParallel(nCores)
fit.lasso  <- cv.glmnet(X[, -y.idx], X[, y.idx],
family = "gaussian",
alpha = 1,
lower.limits = 0,
type.measure = "deviance",
nfolds = 10,
foldid = foldId,
nlambda = 20,
parallel = TRUE)
stopImplicitCluster()
nl1_beta <- coef(fit.lasso, s = "lambda.min")[-1]
# Fit Q1
quantile_beta <- matrix(NA, nrow = length(tau.seq),
ncol = (dim(X)[2] - 1))
quantile_cvm <- numeric()
rownames(quantile_beta) <- paste("tau = ", tau.seq)
# Fit asymmetric LASSO
quantile_beta <- matrix(NA, nrow = length(tau.seq),
ncol = (dim(X)[2] - 1))
quantile_cvm <- numeric()
rownames(quantile_beta) <- paste("tau = ", tau.seq)
for(j in 1:length(tau.seq)) {
#cat(paste(j))
registerDoParallel(nCores)
fit.quantile       <- tune_xrnet(X[, -y.idx], X[, y.idx],
family = "gaussian",
penalty_main = define_penalty("q1", tau.seq[j],
num_penalty = 20),
loss = "deviance",
nfolds = 10,
foldid = foldId,
parallel = TRUE)
stopImplicitCluster()
quantile_beta[j, ] <- coef(fit.quantile)$betas[, , 1]
quantile_cvm[j]    <- min(fit.quantile$cv_mean)
}
q1_beta <- quantile_beta[which.min(quantile_cvm), ]
q1_tau <- which.min(quantile_cvm)
plot(quantile_cvm)
plot(quantile_cvm, type = "b")
plot(quantile_cvm, type = "b", pch = 19)
plot(quantile_cvm, type = "b", pch = 19,
ylab = experssion(Mean~~Cross~~Validation~~Error),
xlab = expression(tau))
plot(quantile_cvm, type = "b", pch = 19,
ylab = expression(Mean~~Cross~~Validation~~Error),
xlab = expression(tau))
plot(quantile_cvm, ylab = "Mean Cross Validation Error",
xlab = expression(tau), xaxt = "n", type = "o", pch = 19)
# Figure 4
plot(quantile_cvm, ylab = "Mean Cross Validation Error",
xlab = expression(tau), xaxt = "n", type = "o", pch = 19)
axis(1, at = seq(2, 19, by = 2), labels = tau.seq[seq(2, 19, by = 2)])
min(fit.nlasso$cvm)
rm(list = ls())
library(data.table)
library(dplyr)
library(hdrm)
library(glmnet)
library(xrnet)
library(doParallel)
library(pROC)
library(VennDiagram)
dat <- readData(bcTCGA)
X   <- cbind(dat$y, dat$X)
X   <- scale(X, scale = TRUE)
colnames(X)[1] <- "BRCA1"
variantId <- colnames(X)
tau.seq <- seq(0.05, 0.95, by = 0.05)
nCores <- detectCores() - 2; nCores
y.idx <- which(variantId == "BRCA1")
set.seed(2021)
train.idx <- sample(1:dim(X)[1], dim(X)[1] * 2/3)
foldId <- sample(1:10, length(train.idx), replace = TRUE)
# Fit LASSO
registerDoParallel(nCores)
fit.lasso  <- cv.glmnet(X[train.idx, -y.idx], X[train.idx, y.idx],
family = "gaussian",
alpha = 1,
type.measure = "deviance",
nfolds = 10,
foldid = foldId,
nlambda = 20,
parallel = TRUE)
stopImplicitCluster()
l1_beta <- coef(fit.lasso, s = "lambda.min")[-1]
registerDoParallel(nCores)
fit.nlasso  <- cv.glmnet(X[train.idx, -y.idx], X[train.idx, y.idx],
family = "gaussian",
alpha = 1,
lower.limits = 0,
type.measure = "deviance",
nfolds = 10,
foldid = foldId,
nlambda = 20,
parallel = TRUE)
stopImplicitCluster()
nl1_beta <- coef(fit.nlasso, s = "lambda.min")[-1]
registerDoParallel(nCores)
fit.nlasso2  <- cv.glmnet(X[train.idx, -y.idx], X[train.idx, y.idx],
family = "gaussian",
alpha = 1,
upper.limits = 0,
type.measure = "deviance",
nfolds = 10,
foldid = foldId,
nlambda = 20,
parallel = TRUE)
stopImplicitCluster()
nl1_beta2 <- coef(fit.nlasso2, s = "lambda.min")[-1]
# Fit Q1
quantile_beta <- matrix(NA, nrow = length(tau.seq),
ncol = (dim(X)[2] - 1))
quantile_cvm <- quantile_sd <- quantile_b0 <- numeric()
rownames(quantile_beta) <- paste("tau = ", tau.seq)
for(j in 1:length(tau.seq)) {
#cat(paste(j))
registerDoParallel(nCores)
fit.quantile       <- tune_xrnet(X[train.idx, -y.idx], X[train.idx, y.idx],
family = "gaussian",
penalty_main = define_penalty("q1", tau.seq[j],
num_penalty = 20),
loss = "deviance",
nfolds = 10,
foldid = foldId,
parallel = TRUE)
stopImplicitCluster()
quantile_b0[j] <- coef(fit.quantile)$beta0
quantile_beta[j, ] <- coef(fit.quantile)$betas[, , 1]
quantile_cvm[j]    <- min(fit.quantile$cv_mean)
quantile_sd[j] <- fit.quantile$cv_sd[which.min(fit.quantile$cv_mean)]
}
q1_beta <- quantile_beta[which.min(quantile_cvm), ]
q1_tau <- which.min(quantile_cvm)
# Figure 4
plot(quantile_cvm, ylab = "Mean Cross Validation Error",
xlab = expression(tau), xaxt = "n", type = "o", pch = 19)
axis(1, at = seq(2, 19, by = 2), labels = tau.seq[seq(2, 19, by = 2)])
min(fit.nlasso$cvm)
# Get linear predictors to calculated predicted R2
l1_eta <- cbind(1, X[-train.idx, -y.idx]) %*% coef(fit.lasso, s = "lambda.min")
nl1_eta <- cbind(1, X[-train.idx, -y.idx]) %*% coef(fit.nlasso, s = "lambda.min")
nl1_eta2 <- cbind(1, X[-train.idx, -y.idx]) %*% coef(fit.nlasso2, s = "lambda.min")
q1_eta <- quantile_b0[which.min(quantile_cvm)] + X[-train.idx, -y.idx] %*% q1_beta
y.test <- X[-train.idx, y.idx]
1 - sum((y.test - l1_eta)^2) / sum((y.test - mean(y.test))^2)
1 - sum((y.test - nl1_eta)^2) / sum((y.test - mean(y.test))^2) # Non-negative LASSO R2
1 - sum((y.test - nl1_eta2)^2) / sum((y.test - mean(y.test))^2) # Non-positive LASSO R2
1 - sum((y.test - q1_eta)^2) / sum((y.test - mean(y.test))^2) # asymLASSO R2
sum(q1_beta > 0) # asymLASSO
sum(q1_beta < 0) # asymLASSO
