Result10_Correlation

res_df       = readRDS("results/equicorr.RDS")
method_list  = c("Mr.ASH","VarBVS","BayesB","Blasso","SuSiE","E-NET","Lasso","Ridge","SCAD","MCP","L0Learn",
                 "Mr.ASH.order","Mr.ASH.init")
method_level = c("Mr.ASH","Mr.ASH.order","Mr.ASH.init","E-NET","Lasso","Ridge",
                 "SCAD","MCP","L0Learn",
                 "VarBVS","BayesB","Blasso","SuSiE")
col          = gg_color_hue(13)[c(1,12,13,2:11)]
rho_list     = c(0.3,0.6,0.9,0.95)
for (i in 1:4) {
res_df[[i+1]]$fit   = rep(method_list, each = 20)
res_df[[i+1]]$fit   = factor(res_df[[i+1]]$fit, levels =  method_level)
}
pp = list()
for (i in 1:4) {
  d       = res_df[[i+1]]
  pp[[i]] = my.box2(d, "fit", "pred", cols = col) +
  theme(axis.line    = element_blank(),
        axis.text.x  = element_text(angle = 45,hjust = 1),
        legend.position = "none") +
  geom_hline(yintercept = mean(d$pred[d$fit == "Mr.ASH.init"]), col = col[3],
             linetype = "dotted", size = 1.5) +
  scale_y_continuous(trans = "log10", breaks = c(1,1.2,1.4,1.6,1.8,2.0)) +
  coord_cartesian(ylim = c(1,1.6))
  subtitle  = ggdraw() + draw_label(paste(paste("rho = ",rho_list[i], sep = ""),"", sep = ""),
                                    fontface  = 'bold', size = 18)
  pp[[i]]   = plot_grid(subtitle, pp[[i]], ncol = 1, rel_heights = c(0.06,0.95))
}
fig_main  = plot_grid(pp[[1]],pp[[2]],pp[[3]],pp[[4]], nrow = 2, rel_widths = c(0.3,0.3,0.3,0.3))
title     = ggdraw() + draw_label("Prediction Error (log-scale)", fontface = 'bold', size = 20) 
subtitle  = ggdraw() + draw_label("Scenario: EquiCorrGauss + SparseNormal, n = 500, p = 2000, pve = 0.5", fontface = 'bold', size = 18)
fig       = plot_grid(title,subtitle,fig_main, ncol = 1, rel_heights = c(0.04,0.06,0.95))
fig

tdat1        = list()
n            = 500
p            = 2000
s            = 20
rho_list     = c(0.3,0.6,0.9,0.95)
method_list  = c("varbvs","bayesb","blasso","susie","enet","lasso","ridge","scad2","mcp2","l0learn")
method_list2 = c("mr.ash", method_list, "mr.ash.order", "mr.ash.init")
method_num   = length(method_list2)
iter_num     = 20
pred         = matrix(0, iter_num, method_num); colnames(pred) = method_list2
time         = matrix(0, iter_num, method_num); colnames(time) = method_list2


for (iter in 1:4) {
rho        = rho_list[iter]
for (i in 1:20) {
  data          = simulate_data(n, p, s = s, seed = i, signal = "normal", rho = rho,
                                design = "equicorrgauss", 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) / 20) - 1)^2)
  pred[i,1]   = fit$rsse / data$sigma / sqrt(n)
  time[i,1]   = fit$t
  
  fit         = fit.mr.ash2(data$X, data$y, data$X.test, data$y.test, seed = i,
                            update.order = lasso.path.order,
                            sa2 = (2^((0:19) / 5) - 1)^2)
  pred[i,j+2] = fit$rsse / data$sigma / sqrt(n)
  time[i,j+2] = fit$t + lasso.time[2]
  
  fit         = fit.mr.ash2(data$X, data$y, data$X.test, data$y.test, seed = i,
                            beta.init = lasso.beta,
                            sa2 = (2^((0:19) / 5) - 1)^2)
  pred[i,j+3] = fit$rsse / data$sigma / sqrt(n)
  time[i,j+3] = fit$t + lasso.time[1]
  
  print(c(pred[i,]))
}
tdat1[[iter]] = data.frame(pred = c(pred), time = c(time), fit = rep(method_list2, each = 20))
}