library(BGGM)
library(ggplot2)
cors <- cor(BGGM::ptsd)
pcors <- corpcor::cor2pcor(cors)
pcors <- ifelse(abs(pcors) < 0.05, 0, pcors)
adj <- ifelse(pcors == 0, 0, 1)
cors <- corpcor::pcor2cor(pcors)



n <- c(100, 250, 500, 1000, 5000, 10000)

res_i <- list()
res_j <- list()

for(j in 1:length(n)){

  for(i in 1:50){

    print(paste("iteration", i, "condition", j))

    Y <- MASS::mvrnorm(n = n[j], rep(0, 20), Sigma = cors)
    bggm_fit <-  explore(Y, prior_sd = 0.40, iter = 5000)
    garbage <- capture.output( bdgraph_fit <- BDgraph::bdgraph(Y) )
    bdgraph_sel <- BDgraph::select(bdgraph_fit)
    sel_50 <- select(bggm_fit, BF_cut = 1.00001)
    sel_3 <- select(bggm_fit, BF_cut = 3)
    bggm_50_res <- BDgraph::compare(adj, sel_50$Adj_10)[6:8,2]
    bggm_3_res <- BDgraph::compare(adj, sel_3$Adj_10)[6:8,2]
    bdgraph_res <- BDgraph::compare(adj, bdgraph_sel)[6:8,2]


res <- data.frame(method = c("BGGM_50", "BGGM_50", "BGGM_50",
                             "BGGM_3", "BGGM_3", "BGGM_3",
                      "BDgraph", "BDgraph", "BDgraph"),
           measure = c("SPC", "SN", "MCC",
                       "SPC", "SN", "MCC",
                       "SPC", "SN", "MCC"),
           score = c(bggm_50_res,
                     bggm_3_res,
                     bdgraph_res),
           N = n[j])

res_i[[i]] <- res

}

  res_j[[j]] <- do.call(rbind.data.frame, res_i)



}

# save(res_j, file = "res_j.Rdata")


windowsFonts(times = windowsFont("Times New Roman"))

dat_res <- do.call(rbind.data.frame, res_j)


# specifcity
dat_res_spc <- reshape::melt(with(subset(dat_res, measure == "SPC"),
     tapply(score, list(method, N), mean)))

dat_res_spc$sd <-as.numeric( with(subset(dat_res, measure == "SPC"),
                       tapply(score, list(method, N), sd))) / sqrt(50)


dat_res_spc$X1 <- factor(dat_res_spc$X1,
                      levels = c("BDgraph", "BGGM_3", "BGGM_50"),
                      labels = c("BDgraph",
                                 "BGGM: BF = 3",
                                 "BGGM: Post.prob = 0.50"))

plot_spc <- ggplot(dat_res_spc,  aes(x = as.factor(X2),
                       y = value,
                      group = X1)) +
  geom_hline(yintercept = 0.95, linetype = "dotted") +
  geom_ribbon(aes(ymax = value + sd,
                  ymin = value - sd,
                  fill = X1),
              color = NA,
              show.legend = F, alpha = 0.25) +
  geom_line(size = 1.5, alpha = 0.85, aes(color = as.factor(X1))) +
  ylab("Specificity") +
  theme_bw(base_family = "times", base_size = 12) +
  theme(legend.position = "top",
      legend.title = element_blank(),
      panel.grid = element_blank()) +
  scale_x_discrete(expand = c(0.01,0.01), breaks = n, labels = n) +
  guides(colour = guide_legend(override.aes = list(alpha = 1))) +
  xlab(expression("Sample Size ("*italic(p)*" = 20)")) +
  scale_color_manual(values = c("#999999", "#0072B2", "#CC79A7")) +
scale_fill_manual(values = c("#999999", "#0072B2", "#CC79A7"))

plot_spc


# sensitivity
dat_res_sn <- reshape::melt(with(subset(dat_res, measure == "SN"),
                                  tapply(score, list(method, N), mean)))
dat_res_sn$sd <-as.numeric( with(subset(dat_res, measure == "SN"),
                                  tapply(score, list(method, N), sd))) / sqrt(50)

dat_res_sn$X1 <- factor(dat_res_sn$X1,
                         levels = c("BDgraph", "BGGM_3", "BGGM_50"),
                         labels = c("BDgraph",
                                    "BGGM: BF = 3",
                                    "BGGM: Post.prob = 0.50"))

plot_sn <- ggplot(dat_res_sn,  aes(x = as.factor(X2),
                                     y = value,
                                     group = X1,
                                     color = as.factor(X1))) +
  geom_hline(yintercept = 0.95, linetype = "dotted") +
  geom_line(size = 1.5, alpha = 0.85) +
  ylab("Specificity") +
  theme_bw(base_family = "times", base_size = 12) +
  theme(legend.position = "top",
        legend.title = element_blank(),
        panel.grid = element_blank()) +
  scale_x_discrete(expand = c(0.01,0.01), breaks = n, labels = n) +
  guides(colour = guide_legend(override.aes = list(alpha = 1))) +
  xlab(expression("Sample Size ("*italic(p)*" = 20)")) +
  scale_color_manual(values = c("#999999", "#0072B2", "#CC79A7"))

plot_sn



# mcc
dat_res_mcc <- reshape::melt(with(subset(dat_res, measure == "MCC"),
                                 tapply(score, list(method, N), mean)))
dat_res_mcc$sd <-as.numeric( with(subset(dat_res, measure == "MCC"),
                                 tapply(score, list(method, N), sd))) / sqrt(50)

dat_res_mcc$X1 <- factor(dat_res_mcc$X1,
                        levels = c("BDgraph", "BGGM_3", "BGGM_50"),
                        labels = c("BDgraph",
                                   "BGGM: BF = 3",
                                   "BGGM: Post.prob = 0.50"))

plot_mcc <- ggplot(dat_res_mcc,  aes(x = as.factor(X2),
                                   y = value,
                                   group = X1,
                                   color = as.factor(X1))) +
  geom_hline(yintercept = 0.95, linetype = "dotted") +
  geom_line(size = 1.5, alpha = 0.85) +
  ylab("Specificity") +
  theme_bw(base_family = "times", base_size = 12) +
  theme(legend.position = "top",
        legend.title = element_blank(),
        panel.grid = element_blank()) +
  scale_x_discrete(expand = c(0.01,0.01), breaks = n, labels = n) +
  guides(colour = guide_legend(override.aes = list(alpha = 1))) +
  xlab(expression("Sample Size ("*italic(p)*" = 20)")) +
  scale_color_manual(values = c("#999999", "#0072B2", "#CC79A7"))

plot_mcc



qgraph::qgraph(pcors)