# Analysis of DP Results
# Created: 09/13/18; Last Modified: 04/23/19

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

# NOTE 042419: 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 120918: The code now contains an embedded check to verify that optimal policies
# do not entail abandonment of health-related behaviours to which people had 
# previously committed. The results of the test are contained in the aptly named list
# "CHECK"; in the event that one or more entries in the list highlight errors, it may
# suggest issues with the generation of the reward and/or transition probability
# matrices or in the solution of the optimization problem that must be addressed 
# before proceeding.

# 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.

# 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/")
MATLAB <- paste0(Model_Path, "/MATLAB/")

GEN_K = list("M", "F")
GEN_K_DF = data.frame(Index = c(1:2), Levels = unlist(GEN_K))
ETH_K = list("C", "AC", "AI")
ETH_K_DF = data.frame(Index = c(1:3), Levels = unlist(ETH_K))
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", 24); VF_List <- vector("list", 24); CHECK <- vector ("list", 24)

# 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 051219: This second map is intended to facilitate horizontal sorting of
# columns in Excel after we export our results.
Map2 <- expand.grid("Education" = as.integer(as.character(EDU_K$Index)), 
                    "Ethnicity" = ETH_K_DF$Index, "Gender" = GEN_K_DF$Index)
Map_List_OP <- vector("list", 24); Map_List_VF <- vector("list", 24)
for (x in 1:nrow(Map2)){
  Map_List_OP[[x]] <- 
    data.frame("Emp_Status" = c(NA, NA, NA, NA), "Age" = c(NA, NA, NA, NA),
               "LTPA" = c(Map2[x, 1], Map2[x, 2], Map2[x, 3], 1),
               "Diet" = c(Map2[x, 1], Map2[x, 2], Map2[x, 3], 2),
               "Meds" = c(Map2[x, 1], Map2[x, 2], Map2[x, 3], 3),
               "LMP" = c(Map2[x, 1], Map2[x, 2], Map2[x, 3], 4))
  Map_List_VF[[x]] <- 
    data.frame("Emp_Status" = c(NA, NA, NA), "Age" = c(NA, NA, NA),
               "Value" = c(Map2[x, 1], Map2[x, 2], Map2[x, 3]))}

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

# Remember that delta, gamma, and theta index level of educational attainment, ethnicity, and
# gender, respectively.
for (delta in seq(4)){
  for (gamma in seq(3)){
    for (theta in seq(2)){
      
      # We use this string to set file paths.
      Case = paste0(theta, gamma, delta)
      # RWD_INT <- load(paste0(MATLAB, "/Baseline_Results/RWD_Matrix_Components/Case_", 
      #                       Case, "/", "RWD_INT_", Case, ".Rda"))
      Ref_States <- data.table(read.table(paste0(RDA, "Ref_States.dat"), sep = ",", 
                                          header = TRUE, stringsAsFactors = F))
      Ref_Actions <- data.table(read.table(paste0(RDA, "Ref_Actions.dat"), sep = ",", 
                                          header = TRUE, stringsAsFactors = F))
      
      # 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, "/Baseline_Results/MATLAB_Output/Optimal_Policies/Case_", 
                          Case, "/", "OptPolicy_", Case, ".dat"), sep = ",", header = TRUE, 
                   stringsAsFactors = F))
      
      # 102418: These next several lines perform a check to ensure the optimal policy was
      # generated correctly, focusing specifically on T2DM states. The rationale for this
      # is that earlier versions produced optimal policies in which people were making 
      # choices which should in principle have fallen outside their choice set (particularly
      # in relation to decisions regarding medication adherence). The code below helped me
      # determine that health-related choices in Years 10 and after did not align with
      # patterns of health-related behaviour to which people had previously committed.
      # IF THE CODE IS WORKING PROPERLY:
      # a) Opt_Pol_Check_A + ... _Check_B + ... _Check_C + ... _Check_D = 228,216
      # b) Opt_Pol_Check_B1 = 0 AND Opt_Pol_Check_D1 = 0
      Opt_Pol_Check_A <- Opt_Pol.00[Duration_1 == 0,]
      Opt_Pol_Check_B <- Opt_Pol.00[Duration_1 < 9 & Duration_1 > 0,]
      Opt_Pol_Check_C <- Opt_Pol.00[Duration_1 == 9,]
      Opt_Pol_Check_D <- Opt_Pol.00[Duration_1 > 9,]
      
      Opt_Pol_Check_B %<>%
        .[,-c("Policy", "ID_Actions_Opt_Pol_01", "ID_States_Opt_Pol_01", "ID_Actions_RefActions",
              "ID_States_RefStates")] %>% 
        .[, c(5, 7, 6, 14, 8, 10, 12, 9, 11, 13, 1:4)]
      Opt_Pol_Check_B %<>%
        rename(Age = Age_1, HS = Health_Stat_1, Duration = Duration_1, ES = Emp_Stat_1,
               H1_LTPA = Per_1_Adherence_LTPA_1, H2_LTPA = Per_2_Adherence_LTPA_1, 
               H1_Diet = Per_1_Adherence_Diet_1, H2_Diet = Per_2_Adherence_Diet_1,
               H1_Meds = Per_1_Adherence_Medication_1, H2_Meds = Per_2_Adherence_Medication_1)
      Opt_Pol_Check_B1 <- 
        Opt_Pol_Check_B[(AD_LTPA!=H1_LTPA | AD_Diet!=H1_Diet | AD_Medication!=H1_Meds),]
      # 102318: In P1 (i.e., Years 1-8), people with T2DM always behave in a manner consistent
      # with their adherence history.
      
      Opt_Pol_Check_D %<>%
        .[,-c("Policy", "ID_Actions_Opt_Pol_01", "ID_States_Opt_Pol_01", "ID_Actions_RefActions",
              "ID_States_RefStates")] %>%
        .[, c(5, 7, 6, 14, 8, 10, 12, 9, 11, 13, 1:4)]
      Opt_Pol_Check_D %<>%
        rename(Age = Age_1, HS = Health_Stat_1, Duration = Duration_1, ES = Emp_Stat_1,
               H1_LTPA = Per_1_Adherence_LTPA_1, H2_LTPA = Per_2_Adherence_LTPA_1, 
               H1_Diet = Per_1_Adherence_Diet_1, H2_Diet = Per_2_Adherence_Diet_1,
               H1_Meds = Per_1_Adherence_Medication_1, H2_Meds = Per_2_Adherence_Medication_1)
      Opt_Pol_Check_D1 <- Opt_Pol_Check_D[(AD_LTPA!=H2_LTPA | AD_Diet!=H2_Diet | AD_Medication!=H2_Meds),]
      Opt_Pol_Check_D2 <- Opt_Pol_Check_D[(AD_LTPA==H2_LTPA & AD_Diet==H2_Diet & AD_Medication==H2_Meds),]
      # 102318: In P2 (i.e., Years), people with T2DM behave in a manner consistent with their
      # adherence history if they chose to adhere to medication in Year 9, and behave in a manner
      # INCONSISTENT with their adherence history otherwise. This means precisely half the
      # contents of the optimal policy for individuals in P2 with T2DM do not align with patterns
      # of behaviour to which they committed in Year 9.
      
      CHECK_Temp <- 
        if((nrow(Opt_Pol_Check_A) + nrow(Opt_Pol_Check_B) + nrow(Opt_Pol_Check_C) + 
            nrow(Opt_Pol_Check_D) == 228216) & 
           nrow(Opt_Pol_Check_B1) == 0 & nrow(Opt_Pol_Check_D1) == 0)
          {"OK"} else {"ERROR"}
      
      Opt_Pol.01 = Opt_Pol.00[, Action:= paste(Per_1_Adherence_LTPA_1, Per_1_Adherence_Diet_1,
                                               Per_1_Adherence_Medication_1, sep = ", ")]
      
      Age_Range <- as.character(seq(from = 25, to = 93, by = 2))
      
      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"))),]
      
      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.01 <- Opt_Pol.01[,c(10:12, 19, 13:18, 5:8, 20)]
      Opt_Pol.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")))
      Opt_Pol.01 <- data.table(Opt_Pol.01)
      
      Opt_Pol.01 <- Opt_Pol.01[, HS_1:= paste(HS_1, Dur_1, Past_AD, sep = "_")]
      Opt_Pol.02 <- Opt_Pol.01[, c("Age_1", "HS_1", "ES_1", "AD_LTPA", "AD_Diet", 
                                   "AD_Medication", "LMP")]
      
      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"))
      colnames(Opt_Pol.02) <- gsub("_NA_NA", "", colnames(Opt_Pol.02))
      
      Opt_Pol.02 %<>% mutate(Gender = GEN_K[[theta]], Ethnicity = ETH_K[[gamma]], 
                             Education = EDU_K[delta, 3])
      
      Opt_Pol.02$Age_1 <- as.character(Opt_Pol.02$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)]
      colnames(Opt_Pol.02) <- c("Gender", "Ethnicity", "Education", "Emp_Status", "Age", 
                                rep(c("LTPA", "Diet", "Meds", "LMP"), 4))
      
      # Now we process the value function.
      Val_Fun.00 = data.table(
        read.table(paste0(MATLAB, "/Baseline_Results/MATLAB_Output/Value_Functions/Case_", 
                          Case, "/", "ValFun_", Case, ".dat"), sep = ",", header = TRUE, 
                   stringsAsFactors = F))
      
      # # * Worked Examples - Part #1 ====
      # 
      # # 052619 - Worked examples
      # load(file = paste(RDA, "Worked_Examples/","WE_T0", LM, "_",
      #                   theta, gamma, delta, ".Rda", sep = ""))
      # load(file = paste(RDA, "Worked_Examples/", "WE_T1", LM, "_",
      #                   theta, gamma, delta, ".Rda", sep = ""))
      # Val_Fun_WE <- data.frame(Val_Fun.00)
      # 
      # WE_T0 <- lapply(RWD_WE_T0, 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))
      # OP_T0 <- vector("list", 3)
      # 
      # for (a in 1:3){OP_T0[[a]] <-
      #   subset(Opt_Pol.00, ID_States_Opt_Pol_01 == WE_T0[[a]]$ID_States[1])}
      # 
      # for (a in 1:3){
      #   WE_T1[[a]] <- lapply(WE_T1[[a]], 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))}
      # 
      # for (i in 1:3){
      #   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[[i]], colNames = T, keepNA = T)
      #   writeData(wb, sheet = 2, WE_T0[[i]], colNames = T, keepNA = T)
      #   for (j in 1:24){writeData(wb, sheet = j + 2, WE_T1[[i]][[j]],
      #                             colNames = T, keepNA = T)}
      #   saveWorkbook(wb, paste0(Model_Path, "/Report_Tables/Worked_Examples/",
      #                           "WE_", i, LM, "_", theta, gamma, delta, ".xlsx"),
      #                overwrite = T)}
      # 
      # # * Worked Examples - Part #2 ====
      # 
      # # 061319: Here we process data for a second set of worked examples,
      # # which, in actuality, do little more than extend Worked Example #2 above
      # # to illustrate the long-term consequences of non-adherence in an employed
      # # 60 year-old Caucasian Male (with a high school diploma or equivalent),
      # #  who has had T2DM for 9 years.
      # load(file = paste(RDA, "/Worked_Examples/","WE_T0_JW", LM, "_",
      #                   theta, gamma, delta, ".Rda", sep = ""))
      # load(file = paste(RDA, "/Worked_Examples/", "WE_T1_JW", LM, "_",
      #                   theta, gamma, delta, ".Rda", sep = ""))
      # Val_Fun_WE_JW <- data.frame(Val_Fun.00)
      # 
      # WE_T0_JW <- lapply(RWD_WE_T0_JW, FUN = function(x)
      #   merge(x, Val_Fun_WE_JW[, c("ID_States_RefStates", "Value_Fun")],
      #         by.x = "ID_States", by.y = "ID_States_RefStates", all.x = TRUE))
      # OP_T0_JW <- vector("list", 2)
      # 
      # for (a in 1:2){OP_T0_JW[[a]] <-
      #   subset(Opt_Pol.00, ID_States_Opt_Pol_01 == WE_T0_JW[[a]]$ID_States[1])}
      # 
      # for (a in 1:2){
      #   WE_T1_JW[[a]] <- lapply(WE_T1_JW[[a]], FUN = function(x)
      #     merge(x, Val_Fun_WE_JW[, c("ID_States_RefStates", "Value_Fun")],
      #           by.x = "ID_States", by.y = "ID_States_RefStates", all.x = TRUE))}
      # 
      # for (i in 1:2){
      #   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_JW[[i]], colNames = T, keepNA = T)
      #   writeData(wb, sheet = 2, WE_T0_JW[[i]], colNames = T, keepNA = T)
      #   for (j in 1:24){writeData(wb, sheet = j + 2, WE_T1_JW[[i]][[j]],
      #                             colNames = T, keepNA = T)}
      #   saveWorkbook(wb, paste0(Model_Path, "/Report_Tables/Worked_Examples/",
      #                           "WE_", i, "_JW", LM, "_", theta, gamma, delta,
      #                           ".xlsx"), overwrite = T)}
      
      # * Processing the Value Function ====
      
      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 = GEN_K[[theta]], 
        Ethnicity = ETH_K[[gamma]],
        Education = EDU_K[delta, 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 
      # statements 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] <- sapply(Val_Fun.02[6:9], 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, gamma, and theta index level of educational attainment, ethnicity, and
      # gender, respectively.
            
      OP_List[[(theta - 1)*(nrow(EDU_K)*length(ETH_K)) + (gamma - 1)*nrow(EDU_K) + 
                 delta]] <- Opt_Pol.02
      VF_List[[(theta - 1)*(nrow(EDU_K)*length(ETH_K)) + (gamma - 1)*nrow(EDU_K) + 
                 delta]] <- Val_Fun.02
      CHECK[[(theta - 1)*(nrow(EDU_K)*length(ETH_K)) + (gamma - 1)*nrow(EDU_K) + 
               delta]] <- CHECK_Temp
      
    }
  }
}

# NOTE 051219: THIS IS NEW
OP_List_H <- OP_List_D0 <- OP_List_D9_A <- OP_List_D9_DA <- OP_List
# https://datawookie.netlify.com/blog/2014/06/concatenating-a-list-of-data-frames/
OP_List_H %<>% 
  lapply(., FUN = function(x) {x <- x[, 6:9]; x}) %>%
  mapply( FUN = function( i, j ) rbind(j[, 3:6], i)
          ,. , Map_List_OP, SIMPLIFY = FALSE ) %>% do.call(cbind, .) %>% 
  cbind(rbind(Map_List_OP[[1]][c(1:4), c(1:2)], OP_List[[1]][, c(4:5)]), .)

OP_List_D0 %<>% 
  lapply(., FUN = function(x) {x <- x[, 10:13]; x}) %>%
  mapply( FUN = function( i, j ) rbind(j[, 3:6], i)
          ,., Map_List_OP, SIMPLIFY = FALSE ) %>% do.call(cbind, .) %>% 
  cbind(rbind(Map_List_OP[[1]][c(1:4), c(1:2)], OP_List[[1]][, c(4:5)]), .)

OP_List_D9_DA %<>% 
  lapply(., FUN = function(x) {x <- x[, 14:17]; x}) %>%
  mapply( FUN = function( i, j ) rbind(j[, 3:6], i)
          ,., Map_List_OP, SIMPLIFY = FALSE ) %>% do.call(cbind, .) %>% 
  cbind(rbind(Map_List_OP[[1]][c(1:4), c(1:2)], OP_List[[1]][, c(4:5)]), .)

OP_List_D9_A %<>% 
  lapply(., FUN = function(x) {x <- x[, 18:21]; x}) %>%
  mapply( FUN = function( i, j ) rbind(j[, 3:6], i)
          ,., Map_List_OP, SIMPLIFY = FALSE ) %>% do.call(cbind, .) %>% 
  cbind(rbind(Map_List_OP[[1]][c(1:4), c(1:2)], OP_List[[1]][, c(4:5)]), .)

VF_List_H <- VF_List_D0 <- VF_List_D9_A <- VF_List_D9_DA <- VF_List

VF_List_H %<>% 
  lapply(., FUN = function(x) {x <- x[, c(4:5, 6)]; x}) %>%
  lapply(., FUN = function(x) {x <- rename(x, Value = "Healthy"); x}) %>%
  mapply( FUN = function( i, j ) rbind(j, i)
          ,. , Map_List_VF, SIMPLIFY = FALSE ) %>% do.call(cbind, .)
VF_List_H <- cbind(VF_List_H[, c(1:2)], 
                   VF_List_H[, grepl("Value", colnames(VF_List_H))])

VF_List_D0 %<>% 
  lapply(., FUN = function(x) {x <- x[, c(4:5, 7)]; x}) %>%
  lapply(., FUN = function(x) {x <- rename(x, Value = "T2DM_0_NA"); x}) %>%
  mapply( FUN = function( i, j ) rbind(j, i)
          ,. , Map_List_VF, SIMPLIFY = FALSE ) %>% do.call(cbind, .)
VF_List_D0 <- cbind(VF_List_D0[, c(1:2)], 
                    VF_List_D0[, grepl("Value", colnames(VF_List_D0))])

VF_List_D9_DA %<>% 
  lapply(., FUN = function(x) {x <- x[, c(4:5, 8)]; x}) %>%
  lapply(., FUN = function(x) {x <- rename(x, Value = "T2DM_9_No"); x}) %>%
  mapply( FUN = function( i, j ) rbind(j, i)
          ,. , Map_List_VF, SIMPLIFY = FALSE ) %>% do.call(cbind, .)
VF_List_D9_DA <- cbind(VF_List_D9_DA[, c(1:2)], 
                       VF_List_D9_DA[, grepl("Value", colnames(VF_List_D9_DA))])

VF_List_D9_A %<>% 
  lapply(., FUN = function(x) {x <- x[, c(4:5, 9)]; x}) %>%
  lapply(., FUN = function(x) {x <- rename(x, Value = "T2DM_9_Yes"); x}) %>%
  mapply( FUN = function( i, j ) rbind(j, i)
          ,. , Map_List_VF, SIMPLIFY = FALSE ) %>% do.call(cbind, .)
VF_List_D9_A <- cbind(VF_List_D9_A[, c(1:2)], 
                      VF_List_D9_A[, grepl("Value", colnames(VF_List_D9_A))])

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

# Create an Excel workbook object and add a worksheet
wb_OP1 = createWorkbook()
for (i in seq(OP_List))
{sht_OP1 = addWorksheet(wb_OP1, paste0(Map[i, 3], ", ", Map[i, 2], ", ", Map[i, 1]))
  dat_OP1 = OP_List[[i]]
  writeData(wb_OP1, sht_OP1, dat_OP1, colNames = T, withFilter = TRUE)
  addStyle(wb_OP1, sht_OP1, style = createStyle(fontSize = 9), rows = 1:71, cols = 1:21, gridExpand = T)}
saveWorkbook(wb_OP1, paste0(Model_Path, "/Report_Tables/Baseline_Results/",
                           "OptPol1_", LM, ".xlsx"), overwrite = T)

wb_OP2 = createWorkbook()
sht_OP2 = list(addWorksheet(wb_OP2, "Healthy"), addWorksheet(wb_OP2, "T2DM_D0"), 
           addWorksheet(wb_OP2, "T2DM_D9_DA"), addWorksheet(wb_OP2, "T2DM_D9_A"))
dat_OP2 = list(OP_List_H, OP_List_D0, OP_List_D9_DA, OP_List_D9_A)
for (i in 1:4){writeData(wb_OP2, sht_OP2[[i]], dat_OP2[[i]], colNames = T, keepNA = T)}
saveWorkbook(wb_OP2, paste0(Model_Path, "/Report_Tables/Baseline_Results/",
                        "OptPol2_", LM, ".xlsx"), overwrite = T)

wb_VF1 = createWorkbook()
for (i in seq(VF_List))
  {sht_VF1 = addWorksheet(wb_VF1, paste0(Map[i, 3], ", ", Map[i, 2], ", ", Map[i, 1]))
  dat_VF1 = VF_List[[i]]
  writeData(wb_VF1, sht_VF1, dat_VF1, colNames = T, withFilter = TRUE)
  addStyle(wb_VF1, sht_VF1, createStyle(fontSize = 9), rows = 1:71, cols = 1:9, gridExpand = T)
  addStyle(wb_VF1, sht_VF1, createStyle(numFmt = "@"), rows = 2:71, cols = 5, gridExpand = T, stack = T)
  addStyle(wb_VF1, sht_VF1, createStyle(numFmt = "#,##0.00"), rows = 2:71, cols = 6:9, gridExpand = T, stack = T)}
saveWorkbook(wb_VF1, paste0(Model_Path, "/Report_Tables/Baseline_Results/",
                           "ValFun1_", LM, ".xlsx"), overwrite = T)

wb_VF2 = createWorkbook()
sht_VF2 = list(addWorksheet(wb_VF2, "Healthy"), addWorksheet(wb_VF2, "T2DM_D0"), 
           addWorksheet(wb_VF2, "T2DM_D9_DA"), addWorksheet(wb_VF2, "T2DM_D9_A"))
dat_VF2 = list(VF_List_H, VF_List_D0, VF_List_D9_DA, VF_List_D9_A)
for (i in 1:4){writeData(wb_VF2, sht_VF2[[i]], dat_VF2[[i]], colNames = T, keepNA = T)}
saveWorkbook(wb_VF2, paste0(Model_Path, "/Report_Tables/Baseline_Results/",
                            "ValFun2_", LM, ".xlsx"), overwrite = T)