# Analysis of DP Results - Robustness Tests
# Created: 06/30/19; Last Modified: 06/30/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)))

# We will eventually require a mechanism 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).
Robustness_Tests_ID <- 2
# NOTE 063019: 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 Robustness_
# Tests_ID below (note that said variable simply indexes a location in this vector).
Robustness_Tests_Types <- 
  c("1_Risk_T2DM", "2_Sleep_Time", "3_No_Savings", "4_Basic_Income")

# 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.
  
# 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, "Robustness_Tests/", 
           Robustness_Tests_Types[Robustness_Tests_ID], 
           "/MATLAB_Output/Optimal_Policies/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))

Opt_Pol.02 %<>% mutate(Gender = Person_Type[1], Ethnicity = Person_Type[2], 
                       Education = Person_Type[3])
Opt_Pol.02b %<>% mutate(Gender = Person_Type[1], Ethnicity = Person_Type[2], 
                        Education = Person_Type[3])

# 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:21, 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("Gender", "Ethnicity", "Education", "Emp_Status", "Age", 
                          rep(c("LTPA", "Diet", "Meds", "LMP"), 4))
Opt_Pol.02b %<>%
  arrange(., Age_1) %>%
  .[,c(19:21, 1:2, 3, 7, 11, 15, 4, 8, 12, 16, 5, 9, 13, 17, 6, 10, 14, 18)]
colnames(Opt_Pol.02b) <- c("Gender", "Ethnicity", "Education", "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, "Robustness_Tests/", 
                     Robustness_Tests_Types[Robustness_Tests_ID], 
                     "/MATLAB_Output/Value_Functions/ValFun_", Case, ".dat"), 
              sep = ",", header = T, stringsAsFactors = F))

# * Worked Examples ====

# 070119 - Worked examples
load(file = paste0(RDA, "Worked_Examples/","WE_T0", LM, "_RT_",
                  Case, ".Rda"))
load(file = paste0(RDA, "Worked_Examples/", "WE_T1", LM, "_RT_",
                  Case, ".Rda"))
Val_Fun_WE <- data.frame(Val_Fun.00)

WE_T0 <- merge(RWD_WE_T0, Val_Fun_WE[, c("ID_States_RefStates", "Value_Fun")],
        by.x = "ID_States", by.y = "ID_States_RefStates", all.x = TRUE)

OP_T0 <- subset(Opt_Pol.00, ID_States_Opt_Pol_01 == WE_T0$ID_States[1])

WE_T1 <- lapply(WE_T1, FUN = function(x)
  merge(x, Val_Fun_WE[, c("ID_States_RefStates", "Value_Fun")],
        by.x = "ID_States", by.y = "ID_States_RefStates", all.x = TRUE))

wb = createWorkbook()
addWorksheet(wb, "OP_T0"); addWorksheet(wb, "RWD_T0")
for (j in 1:24){addWorksheet(wb, paste0("RWD_T1_", j))}
writeData(wb, sheet = 1, OP_T0, colNames = T, keepNA = T)
writeData(wb, sheet = 2, WE_T0, colNames = T, keepNA = T)
for (j in 1:24){writeData(wb, sheet = j + 2, WE_T1[[j]],
                          colNames = T, keepNA = T)}
saveWorkbook(wb, paste0(Model_Path, "/Report_Tables/Worked_Examples/",
                        "WE", LM, "_RT_", Case, ".xlsx"), overwrite = T)

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(Gender = Person_Type[1], Ethnicity = Person_Type[2], 
                       Education = Person_Type[3])

# https://stackoverflow.com/questions/20466800/non-numeric-argument-to-mathematical-function
Val_Fun.02 <- data.frame(Val_Fun.02[,c(7:9, 2:1, 3:6)])
colnames(Val_Fun.02)[4:5] <- 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[6:9] <- 
  # 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[6:9], function(x) format(round(x, 0), big.mark = ",", justify = "r"))

setnames(Val_Fun.02, 6:9, c("Healthy", "Newly-Diagnosed T2DM",
                             "T2DM, Poor Adherence", "T2DM, Good Adherence"))

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

wb_OP = createWorkbook()
sht_OP = addWorksheet(wb_OP, "Sheet_1")
dat_OP = Opt_Pol.02
# 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/Robustness_Tests/",
                    Robustness_Tests_Types[Robustness_Tests_ID], 
                    "/OptPolicy_", Case, "_", 
                    Robustness_Tests_Types[Robustness_Tests_ID], ".xlsx"),
             overwrite = T)

wb_VF = createWorkbook()
sht_VF = addWorksheet(wb_VF, "Sheet_1")
dat_VF = Val_Fun.02
# 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/Robustness_Tests/",
                    Robustness_Tests_Types[Robustness_Tests_ID], 
                    "/ValFun_", Case, "_", 
                    Robustness_Tests_Types[Robustness_Tests_ID], ".xlsx"), 
             overwrite = T)