# Validation of health state transition probabilities - Healthy to T2DM

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

# 042619: I want to verify that a subset of health state transition
# probabilities used in the decision model are used correctly. In
# particular, I want to look at the likelihood of transitioning from being
# healthy to being diagnosed with T2DM.
# - We already have the TP matrices - we'll have generated them in setting 
# up the optimization problem. We just need to specify where they are
# located. This, in turn, varies according to person type, which we will
# thus also need to specify. To simply things, we can filter the TP
# matrices to focus on a subset of state transitions, which should reduce
# the size of the resultant data structures significantly.
# - We also have the values we wish to check against, which derive from the
# DPoRT calculations.

# 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/")
# Paths for components of the transition probability and reward matrices
path_TP <- paste0(MATLAB, "Baseline_Results/TP_Matrix_Components/Case_")
path_RWD <- paste0(MATLAB, "Baseline_Results/RWD_Matrix_Components/Case_")

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

# We are interested in females (gender is given by theta) of Asian-Indian
# extraction (ethnicity is given by gamma) with a university certificate,
# diploma or degree at bachelor level or above (education is given by delta).
delta = 4; gamma = 3; theta = 2
Gender.U = theta; Ethnicity.U = gamma; Education.U = delta
Case <- paste0(theta, gamma, delta);

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)

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

# DPoRT Data --------------------------------------------------------------

Pr_Risk.01 <- Pr_Risk[,c(1:4, 25:40)]
Pr_Risk.01 <- Pr_Risk.01 %>%
  mutate("Chars" = paste0(Ethnicity, "_", Gender, "_", Caloric_Intake))
Pr_Risk.01 <- Pr_Risk.01[,c(2,21,5:20)]
Pr_Risk.02 <- melt(Pr_Risk.01, id = c(1:2), measure = c(3:ncol(Pr_Risk.01)))
Pr_Risk.02 <- dcast(Pr_Risk.02, Age_Group + variable ~ Chars, value.var = "value")
Pr_Risk.02 <- Pr_Risk.02[,c(1:2, 14:11, 6:3, 10:7)]
Pr_Risk.02$variable %<>%
  gsub("Low_Active", "Low.Active", .) %>%
  gsub("Very_Active", "Very.Active", .)
Pr_Risk.02 <- separate(Pr_Risk.02, variable, into = c("LTPA", "PS", "TQ", "Risk"), sep = "_")
Pr_Risk.02$LTPA %<>%
  gsub("Low.Active", "Low_Active", .) %>%
  gsub("Very.Active", "Very_Active", .)

# http://www.cookbook-r.com/Manipulating_data/Changing_the_order_of_levels_of_a_factor/
Pr_Risk.02$LTPA <- factor(Pr_Risk.02$LTPA, 
                          levels = c("Sedentary", "Low_Active", "Active", "Very_Active"))
Pr_Risk.02 <- arrange(Pr_Risk.02, TQ, PS, Age_Group, LTPA)
Pr_Risk.02 <- Pr_Risk.02 %>%
  mutate("variable" = paste(LTPA, PS, TQ, Risk, sep = "_"))
Pr_Risk.02 <- Pr_Risk.02[,c(1,18,6:17)]

Healthy_T2DM.01 <- 
  subset(Pr_Risk.02, select = c("Age_Group","variable",
                                paste(ETH_K[[Ethnicity.U]], GEN_K[[Gender.U]], "Recommended", sep = "_"),
                                paste(ETH_K[[Ethnicity.U]], GEN_K[[Gender.U]], "Excessive", sep = "_")))

Healthy_T2DM.01 <- Healthy_T2DM.01 %>%
  filter((Age_Group == "Age_19_44" & variable == 
            paste("Active", EDU_K[Education.U, 4], 
                  IsTopQuin[match("Age_19_44",IsTopQuin$Age_Group),4], "Risk", sep = "_")) | 
           (Age_Group == "Age_19_44" & variable == 
              paste("Sedentary", EDU_K[Education.U, 4], 
                    IsTopQuin[match("Age_19_44",IsTopQuin$Age_Group),4], "Risk", sep = "_")) |
           (Age_Group == "Age_45_64" & variable == 
              paste("Active", EDU_K[Education.U, 4], 
                    IsTopQuin[match("Age_45_64",IsTopQuin$Age_Group),4], "Risk", sep = "_")) | 
           (Age_Group == "Age_45_64" & variable == 
              paste("Sedentary", EDU_K[Education.U, 4], 
                    IsTopQuin[match("Age_45_64",IsTopQuin$Age_Group),4], "Risk", sep = "_")) | 
           (Age_Group == "Age_65_95" & variable == 
              paste("Active", EDU_K[Education.U, 4], 
                    IsTopQuin[match("Age_65_95",IsTopQuin$Age_Group),4], "Risk", sep = "_")) | 
           (Age_Group == "Age_65_95" & variable == 
              paste("Sedentary", EDU_K[Education.U, 4], 
                    IsTopQuin[match("Age_65_95",IsTopQuin$Age_Group),4], "Risk", sep = "_")))
Healthy_T2DM.01$Age_Group <- as.character(Healthy_T2DM.01$Age_Group)

for (i in 1:length(Diet)){
  for (j in 1:length(Activity)){
    Healthy_T2DM.02 <-
      # NOTE 061018: Here we combine several programming techniques to significantly
      # condense what would have otherwise been a rather lengthy statement. In essense,
      # we'd like to perform a sequence of left joins on a character field that reflects
      # three distinct individual characteristics, namely: whether they are active or
      # sedentary; whether they attended post-secondary education; and whether they are in
      # the top income quintile. This would be simple enough in itself. However, assembling
      # the strings using the PASTE function is problematic, because the tendency for 
      # employment income to change over an individual's career means they may or may not
      # be in the top quintile at any given point in time, which, in turn, means that the
      # string we are looking for varies even for individuals who are otherwise the same.
      # For example, an active male in the 45-64 age range who engages in moderate caloric
    # consumption would be identified by the string "Active_PS_TQ_Risk", while the same
    # male in the 65+ age range would be identified by the string "Active_PS_NoTQ_Risk".
    # 
    # One way to deal with this is to use the GREP function in combination with a wild-
    # card character (where ".*" denotes a string of arbitrary length), but we would
    # also like to be able to use standard subsetting techniques to account for personal
    # characteristics. To get the best of both worlds, we nest a PASTE function inside of
    # a GREP function, expressing the wildcard as just another substring. Once the string
    # is concatenated, however, R seems to understand what we were intending, and
    # identifies the appropriate rows (see:
    # https://stackoverflow.com/questions/26813667/how-to-use-grep-to-find-exact-match).
    # 
    # Going one step further, we insert elements of the substrings relating to LTPA
    # participation and diet into matrices, and then employ a for-loop to obtain all
    # possible combinations of the matrix elements, reducing the amount of coding
    # required by nearly three-quarters (I say "nearly", because there is an up-front
    # investment necessary to set up the for-loop).
    left_join(Healthy_T2DM.02, Healthy_T2DM.01[grep(paste(
      Activity[j], EDU_K[Education.U, 4],".*", "Risk", sep = "_"), Healthy_T2DM.01$variable),
      c("Age_Group", paste(ETH_K[[Ethnicity.U]], GEN_K[[Gender.U]], Diet[i], sep = "_"))],
      by = "Age_Group")
  }
}

colnames(Healthy_T2DM.02)[3:6] <- c("LTPA_A>Diet_A", "LTPA_NA>Diet_A", "LTPA_A>Diet_NA",
                                    "LTPA_NA>Diet_NA")

Healthy_T2DM.02 <- Healthy_T2DM.02 %>% 
  melt(., id.vars = c("Age_1", "Age_Group"),
       variable.names = c("LTPA_A.Diet_A","LTPA_NA.Diet_A","LTPA_A.Diet_NA",
                          "LTPA_NA.Diet_NA"),
       value.name = "p_healthy_T2DM") %>%
  mutate("Health_Stat_1" = "Healthy") %>%
  mutate("Health_Stat_2" = "T2DM") %>%
  mutate("p_healthy_healthy" = 1 - p_healthy_T2DM)

Healthy_T2DM.02$variable <- as.character(Healthy_T2DM.02$variable)

Healthy_T2DM.02 <- separate(Healthy_T2DM.02, col = "variable", 
                            into = c("LTPA", "Diet"), sep = ">", remove = TRUE)
Healthy_T2DM.02$p_healthy_T2DM[Healthy_T2DM.02$Age_1 == 25] <- 0

# SLID Data ---------------------------------------------------------------

Emp.01 <- SLID_2011_S6 %>%
  # Interestingly, despite the fact these variables are categorized as integers, it is
  # still necessary to apply quotation marks in referring to them for the filter to
  # operate correctly.
  filter(ecsex99==as.character(Gender.U), educ_at==Education_SLID) %>%
  subset(., select = c(age_cat, ml01v28, EmpPr_01.01, UnempPr_01.01))

Emp.01$ml01v28 <- as.character(Emp.01$ml01v28)
Emp.01$ml01v28[Emp.01$ml01v28=="10"] <- "E"
Emp.01$ml01v28[Emp.01$ml01v28=="20"] <- "U"
Emp.01 <- Emp.01 %>% rename("Age_Group" = age_cat, "Start_Emp" = ml01v28)
Emp.02 <- melt(Emp.01, id.vars = c("Age_Group", "Start_Emp"), 
               variable.name = "End_Emp", value.name = "pEmpStatChange")
Emp.02$End_Emp <- as.character(Emp.02$End_Emp)
Emp.02$End_Emp[Emp.02$End_Emp=="EmpPr_01.01"] <- "E"
Emp.02$End_Emp[Emp.02$End_Emp=="UnempPr_01.01"] <- "U"

# Life Table Data ---------------------------------------------------------

TP.02 <- 
  TP.01[TP.01$Age %in% c(25:94) & TP.01$Gender == GEN_K[theta] & 
          TP.01$Ethnicity == ETH_K[gamma] & TP.01$Employment == "Employed",]
TP.02 %<>% rename(., "pMortality" = "Probability")

# TP Matrices -------------------------------------------------------------

# Download Ref_States and Ref_Actions
Path_Ref <- "C:/Users/starg/Documents/Dissertation/Diabetes_Model_042019/R_Objects/"
Ref_States <- read.table(paste0(Path_Ref, "Ref_States.dat"), header = T, sep = ",")
Ref_Actions <- read.table(paste0(Path_Ref, "Ref_Actions.dat"), header = T, sep = ",")

Path_TPs <- paste0("C:/Users/starg/Documents/Dissertation/Diabetes_Model_042019",
                   "/MATLAB/Baseline_Results/TP_Matrix_Components/")
Ref_States.01 <- select(Ref_States, c("ID_States", "Health_Stat_1"))

# Need to load all of the pertinent files into a list.
TP_list <- vector("list", 24)
for (i in 1:24){
  TP_list[[i]] <- 
    read_xlsx(paste0(Path_TPs, "Case_", Case, "/TP_", i, "_", Case, ".xlsx"), 
              col_names = c("Start", "End", "Prob"), col_types = c("numeric"))
}

TP_list.01 <- TP_list;
TP_list.01 <-
  lapply(TP_list, FUN = function(x)
    x <- x %>% 
      left_join(., Ref_States[, c("ID_States", "Age_1", "Health_Stat_1", "Emp_Stat_1")], 
                by = c("Start" = "ID_States")) %>% 
      filter(., Health_Stat_1 == "Healthy") %>%
      left_join(., Ref_States[, c("ID_States","Age_1", "Health_Stat_1", "Emp_Stat_1")], 
                by = c("End" = "ID_States")) %>%
      arrange(., Start, End))
names(TP_list.01) <- c(1:24)

for (i in 1:24){TP_list.01[[i]] %<>% mutate(., "ID_Actions" = i)}

TP_list.01 %<>%
  lapply(., FUN = function(x)
    x <- x %>% 
      left_join(., Ref_Actions, by = "ID_Actions") %>%
      arrange(., Start, End))

TP_list.02 <- TP_list.01 %>%
  lapply(., FUN = function(x)
    {x[c(5:6, 8:9, 11:14)] <- sapply(x[c(5:6, 8:9, 11:14)], as.character); x}) 

TP_list.02 %<>% 
  lapply(., FUN = function(x) 
    x <- x %>% left_join(., TP.02[, c("Age", "pMortality")], 
                         by = c("Age_1.x" = "Age")) %>%
      mutate(Age_Group = ifelse(Age_1.x < 55, 1, 2)) %>%
      left_join(., Emp.02, by = c("Age_Group", 
                                        "Emp_Stat_1.x" = "Start_Emp",
                                        "Emp_Stat_1.y" = "End_Emp")))

TP_list.03 <- do.call("rbind", TP_list.02)

wb = createWorkbook()
for (i in seq(TP_list.02))
{sht = addWorksheet(wb, paste0("Action", names(TP_list.02[i])))
dat = TP_list.02[[i]]
writeData(wb, sht, dat, colNames = T)}
saveWorkbook(wb, paste0("C:/Users/starg/Desktop/", "TP_042619", ".xlsx"), 
             overwrite = T)

