Result1_Ridge

res_df       = readRDS("results/ridge_pve0.5.RDS")
p_list       = c(50,100,200,500,1000,2000)
method_list  = c("Mr.ASH","VarBVS","BayesB","Blasso","SuSiE","E-NET","Lasso","Ridge","SCAD","MCP","L0Learn","Ridge.opt")
col          = gg_color_hue(13)[1:11]
sdat = data.frame()
for (i in 1:6) {
  sdat       = rbind(sdat, data.frame(pred = colMeans(matrix(res_df[[i]]$pred,20,12)),
                                      time = colMeans(matrix(res_df[[i]]$time,20,12)),
                                      p = p_list[i],
                                      fit = method_list))
}
sdat$fit    = factor(sdat$fit, levels =  c("Mr.ASH","E-NET","Lasso","Ridge",
                                           "SCAD","MCP","L0Learn",
                                           "VarBVS","BayesB","Blasso","SuSiE",
                                           "Ridge.opt"))
sdat1       = sdat[sdat$fit %in% c("Mr.ASH","E-NET","Lasso","Ridge","Ridge.opt"),]
sdat2       = sdat[sdat$fit %in% c("Mr.ASH","E-NET","SCAD","MCP","L0Learn","Ridge.opt"),]
sdat3       = sdat[sdat$fit %in% c("Mr.ASH","VarBVS","BayesB","Blasso","SuSiE","Ridge.opt"),]

p1 = ggplot(sdat1) + geom_line(aes(x = p, y = pred, color = fit)) +
  geom_point(aes(x = p, y = pred, color = fit, shape = fit), size = 2.5) +
  theme_cowplot(font_size = 14) +
  scale_x_continuous(trans = "log10", breaks = p_list) +
  labs(y = "predictior error (rmse / sigma)", x = "number of coefficients (p)") +
  theme(axis.line = element_blank(),
        plot.title = element_text(hjust = 0.5)) +
  scale_color_manual(values = c(col[c(1,2,3,4)],"gray50")) +
  scale_shape_manual(values = c(19,17,24,25,15)) +
  scale_y_continuous(trans = "log10", limits = c(1.04,1.46), breaks = c(1.1,1.2,1.3,1.4))
p2 = ggplot(sdat2) + geom_line(aes(x = p, y = pred, color = fit)) +
  geom_point(aes(x = p, y = pred, color = fit, shape = fit), size = 2.5) +
  theme_cowplot(font_size = 14) +
  scale_x_continuous(trans = "log10", breaks = p_list) +
  labs(y = "", x = "number of coefficients (p)") +
  theme(axis.line = element_blank(),
        plot.title = element_text(hjust = 0.5)) +
  scale_color_manual(values = c(col[c(1,2,5,6,7)],"gray50")) +
  scale_shape_manual(values = c(19,17,9,3,11,15)) +
  scale_y_continuous(trans = "log10", limits = c(1.04,1.46), breaks = c(1.1,1.2,1.3,1.4))
p3 = ggplot(sdat3) + geom_line(aes(x = p, y = pred, color = fit)) +
  geom_point(aes(x = p, y = pred, color = fit, shape = fit), size = 2.5) +
  theme_cowplot(font_size = 14) +
  scale_x_continuous(trans = "log10", breaks = p_list) +
  labs(y = "", x = "number of coefficients (p)") +
  theme(axis.line = element_blank(),
        plot.title = element_text(hjust = 0.5)) +
  scale_color_manual(values = c(col[c(1,8,9,10,11)],"gray50")) +
  scale_shape_manual(values = c(19,4,5,7,8,15)) +
  scale_y_continuous(trans = "log10", limits = c(1.04,1.46), breaks = c(1.1,1.2,1.3,1.4))
fig_main  = plot_grid(p1,p2,p3, nrow = 1, rel_widths = c(0.35,0.35,0.3))
title     = ggdraw() + draw_label("Prediction Error (log-scale)", fontface = 'bold', size = 20) 
subtitle  = ggdraw() + draw_label("Scenario: IndepGauss + Normal, n = 500, p = 50,100,200,500,1000,2000, pve = 0.5", fontface  = 'bold', size = 18) 
fig       = plot_grid(title,subtitle,fig_main, ncol = 1, rel_heights = c(0.1,0.06,0.95))
fig

p4  = my.box2(sdat[sdat$fit != "Ridge.opt",], "fit", "time",
             gg_color_hue(13)[1:11]) +
  theme(axis.line    = element_blank(),
        axis.text.x  = element_text(angle = 45,hjust = 1),
        legend.position = "none") +
  scale_y_continuous(trans = "log10")
title     = ggdraw() + draw_label("Computation time (log-scale)", fontface = 'bold', size = 20) 
subtitle  = ggdraw() + draw_label("Scenario: IndepGauss + Normal, n = 500, p = 50,100,200,500,1000,2000, pve = 0.5", fontface  = 'bold', size = 18) +
  labs(y = "computation time (sec)")
p0        = ggplot() + geom_blank() + theme_cowplot() + theme(axis.line = element_blank())
fig_main  = plot_grid(p0,p4,p0, nrow = 1, rel_widths = c(0.3,0.6,0.3))
fig       = plot_grid(title,subtitle,fig_main, ncol = 1, rel_heights = c(0.1,0.06,0.95))
fig

library(Matrix); library(ggplot2); library(cowplot); library(susieR); library(BGLR);
library(glmnet); library(mr.ash); library(ncvreg); library(L0Learn); library(varbvs);
standardize = FALSE
source('code/method_wrapper.R')
source('code/sim_wrapper.R')
tdat1        = list()
n            = 500
p_range      = c(50,100,200,500,1000,2000)
method_list  = c("varbvs","bayesb","blasso","susie","enet","lasso","ridge","scad","mcp","l0learn")
method_list2 = c("mr.ash", method_list, "ridge.opt")
method_num   = length(method_list) + 2
iter_num     = 20
pred         = matrix(0, iter_num, method_num);
time         = matrix(0, iter_num, method_num);
colnames(pred) <- colnames(time) <- method_list2

for (iter in 1:6) {
p               = p_range[iter]
for (i in 1:20) {
  data          = simulate_data(n, p, s = p, seed = i, signal = "normal", pve = 0.5)
 
  for (j in 1:length(method_list)) {
    fit.method    = get(paste("fit.",method_list[j],sep = ""))
    fit           = fit.method(data$X, data$y, data$X.test, data$y.test, seed = i)
    pred[i,j+1]   = fit$rsse / data$sigma / sqrt(n)
    time[i,j+1]   = fit$t
    
    if (method_list[j] == "lasso") {
      lasso.path.order = mr.ash:::path.order(fit$fit$glmnet.fit)
      lasso.beta       = as.vector(coef(fit$fit))[-1]
      lasso.time       = c(fit$t, fit$t2)
    }
  }
 
  fit         = fit.mr.ash(data$X, data$y, data$X.test, data$y.test, seed = i,
                           sa2 = (2^((0:19) / 5 / 2^iter) - 1)^2)
  pred[i,1]   = fit$rsse / data$sigma / sqrt(n)
  time[i,1]   = fit$t
  
  fit         = fit.ridge.opt(data$X, data$y, data$X.test, data$y.test, data$sigma, seed = i)
  pred[i,method_num] = fit$rsse / data$sigma / sqrt(n)
  time[i,method_num] = -Inf
  cat(i," ")
}
cat("\n")
print(c(colMeans(pred)))
tdat1[[iter]] = data.frame(pred = c(pred), time = c(time),
                           fit = rep(method_list2, each = 20))
}