# Analysis of DP Results - Sensitivity Analysis
# Created: 09/20/18; Last Modified: 04/24/19

# Objective/Rationale -----------------------------------------------------

# NOTE 042419: Extended the sensitivity analysis to allow for examination of the impact
# of changes to the price of consumption and leisure (i.e., wages), as well as to the
# parameters of the CES utility function. In addition, by setting the default (RSTUDIO) 
# working directory to the DIABETES_MODEL subfolder, I have simplified the process of 
# detecting the location of the active script. 
# 
# NOTE 111818: The code now detects the path where the script is located, eliminating 
# the need for the user to change this manually upon downloading the script for the 
# first time.
# 
# This script facilitates analysis of the solution to the dynamic programming 
# problem. In particular, it is intended to support the sensitivity analysis
# by aggregating results across several runs of the decision model.

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

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

# Call required libraries
library(data.table); library(dplyr); library(httr); library(magrittr); library(purrr)
library(readxl); library(reshape2); library(rvest); library(stringr); library(tibble)
library(tidyr); library(openxlsx); library(xtable); library(here)

Model_Path <- here()

# String for date last modified
LM <- "_042019"
# Filename for output
Filename <- "RWD&TP_Matrices_Construction"
# Paths for output
RDA <- paste0(Model_Path, "/R_Objects/")
# Paths for the optimal policies and value functions generated in MATLAB.
MATLAB <- paste0(Model_Path, "/MATLAB/")

# We want to generate a series of lists offering a short-hand description of the
# characteristics of the decision-maker. We use this later in the table captions.
# Even though we are not expecting to vary the attributes of the decision-maker
# as part of the sensitivity analysis, it is important to always be cognizant of
# what they are.
GEN_K = list("M", "F")
ETH_K = list("C", "AC", "AI")
EDU_K <- 
  data.frame(Index = as.character(seq(1:4)), 
             Levels = c("1. No certificate, diploma or degree", 
                        "2. High school diploma or equivalent", 
                        "3. Postsecondary certificate or diploma below bachelor level", 
                        "4. University certificate, diploma or degree at bachelor level or above"),
             Levels2 = c("< HS", "HS", "Below Bach.", "Bach. or above"),
             Pres_Code = c("Educ_1", "Educ_2", "Educ_3", "Educ_4"),
             Attended_PS = c(rep("NoPS", 2), rep("PS", 2)))

# Create lists to hold the summary tables we will produce below. We will subsequently use this
# construct to combine the results from many different scenarios, enabling us to compare how
# the optimal policies and value functions vary with the decision-maker's characteristics and
# key parameter values.
OP_List <- vector("list", 5); OP_List_B <- vector("list", 5); VF_List <- vector("list", 5)

# We will eventually require a mechanisms that maps results to the lists we just created. To
# this end, we create a matrix containing all possible combinations of personal attributes.
Map <- expand.grid("Education" = EDU_K$Levels2[1:4], "Ethnicity" = ETH_K, "Gender" = GEN_K)

# NOTE 092018: This next line is very important, because it determines the characteristics of
# the decision-maker for the purposes of the decision analysis.
Person_Type <- c(2, 3, 4)
# We use this string to set file paths, since the filenames for the optimal policies and
# value functions include strings denoting the decision-maker's characteristics.
Case = paste0(Person_Type[1], Person_Type[2], Person_Type[3])
# NOTE 042419: Updated sensitivity analysis types.
# NOTE 092018: This line indicates the parameter we are currently considering, and will need
# to be modified every time we examine a new one by changing the value of Sens_Analysis_ID
# below (note that said variable simply indexes a location in this vector).
Sens_Analysis_ID <- 4
Sens_Analysis_Types <- 
  c("1_Functional_Form", "2_Alpha", "3_Discount_Factor", "4_Price_Consump", "5_Wage", 
    "6_ARPD", "7_Savings_Percentage", "8_Dietary_Costs")
# Sens_Analysis_Table indicates the levels of all of the parameters included in the
# sensitivity analysis. Since this script presents results rather than generating them, the
# function of the table is primarily to faciliate appropriate labelling of the table rows.
Sens_Analysis_Table <- 
  data.table(SA_Type = Sens_Analysis_Types,
             Level_1 = c(-10.00, 0.00, 0.00, 0.50, 0.50, 0.00, 0.00, 0.00),
             Level_2 = c(-1.00, 0.25, 0.50, 0.75, 0.75, 0.01, 0.05, 0.50), 
             Level_3 = c(+0.00, 0.50, 0.75, 1.00, 1.00, 0.025, 0.10, 1.00),
             Level_4 = c(+0.50, 0.75, 0.95, 1.50, 1.50, 0.05, 0.15, 1.50), 
             Level_5 = c(+1.00, 1.00, 1.00, 2.00, 2.00, 0.10, 0.20, 2.00))

Table_1 <- data.frame("Sens_Analysis_ID" = c(1:8), "SA_Type" = Sens_Analysis_Types,
                      "Start" = rep(1, 8), "End" = c(1, 1, 1, 3, 3, 1, 1, 1), 
                      stringsAsFactors = F)
Table_2 <- data.frame("Rho" = c(1:3), "Value" = c(-10.00, 0, 1.00))
Table_2$Value <- as.character(Table_2$Value)

# Processing of Optimal Policies and Value Functions ----------------------

# Delta refers here to the scenario under consideration, and corresponds to one of the
# columns in Sens_Analysis_Table. Note that this interpretation of the value of DELTA
# differs from the script developed to present results for the baseline analysis.
# NOTE 091618: We do not necessarily need to use contiguous integers in the FOR loops.
for (a in Table_1$Start[Sens_Analysis_ID]:Table_1$End[Sens_Analysis_ID]){
for (delta in seq(5)){
  
Sens_Anal_String_1 <- 
  if(Sens_Analysis_ID == 4 | Sens_Analysis_ID == 5) {
  paste0("/", "Rho_", Table_2$Value[a])
  } else {""}

Sens_Anal_String_2 <- 
  if(Sens_Analysis_ID == 4 | Sens_Analysis_ID == 5) {
    paste0("_", "Rho_", Table_2$Value[a])
  } else {""}
  
  # MATLAB output its results as a series of DAT files, which we now import into R. To
  # begin with we categorize these as data tables, so that we can play around with some
  # of the innovative functionality associated with this package and class of data
  # structure ...
  # https://stackoverflow.com/questions/11664075/import-dat-file-into-r
  Opt_Pol.00 <- 
    data.table(read.table(
      paste0(MATLAB, "Sensitivity_Analysis/", 
             Sens_Analysis_Types[Sens_Analysis_ID], 
             "/MATLAB_Output/Optimal_Policies/Case_", delta, Sens_Anal_String_1,
             "/OptPolicy_", Case, ".dat"), sep = ",", header = T, 
      stringsAsFactors = F))
      
      # ... like this, for example - there is no need even to apply DPLYR's MUTATE function
      # to create a new variable.
      Opt_Pol.01 <- Opt_Pol.00[, Action:=paste(Per_1_Adherence_LTPA_1, Per_1_Adherence_Diet_1,
                                               Per_1_Adherence_Medication_1, sep = ", ")]
      
      # There's way too much information in the dataframe to present to the reader - 10
      # year intervals is probably sufficient to convey how age affects health-related
      # decision-making. Begin by defining the ages we'd like to include.
      Age_Range <- as.character(seq(25, 85, by = 5))
      # Filter Opt_Pol.01 so it includes only the ages we just listed, individuals who
      # were recently diagnosed OR have had T2DM for 9 years (and can therefore modify
      # their health-related behaviours at this point if they wish to do so), and
      # individual who were either fully compliant with LTPA, diet and medication in
      # their first decade with T2DM, or not compliant at all.
      Opt_Pol.01b <-
        Opt_Pol.01[is.na(Age_1) == F & ((Duration_1 %in% c(NA, "0")) |
                      (Duration_1 == 9 & Action %in% c("A, A, A", "DA, DA, DA"))),]
      Opt_Pol.01 <- 
        Opt_Pol.01[Age_1 %in% Age_Range &
                     ((Duration_1 %in% c(NA, "0")) | 
                        (Duration_1 == 9 & Action %in% c("A, A, A", "DA, DA, DA"))),]
      
      # Now rename some variables. The existing labels are far too long, and make it
      # difficult to fit all the pertinent variables onto a single screen/page.
      Opt_Pol.01 <- Opt_Pol.01 %>%
        rename(OP_Actions = ID_Actions_Opt_Pol_01, OP_States = ID_States_Opt_Pol_01, 
               RefActions = ID_Actions_RefActions, RefStates = ID_States_RefStates,
               Dur_1 = Duration_1, HS_1 = Health_Stat_1, LTPA_1 = Per_1_Adherence_LTPA_1,
               LTPA_2 = Per_2_Adherence_LTPA_1, Diet_1 = Per_1_Adherence_Diet_1,
               Diet_2 = Per_2_Adherence_Diet_1, Meds_1 = Per_1_Adherence_Medication_1,
               Meds_2 = Per_2_Adherence_Medication_1, ES_1 = Emp_Stat_1)
      
      Opt_Pol.01b <- Opt_Pol.01b %>%
        rename(OP_Actions = ID_Actions_Opt_Pol_01, OP_States = ID_States_Opt_Pol_01, 
               RefActions = ID_Actions_RefActions, RefStates = ID_States_RefStates,
               Dur_1 = Duration_1, HS_1 = Health_Stat_1, LTPA_1 = Per_1_Adherence_LTPA_1,
               LTPA_2 = Per_2_Adherence_LTPA_1, Diet_1 = Per_1_Adherence_Diet_1,
               Diet_2 = Per_2_Adherence_Diet_1, Meds_1 = Per_1_Adherence_Medication_1,
               Meds_2 = Per_2_Adherence_Medication_1, ES_1 = Emp_Stat_1)
      
      # Rerrange the order of the columns.
      Opt_Pol.01 <- Opt_Pol.01[,c(10:12, 19, 13:18, 5:8, 20)]
      Opt_Pol.01b <- Opt_Pol.01b[,c(10:12, 19, 13:18, 5:8, 20)]
      
      Opt_Pol.01 %<>% 
        # Now sort the rows by health status, THEN T2DM duration, THEN employment status,
        # THEN age.
        arrange(., HS_1, Dur_1, ES_1, Age_1) %>%
        # Create a new variable that indicates whether or not an individual who has had
        # T2DM for at least nine years was completely adherent to LTPA, healthy eating
        # habits and prescribed courses of pharmacotherapy in the past.
        mutate(Past_AD = ifelse(Dur_1 != "9", NA, ifelse(Action == "A, A, A", "Yes", "No")))
      
      Opt_Pol.01b %<>% 
        arrange(., HS_1, Dur_1, ES_1, Age_1) %>%
        mutate(Past_AD = ifelse(Dur_1 != "9", NA, ifelse(Action == "A, A, A", "Yes", "No")))
      
      # Convert Opt_Pol.01 back to a data table; for some reason, certain operations appear
      # to cause a data structure's type to change.
      Opt_Pol.01 <- data.table(Opt_Pol.01); Opt_Pol.01b <- data.table(Opt_Pol.01b)
      # Create a new variable that integrates health status, T2DM duration, and past
      # adherence to health-related behaviours. Note that if T2DM duration = NA, the string
      # "NA" will actual appear in the new column (i.e., concatenated with the other 
      # elements just described). 
      Opt_Pol.01 <- Opt_Pol.01[, HS_1:= paste(HS_1, Dur_1, Past_AD, sep = "_")]
      Opt_Pol.01b <- Opt_Pol.01b[, HS_1:= paste(HS_1, Dur_1, Past_AD, sep = "_")]
      # Now select a subset of the original variables, and pull these into a new data table.
      Opt_Pol.02 <- Opt_Pol.01[, c("Age_1", "HS_1", "ES_1", "AD_LTPA", "AD_Diet", 
                                   "AD_Medication", "LMP")]
      Opt_Pol.02b <- Opt_Pol.01b[, c("Age_1", "HS_1", "ES_1", "AD_LTPA", "AD_Diet", 
                                   "AD_Medication", "LMP")]
      
      # Now we want to flatten out the table so that it more closely resembles the way it
      # will ultimately appear in the thesis. We put age in the rows, and everything else
      # in the columns.
      Opt_Pol.02 <- 
        dcast.data.table(Opt_Pol.02, ES_1 + Age_1 ~ HS_1, 
                         value.var = c("AD_LTPA", "AD_Diet", "AD_Medication", "LMP"))
      Opt_Pol.02b <- 
        dcast.data.table(Opt_Pol.02b, ES_1 + Age_1 ~ HS_1, 
                         value.var = c("AD_LTPA", "AD_Diet", "AD_Medication", "LMP"))
      # Eliminate the string "_NA_NA" from the column names, wherever it appears (because
      # it looks silly).
      colnames(Opt_Pol.02) <- gsub("_NA_NA", "", colnames(Opt_Pol.02))
      colnames(Opt_Pol.02b) <- gsub("_NA_NA", "", colnames(Opt_Pol.02b))
      
      # Create new columns to denote the gender, ethnicity, and education of the
      # decision-maker. It turns out we won't actually use these columns, but eliminating
      # them would have required significant changes to the code below (and I might
      # identify a use for them later).
      Sens_Analysis_Table <- data.frame(Sens_Analysis_Table)
      Opt_Pol.02 <- Opt_Pol.02 %>% 
        mutate(Parameter = paste0(Sens_Analysis_Table[Sens_Analysis_ID, "SA_Type"], "; ",
                                  Sens_Analysis_Table[Sens_Analysis_ID, c(delta + 1)]))
      Opt_Pol.02b <- Opt_Pol.02b %>% 
        mutate(Parameter = paste0(Sens_Analysis_Table[Sens_Analysis_ID, "SA_Type"], "; ",
                                  Sens_Analysis_Table[Sens_Analysis_ID, c(delta + 1)]))
      
      # Rearrange the columns, starting with the ones we just created.
      Opt_Pol.02$Age_1 <- as.character(Opt_Pol.02$Age_1)
      Opt_Pol.02b$Age_1 <- as.character(Opt_Pol.02b$Age_1)
      
      Opt_Pol.02 %<>%
        arrange(., Age_1) %>%
        .[,c(19, 1:2, 3, 7, 11, 15, 4, 8, 12, 16, 5, 9, 13, 17, 6, 10, 14, 18)]
      Opt_Pol.02b %<>%
        arrange(., Age_1) %>%
        .[,c(19, 1:2, 3, 7, 11, 15, 4, 8, 12, 16, 5, 9, 13, 17, 6, 10, 14, 18)]
      # Now rename the columns. The first three column names only appear once. The other
      # names, however, should repeat. Rather than manually naming each column, we employ
      # the REP function, looping over the same four names four times (i.e., because
      # there are two possible employment statuses and four possible health statuses 
      # [i.e., healthy, newly-diagnosed with T2DM, established T2DM with previous good
      # adherence, and established T2DM with previous poor adherence]).
      colnames(Opt_Pol.02) <- c("Parameter", "Emp_Status", "Age", 
                                rep(c("LTPA", "Diet", "Meds", "LMP"), 4))
      colnames(Opt_Pol.02b) <- c("Parameter", "Emp_Status", "Age", 
                                rep(c("LTPA", "Diet", "Meds", "LMP"), 4))
      
      # Now we process the value function. We apply more or less the same steps as we did
      # for the optimal policies. The only difference is that there are fewer columns to
      # work with here.
      Val_Fun.00 <- data.table(
        read.table(paste0(MATLAB, "Sensitivity_Analysis/", Sens_Analysis_Types[Sens_Analysis_ID], 
                          "/MATLAB_Output/Value_Functions/Case_", delta, Sens_Anal_String_1,
                          "/ValFun_", Case, ".dat"), 
                   sep = ",", header = T, stringsAsFactors = F))
      
      Val_Fun.01 = Val_Fun.00[, Action:= paste(Per_1_Adherence_LTPA_1, Per_1_Adherence_Diet_1,
                                               Per_1_Adherence_Medication_1, sep = ", ")]
      
      Val_Fun.01 <- Val_Fun.01[Age_1 %in% Age_Range & 
                                 ((Duration_1 %in% c(NA, "0")) | 
                                    (Duration_1 == 9 & Action %in% c("A, A, A", "DA, DA, DA"))),]
      
      Val_Fun.01 <- Val_Fun.01 %>%
        rename(Dur_1 = Duration_1, HS_1 = Health_Stat_1, LTPA_1 = Per_1_Adherence_LTPA_1,
               LTPA_2 = Per_2_Adherence_LTPA_1, Diet_1 = Per_1_Adherence_Diet_1,
               Diet_2 = Per_2_Adherence_Diet_1, Meds_1 = Per_1_Adherence_Medication_1,
               Meds_2 = Per_2_Adherence_Medication_1, ES_1 = Emp_Stat_1)
      
      Val_Fun.01 <- Val_Fun.01[,c(4:6, 13, 7:12, 1, 14)]
      Val_Fun.01 %<>% 
        arrange(., HS_1, Dur_1, ES_1, Age_1) %>%
        mutate(Past_AD = ifelse(Dur_1 != "9", NA, ifelse(Action == "A, A, A", "Yes", "No")))
      
      Val_Fun.01 <- data.table(Val_Fun.01)
      Val_Fun.01 <- Val_Fun.01[, HS_1:= paste(HS_1, Dur_1, Past_AD, sep = "_")]
      
      Val_Fun.02 <- Val_Fun.01[, c("Age_1", "HS_1", "ES_1", "Value_Fun")]
      
      Val_Fun.02 <- dcast.data.table(Val_Fun.02, Age_1 + ES_1 ~ HS_1, 
                                     value.var = c("Value_Fun"))
      colnames(Val_Fun.02) <- gsub("_NA_NA", "", colnames(Val_Fun.02))
      # http://viktoriawagner.weebly.com/uploads/1/4/5/2/14527152/r-big-vegetation-data---iavs-23-nov-2015.pdf
      
      Val_Fun.02 %<>% mutate(Parameter = paste0(
        Sens_Analysis_Table[Sens_Analysis_ID, "SA_Type"], "; ",
        Sens_Analysis_Table[Sens_Analysis_ID, delta + 1]))
      
      # https://stackoverflow.com/questions/20466800/non-numeric-argument-to-mathematical-function
      Val_Fun.02 <- data.frame(Val_Fun.02[,c(7, 2:1, 3:6)])
      colnames(Val_Fun.02)[2:3] <- c("Emp_Status", "Age")
      # NOTE 091618: This function did NOT work when the structure was defined as a data table.
      # Unfortunately, it appears not all function are compatible across both types - despite 
      # suggestions to the contrary.
      # https://www.dummies.com/programming/r/how-to-format-numbers-in-r/
      # https://stackoverflow.com/questions/22070777/represent-numeric-value-with-typical-dollar-amount-format
      # https://stackoverflow.com/questions/3443687/formatting-decimal-places-in-r
      Val_Fun.02[4:7] <- 
        # The rationale for formatting the columns of Val_Fun.02 here is that the formatting
        # ultimately carries over into the tables we will insert into the report. We eliminate
        # all decimal places, using commas to separate every three decimal places when the
        # number consists of four digits or more. We also right-align the values in each cell.
        sapply(Val_Fun.02[4:7], function(x) format(round(x, 0), big.mark = ",", justify = "r"))

      # NOTE 091618: Need to assign the tabulation results to a slot in the lists we prepared
      # during the initialization stage.
      # Remember that delta denotes the scenario under consideration with respect to this
      # component of the sensitivity analysis.
      OP_List[[delta]] <- Opt_Pol.02
      OP_List_B[[delta]] <- Opt_Pol.02b
      VF_List[[delta]] <- Val_Fun.02
}

# Reporting Results -------------------------------------------------------

OP_Table.00 <- data.table(Reduce(function(x,y) rbind(x,y), OP_List, accumulate = FALSE))
OP_Table_B.00 <- data.table(Reduce(function(x,y) rbind(x,y), OP_List_B, accumulate = FALSE))
VF_Table.00 <- data.table(Reduce(function(x,y) rbind(x,y), VF_List, accumulate = FALSE))

OP_Table.00 %<>% .[order(Age, Parameter),]
OP_Table_B.00 %<>% .[order(Age, Parameter),]
VF_Table.00 %<>% .[order(Age, Parameter),]
setnames(VF_Table.00, 3:6, c("Healthy", "Newly-Diagnosed T2DM",
                             "T2DM, Poor Adherence", "T2DM, Good Adherence"))

wb_OP = createWorkbook()
sht_OP = addWorksheet(wb_OP, "Sheet_1")
dat_OP = OP_Table.00
# Add data to the worksheet we just created
writeData(wb_OP, sht_OP, dat_OP, colNames = T)
# Add the percent style to the desired cells
saveWorkbook(wb_OP, paste0(Model_Path, "/Report_Tables/Sensitivity_Analysis/",
                        Sens_Analysis_Types[Sens_Analysis_ID], "/OptPolicy_",
                        Case, "_", Sens_Analysis_Types[Sens_Analysis_ID], 
                        Sens_Anal_String_2,".xlsx"), 
             overwrite = T)

wb_VF = createWorkbook()
sht_VF = addWorksheet(wb_VF, "Sheet_1")
dat_VF = VF_Table.00
# Add data to the worksheet we just created
writeData(wb_VF, sht_VF, dat_VF, colNames = T)
# Add the percent style to the desired cells
saveWorkbook(wb_VF, paste0(Model_Path, "/Report_Tables/Sensitivity_Analysis/",
                           Sens_Analysis_Types[Sens_Analysis_ID], "/ValFun_",
                           Case, "_", Sens_Analysis_Types[Sens_Analysis_ID], 
                           Sens_Anal_String_2, ".xlsx"), 
             overwrite = T)
}