# Control Module for the Patient Decision Model - Robustness Tests
# Created: 09/22/18; Last Modified: 06/29/19

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

# NOTE 042019: 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 passes parameters to RWD&TP_Matrices_Construction_090618.R,
# before executing said script. This is done primarily to automate the
# generation of TP and RWD matrix components for individuals with differing
# personal characteristics. It also serves as a template for the script
# operationalizing the sensitivity analysis.

# 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(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/")

# Load RDA files derived from other scripts, beginning with the UKPDS-OM2 simulation
# results
filenames <- c(paste0(RDA, "UKPDS_OM2_Results_Processing", LM, ".Rda"), 
               # Load life tables for probabilities relating to risk of mortality.
               paste0(RDA, "Life_Tables", LM, ".Rda"), 
               # Load vector of life expectancies (which are derived from the life tables).
               paste0(RDA, "Life_Exp", LM, ".Rda"),
               # Load health state transition probabilities derived from application of the DPoRT
               # risk equations.
               paste0(RDA, "DPoRT_Calculations", LM, ".Rda"),
               # Loading our calculations from the analysis of the 2011 SLID data.
               paste0(RDA, "SLID_2011_Calculations", LM, ".Rda"),
               # Load estimates of gross individual employment income derived from analysis of 2011
               # NHS data
               paste0(RDA, "NHS_2011_Calculations", LM, "_1", ".Rda"),
               # Load tax rates used in analysis of 2011 NHS data
               paste0(RDA, "NHS_2011_Calculations", LM, "_2", ".Rda"),
               # Load presenteeism calculations
               paste0(RDA, "Presenteeism_Calculations", LM, ".Rda"))

# https://stackoverflow.com/questions/3764292/loading-many-files-at-once
lapply(filenames, load, environment())

# * Setting up for robustness tests ====

# NOTE 092018: This next line is VERY important, because it determines the characteristics of
# the decision-maker for the purposes of the decision analysis, AS WELL AS the the parameter
# being varied.
Person_Type <- c(2, 3, 4); Robustness_Tests_ID <- 2
# These integer values index positions in the three lists constructed above. BY
# MODIFYING THESE THREE VALUES, IT IS POSSIBLE TO GENERATE TRANSITION PROBABILITY AND
# REWARD MATRICES FOR ALL COMBINATIONS OF FIXED PERSONAL CHARACTERISTICS.
Gender.U = Person_Type[1]; Ethnicity.U = Person_Type[2]; Education.U = Person_Type[3]

# 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")

# * Parameters relating to retirement income and contributions ====

# For these calculations, we use the average of maximum pensionable earnings for the
# last 5 years ($53,480; 2013 - $51,100, 2014 - $52,500, 2015 - $53,600, 
# 2016 - $54,900, 2017 - $55,300):
# http://www.drpensions.ca/cpp-rate-table.html
Max_Pen_Earn_K <- c(53480)
# We also generate a separate variable for 2017 MPE.
Max_Pen_Earn_2017_K <- c(55300)
# The sources for the remainder of these parameters are cited in the thesis. The
# OAS amounts are effective October 2017.
OAS_Amt_K <- 585.49; OAS_Threshold_K <- 72809.00
GIS_Amt_K <- 874.48; GIS_Threshold_K <- 3500.00
CPP_Contrib_Rate_K <- 0.0495; Max_CPP_Contrib_K <- 2564.10; Basic_Exemption_K <- 3500.00
EI_Contrib_Rate_K <- 0.0163; Max_EI_Contrib_K <- 836.19
# These parameters relate to the individual's private pension and represent,
# respectively, the rate at which savings accumulate, the proportion of income (either
# pre-tax or post-tax) saved in any given year, and the year in which the individual
# begins withdrawing income from their private retirement fund.
Rate_K <- 0.03; Savings_Per_K <- 0.10; Year_Ret_K <- 65

# * Parameters relating to household expenditures ====

# Annual pharmaceutical expenditures, based on results of the pharmacist consultation.
# https://stackoverflow.com/questions/7070173/r-friendly-way-to-convert-r-data-frame-column-to-a-vector
RxExp.Early.Annual_K <- read_xlsx(path = paste0(Model_Path, "/Raw_Data/",
                                                "Pharmaceutical_Costs_V2.xlsx"), 
                                  sheet = "ND", range = "I7", col_names = FALSE)[[1]]
RxExp.Adv.WC.Annual_K <- read_xlsx(path = paste0(Model_Path, "/Raw_Data/",
                                                 "Pharmaceutical_Costs_V2.xlsx"), 
                                   sheet = "WM", range = "I12", col_names = FALSE)[[1]] 
RxExp.Adv.PC.Annual_K <- read_xlsx(path = paste0(Model_Path, "/Raw_Data/",
                                                 "Pharmaceutical_Costs_V2.xlsx"), 
                                   sheet = "PM", range = "I15", col_names = FALSE)[[1]] 
Min_Ded_K <- 100

# Now we load the Pharmacare rates for the 2018/19 benefit year, which were collected
# from the province's website: http://www.gov.mb.ca/health/pharmacare/estimator.html

url <- 'https://www.gov.mb.ca/health/pharmacare/estimator.html'
webpage <- read_html(url)

Pharm.00 <- html_nodes(webpage, css = "p+ div tr+ tr td")
Pharm.01 <- as.matrix(html_text(Pharm.00))

Pharm.01_Income <- Pharm.01[1]
Pharm.01_Income <- unlist(strsplit(Pharm.01_Income, split = "\r\n"))
# https://stackoverflow.com/questions/5992082/how-to-remove-all-whitespace-from-a-string
# https://www.rforexcelusers.com/remove-currency-dollar-sign-r/
# https://stackoverflow.com/questions/33949945/replace-multiple-strings-in-one-gsub-or-chartr-statement-in-r
# https://stackoverflow.com/questions/11129107/removing-non-alpanumeric-characters-from-an-ordered-collection-of-objects-list

# Remove all non-numeric characters, excepting dashes
Pharm.01_Income %<>%
  gsub(pattern = "[^0-9-]", replacement = "", .) %>%
  data.frame(., stringsAsFactors = FALSE) %>%
  rename("Income" = ".") %>%
  separate(., col = "Income", into = c("Lower.Limit", "Upper.Limit"), sep = "-")

# https://stackoverflow.com/questions/15629885/replace-na-in-column-with-value-in-adjacent-column
Pharm.01_Income[is.na(Pharm.01_Income)] <- 1000000

Pharm.01_DedPer <- Pharm.01[2]
Pharm.01_DedPer <- unlist(strsplit(Pharm.01_DedPer, split = "%\r\n"))
Pharm.01_DedPer %<>%
  gsub(pattern = "[^0-9.]", replacement = " ", .) %>%
  data.frame(., stringsAsFactors = FALSE) %>%
  rename("Rate" = ".")

Pharmacare_K <- cbind(Pharm.01_Income, Pharm.01_DedPer)
# https://stackoverflow.com/questions/18503177/r-apply-function-on-specific-dataframe-columns
Pharmacare_K[1:3] <- sapply(Pharmacare_K[1:3], as.numeric)
Pharmacare_K$Rate <- Pharmacare_K$Rate/100

# Dietary expenditures
Age_Lookup_K <- data.frame("Age" = seq(from = 25, to = 95))
Age_Lookup_K <- Age_Lookup_K %>% 
  mutate(Age_Group = case_when(Age < 45 ~ "Age_19_44",
                               Age < 65 ~ "Age_45_64",
                               Age >= 65 ~ "Age_65_95"))

Diet_Cost <- data.table(read_xlsx(path = paste0(Model_Path, "/Raw_Data/",
                                                "Meal_Plan_MS_v9.xlsm"), 
                                  sheet = "EER", range = "A21:H33", col_names = TRUE))
Diet_Cost <- Diet_Cost %>% rename(Diet_Type = "Diet Type")
# https://stackoverflow.com/questions/9202413/how-do-you-delete-a-column-by-name-in-data-table
# https://stackoverflow.com/questions/16638484/remove-multiple-columns-from-data-table
Diet_Cost.H <- Diet_Cost[grep("AHEI Score: 102.5", Diet_Type),]
Diet_Cost.UH <- Diet_Cost[grep("AHEI Score: 2.5", Diet_Type),]
Diet_Cost.H <- Diet_Cost.H[, c("Diet_Type") := NULL]
Diet_Cost.UH <- Diet_Cost.UH[, c("Diet_Type") := NULL]

Diet_Cost.H <- melt(Diet_Cost.H, id.vars = c("Age_Group"), variable.name = c("Chars"),
                    value.name = "Cost")
Diet_Cost.H <- separate(data = Diet_Cost.H, col = Chars, into = c("Ethnicity", "Gender")
                        , sep = "_")
Diet_Cost.UH <- melt(Diet_Cost.UH, id.vars = c("Age_Group"), variable.name = c("Chars"),
                     value.name = "Cost")
Diet_Cost.UH <- separate(data = Diet_Cost.UH, col = Chars, into = c("Ethnicity", "Gender")
                         , sep = "_")

Diet_Cost.H.01_K <- Diet_Cost.UH.01_K <- 
  data.frame(expand.grid("Age" = seq(from = 25, to = 95), "Gender" = c("Male", "Female"), 
                         Ethnicity = c("Caucasian", "Afro-Caribbean", "Asian-Indian"), 
                         stringsAsFactors = F))
Diet_Cost.H.01_K <- Diet_Cost.H.01_K %>%
  mutate(Age_Group = Age_Lookup_K[findInterval(Age,Age_Lookup_K$Age),2])
Diet_Cost.H.01_K <- merge(Diet_Cost.H.01_K, Diet_Cost.H, 
                          by = c("Age_Group", "Gender", "Ethnicity"), all.x = T)
Diet_Cost.UH.01_K <- Diet_Cost.UH.01_K %>%
  mutate(Age_Group = Age_Lookup_K[findInterval(Age,Age_Lookup_K$Age),2])
Diet_Cost.UH.01_K <- merge(Diet_Cost.UH.01_K, Diet_Cost.UH, 
                           by = c("Age_Group", "Gender", "Ethnicity"), all.x = T)

# LTPA expenditures
LTPA_Cost_List_K <- c(0,0)

# * Parameters relating to time allocation ====

# Time allocated to sleep each year
Time_Req_Sleep_K <- 8*365
# Time available for any other activity
Time_Available_K <- 24*365 - Time_Req_Sleep_K
# Time allocated to LTPA (for those adherent to it)
Time_Req_LTPA_K <- (150/60)*52
# Time allocated to food preparation
Time_Req_Diet_H_K <- 90*365/60; Time_Req_Diet_UH_K <- 30*365/60
# Time allocated to job search (for individuals younger than 65 who are
# currently unemployed)
Time_Req_Search_K <- 9.5*52

# * Other miscellaneous parameters ====

# 052519: Penalty for non-admissible actions
Penalty_K <- -1000000
# 042019: Price of consumption, the multiplier applied to the wage, and Rho, which is 
# a key parameter of the CES utility function, and which relates to the elasticity of
# substitution in the following way: sigma = (1/(1-Rho)).
P_c_K <- 1; Mult_w_K <- 1; Rho_K <- 0;
# Impact of healthy, prediabetic, and T2DM states on probability of sustaining and/or
# recovering employment
Emp_Diff_Healthy_K = 0; Emp_Diff_Pred_K = 0; Emp_Diff_T2DM_K = 0
# Starting age for the decision problems
Start_Age_K = 24
# Monthly EIA for individuals with no dependents in the General Assistance category: 
# https://www.gov.mb.ca/fs/eia/eia_general.html
EIA_General_K = 12*806
# Monthly EIA for individuals with no dependents who is considered to have a disability:
# https://www.gov.mb.ca/fs/eia/EIA_Disability_K_K.html
EIA_Disability_K = 12*1027
Annual_Work_Hrs.FT_K <- 2000; Annual_Work_Hrs.PT_K <- 1000
Pres_Share_K <- 0.5
Paid_SDs.FT_K <- 7; Paid_SDs.PT_K <- Paid_SDs.FT_K/2
# "ARPD" = "age-related productivity decline"
ARPD_K <- 0.025
# Parameters for the utility function
Alpha_K <- 0.5

# * Personal characteristics ====
GEN_K = list("Male", "Female")
ETH_K = list("Caucasian", "Afro-Caribbean", "Asian-Indian")

EDU_K <- 
  data.frame(Index = as.character(seq(1:5)), 
             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",
                        "5. Postsecondary certificate, diploma or degree"), 
             Pres_Code = c("Educ_1", "Educ_2", "Educ_3", "Educ_4", "Educ_5"),
             # 081718: The next line effectively defines PSE for the purposes of our
             # calculations. The paper in which the DPoRT 2.0 Risk Equations were published
             # included a variable that denoted any PS attendance. Reviewing the data dictionary
             # and derived variable guide for the CCHS 2012 confirms the authors' interpretation
             # was correct.
             # 
             # However, the NHS data indicates only the individual's level of educational
             # attainment (i.e., whether they actually progressed far enough in their education 
             # to receive some kind of credential). We bridge the gap by assuming a person must
             # have obtained a PS certificate or diploma below bachelor level to satisfy the 
             # definition. This is conservative in that it potentially excludes many individuals
             # who attended PS educational institutions but did not progress far enough to 
             # receive a credential. However, it avoids the risk of including a potentially 
             # large segment of individuals who did not attend PSE after completing high school 
             # (NOTE: the next line of code been revised - an earlier version of the code 
             # included everyone with at least a HS education in the definition of "some PSE").
             Attended_PS = c(rep("NoPS", 2), rep("PS", 3)))

Age_Align_K <- data.frame("Age_Cat" = c("25_34", "35_44", "45_54", "55_64", "65_"),
                          "Age_Group" = c("Age_19_44", "Age_19_44", "Age_45_64", "Age_45_64",
                                          "Age_65_95"), stringsAsFactors = FALSE)

rm(Diet_Cost, Diet_Cost.H, Diet_Cost.UH, Pharm.00, Pharm.01, Pharm.01_DedPer, Pharm.01_Income,
   webpage, url)

# Generation of Reward and Transition Probability Matrices ----------------

# Paths for components of the transition probability matrices (i.e., for the output generated 
# by this script).
path_TP <- paste0(MATLAB, "Robustness_Tests/", 
                   Robustness_Tests_Types[Robustness_Tests_ID],
                   "/TP_Matrix_Components/")
# Paths for components of the reward matrices (i.e., for the output generated by this script).
path_RWD <- paste0(MATLAB, "Robustness_Tests/", 
                   Robustness_Tests_Types[Robustness_Tests_ID],
                   "/RWD_Matrix_Components/")

# This next variable is included specifically to faciliate integration of 2011 SLID
# data. Anyone who has at minimum a high school diploma or equivalent is assigned a
# value of 1 (note that this does not align precisely with the 2011 SLID codebook,
# since many individuals completing secondary school will not go on to participate in
# ANY post-secondary education).
Education_SLID = as.character(ifelse(Education.U > 1, 1, 0))

IsTopQuin <- Emp_Income[Emp_Income$Sex == GEN_K[[Gender.U]] & Emp_Income$Highest_Degree ==
                          EDU_K[Education.U, 2],c(1:3, 5)]
IsTopQuin <- merge(IsTopQuin, Age_Align_K, by = "Age_Cat", all.x = TRUE)
IsTopQuin$IsTopQuin <- ifelse(IsTopQuin$IsTopQuin==FALSE, "NoTQ", "TQ")

# NOTE 092318: Now we call the appropriate script to initiate the generation of the reward
# and transition probability matrices (note that at present, this is slighly different
# than the script used for the baseline analysis - for reasons we discuss below.
source(paste0(Model_Path, "/RWD&TP_Matrices_Construction", LM, "_RT", ".R"))

# NOTE 092318: For each model run, we save the reward and transition probability matrix
# components in the appropriate place.
for (i in seq(TP.03))
{fwrite(x = TP.03[[i]], sep = ",", col.names = FALSE, file = 
          paste0(path_TP, "TP_", i, "_",
                 Gender.U, Ethnicity.U, Education.U, ".csv"))}
for (i in seq(RWD.03))
{fwrite(x = RWD.03[[i]], sep = ",", col.names = FALSE, file = 
          paste0(path_RWD, "RWD_", i, "_",
                 Gender.U, Ethnicity.U, Education.U, ".csv"))}