# R Code for Analysis of Canadian Chronic Disease Surveillance System (CCDSS) Data
# Created: 12/22/18; Last Modified: 04/20/19

# Introduction & Rationale ------------------------------------------------

# This code is designed to analyze and prepare visualizations based on CCDSS data
# relating to diabetes, for the purposes of validating output generated by the
# simulation model.

# 012419: Updated code to incorporate results from the third scenario.

# Initialization ----------------------------------------------------------

# The next line of code clears the variable list, if required.
rm(list=ls())

# Call requisite libraries
library(dplyr); library(tidyr); library(reshape2); library(magrittr)
library(data.table); library(ggplot2); library(openxlsx); library(here)

Model_Path <- here()
# String for date last modified
LM <- "_042019"
# Filename for output
Filename <- "CCDSS_Calculations"
# Paths for output
Path <- paste0("C:/Users/starg/Documents/PhD Dissertation/Combined Document/",
               "CCDSS Data/Processed/PHAC_Infobase_CCDSS_")

# Preparing the Data ------------------------------------------------------

# Read PHAC Infobase CCDSS CSV file (NOTE: There does not currently appear to be
# any way to download the file directly from the Infobase at this point in time).
CCDSS.00 <- vector("list", length = 2)
Data_Types <- c("Inc_CR", "Prev_CR")
names(CCDSS.00) <- Data_Types

CCDSS.00[["Inc_CR"]] <- fread(paste0(Path, "Incidence_Crude", ".csv"))
CCDSS.00[["Prev_CR"]] <- fread(paste0(Path, "Prevalence_Crude", ".csv"))
CCDSS.00[["Inc_CR"]] <- CCDSS.00[["Inc_CR"]] %>%
  rename("Age_Group" = `Age Group`, "Lower_CI" = "Rate_95%_CI_Lower",
         "Upper_CI" = "Rate_95%_CI_Upper")
CCDSS.00[["Prev_CR"]] <- CCDSS.00[["Prev_CR"]] %>%
  rename("Age_Group" = `Age Group`, "Lower_CI" = "Rate_95%_CI_Lower",
         "Upper_CI" = "Rate_95%_CI_Upper")

CCDSS.00[["Inc_CR"]]$Age_Group <- gsub("'", "", CCDSS.00[["Inc_CR"]]$Age_Group)
CCDSS.00[["Prev_CR"]]$Age_Group <- gsub("'", "", CCDSS.00[["Prev_CR"]]$Age_Group)

CCDSS.01 <- CCDSS.00 %>% 
  lapply(., FUN = function(x) 
  {x <- x %>%
    filter(Area == "Manitoba", Year == "2015", 
           Age_Group %in% c("20-34", "35-49", "50-64", "65-79", "80+"),
           Sex == "Both sexes") %>%
    select(Age_Group, Rate, Lower_CI, Upper_CI); x})

CCDSS.01$Inc_CR <- mutate(CCDSS.01$Inc_CR, Metric = "Incidence")
CCDSS.01$Prev_CR <- mutate(CCDSS.01$Prev_CR, Metric = "Prevalence")

CCDSS.01$Inc_CR$Age_Group <- as.factor(CCDSS.01$Inc_CR$Age_Group)
CCDSS.01$Prev_CR$Age_Group <- as.factor(CCDSS.01$Prev_CR$Age_Group)
CCDSS.01$Inc_CR[2:4] %<>%
  sapply(., FUN = function(x) sub(",", "", x)) %>%
  sapply(., as.numeric)

Age_Groups <- data.frame(Age = c(26:95))
Age_Groups %<>%
  mutate("Age_Group" =
           case_when(Age < 35 ~ "20-34",
                     Age < 50 ~ "35-49",
                     Age < 65 ~ "50-64",
                     Age < 80 ~ "65-79",
                     TRUE ~ "80+"))

# Summarizing Prevalence Data -------------------------------------------------

Scenario_1P <- read.xlsx(xlsxFile = 
                          paste0(Model_Path, "/MATLAB/Simulation_Results/",
                                 "Tables/", "Sim1_Table.xlsx"), sheet = 1, startRow = 1, 
                        colNames = T, cols = c(1:3))
Scenario_1P$Age <- Scenario_1P$Elapsed_Time + 25
Scenario_1P <- merge(Scenario_1P, Age_Groups, by = "Age", sort = FALSE)

Scenario_2P <- read.xlsx(xlsxFile = 
                          paste0(Model_Path, "/MATLAB/Simulation_Results/",
                                 "Tables/", "Sim2_Table.xlsx"), sheet = 1, startRow = 1, 
                        colNames = T, cols = c(1:3))
Scenario_2P$Age <- Scenario_2P$Elapsed_Time + 25
Scenario_2P <- merge(Scenario_2P, Age_Groups, by = "Age", sort = FALSE)

Sim_Prev.00 <- cbind(Scenario_1P[c("Age", "Age_Group", "Prev_Opt", "Prev_Alt_A")], 
                  Scenario_2P["Prev_Alt_DA"])
Sim_Prev.00 <- filter(Sim_Prev.00, Age != 95)
Sim_Prev.00[3:5] <- Sim_Prev.00[3:5]*100

Sim_Prev.00 %<>%
  group_by(., Age_Group) %>%
  summarise(., Optimal = mean(Prev_Opt),
            Alternate_A = mean(Prev_Alt_A),
            Alternate_DA = mean(Prev_Alt_DA))

Sim_Prev.00 <- merge(Sim_Prev.00, CCDSS.01$Prev_CR[1:4], by = "Age_Group", sort = FALSE)
Sim_Prev.00 <- Sim_Prev.00[, c(1, 5:7, 2:4)]
Sim_Prev.00 <- rename(Sim_Prev.00, "CCDSS" = "Rate", "CCDSS_LowerCI" = "Lower_CI",
                   "CCDSS_UpperCI" = "Upper_CI")

Sim_Prev.01 <- 
  Sim_Prev.00[c("Age_Group", "CCDSS", "Optimal", "Alternate_A", "Alternate_DA")]
Sim_Prev.01 <- rename(Sim_Prev.01, "Alternate (Adhere)" = "Alternate_A",
                      "Alternate (Don't Adhere)" = "Alternate_DA")

Sim_Prev.01 <- melt(Sim_Prev.01, id.vars = "Age_Group", variable.name = "Policy", 
                    value.name = "Prevalence")

# Summarizing Incidence Data ----------------------------------------------

Scenario_1I <- read.xlsx(xlsxFile = 
                           paste0(Model_Path, "/MATLAB/Simulation_Results/",
                                  "Tables/", "Sim1_Table.xlsx"), sheet = 1, startRow = 1, 
                         colNames = T, cols = c(1, 4:5))
Scenario_1I$Age <- Scenario_1I$Elapsed_Time + 25
Scenario_1I <- merge(Scenario_1I, Age_Groups, by = "Age", sort = FALSE)
Scenario_1I[3:4] <- sapply(Scenario_1I[3:4], FUN = function(x) 100*x)

Scenario_2I <- read.xlsx(xlsxFile = 
                           paste0(Model_Path, "/MATLAB/Simulation_Results/",
                                  "Tables/", "Sim2_Table.xlsx"), sheet = 1, startRow = 1, 
                         colNames = T, cols = c(1, 4:5))
Scenario_2I$Age <- Scenario_2I$Elapsed_Time + 25
Scenario_2I <- merge(Scenario_2I, Age_Groups, by = "Age", sort = FALSE)
Scenario_2I[3:4] <- sapply(Scenario_2I[3:4], FUN = function(x) 100*x)

Sim_Inc.00 <- cbind(Scenario_1I[c("Age", "Age_Group", "Inc_Opt", "Inc_Alt_A")], 
                     Scenario_2I["Inc_Alt_DA"])
Sim_Inc.00 <- filter(Sim_Inc.00, Age != 95)

Sim_Inc.00 %<>%
  group_by(., Age_Group) %>%
  summarise(., Optimal = mean(Inc_Opt),
            Alternate_A = mean(Inc_Alt_A),
            Alternate_DA = mean(Inc_Alt_DA))

Sim_Inc.00 <- merge(Sim_Inc.00, CCDSS.01$Inc_CR[1:4], by = "Age_Group", sort = FALSE)
Sim_Inc.00 <- Sim_Inc.00[, c(1, 5:7, 2:4)]
Sim_Inc.00 <- rename(Sim_Inc.00, "CCDSS" = "Rate", "CCDSS_LowerCI" = "Lower_CI",
                      "CCDSS_UpperCI" = "Upper_CI")

Sim_Inc.01 <- 
  Sim_Inc.00[c("Age_Group", "CCDSS", "Optimal", 
               "Alternate_A", "Alternate_DA")]
Sim_Inc.01 <- rename(Sim_Inc.01, "Alternate (Adhere)" = "Alternate_A",
                      "Alternate (Don't Adhere)" = "Alternate_DA")
Sim_Inc.01 <- melt(Sim_Inc.01, id.vars = "Age_Group", variable.name = "Policy", 
                    value.name = "Incidence")

# Visualizing Results --------------------------------------------------------

# References about creating bar charts using ggplot2:
# http://www.r-graph-gallery.com/218-basic-barplots-with-ggplot2/
# http://www.sthda.com/english/wiki/ggplot2-axis-scales-and-transformations
Prev_Validation <-
  ggplot(Sim_Prev.01, 
         aes(x = Age_Group, y = Prevalence, group = Policy,
             fill = Policy, label = sprintf("%0.1f", round(Prevalence, digits = 1)))) +
  # Altering spacing between bars:
  # https://stackoverflow.com/questions/25887066/remove-space-between-bars-ggplot2
  geom_bar(position = position_dodge(width = 0.9), stat = "identity") +
  # Placing labels in bar chart:
  # https://stackoverflow.com/questions/6644997/showing-data-values-on-stacked-bar-chart-in-ggplot2
  # https://stackoverflow.com/questions/40211451/geom-text-how-to-position-the-text-on-bar-as-i-want
  geom_text(
    aes(x = Age_Group, y = Prevalence, group = Policy,
        # Changing number of decimal points displayed:
        # https://stackoverflow.com/questions/45799000/how-to-round-percentage-to-2-decimal-places-in-ggplot2
        # https://stackoverflow.com/questions/38369855/how-to-put-exact-number-of-decimal-places-on-label-ggplot-bar-chart
        label = sprintf("%0.1f", round(Prevalence, digits = 1))), 
    position = position_dodge(width = 0.9), vjust = -0.5, size = 3) +
  labs(y = "Prevalence (%)", x = "Age Group") +
  scale_fill_grey(start = 0.75, end = 0.25) + theme_classic() + 
  # Placement of legend:
  # http://www.sthda.com/english/wiki/ggplot2-legend-easy-steps-to-change-the-position-and-the-appearance-of-a-graph-legend-in-r-software
  theme(legend.position = "bottom") + 
  # Changing axis limits:
  # https://stackoverflow.com/questions/36669095/y-axis-wont-start-at-0-in-ggplot
  scale_y_continuous(expand = c(0,0), limits = c(0, 100))

Inc_Validation <-
  ggplot(Sim_Inc.01, 
         aes(x = Age_Group, y = Incidence, group = Policy,
             fill = Policy, label = sprintf("%0.1f", round(Incidence, digits = 1)))) +
  geom_bar(position = position_dodge(width = 0.9), stat = "identity") +
  geom_text(
    aes(x = Age_Group, y = Incidence, group = Policy,
        label = sprintf("%0.1f", round(Incidence, digits = 1))), 
    position = position_dodge(width = 0.9), vjust = -0.5, size = 2) +
  labs(y = "Incidence (%)", x = "Age Group") +
  scale_fill_grey(start = 0.75, end = 0.25) + theme_classic() + 
  theme(legend.position = "bottom") + 
  scale_y_continuous(expand = c(0,0), limits = c(0, 10000))

ggsave(filename = paste0("C:/Users/starg/Documents/PhD Dissertation/Combined Document/", 
                         "Prev_Validation", LM, ".jpg"), 
       plot = Prev_Validation, width = 15.7, height = 10, units = "cm")
ggsave(filename = paste0("C:/Users/starg/Documents/PhD Dissertation/Combined Document/", 
                         "Inc_Validation", LM, ".jpg"), 
       plot = Inc_Validation, width = 15.7, height = 10, units = "cm")