library(tidyverse)


#Figure 3 in Simulation 1


plot_fig <- function(link, x_name, x_numbers, to_do,title){
  
  load(link)
  data_ <- data.frame(Numbers = x_numbers,
                      Bonferroni = colMeans(mat1),
                      Fisher = colMeans(mat3),
                      Ours = colMeans(mat4))
  
  if(to_do == "Type_1") {
    p <- ggplot(data_, aes(x=Numbers))+
      geom_hline(yintercept = 0.10)+
      scale_y_continuous(breaks=c(0, 0.10, 0.30, 0.50), limits = c(0.0,0.50))+
      geom_line(aes(y = Bonferroni,color = "Bonferroni"),size=2) +
      geom_line(aes(y = Fisher, color =  "MCRTs+F"),size=2) +
      geom_line(aes(y = Ours,color =  "MCRTs+Z"),size=2) +
      geom_point(aes(y = Bonferroni,color = "Bonferroni",shape = "Bonferroni"),size=5)+
      geom_point(aes(y = Fisher,color =  "MCRTs+F",shape =  "MCRTs+F"),size=5)+
      geom_point(aes(y = Ours,color = "MCRTs+Z",shape = "MCRTs+Z"),size=5)+
      #facet_grid(. ~ title)+
      scale_color_manual(name = NULL, values = c("Bonferroni" = "#999999", "MCRTs+F" = "#009E73", "MCRTs+Z" = "#D55E00"))+
      scale_shape_manual(name = NULL, values = c("Bonferroni" =15, "MCRTs+F" = 19, "MCRTs+Z" = 17))+
      xlab(x_name)+
      ylab("Type I error rate")
    p <- p + theme_bw(base_line_size =0.5,base_rect_size = 1.0)
    p <-p + theme(
      legend.position="top",
      text = element_text(size=32, face="bold"),
      axis.title.x = element_text(color="black", size=32, face="bold",vjust = -2),
      axis.title.y = element_text(color="black", size=32, face="bold",vjust = 5),
      plot.margin=unit(c(1.0,1.0,1.0,1.0),"cm"))
  } else if (to_do == "Power"){
    
    p <- ggplot(data_, aes(x=Numbers))+
      geom_line(aes(y = Bonferroni,color = "Bonferroni"),size=2) +
      geom_line(aes(y = Fisher, color =  "MCRTs+F"),size=2) +
      geom_line(aes(y = Ours,color =  "MCRTs+Z"),size=2) +
      geom_point(aes(y = Bonferroni,color = "Bonferroni",shape = "Bonferroni"),size=5)+
      geom_point(aes(y = Fisher,color =  "MCRTs+F",shape =  "MCRTs+F"),size=5)+
      geom_point(aes(y = Ours,color = "MCRTs+Z",shape = "MCRTs+Z"),size=5)+
      scale_color_manual(name = NULL, values = c("Bonferroni" = "#999999", "MCRTs+F" = "#009E73", "MCRTs+Z" = "#D55E00"))+
      scale_shape_manual(name = NULL, values = c("Bonferroni" =15, "MCRTs+F" = 19, "MCRTs+Z" = 17))+
      xlab(x_name)+
      ylab("Power")
    p <- p + theme_bw(base_line_size =0.5,base_rect_size = 1.0)
    p <-p + theme(
      legend.position="top",
      text = element_text(size=32, face="bold"),
      axis.title.x = element_text(color="black", size=32, face="bold",vjust = -2),
      axis.title.y = element_text(color="black", size=32, face="bold",vjust = 5),
      plot.margin=unit(c(1.0,1.0,1.0,1.0),"cm"))
  } else{
    stop("Error")
  }
  ggsave(p,  width=6,height=6, dpi = 300, filename =paste("figures/sim_1/",x_name,"_",to_do,".png"))
  return(p)
}

p <- plot_fig(link="results/sim_1_units_1.RData",  x_name="Number of units", x_numbers=c(100,200,300,400,500),to_do = "Type_1",title)
p <- plot_fig(link="results/sim_1_units_2.RData",  x_name="Number of units", x_numbers=c(100,200,300,400,500),to_do = "Power",title)

p <- plot_fig(link="results/sim_1_steps_1.RData",  x_name="Number of time steps", x_numbers=c(4,6,8,10,12),to_do = "Type_1",title)
p <- plot_fig(link="results/sim_1_steps_2.RData",  x_name="Number of time steps", x_numbers=c(4,6,8,10,12),to_do = "Power",title)

p <- plot_fig(link="results/sim_1_lags_1.RData",  x_name="Time lag", x_numbers=c(0,1,2,3,4),to_do = "Type_1",title)
p <- plot_fig(link="results/sim_1_lags_2.RData",  x_name="Time lag", x_numbers=c(0,1,2,3,4),to_do = "Power",title)

p <- plot_fig(link="results/sim_1_sizes_2.RData",  x_name="Effect size", x_numbers=c(0.01,0.02,0.03,0.04,0.05),to_do = "Power",title)


