# ######################################################
#   AnalyzeExpertLens_Test.R
#   Pane (adapted from Smart) January 04 2021
#   Creates OPTIC Delphi Visualization Examples
#   data input: Excel file with data (cleaned)
#
#   data output: graphics
######################################################
#### Load R packages ####
library(readxl) # this will allow us to read in excel files
library(tidyverse) # used for data cleaning
library(grid) # used to insert additional text on graphic

#### Load Functions ####
clean_panel = function(data, panel){
  ####################################
  # Load Function to clean panel data
  # Inputs:
  # data - is the data.frame
  # panel - is the panel you want to produce a figure for. 
  # panel cont. - Possible options for panel are "A" or "B"
  ####################################

  # drop empty rows
  data = data %>%
    filter(User != "") # dplyr::filter() keeps everything not equal to the condition
  
  # label outcomes
  if(panel == "A"){
    data = data %>%
      # mutate will create a new variable
      # case_when is a robust version of ifelse()
      mutate(Outcome = case_when(Outcome == 1 ~ "Naloxone Distribution",
                                 Outcome == 2 ~ "OUD Prevalence",
                                 Outcome == 3 ~ "Nonfatal Overdose",
                                 Outcome == 4 ~ "Fatal Overdose",
                                 TRUE ~ NA_character_) %>% 
               # create a labeled factor variable
               factor(., levels = c("Naloxone Distribution", "OUD Prevalence", "Nonfatal Overdose",
                                    "Fatal Overdose")))
  }
  if(panel == "B"){
    data = data %>%
      mutate(Outcome = case_when(Outcome == 1 ~ "Acceptability",
                                 Outcome == 2 ~ "Feasibility",
                                 Outcome == 3 ~ "Affordability",
                                 Outcome == 4 ~ "Equity",
                                 TRUE ~ NA_character_) %>% 
               factor(., levels = c("Acceptability", "Feasibility", "Affordability",
                                    "Equity")))
  }
  
  #Label Policy Variables
  data = data %>%
    mutate(Policy = case_when(Policy == 6  ~ "Liability protections (Rx)",
                              Policy == 7 ~ "Provide educ/training (Rx)",
                              Policy == 8 ~ "Third Party Rx",
                              Policy == 9 ~ "Rx Mandate (Opioid Dose Only)",
                              Policy == 10 ~ "Rx Mandate (OUD, OD history)",
                              Policy == 11 ~ "Standing Order (not statewide)",
                              Policy == 12 ~ "Statewide Order",
                              Policy == 13 ~ "Pharmacist Rx Authority",
                              (Policy == 14 & panel == "A") |
                                (Policy == 15 & panel == "B") ~ "OTC Supply",
                              (Policy == 15 & panel == "A") |
                                (Policy == 16 & panel == "B") ~ "Liability Protections (disp)",
                              (Policy == 16 & panel == "A") |
                                (Policy == 17 & panel == "B") ~ "Provide educ/training (disp)",
                              (Policy == 17 & panel == "A") |
                                (Policy == 18 & panel == "B") ~ "Insurance Coverage",
                              (Policy == 18 & panel == "A") |
                                (Policy == 20 & panel == "B") ~ "State Subsidized Insurance",
                              (Policy == 19 & panel == "A") |
                                (Policy == 21 & panel == "B") ~ "Statewide 'Free Naloxone'",
                              (Policy == 20 & panel == "A") |
                                (Policy == 22 & panel == "B") ~ "Liability Protections (admin)",
                              TRUE ~ NA_character_) %>%
             # If you want to change the order of the labels on the graphics,
             # then you should just reorder the following:
             factor(., levels = c("Statewide 'Free Naloxone'",
                                  "State Subsidized Insurance", 
                                  "Insurance Coverage",
                                  "Liability Protections (admin)",
                                  "OTC Supply",
                                  "Pharmacist Rx Authority",
                                  "Statewide Order",
                                  "Standing Order (not statewide)",
                                  "Provide educ/training (disp)",
                                  "Liability Protections (disp)",
                                  "Rx Mandate (OUD, OD history)",
                                  "Rx Mandate (Opioid Dose Only)",
                                  "Third Party Rx",
                                  "Provide educ/training (Rx)",
                                  "Liability protections (Rx)")))
  
  # Sort Data
  data = data %>%
    # sort data
    arrange(Policy, Outcome)
  
  data = data %>%
    gather(., rating_num, rating, R1_Rating:R3_Rating, factor_key = T)
  
  return(data)
}

get_summary = function(data, policy, outcome){
  ####################################
  # Load Function to get summary statistics
  # Inputs:
  # data - is the data.frame
  # policy - is the column name in data representing the relevant policy.
  # outcome - is the column name in data representing the relevant outcome.
  ####################################
  
  # convert variables so can use in dplyr
  policy = enquo(policy)
  outcome = enquo(outcome)

  summary = data %>%
    # replace zeroes as missing rates
    mutate(new_rating = case_when(rating == 0 ~ NA_real_,
                              TRUE ~ rating)) %>%
    # group data
    group_by(!!policy, !!outcome, rating_num) %>%
    # get summaries by group (see group above)
    summarize(median = median(new_rating, na.rm=T),
              LowerIPR = quantile(new_rating, probs=.3, na.rm=T),
              UpperIPR = quantile(new_rating, probs=.7, na.rm=T)) %>%
    # add columns to data.frame
    mutate(IPRCP = (LowerIPR + UpperIPR)/2,
           AssymIndex = abs(5 - IPRCP),
           IPRAS = 2.35 + (1.5 * AssymIndex),
           IPR = UpperIPR - LowerIPR,
           DisagIndex = IPR / IPRAS,
           Disagreement = DisagIndex > 1)
  
  return(summary)
}

graphic_panel = function(summary, panel, rating){
  ####################################
  # Load function to get Panel A graphic
  # Inputs:
  # summary - is the output from get_summary()
  # panel - is the panel you want to produce a figure for. 
  # panel cont. - Possible options for panel are "A" or "B"
  # rating - The rating you want your figures for ("both", "R1_Rating", or "R3_Rating")
  ####################################
  if(rating %in% c("R1_Rating", "R3_Rating", "both")){
    if(panel == "A"){
      scale_y = scale_y_continuous(
        name = "",  #"Decrease     Increase             Decrease     Increase"
        breaks=c(1, 3, 5, 7, 9),
        limits = c(1,9))
    }
    if(panel == "B"){
      scale_y = scale_y_continuous(
        name = "", #"Low           Moderate     High          Low   Moderate   High"
        breaks=c(1, 3, 5, 7, 9),
        limits = c(1,9))  
    }
    if(rating == "both"){
      my_plot = ggplot(summary, aes(x = Policy, 
                                    y = median, 
                                    group = forcats::fct_rev(rating_num), 
                                    colour = forcats::fct_rev(rating_num))) +
        geom_point(position = position_dodge(.95)) + coord_flip() +
        facet_wrap(~ Outcome, nrow=1) +
        theme_bw() + 
        scale_y + xlab("") +
        geom_hline(yintercept = 5) +
        geom_errorbar(aes(ymin=LowerIPR, ymax=UpperIPR), width=.2, 
                      position = position_dodge(.95)) +
        theme(legend.position = "bottom") +
        scale_color_manual(name = "", labels = c("R3_Rating"="Round 3 Rating: Median (30th-70th percentile)",
                                                 "R1_Rating"="Round 1 Rating: Median (30th-70th percentile)"),
                           values=c("#F81C00", "#1A476F"),
                           guide = guide_legend(reverse = TRUE)) +
        geom_vline(xintercept=4.5, colour="#808080", linetype="dashed") +
        geom_vline(xintercept=10.5, colour="#808080", linetype="dashed")
    } else{
      summary = summary %>% filter(rating_num == rating)
      my_plot = ggplot(summary, aes(x = Policy, y = median, colour = rating_num)) +
        geom_point(position = position_dodge(.95)) + coord_flip() +
        facet_wrap(~ Outcome, nrow=1) +
        theme_bw() + 
        scale_y + xlab("") +
        geom_hline(yintercept = 5) +
        geom_errorbar(aes(ymin=LowerIPR, ymax=UpperIPR), width=.2, 
                      position = position_dodge(.95))  #position_dodge(.9)
      if(rating == "R1_Rating"){
        my_plot = my_plot + theme(legend.position = "bottom") +
          scale_color_manual(name = "", 
                             labels = c("Round 1 Rating: Median (30th-70th percentile)"),
                             values = "#1A476F") +
          geom_vline(xintercept=4.5, colour="#808080", linetype="dashed") +
          geom_vline(xintercept=10.5, colour="#808080", linetype="dashed")
      } else if(rating == "R3_Rating"){
        my_plot = my_plot + theme(legend.position = "bottom") +
          scale_color_manual(name = "", 
                             labels = c("Round 3 Rating: Median (30th-70th percentile)"),
                             values = "#F81C00") +
          geom_vline(xintercept=4.5, colour="#808080", linetype="dashed") +
          geom_vline(xintercept=10.5, colour="#808080", linetype="dashed")
      }
      
    }
  } else{
    stop("Please select a valid rating.")
  }
  return(my_plot)
}

# Load function to save graphic to pre-specified directory
save_graphic = function(file, graphic, panel, rating){
  ####################################
  # Load function to get Panel A graphic
  # Inputs:
  # file - is the location/name of your pdf graphic file. 
  # graphic - is the object created from graphic_panel()
  # panel - is the panel you want to produce a figure for. 
  # panel cont. - Possible options for panel are "A" or "B"
  # rating - The rating you want your figures for ("both", "R1_Rating", or "R3_Rating")
  ####################################
  pdf(file = file, width=10, height=7)
  if(panel == "A"){
    print(graphic)
    grid.text("Increase", x = unit(0.37, "npc"), y = unit(0.1, "npc"), rot=0,
              gp = gpar(fontsize = 7, fontface = "bold"))
    grid.text("No Effect", x = unit(0.305, "npc"), y = unit(0.1, "npc"), rot=0,
              gp = gpar(fontsize = 7, fontface = "bold"))
    grid.text("Reduction", x = unit(0.24, "npc"), y = unit(0.1, "npc"), rot=0,
              gp = gpar(fontsize = 7, fontface = "bold"))
    
    grid.text("Increase", x = unit(0.57, "npc"), y = unit(0.1, "npc"), rot=0,
              gp = gpar(fontsize = 7, fontface = "bold"))
    grid.text("No Effect", x = unit(0.505, "npc"), y = unit(0.1, "npc"), rot=0,
              gp = gpar(fontsize = 7, fontface = "bold"))
    grid.text("Reduction", x = unit(0.44, "npc"), y = unit(0.1, "npc"), rot=0,
              gp = gpar(fontsize = 7, fontface = "bold"))
    
    grid.text("Increase", x = unit(0.765, "npc"), y = unit(0.1, "npc"), rot=0,
              gp = gpar(fontsize = 7, fontface = "bold"))
    grid.text("No Effect", x = unit(0.70, "npc"), y = unit(0.1, "npc"), rot=0,
              gp = gpar(fontsize = 7, fontface = "bold"))
    grid.text("Reduction", x = unit(0.635, "npc"), y = unit(0.1, "npc"), rot=0,
              gp = gpar(fontsize = 7, fontface = "bold"))
    
    grid.text("Increase", x = unit(0.96, "npc"), y = unit(0.1, "npc"), rot=0,
              gp = gpar(fontsize = 7, fontface = "bold"))
    grid.text("No Effect", x = unit(0.895, "npc"), y = unit(0.1, "npc"), rot=0,
              gp = gpar(fontsize = 7, fontface = "bold"))
    grid.text("Reduction", x = unit(0.83, "npc"), y = unit(0.1, "npc"), rot=0,
              gp = gpar(fontsize = 7, fontface = "bold"))

  }
  if(panel == "B"){
    print(graphic + geom_rect(aes(xmin = 0, xmax=16, ymin = 4 , ymax = 6), 
                              fill = "grey", color=alpha("grey",0), alpha=0.04), show.legend=F )
    grid.text("High", x = unit(0.37, "npc"), y = unit(0.1, "npc"), rot=0,
              gp = gpar(fontsize = 9, fontface = "bold"))
    grid.text("Moderate", x = unit(0.305, "npc"), y = unit(0.1, "npc"), rot=0,
              gp = gpar(fontsize = 9, fontface = "bold"))
    grid.text("Low", x = unit(0.24, "npc"), y = unit(0.1, "npc"), rot=0,
              gp = gpar(fontsize = 9, fontface = "bold"))
    
    grid.text("High", x = unit(0.57, "npc"), y = unit(0.1, "npc"), rot=0,
              gp = gpar(fontsize = 9, fontface = "bold"))
    grid.text("Moderate", x = unit(0.505, "npc"), y = unit(0.1, "npc"), rot=0,
              gp = gpar(fontsize = 9, fontface = "bold"))
    grid.text("Low", x = unit(0.44, "npc"), y = unit(0.1, "npc"), rot=0,
              gp = gpar(fontsize = 9, fontface = "bold"))
    
    grid.text("High", x = unit(0.765, "npc"), y = unit(0.1, "npc"), rot=0,
              gp = gpar(fontsize = 9, fontface = "bold"))
    grid.text("Moderate", x = unit(0.7, "npc"), y = unit(0.1, "npc"), rot=0,
              gp = gpar(fontsize = 9, fontface = "bold"))
    grid.text("Low", x = unit(0.635, "npc"), y = unit(0.1, "npc"), rot=0,
              gp = gpar(fontsize = 9, fontface = "bold"))
    
    grid.text("High", x = unit(0.96, "npc"), y = unit(0.1, "npc"), rot=0,
              gp = gpar(fontsize = 9, fontface = "bold"))
    grid.text("Moderate", x = unit(0.895, "npc"), y = unit(0.1, "npc"), rot=0,
              gp = gpar(fontsize = 9, fontface = "bold"))
    grid.text("Low", x = unit(0.83, "npc"), y = unit(0.1, "npc"), rot=0,
              gp = gpar(fontsize = 9, fontface = "bold"))
  }
  dev.off()
}

#### End of functions ####

#### Panel A Analysis ####
# Read in Panel A data - specify sheet.
dta = readxl::read_excel("data/NaloxoneResults_Restructure.xlsx",  sheet = "A_Forstata")

# Clean Panel A data
dta = clean_panel(data = dta, panel = "A") #Note: panel must be A or B

# Get Summary Statistics
summary_stats_PanelA = get_summary(data = dta, policy = Policy, outcome = Outcome)

# Create Panel A Graphic
plot_PanelA_R1 = graphic_panel(summary = summary_stats_PanelA, panel = "A", rating = "R1_Rating") 
plot_PanelA_R3 = graphic_panel(summary = summary_stats_PanelA, panel = "A", rating = "R3_Rating") 
plot_PanelA_both = graphic_panel(summary = summary_stats_PanelA, panel = "A", rating = "both") 

# Save Panel A Graphic
save_graphic(file = "output/PanelA_Round1_Ratings.pdf", graphic = plot_PanelA_R1, 
             panel = "A", rating="R1_Rating")
save_graphic(file = "output/PanelA_Round3_Ratings.pdf", graphic = plot_PanelA_R3, 
             panel = "A", rating="R3_Rating")
save_graphic(file = "output/PanelA_BothRounds_Ratings.pdf", graphic = plot_PanelA_both, 
             panel = "A", rating="both")

#### Panel B Analysis ####
# Read in Panel B data - specify sheet.
dta_B = readxl::read_excel("data/NaloxoneResults_Restructure.xlsx", sheet = "B_Forstata")

# Clean Panel B data
dta_B = clean_panel(data = dta_B, panel = "B") #Note: panel must be A or B

# Get Summary Statistics
summary_stats_PanelB = get_summary(data = dta_B, policy = Policy, outcome = Outcome)

# Create Panel B Graphic
plot_PanelB_R1 = graphic_panel(summary = summary_stats_PanelB, panel = "B", rating = "R1_Rating") 
plot_PanelB_R3 = graphic_panel(summary = summary_stats_PanelB, panel = "B", rating = "R3_Rating") 
plot_PanelB_both = graphic_panel(summary = summary_stats_PanelB, panel = "B", rating = "both") 

# Save Panel A Graphic
save_graphic(file = "output/PanelB_Round1_Ratings.pdf", graphic = plot_PanelB_R1, 
             panel = "B", rating="R1_Rating")
save_graphic(file = "output/PanelB_Round3_Ratings.pdf", graphic = plot_PanelB_R3, 
             panel = "B", rating="R3_Rating")
save_graphic(file = "output/PanelB_BothRounds_Ratings.pdf", graphic = plot_PanelB_both, 
             panel = "B", rating="both")

