# This code can be used to generate the manuscript figures  

library(tidyverse)
library(ggh4x)
library(glue)

# Reading, combining, and writing the 48 combinations of \phi_1 and \theta files
out=c(1:48)
x=c()
for (i in c(1:48)){
  mypath <- file.path("/Users/Journal of Data Science/results",paste("senario-", out[i], ".csv", sep = ""))
  y=read.csv(mypath, header = T)
  x=rbind(x,y)
}
write.csv(x,file="/Users/Journal of Data Science/results/Summary_norm_n_36_var_1_ci.csv", row.names = F)

# Recall stored files
res1=read.csv(file="/Users/simulation_results/Summary_norm_n_31_ci.csv", header = T)
res2=read.csv(file="/Users/simulation_results/Summary_norm_n_31_var_0.25_ci.csv", header = T)
res3=read.csv(file="/Users/simulation_results/Summary_norm_n_31_var_0.5_ci.csv", header = T)
res4=read.csv(file="/Users/simulation_results/Summary_snorm_n_31_ci.csv", header = T)
res5=read.csv(file="/Users/simulation_results/Summary_t_n_31_ci.csv", header = T)

res6=read.csv(file="/Users/simulation_results/Summary_norm_n_63_ci.csv", header = T)
res7=read.csv(file="/Users/simulation_results/Summary_norm_n_63_var_0.25_ci.csv", header = T)
res8=read.csv(file="/Users/simulation_results/Summary_norm_n_63_var_0.5_ci.csv", header = T)
res9=read.csv(file="/Users/simulation_results/Summary_snorm_n_63_ci.csv", header = T)
res10=read.csv(file="/Users/simulation_results/Summary_t_n_63_ci.csv", header = T)

res11=read.csv(file="/Users/simulation_results/Summary_norm_n_127_ci.csv", header = T)
res12=read.csv(file="/Users/simulation_results/Summary_norm_n_127_var_0.5_ci.csv", header = T)
res13=read.csv(file="/Users/simulation_results/Summary_norm_n_127_var_0.25_ci.csv", header = T)
res14=read.csv(file="/Users/simulation_results/Summary_snorm_n_127_ci.csv", header = T)
res15=read.csv(file="/Users/simulation_results/Summary_t_n_127_ci.csv", header = T)


xall <- res1 %>% # Change the file name to plot new results
  gather(key=EST, value = value, c(asy.fm.ind.t, asy.ind.t, asy.boot1.ind.t, asy.boot2f.ind.t, perc.ind.t, perc.boot1.ind.t, perc.boot2f.ind.t,percbc.ind.t,bca.ind.t,
                                   asy.fm.w, asy.w, asy.boot1.w, asy.boot2f.w, perc.w, perc.boot1.w, perc.boot2f.w,percbc.w,bca.w,
                                   asy.fm.ind.l, asy.ind.l, asy.boot1.ind.l, asy.boot2f.ind.l, perc.ind.l, perc.boot1.ind.l, perc.boot2f.ind.l,percbc.ind.l,bca.ind.l,
                                   asy.fm.ind.u, asy.ind.u, asy.boot1.ind.u, asy.boot2f.ind.u, perc.ind.u, perc.boot1.ind.u, perc.boot2f.ind.u,percbc.ind.u,bca.ind.u)) %>% 
  mutate(measure = c(rep("Covarge",9*48),rep("Width",9*48),rep("SL.L",9*48),rep("SL.U",9*48))) %>% 
  mutate(EST.Short = case_when(
    EST%in%c("asy.fm.ind.t","asy.fm.w","asy.fm.ind.l","asy.fm.ind.u") ~ "ASY",
    EST%in%c("asy.ind.t","asy.w","asy.ind.l","asy.ind.u") ~ "NORMB",
    EST%in%c("asy.boot1.ind.t","asy.boot1.w","asy.boot1.ind.l","asy.boot1.ind.u") ~ "NORMB-SBC",
    EST%in%c("asy.boot2f.ind.t","asy.boot2f.w","asy.boot2f.ind.l","asy.boot2f.ind.u") ~ "NORM-FDBC",
    EST%in%c("perc.ind.t","perc.w","perc.ind.l","perc.ind.u") ~ "PERC",
    EST%in%c("perc.boot1.ind.t","perc.boot1.w","perc.boot1.ind.l","perc.boot1.ind.u") ~ "PERC-SBC",
    EST%in%c("perc.boot2f.ind.t","perc.boot2f.w","perc.boot2f.ind.l","perc.boot2f.ind.u") ~ "PERC-FDBC",
    EST%in%c("bca.ind.t","bca.w","bca.ind.l","bca.ind.u") ~ "BCa",
    TRUE ~ "percbc"))


# 4 by 6 plots
xall %>% 
  mutate(measure = factor(measure, levels = c("Covarge", "Width","SL.L","SL.U"),
                          labels = c("Covarge", "Width","SL.L","SL.U"))) %>% 
  filter(!EST.Short%in%c("percbc")) %>%
  ggplot(aes(x=ar.coef, y=value, col=EST.Short)) +
  geom_point(size=0.5) +
  geom_line() +
  labs(x=expression(phi), y="Performance Measure") +
  scale_color_discrete(name="Method") +
  facet_grid(measure~forcats::fct_inorder(glue('theta*" = {r1}"')), scales = "free",labeller = label_parsed ) +
  facetted_pos_scales( y = list(
    measure == "Covarge" ~ scale_y_continuous(limits = c(0.8, 1)),
    measure == "Width" ~ scale_y_continuous(limits = c(0.05, 0.5)),
    measure == "SL.L" ~ scale_y_continuous(limits = c(0, 0.065)),
    measure == "SL.U" ~ scale_y_continuous(limits = c(0, 0.16)))) +
  #facet_grid(measure~'theta*" = "' r1, scales = "free_x",labeller = label_parsed ) +
  theme(panel.background = element_rect(fill = "white",
                                        colour = "lightgray",
                                        size = 0.5, linetype = "solid"),
        panel.grid.major = element_line(size = 0.05,linetype = "solid", colour = "gray"), 
        panel.grid.minor = element_line(size = 0.025, linetype = "solid",colour = "gray"),
        legend.position = "bottom")
