# Utility for Construction of Reward and Transition Probability Matrices
# Created: 04/09/18; Last Modified: 04/20/19

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

# This script generates the transition probability and reward matrices that lie at
# heart of the decision model.

# IMPORTANT: THIS CODE RELIES HEAVILY ON A NUMBER OF RDA FILES IMPORTED FROM OTHER 
# R SCRIPTS. THE SCRIPTS USED TO GENERATE THOSE FILES MUST THEREFORE BE IMPLEMENTED
# FIRST. REFER TO THE README FILE FOR MORE INFORMATION.

# NOTE 062419: Replaced use of the MERGE function (part of the BASE package)
# with the LEFT_JOIN function (part of the DPLYR package). This accelerates
# merges significantly, and probably shaves ~15 minutes off execution time
# for each loop (a total of 6 hours across all loops). In addition, fixed an 
# error that referred to a non-existent variable "Health_Stat", rather than 
# "Health_Stat_1", which may have affected calculation of utility for
# deceased individuals.

# NOTE 102118: Fixed an error relating to calculation of absenteeism after having
# had T2DM for 10 years.

# NOTE 100418: Updated script to generate and store efficiency loss calculations,
# which can then be used to carry out simulations examining the impact of a range
# of health behaviors and interventions on indirect costs associated with T2DM.

# NOTE 090618: Updated script to incorporate most recent Pharmacare deductible
# rates from the Manitoba Health, Seniors and Active Living website:
# http://www.gov.mb.ca/health/pharmacare/estimator.html

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

# Characterize the choice and state spaces ----

# * Process UKPDS simulation data ====

Output6_TEST <- Output6 %>%
  filter(Gender == GEN_K[[Gender.U]], Ethnicity == ETH_K[[Ethnicity.U]]) 

# %>%
# rename("Rate" = "Rate.F", "EF_Survival" = "EF_Survival.F", 
#        "Probability" = "Probability.F", "Expected_Utility" = "Expected_Utility.F",
#        "Annual_Absences" = "Annual_Absences.F", 
#        "Adjusted_Utility" = "Adjusted_Utility.F")

Output6_TEST$Age <- as.integer(Output6_TEST$Age)

Output6_TEST <- Output6_TEST %>% rename("Age_at_Diagnosis" = "Age") 
Output6_TEST$Age_at_Diagnosis <- as.integer(Output6_TEST$Age_at_Diagnosis)
Output6_TEST$Duration <- as.integer(Output6_TEST$Duration)

# We note that "age" as it is captured in the UKPDS results is believed to reflect (this
# must be verified) the age at which the individual was DIAGNOSED, implying that their
# actual age is the sum of age at diagnosis and duration of T2DM (upon first glance, the
# data would suggest that there are 27 year-old individuals who have had T2DM for 50 or more
# years, which obviously doesn't make much logical sense).
# 
# # A lot of the UKPDS simulation results involve individuals who are at advanced ages when
# diagnosed. Given the limits to the human lifespan, while the UKPDS might obtain results
# on individuals who were diagnosed at age 90 and have had T2DM for 60 years, such results
# are of questionable value and can be discarded. We note that when we do so, the number of
# records in the data drops by more than half.
Output6_TEST <- Output6_TEST %>% 
  mutate(Age = Age_at_Diagnosis + Duration) %>%
  filter(Age < 96)

# The next step relates to individuals recently diagnosed with T2DM (i.e., less than a 
# decade).
# 102118: In carrying out simulations based on the results of the dynamic optimization
# exercise, I found that no losses attributable to absenteeism were recorded for year
# 10. The issue appears to be related to these three lines of code, which stipulate
# that adherence levels should be set equal to NA for T2DM durations < 11 years. 
# The intent here was to recognize that individual choices should only become
# part of their adherence "history" once certain milestones were passed. In this case,
# the model assumes people can revisit prior health decisions upon being diagnosed, and
# again after having had T2DM for 9 years. As regards the latter, we assume that
# whatever choices were made in Year 9 become part of their permanent "record" beginning
# in Year 10. 
# 102218: On first glance, one might assume that values in the columns referenced below 
# should be set to NA up to Year 9, NOT Year 10. HOWEVER, this was in fact correct, and
# changing the value from 11 to 10 caused issues - in particular, the number of model
# states increased from 228,357 to 229,168. This occurred because the number of states
# associated with duration == 9 increased by almost a factor of magnitude. Upon further
# investigation (with reference to manipulations undertaken with RWD_T2DM), I determined 
# that at one point in the code, both duration and age are decremented by one (i.e., 11
# becomes 10, 10 becomes 9, etc.). If we set the value below to 10, it therefore means
# that NA values are not ultimately applied to columns recording second-phase adherence
# levels, meaning there are more unique duration 9 states than there ought to be. As a 
# result, when we eliminate duplication states, we end up retaining some that we had
# really intended to drop. Hence - too many states.
Output6_TEST$Per_2_Adherence_LTPA[Output6_TEST$Duration < 11] <- NA
Output6_TEST$Per_2_Adherence_Diet[Output6_TEST$Duration < 11] <- NA
Output6_TEST$Per_2_Adherence_Medication[Output6_TEST$Duration < 11] <- NA
# So we ran up against a problem: we would like to get rid of duplicate observations for
# individuals who have had T2DM for less than a decade, but found that because the values
# derived from the simulations were not exactly equal (even when all categorical variables
# were), we ended up with multiple observations for some combinations of health behaviours.
# One way to address this would be to simply utilize the first value (the theory being that
# the values are different only because the simulation results are subject to variability).
# The next line of code does precisely this. In words, it specifies that we should keep only
# rows that are not duplicates (i.e., of elements with smaller subscripts). But we need to
# go further, and indicate which columns should be utilized as the basis for the 
# determination that a row is or is not unique. The term in parentheses (i.e., 
# "Output6_TEST[1:10]") accomplishes this nicely.
# https://stackoverflow.com/questions/7790732/unique-for-more-than-one-variable
Output6.02 <- Output6_TEST[!duplicated(Output6_TEST[1:10]),]

Emp_Stat_1 <- data.frame(c("Employed", "Unemployed"))
Output6.02 <- Output6.02 %>% mutate("Emp_Stat_1" = NA)

Output6.02 <- Output6.02[rep(row.names(Output6.02),2),] 
Output6.02[c(1:(nrow(Output6.02)/2)),"Emp_Stat_1"] <- "Employed"
Output6.02[c((nrow(Output6.02)/2)+1):nrow(Output6.02),"Emp_Stat_1"] <- "Unemployed"

if (SensAnal == 1) {
  
} else {
  
  # * Identify all possible states (i.e., state space) ====
  
  Factors <- 
    list(Gender = GEN_K[[Gender.U]], Ethnicity = ETH_K[[Ethnicity.U]], 
         Age_at_Diagnosis = NA, Age_1 = (seq(1:70)+Start_Age_K), Duration_1 = NA, 
         Health_Stat_1 = "Healthy", Per_1_Adherence_LTPA_1 = NA, Per_2_Adherence_LTPA_1 = NA, 
         Per_1_Adherence_Diet_1 = NA, Per_2_Adherence_Diet_1 = NA, 
         Per_1_Adherence_Medication_1 = NA, Per_2_Adherence_Medication_1 = NA, 
         Emp_Stat_1 = c("Employed", "Unemployed"))
  
  TP_Healthy.01 <- expand.grid(Factors, stringsAsFactors = FALSE)
  
  TP_T2DM <- Output6.02 %>%
    subset(., select = -c(11:16)) %>%
    add_column(., "Health_Stat_1" = "T2DM", .before = "Age_at_Diagnosis") %>%
    rename("Age_1" = "Age", "Duration_1" = "Duration",
           "Per_1_Adherence_LTPA_1" = "Per_1_Adherence_LTPA", 
           "Per_2_Adherence_LTPA_1" = "Per_2_Adherence_LTPA", 
           "Per_1_Adherence_Diet_1" = "Per_1_Adherence_Diet", 
           "Per_2_Adherence_Diet_1" = "Per_2_Adherence_Diet", 
           "Per_1_Adherence_Medication_1" = "Per_1_Adherence_Medication", 
           "Per_2_Adherence_Medication_1" = "Per_2_Adherence_Medication")
  
  TP_T2DM <- TP_T2DM[,c(1:4, 12, 5:11, 13)]
  TP_T2DM[TP_T2DM == "Not Applicable"] <- NA
  TP_T2DM[,"Duration_1"] <- TP_T2DM$Duration_1 - 1 
  TP_T2DM[,"Age_1"] <- TP_T2DM$Age_1 - 1
  TP_T2DM[TP_T2DM$Duration_1 == 0, c("Per_1_Adherence_LTPA_1", "Per_1_Adherence_Diet_1",
                                     "Per_1_Adherence_Medication_1")] <- NA
  TP_T2DM <- unique(TP_T2DM)
  
  TP_Death <- 
    c(Gender = GEN_K[[Gender.U]], Ethnicity = ETH_K[[Ethnicity.U]], 
      Age_at_Diagnosis = NA, Age_1 = NA, Duration_1 = NA, Health_Stat_1 = "Dead",
      Per_1_Adherence_LTPA_1 = NA, Per_2_Adherence_LTPA_1 = NA, 
      Per_1_Adherence_Diet_1 = NA, Per_2_Adherence_Diet_1 = NA, 
      Per_1_Adherence_Medication_1 = NA, Per_2_Adherence_Medication_1 = NA, 
      Emp_Stat_1 = NA)
  
  # This is the list of possible states the individual can occupy
  State_List_TP <- rbind(TP_Healthy.01, TP_T2DM, TP_Death)
  State_List_TP.02a <- State_List_TP[, c(4:13)]
  rownames(State_List_TP.02a) <- seq(1, nrow(State_List_TP.02a))
  State_List_TP.02a <- rownames_to_column(State_List_TP.02a, var = "ID")
  State_List_TP.02b <- State_List_TP.02a %>%
    rename("Age_2" = "Age_1", "Duration_2" = "Duration_1", 
           "Health_Stat_2" = "Health_Stat_1", "Emp_Stat_2" = "Emp_Stat_1",
           "Per_1_Adherence_LTPA_2" = "Per_1_Adherence_LTPA_1", 
           "Per_2_Adherence_LTPA_2" = "Per_2_Adherence_LTPA_1", 
           "Per_1_Adherence_Diet_2" = "Per_1_Adherence_Diet_1", 
           "Per_2_Adherence_Diet_2" = "Per_2_Adherence_Diet_1", 
           "Per_1_Adherence_Medication_2" = "Per_1_Adherence_Medication_1", 
           "Per_2_Adherence_Medication_2" = "Per_2_Adherence_Medication_1")
  
  # * Identify all possible decisions (i.e., choice set) ####
  
  # Now let's enumerate all possible choices the individual can make.
  Decision_Combos <- 
    list(AD_LTPA <- c("Adhere", "Don't Adhere"), AD_Diet <- c("Adhere", "Don't Adhere"),
         AD_Medication <- c("Adhere", "Don't Adhere"), LMP <- c("Full Time", "Part Time",
                                                                "Retire"))
  
  Possible_Decisions <- expand.grid(Decision_Combos)
  Possible_Decisions <- rownames_to_column(Possible_Decisions)
  
  names(Possible_Decisions) <- c("ID_Actions", "AD_LTPA", "AD_Diet", "AD_Medication", "LMP")
  # http://www.sthda.com/english/wiki/writing-data-from-r-to-txt-csv-files-r-base-functions
  # 051218: Learned valuable lesson regarding GSUB. One MUST specify the object to which
  # the function is being applied in the third function slot; otherwise, R appears to
  # misinterpret what is to replace what and returns an error message (i.e., "argument 
  # 'pattern' has length > 1 and only the first element will be used").
  # https://stackoverflow.com/questions/29271549/replace-all-occurrences-of-a-string-in-a-data-frame
  Ref_Actions <- data.frame(lapply(Possible_Decisions, function(x){x <- x %>%
    gsub("Don't Adhere", "DA",.) %>% gsub("Adhere", "A",.) %>% 
    gsub("Full Time", "FT", .) %>% gsub("Part Time", "PT", .) %>%
    gsub("Retire", "RT", .)}), stringsAsFactors = FALSE)
  
  write.table(Ref_Actions, 
              file = paste(RDA, "Ref_Actions", ".dat", sep = ""), quote = FALSE, 
              sep = ",", row.names = FALSE)
  
  Decision_List <- vector("list", 24)
  
  for (i in 1:24){
    Decision_List[[i]] = Possible_Decisions[i, ]
  }
  
  # Transition Probability Matrix ----
  
  # * Set Up Data Structures for Transition Probability Matrix ====
  
  # * * Individuals who are initially healthy ####
  Healthy_Trans_A.01 <-   
    list(Health_Stat_1 = "Healthy", Health_Stat_2 = "Healthy", Age_2 = NA, 
         Duration_2 = NA, Per_1_Adherence_LTPA_2 = NA, Per_2_Adherence_LTPA_2 = NA, 
         Per_1_Adherence_Diet_2 = NA, Per_2_Adherence_Diet_2 = NA, 
         Per_1_Adherence_Medication_2 = NA, Per_2_Adherence_Medication_2 = NA, 
         Emp_Stat_2 = c("Employed", "Unemployed"))
  
  Healthy_Trans_A.02 <- expand.grid(Healthy_Trans_A.01, stringsAsFactors = FALSE)
  
  Healthy_Trans_B.01 <-   
    list(Health_Stat_1 = "Healthy", Health_Stat_2 = "T2DM", Age_2 = NA, Duration_2 = 0, 
         Per_1_Adherence_LTPA_2 = NA, Per_2_Adherence_LTPA_2 = NA, 
         Per_1_Adherence_Diet_2 = NA, Per_2_Adherence_Diet_2 = NA, 
         Per_1_Adherence_Medication_2 = NA, Per_2_Adherence_Medication_2 = NA, 
         Emp_Stat_2 = c("Employed", "Unemployed"))
  
  Healthy_Trans_B.02 <- expand.grid(Healthy_Trans_B.01, stringsAsFactors = FALSE)
  
  Healthy_Trans_C.01 <- 
    list(Health_Stat_1 = "Healthy", Health_Stat_2 = "Dead", Age_2 = NA, Duration_2 = NA, 
         Per_1_Adherence_LTPA_2 = NA, Per_2_Adherence_LTPA_2 = NA, 
         Per_1_Adherence_Diet_2 = NA, Per_2_Adherence_Diet_2 = NA, 
         Per_1_Adherence_Medication_2 = NA, Per_2_Adherence_Medication_2 = NA, 
         Emp_Stat_2 = NA)
  
  Healthy_Trans_C.02 <- expand.grid(Healthy_Trans_C.01, stringsAsFactors = FALSE)
  
  Healthy_Trans.01 <- rbind(Healthy_Trans_A.02,Healthy_Trans_B.02,Healthy_Trans_C.02)
  
  TP_Healthy.02 <- left_join(TP_Healthy.01, Healthy_Trans.01, by = "Health_Stat_1" )
  TP_Healthy.02[,"Age_2"] <- TP_Healthy.02$Age_1 + 1
  # Only one transition is possible at age 94, so filter out all other conceivable
  # state changes. Additionally, people are assumed not to begin transitioning to T2DM
  # until after age 26, so filter out transitions to T2DM at age 25.
  TP_Healthy.02 <- TP_Healthy.02 %>%
    filter(., !(Age_1==94 & !Health_Stat_2=="Dead")) %>%
    filter(., !(Age_1==25 & Health_Stat_2=="T2DM"))
  
  # * * Individuals who are initially diabetic ####
  
  # For the most part, only two transitions are possible from such states - continuing 
  # on with T2DM, or dying.
  
  # This is a little more complicated than transitions initiated while in the healthy
  # state, so we split it into three pieces, beginning with Years 0 to 9.
  
  T2DM_Trans_0_9.01 <-   
    list(Health_Stat_1 = "T2DM", Health_Stat_2 = "T2DM", Age_2 = NA, 
         Duration_2 = NA, Per_1_Adherence_LTPA_2 = c("Adherent","Non-Adherent"), 
         Per_2_Adherence_LTPA_2 = NA, 
         Per_1_Adherence_Diet_2 = c("Adherent","Non-Adherent"), 
         Per_2_Adherence_Diet_2 = NA, 
         Per_1_Adherence_Medication_2 = c("Adherent","Non-Adherent"), 
         Per_2_Adherence_Medication_2 = NA, 
         Emp_Stat_2 = c("Employed", "Unemployed"))
  
  T2DM_Trans_0_9.02 <- expand.grid(T2DM_Trans_0_9.01, stringsAsFactors = FALSE)
  
  # As just noted, we are focused on Years 0 to 9 right now, so filter out everything
  # else from the original list of T2DM states.
  TP_T2DM_0_9.01 <- TP_T2DM %>%
    filter(Duration_1 < 9)
  
  TP_T2DM_0_9.02 <- left_join(TP_T2DM_0_9.01, T2DM_Trans_0_9.02, by = "Health_Stat_1")
  # Filter so as to retain only rows that reflect continuity in terms of adherence to
  # health-related behaviours (i.e., if the individual exercised last year, they should
  # exercise this year as well). There is, of course, one exception (i.e., Duration = 
  # 10), which we will get to shortly.
  TP_T2DM_0_9.02 <- TP_T2DM_0_9.02 %>%
    filter(Duration_1 == 0 | Duration_1 > 0 & 
             (Per_1_Adherence_LTPA_1 == Per_1_Adherence_LTPA_2) & 
             (Per_1_Adherence_Diet_1 == Per_1_Adherence_Diet_2) &
             (Per_1_Adherence_Medication_1 == Per_1_Adherence_Medication_2))
  
  # Increment age and T2DM duration by 1.
  TP_T2DM_0_9.02[,"Age_2"] <- TP_T2DM_0_9.02$Age_1 + 1
  TP_T2DM_0_9.02[,"Duration_2"] <- TP_T2DM_0_9.02$Duration_1 + 1
  
  T2DM_Trans_9_10.01 <- T2DM_Trans_10_94.01 <-   
    list(Health_Stat_1 = "T2DM", Health_Stat_2 = "T2DM", Age_2 = NA, 
         Duration_2 = NA, Per_1_Adherence_LTPA_2 = c("Adherent","Non-Adherent"), 
         Per_2_Adherence_LTPA_2 = c("Adherent","Non-Adherent"), 
         Per_1_Adherence_Diet_2 = c("Adherent","Non-Adherent"), 
         Per_2_Adherence_Diet_2 = c("Adherent","Non-Adherent"), 
         Per_1_Adherence_Medication_2 = c("Adherent","Non-Adherent"), 
         Per_2_Adherence_Medication_2 = c("Adherent","Non-Adherent"), 
         Emp_Stat_2 = c("Employed", "Unemployed"))
  
  T2DM_Trans_9_10.02 <- T2DM_Trans_10_94.02 <- 
    expand.grid(T2DM_Trans_10_94.01, stringsAsFactors = FALSE)
  
  # Analogous to the first part of this process, we are focused on Years 10 to 94 right 
  # now, so filter out everything else from the original list of T2DM states.
  TP_T2DM_10_94.01 <- TP_T2DM %>%
    filter(Duration_1 > 9 & Duration_1 < 68)
  
  TP_T2DM_10_94.02 <- left_join(TP_T2DM_10_94.01, T2DM_Trans_10_94.02, 
                                by = "Health_Stat_1")
  
  TP_T2DM_10_94.02 <- TP_T2DM_10_94.02 %>%
    filter((Per_1_Adherence_LTPA_1 == Per_1_Adherence_LTPA_2) &
             (Per_1_Adherence_Diet_1 == Per_1_Adherence_Diet_2) &
             (Per_1_Adherence_Medication_1 == Per_1_Adherence_Medication_2) &
             (Per_2_Adherence_LTPA_1 == Per_2_Adherence_LTPA_2) & 
             (Per_2_Adherence_Diet_1 == Per_2_Adherence_Diet_2) &
             (Per_2_Adherence_Medication_1 == Per_2_Adherence_Medication_2))
  
  # Increment age and T2DM duration by 1.
  TP_T2DM_10_94.02[,"Age_2"] <- TP_T2DM_10_94.02$Age_1 + 1
  TP_T2DM_10_94.02[,"Duration_2"] <- TP_T2DM_10_94.02$Duration_1 + 1
  
  # Finally, we look at transitions occurring between the 9th and 10th years of T2DM.
  TP_T2DM_9_10.01 <- TP_T2DM %>% filter(Duration_1 == 9)
  
  TP_T2DM_9_10.02 <- left_join(TP_T2DM_9_10.01, T2DM_Trans_9_10.02, 
                               by = "Health_Stat_1")
  
  TP_T2DM_9_10.02 <- TP_T2DM_9_10.02 %>%
    filter((Per_1_Adherence_LTPA_1 == Per_1_Adherence_LTPA_2) &
             (Per_1_Adherence_Diet_1 == Per_1_Adherence_Diet_2) &
             (Per_1_Adherence_Medication_1 == Per_1_Adherence_Medication_2))
  
  TP_T2DM_9_10.02[,"Age_2"] <- TP_T2DM_9_10.02$Age_1 + 1
  TP_T2DM_9_10.02[,"Duration_2"] <- TP_T2DM_9_10.02$Duration_1 + 1
  
  T2DM_Trans_Dead.01 <- 
    list(Health_Stat_1 = "T2DM", Health_Stat_2 = "Dead", Age_2 = NA, Duration_2 = NA, 
         Per_1_Adherence_LTPA_2 = NA, Per_2_Adherence_LTPA_2 = NA, 
         Per_1_Adherence_Diet_2 = NA, Per_2_Adherence_Diet_2 = NA, 
         Per_1_Adherence_Medication_2 = NA, Per_2_Adherence_Medication_2 = NA, 
         Emp_Stat_2 = NA)
  
  T2DM_Trans_Dead.02 <- expand.grid(T2DM_Trans_Dead.01, stringsAsFactors = FALSE)
  
  TP_T2DM_Dead.02 <- left_join(TP_T2DM, T2DM_Trans_Dead.02, by = "Health_Stat_1")
  
  TP_T2DM.02 <- rbind(TP_T2DM_0_9.02, TP_T2DM_9_10.02, TP_T2DM_10_94.02, 
                      TP_T2DM_Dead.02)
  
  # As with individuals beginning a period in the healthy state, people aged 94 with
  # T2DM are assumed to die before their next birthday, so we filter out all other
  # conceivable transitions.
  TP_T2DM.02 <- TP_T2DM.02 %>% filter(., !(Age_1==94 & !Health_Stat_2=="Dead"))
  
  # * * Transitions from death ####
  
  # This is the simplest case, since there is only one possible transition from death,
  # and that is back to death.
  TP_Dead.02 <-   
    t(c(Gender = GEN_K[[Gender.U]], Ethnicity = ETH_K[[Ethnicity.U]], Health_Stat_1 = "Dead", 
        Age_at_Diagnosis = NA, Age_1 = NA, Duration_1 = NA, Per_1_Adherence_LTPA_1 = NA, 
        Per_2_Adherence_LTPA_1 = NA, Per_1_Adherence_Diet_1 = NA, 
        Per_2_Adherence_Diet_1 = NA,Per_1_Adherence_Medication_1 = NA, 
        Per_2_Adherence_Medication_1 = NA, Emp_Stat_1 = NA,
        Health_Stat_2 = "Dead", Age_2 = NA, Duration_2 = NA, 
        Per_1_Adherence_LTPA_2 = NA, Per_2_Adherence_LTPA_2 = NA, 
        Per_1_Adherence_Diet_2 = NA, Per_2_Adherence_Diet_2 = NA, 
        Per_1_Adherence_Medication_2 = NA, Per_2_Adherence_Medication_2 = NA, 
        Emp_Stat_2 = NA))
  
  # * * Compile complete list of all possible transitions ####
  
  Transition_List.01 <- rbind(TP_Healthy.02, TP_T2DM.02, TP_Dead.02)
  Transition_List.02 <- Transition_List.01[, c(3:23)]
  
  # This list will ultimately hold the rows of the transition probability matrix.
  TP <- vector("list", 24)
  # https://www.datamentor.io/r-programming/examples/odd-even
  
  # * Populating the matrix ====
  
  # * * Likelihood of mortality ####
  
  # Start with mortality risk. Our reference case is an Asian-Indian Female, so let's
  # filter out everything else using DPLYR (noting that the life tables do not differentiate
  # by ethnicity anyway). In addition, we are interested only in individuals between the ages
  # of 25 and 95.
  Mortality.01 <- TP.01
  Mortality.01$Age <- as.integer(Mortality.01$Age)
  
  Mortality.01 <- Mortality.01 %>%
    filter(Gender == GEN_K[[Gender.U]], Ethnicity == ETH_K[[Ethnicity.U]], 
           Age >= 25 & Age <= 94)
  
  Mortality.02 <- unique(Mortality.01[,-c(2:4)])
  Mortality.02 <- Mortality.02 %>%
    rename(Age_1 = Age) %>%
    rename(p_live = Probability) %>%
    mutate(p_mort = 1 - p_live) %>%
    mutate("Health_Stat_2a" = "Dead") %>%
    mutate("Health_Stat_2b" = "T2DM") %>%
    mutate("Health_Stat_2c" = "Healthy")
  
  Mortality.02[Mortality.02$Age==94,"p_live"] = 0.0000
  Mortality.02[Mortality.02$Age==94,"p_mort"] = 1.0000
  Mortality.03a = subset(Mortality.02, 
                         select = c("Age_1","Health_Stat_2a","p_mort"))
  Mortality.03b = subset(Mortality.02, 
                         select = c("Age_1","Health_Stat_2b","p_live"))
  Mortality.03c = subset(Mortality.02, 
                         select = c("Age_1","Health_Stat_2c","p_live"))
  
  Mortality.03 <- bind_rows(Mortality.03a, Mortality.03b, Mortality.03c)
  Mortality.03 <- Mortality.03 %>%
    mutate(p_event = coalesce(p_mort, p_live)) %>%
    mutate(Health_Stat_2 = coalesce(Health_Stat_2a, Health_Stat_2b, Health_Stat_2c)) %>%
    subset(., select = -c(p_mort, p_live, Health_Stat_2a, Health_Stat_2b, Health_Stat_2c))
  
  TP_Healthy_Dead <- Mortality.03a[,c(1,3)]
  TP_Healthy_Dead <- TP_Healthy_Dead %>%
    mutate(., "Health_Stat_1" = "Healthy", "Health_Stat_2" = "Dead") %>%
    rename(., "p_mort.healthy" = "p_mort")
  TP_Healthy_Live <- TP_Healthy_Dead %>% mutate(., "p_live.healthy" = 
                                                  1 - p_mort.healthy)
  TP_Healthy_Live <- TP_Healthy_Live[, -c(2,4)]
  # These next four lines may not actually be necessary. EDIT: In fact, they ARE
  # necessary, primarily because we need to prevent merging on rows in which the
  # person experiences mortality by the end of the year.
  TP_Healthy_Live <- TP_Healthy_Live[rep(row.names(TP_Healthy_Live), 2), ]
  TP_Healthy_Live$Health_Stat_2 <- NA
  TP_Healthy_Live[c(1:nrow(TP_Healthy_Live)/2), "Health_Stat_2"] <- "Healthy"
  TP_Healthy_Live[c((nrow(TP_Healthy_Live)/2+1):nrow(TP_Healthy_Live)), 
                  "Health_Stat_2"] <- "T2DM"
  
  # This integrates the estimated annual likelihood of mortality from the life tables
  # into the TP matrix. We call LAPPLY twice. The first time, we merge data pertaining
  # to the likelihood of mortality, whereas the second, we carry out an analogous
  # operation using data pertaining to the likelihood of remaining alive. It is
  # important to note that this process works because we deliberately created DF
  # columns for TP_Healthy_Dead and TP_Healthy_Live that specify initial and final
  # health status, enabling us to link the corresponding observations to ONLY those
  # rows of each element of the TP list with individuals experiencing these transitions.
  Transition_List.02 <- Transition_List.02 %>% 
    merge(., TP_Healthy_Dead, 
          by = c("Health_Stat_1", "Health_Stat_2", "Age_1"), all.x = TRUE) %>%
    merge(., TP_Healthy_Live, 
          by = c("Health_Stat_1", "Health_Stat_2", "Age_1"), all.x = TRUE)
  
  # 050518: Originally, assigning ID numbers to each entry in the TP matrix was
  # complicated by the fact that integer values were being assigned to AGE_2 for
  # healthy people who die in the period under consideration. Since there is no
  # age assigned to this particular combination of statuses (i.e., Health_Stat_2 =
  # "Dead" and "Age_2" = Age_1 + 1) in STATE_LIST_TP.02B, the final assignment of IDs
  # was returning NA values for these transitions.
  Transition_List.02$Age_2[Transition_List.02$Health_Stat_1 == "Healthy" & 
                             Transition_List.02$Health_Stat_2 == "Dead"] <- NA
  
  # For the first year of T2DM, we assume that the risk of mortality is the same as it
  # would be for individuals who are healthy (I would prefer not to use the UKPDS for
  # this, as we are implicitly assuming that people have not yet selected a pattern of
  # health-related behaviour to manage their T2DM). Since the likelihood of an event
  # depends not only on the state the individual occupies, but also the choices they
  # make in any given period, we could potentially still use the UKPDS for this, but
  # for the moment we stick with this approach below.
  TP_T2DM_Dead_0_1 <- TP_Healthy_Dead
  TP_T2DM_Dead_0_1$Health_Stat_1 <- "T2DM" 
  TP_T2DM_Dead_0_1 <- TP_T2DM_Dead_0_1 %>%
    mutate(., "Duration_1" = 0) %>%
    filter(Age_1 > 26) %>%
    rename(., "p_mort.T2DM_0_1" = "p_mort.healthy")
  
  TP_T2DM_Live_0_1 <- TP_Healthy_Live
  row.names(TP_T2DM_Live_0_1) <- seq(1:140)
  TP_T2DM_Live_0_1 <- TP_T2DM_Live_0_1[c(1:70),] 
  TP_T2DM_Live_0_1 <- TP_T2DM_Live_0_1 %>% 
    mutate(., "Duration_1" = 0) %>%
    filter(Age_1 > 26) %>%
    rename("p_live.T2DM_0_1" = "p_live.healthy") 
  TP_T2DM_Live_0_1$Health_Stat_2 <- TP_T2DM_Live_0_1$Health_Stat_1 <- "T2DM"
  
  # We integrate the probabilities in chunks. We begin with Year 0 probabilities, which
  # as already noted, are assumed to be the same as for individuals who are initially
  # healthy.
  Transition_List.02 <- Transition_List.02 %>% 
    # Here we apply effectively the same approach as we did for individuals who are
    # initially healthy. There is one important difference, namely the requirement to
    # specify duration of T2DM (i.e., because we need to limit the merge to a subset of
    # records relating to individuals who have just been diagnosed).
    merge(., TP_T2DM_Dead_0_1, 
          by = c("Health_Stat_1", "Health_Stat_2", "Age_1", "Duration_1"),
          all.x = TRUE) %>%
    merge(., TP_T2DM_Live_0_1, 
          by = c("Health_Stat_1", "Health_Stat_2", "Age_1", "Duration_1"),
          all.x = TRUE) 
  
  # An important question to consider is why we seem to gloss over the assignment of
  # probabilities for transitions occurring between the 9th and 10th years living with
  # T2DM. There are two components to the answer. First, as regards mortality, we do
  # not need to worry about the choices the individual makes.
  TP_T2DM_Dead_1_94 <- Output6.02 %>%
    subset(., select = c(4:10, 13, 17)) %>%
    rename(., "Age_1" = "Age", "Duration_1" = "Duration", 
           "Per_1_Adherence_LTPA_1" = "Per_1_Adherence_LTPA", 
           "Per_2_Adherence_LTPA_1" = "Per_2_Adherence_LTPA", 
           "Per_1_Adherence_Diet_1" = "Per_1_Adherence_Diet", 
           "Per_2_Adherence_Diet_1" = "Per_2_Adherence_Diet", 
           "Per_1_Adherence_Medication_1" = "Per_1_Adherence_Medication", 
           "Per_2_Adherence_Medication_1" = "Per_2_Adherence_Medication") %>%
    mutate(., "Health_Stat_1" = "T2DM") %>%
    mutate(., "Health_Stat_2" = "Dead") %>%
    rename(., "p_mort.T2DM_1_94" = "Probability")
  TP_T2DM_Dead_1_94[TP_T2DM_Dead_1_94 == "Not Applicable"] <- NA
  TP_T2DM_Dead_1_94$Duration_1 <- TP_T2DM_Dead_1_94$Duration_1 - 1
  TP_T2DM_Dead_1_94$Age_1 <- TP_T2DM_Dead_1_94$Age_1 - 1
  TP_T2DM_Dead_1_94 <- unique(TP_T2DM_Dead_1_94)
  # We filter out cases where duration of T2DM < 1 year, because we just addressed this
  # separately.
  TP_T2DM_Dead_1_94 <- filter(TP_T2DM_Dead_1_94, Duration_1 > 0)
  TP_T2DM_Dead_1_94$p_mort.T2DM_1_94[TP_T2DM_Dead_1_94$Age_1==94] <- 1
  
  TP_T2DM_Live_1_94 <- TP_T2DM_Dead_1_94 %>% 
    mutate(., "p_live.T2DM_1_94" = 1 - p_mort.T2DM_1_94)
  TP_T2DM_Live_1_94$Health_Stat_2 <- "T2DM"
  TP_T2DM_Live_1_94 <- subset(TP_T2DM_Live_1_94, select = -c(8))
  TP_T2DM_Live_1_94 <- TP_T2DM_Live_1_94[, c(1:7,11,8:10)]
  
  # Now we integrate probabilities for individuals who have had T2DM for >= 1 year.
  # Note that we have to go through the labourious process of indicating which columns
  # should be included from the DF with which we will be merging, because we have
  # deliberately chosen NOT to include the probability of living through any given
  # year (we will use this information below).
  Transition_List.02 <- Transition_List.02 %>%
    merge(., TP_T2DM_Dead_1_94, 
          by = c("Health_Stat_1", "Health_Stat_2", "Age_1", "Duration_1", 
                 "Per_1_Adherence_LTPA_1", "Per_2_Adherence_LTPA_1", 
                 "Per_1_Adherence_Diet_1", "Per_2_Adherence_Diet_1",
                 "Per_1_Adherence_Medication_1", "Per_2_Adherence_Medication_1"),
          all.x = TRUE) %>% 
    merge(., TP_T2DM_Live_1_94, 
          by = c("Health_Stat_1", "Health_Stat_2", "Age_1", "Duration_1", 
                 "Per_1_Adherence_LTPA_1", "Per_2_Adherence_LTPA_1", 
                 "Per_1_Adherence_Diet_1", "Per_2_Adherence_Diet_1",
                 "Per_1_Adherence_Medication_1", "Per_2_Adherence_Medication_1"),
          all.x = TRUE)
  
  # NOTE: May want to double-check here that there is no more than one probability in
  # each row.
  Transition_List.02 <- Transition_List.02 %>% 
    # "DS" stands for "demographic status" (i.e., whether the person is alive or dead).
    mutate("p_DS_change" = coalesce(p_mort.healthy, p_live.healthy, p_mort.T2DM_0_1, 
                                    p_live.T2DM_0_1, p_mort.T2DM_1_94, p_live.T2DM_1_94))
  
  # People who have died will inevitably remain in that state.
  Transition_List.02$p_DS_change[Transition_List.02$Health_Stat_1 == "Dead"] <- 1
  Transition_List.02 <- subset(Transition_List.02, select = -c(22:27))
  
  for (i in 1:length(TP)) {
    TP[[i]] = Transition_List.02
  }
  
  # * * Likelihood of other transitions ####
  
  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)
  
  Healthy_T2DM.02 <- data.frame("Age_1" = as.integer(seq(1:69)+24))
  
  Healthy_T2DM.02 <- Healthy_T2DM.02 %>%
    mutate("Age_Group" =
             case_when(Age_1 < 45 ~ "Age_19_44",
                       Age_1 < 65 ~ "Age_45_64",
                       Age_1 >= 65 ~ "Age_65_95")) 
  
  Activity <- c("Active", "Sedentary"); Diet <- c("Recommended", "Excessive")
  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
  
  Healthy_Healthy.02 <- subset(Healthy_T2DM.02, select = c(1,6:8))
  Healthy_Healthy.02[, "Health_Stat_2"] <- "Healthy"
  Healthy_Healthy.02$p_healthy_healthy[Healthy_Healthy.02$Age_1 == 25] <- 1
  
  Healthy_T2DM.02 <- subset(Healthy_T2DM.02, select = c(1,5:7))
  
  Healthy_T2DM_List <- vector("list", 4)
  for (i in 1:length(Healthy_T2DM_List)) {
    Healthy_T2DM_List[[i]] = Healthy_T2DM.02[c(((i-1)*69+1):(i*69)),]
  }
  Healthy_T2DM_List <- rep(Healthy_T2DM_List, 6)
  
  Healthy_Healthy_List <- vector("list", 4)
  for (i in 1:length(Healthy_Healthy_List)) {
    Healthy_Healthy_List[[i]] = Healthy_Healthy.02[c(((i-1)*69+1):(i*69)),]
  }
  Healthy_Healthy_List <- rep(Healthy_Healthy_List, 6)
  
  Healthy_T2DM_List %<>% 
    lapply(., function(x) {x[,1] <- as.character(x[,1]); x})
  Healthy_Healthy_List %<>% 
    lapply(., function(x) {x[,1] <- as.character(x[,1]); x})
  
  TP <- TP %>%
    mapply( FUN = function( i, j ) 
      left_join(x = i, y = j, by = c("Age_1", "Health_Stat_1", "Health_Stat_2"))
      , .
      , Healthy_T2DM_List
      , SIMPLIFY = FALSE ) %>%
    mapply( FUN = function( i, j ) 
      left_join(x = i, y = j, by = c("Age_1", "Health_Stat_1", "Health_Stat_2"))
      , .
      , Healthy_Healthy_List
      , SIMPLIFY = FALSE )
 
  Yr_0_Choices.01 <- list(Duration_1 = 0, Health_Stat_2 = "T2DM", 
                          Per_1_Adherence_LTPA_1 = c("Adherent", "Non-Adherent"), 
                          Per_1_Adherence_Diet_1 = c("Adherent", "Non-Adherent"), 
                          Per_1_Adherence_Medication_1 = c("Adherent", "Non-Adherent"))
  Yr_0_Choices.01 <- expand.grid(Yr_0_Choices.01)
  Yr_0_Choices.01$p_Yr_0_Choice <- 0
  Yr_0_Choices.02 <- vector("list", 8)
  for (i in 1:length(Yr_0_Choices.02)) {
    Yr_0_Choices.02[[i]] <- Yr_0_Choices.01
  }
  
  Yr_0_Choices.02 <- Yr_0_Choices.02 %>%
    mapply( FUN = function( i, j )
    {j[i,6] = 1; j} 
    , seq_along( along.with = . )
    , .
    , SIMPLIFY = FALSE ) 
  
  Yr_0_Choices.02 <- rep(Yr_0_Choices.02, 3)
  
  Yr_9_Choices.02 <- Yr_0_Choices.02 %>% lapply(., FUN = function(x)
  {x <- rename(x, "Per_2_Adherence_LTPA_1" = "Per_1_Adherence_LTPA_1",
               "Per_2_Adherence_Diet_1" = "Per_1_Adherence_Diet_1",
               "Per_2_Adherence_Medication_1" = "Per_1_Adherence_Medication_1", 
               "p_Yr_9_Choice" = "p_Yr_0_Choice"); x})
  
  Yr_9_Choices.02 <- Yr_9_Choices.02 %>% lapply(., FUN = function(x) {x$Duration_1 <- 9; x})
  
  Yr_0_Choices.02 %<>% 
    lapply(., function(x) {x[,1] <- as.character(x[,1]); x}) %>%
    lapply(., function(x) {x[,2:5] <- sapply(x[,2:5], as.character); x})
  Yr_9_Choices.02 %<>% 
    lapply(., function(x) {x[,1] <- as.character(x[,1]); x}) %>%
    lapply(., function(x) {x[,2:5] <- sapply(x[,2:5], as.character); x})
  
  TP <- TP %>%
    mapply( FUN = function( i, j ) 
      left_join(x = i, y = j, 
                by = c("Duration_1","Health_Stat_2", 
                       "Per_1_Adherence_LTPA_2" = "Per_1_Adherence_LTPA_1",
                       "Per_1_Adherence_Diet_2" = "Per_1_Adherence_Diet_1", 
                       "Per_1_Adherence_Medication_2" = "Per_1_Adherence_Medication_1"))
      , .
      , Yr_0_Choices.02
      , SIMPLIFY = FALSE ) %>%
    mapply( FUN = function( i, j ) 
      left_join(x = i, y = j, 
                by = c("Duration_1","Health_Stat_2", 
                       "Per_2_Adherence_LTPA_2" = "Per_2_Adherence_LTPA_1",
                       "Per_2_Adherence_Diet_2" = "Per_2_Adherence_Diet_1", 
                       "Per_2_Adherence_Medication_2" = "Per_2_Adherence_Medication_1"))
      , .
      , Yr_9_Choices.02
      , SIMPLIFY = FALSE )
  
  # * * Employment State Transitions ####
  
  # As noted in the dissertation, the available evidence does not bear out the
  # hypothesis that prediabetes and T2DM affect the likelihood of sustaining or
  # recovering employment (i.e., independent of the action taken by the individual
  # themselves), so these are proxy values we can adjust as required (i.e., in the
  # future, the dataframe Emp_Diff might be loaded from an external source, as
  # opposed to being derived from the structure of the SLID_2011_S6 dataframe).
  Emp_Diff <- SLID_2011_S6[,c(1:4)] %>%
    mutate(d_Healthy.EmpPr_01.01 = Emp_Diff_Healthy_K, 
           d_Pred.EmpPr_01.01 = Emp_Diff_Pred_K, 
           d_T2DM.EmpPr_01.01 = Emp_Diff_T2DM_K)
  
  # Link the proxy or empirical results regarding employment probability differentials
  # to the SLID data.
  SLID_2011_S6.01 <- merge(SLID_2011_S6, Emp_Diff, 
                           by = c("ecsex99", "age_cat", "educ_at", "ml01v28"), 
                           all.x = TRUE)
  
  Emp.01 <- SLID_2011_S6.01 %>%
    # 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, 
                         d_Healthy.EmpPr_01.01, d_Pred.EmpPr_01.01,
                         d_T2DM.EmpPr_01.01))
  
  Emp.01$ml01v28 <- as.character(Emp.01$ml01v28)
  Emp.01$ml01v28[Emp.01$ml01v28=="10"] <- "Employed"
  Emp.01$ml01v28[Emp.01$ml01v28=="20"] <- "Unemployed"
  Emp.01 <- Emp.01 %>% rename("Age_Group" = age_cat)
  
  # The strategy here draws upon the acknowledgement that there are essentially 24
  # distinct probabilities to consider (2 age categories x 2 initial employment states 
  # x 2 final employment states x 3 initial health states). EMP.01, however, consists
  # of 5 observations (2 probabilities of transitioning to employment and unemployment, 
  # respectively, and 3 sets of "adjustment coefficients" that reflect the impact of
  # occupying particular health states [i.e., healthy, prediabetes and T2DM on the 
  # former]) on 4 rows (2 age categories x 2 initial employment states).
  # 
  # We need to restructure this in such a way that there are only two observations in
  # each row - one probability, and one adjustment coefficient. Since each row should
  # reflect a different combination of final employment state and initial health state
  # (i.e., in addition to age category and initial employment state), this suggests
  # that the table we want will consist of 24 rows (i.e., 4 starting rows x 2 final
  # employment states x 3 initial health states). We begin by creating a list that will
  # hold the six sub-tables we are about to construct, that will ultimately be merged
  # to produce what we are looking for.
  Emp_List.01 <- vector("list", 6)
  # To begin with, each list entry holds the same data structure.
  for (i in 1:length(Emp_List.01)){
    Emp_List.01[[i]] <- Emp.01
  }
  
  # Each of the six sub-tables should draw on a different combination of columns from
  # the original table, so we create a second list in which each of these combinations
  # is specified.
  Emp_List.02 <- list(c(3,5), c(3,6), c(3,7), c(4,5), c(4,6), c(4,7))
  
  # Now we reference EMP_LIST.02 to pick off the columns we want from each entry of
  # EMP_LIST.01, which, as already noted, are copies of the original table.
  Emp_List.03 <- Emp_List.01 %>%
    mapply( FUN = function( i, j )
      subset(i, select = c(1:2, j))
      , .
      , Emp_List.02
      , SIMPLIFY = FALSE )
  
  # Here we use the labels associated with columns in each sub-table to generate new
  # columns. The existing column names tell us what data should be in the new columns,
  # so we apply the GREPL function to ascertain the final employment status and initial
  # health status associated with each one.
  Emp_List.03 <- lapply(Emp_List.03, FUN = function(x) {x <- x %>% 
    mutate("Emp_Stat_2" = ifelse(grepl("Unemp", names(x[3]), fixed = TRUE)
                                 , "Unemployed"
                                 , "Employed")) %>%
    mutate("Health_Stat_1" = case_when(
      grepl("Healthy", names(x[4]), fixed = TRUE) == TRUE ~ "Healthy",
      grepl("Pred", names(x[4]), fixed = TRUE) == TRUE ~ "Prediabetes",
      grepl("T2DM", names(x[4]), fixed = TRUE) == TRUE ~ "T2DM")); x})
  
  # Specify the labels to apply to each list entry (this is done because it turns out
  # to be laborious to rename the sub-tables individually, which is necessitated by
  # the fact that each sub-table begins with a slightly different combination of
  # labels).
  Emp_List_Names <- c("Age_Group", "Emp_Stat_1", "p_ES_change_OG", "p_ES_adjustment",
                      "Emp_Stat_2", "Health_Stat_1")
  
  # Now apply the labels we just produced.
  Emp_List.03 <- lapply(Emp_List.03, FUN = function(x) {names(x) <- Emp_List_Names; x})
  
  # Employ the REDUCE function to recursively bind together the six sub-tables into 
  # one.
  Emp.01 <- Reduce(function(x,y) rbind(x,y), Emp_List.03, accumulate = FALSE)
  
  # Generate a new dataframe. This one distinguishes uses integer values to denote age,
  # rather than age categories; to align with the SLID data, however, we create a new
  # column that indicates which SLID age category each individual would fall into,
  # enabling us to join on the SLID results.
  Emp.02 <- data.frame("Age_1" = as.integer(seq(25:94)+24))
  Emp.02 <- Emp.02 %>%
    mutate("Age_Group" =
             case_when(Age_1 < 55 ~ 1,
                       Age_1 >= 55 ~ 2)) %>%
    left_join(., Emp.01[,c("Age_Group","Emp_Stat_1","p_ES_change_OG", 
                           "p_ES_adjustment", "Emp_Stat_2", "Health_Stat_1")], 
              by = "Age_Group") %>%
    # This line of code is important. Here we create a new variable that applies the
    # "adjustment coefficients" to the original employment state transition
    # probabilities. Since the sum of the probabilities of transitioning from one state
    # to any other state must equal one, any amount added to any one of these
    # probabilities must be deducted from one or more of the others to preserve the
    # overall sum. Since an individual can only transition to one of two possible
    # employment states, the implication is that adjustment coefficients with the same
    # magnitude but opposite signs must be applied to the corresponding transition
    # probabilities. By convention, we assume the coefficients reflect the reduced
    # likelihood of sustaining or recovering employment, and, as such, are negative;
    # this means that the coefficient applied to the likelihood of losing one's
  # position or of failing to exit unemployment will be positive.
  mutate("p_ES_change" = ifelse(Emp_Stat_1 == "Employed", 
                                p_ES_change_OG + p_ES_adjustment,
                                p_ES_change_OG - p_ES_adjustment))
  
  # Individuals who are unemployed when they reach the standard age of retirement or
  # who subsequently become unemployed are assumed to retire. As such, the probability
  # of remaining unemployed in all subsequent years is one, while the probability of
  # obtaining new work is zero.
  Emp.02$p_ES_change <- 
    ifelse((Emp.02$Age_1 > 64 & Emp.02$Emp_Stat_1 == "Unemployed" & 
              Emp.02$Emp_Stat_2 == "Unemployed"), 1, Emp.02$p_ES_change)
  Emp.02$p_ES_change <- 
    ifelse((Emp.02$Age_1 > 64 & Emp.02$Emp_Stat_1 == "Unemployed" & 
              Emp.02$Emp_Stat_2 == "Employed"), 0, Emp.02$p_ES_change)
  
  Emp_Trans.01 <- vector("list", 3)
  
  for (i in 1:length(Emp_Trans.01)){
    Emp_Trans.01[[i]] <- Emp.02[, c(1, 3, 6:8)]
  }
  
  # If the individual retires, they are assumed to exit the labour market permanently.
  # Recalling that the final third of the entries in the list of transition probability
  # components is associated with the retirement decision, individuals who are 65 or
  # older when the decision is made are assumed to transition to unemployment in the
  # subsequent period with a probability of one - which, of course, implies that people
  # who have made this choice are certain NOT to sustain employment (i.e., the 
  # probability of remaining in this state is zero). Note that the case where people
  # are at or beyond the standard age of retirement and are unemployed has already been
  # dealt with above.
  Emp_Trans.02 <- Emp_Trans.01 %>%
    mapply( FUN = function( i, j )
      if( i == 3 )
      {j[(j$Age_1 > 64 & j$Emp_Stat_1 == "Employed" & 
            j$Emp_Stat_2 == "Unemployed"),5] = 1; j} 
      else
      {j = j}
      , seq_along( along.with = . )
      , .
      , SIMPLIFY = FALSE ) %>%
    mapply( FUN = function( i, j )
      if( i == 3 )
      {j[(j$Age_1 > 64 & j$Emp_Stat_1 == "Employed" & 
            j$Emp_Stat_2 == "Employed"),5] = 0; j} 
      else
      {j = j}
      , seq_along( along.with = . )
      , .
      , SIMPLIFY = FALSE )
  
  Emp_Trans.02_01 <- list(Emp_Trans.02[[1]]); Emp_Trans.02_01 <- rep(Emp_Trans.02_01, 8)
  Emp_Trans.02_02 <- list(Emp_Trans.02[[2]]); Emp_Trans.02_02 <- rep(Emp_Trans.02_02, 8)
  Emp_Trans.02_03 <- list(Emp_Trans.02[[3]]); Emp_Trans.02_03 <- rep(Emp_Trans.02_03, 8)
  Emp_Trans.02 <- c(Emp_Trans.02_01, Emp_Trans.02_02, Emp_Trans.02_03)
  
  Emp_Trans.02 %<>% 
    lapply(., function(x) {x$Age_1 <- as.character(x$Age_1); x})
  
  TP <- TP %>%
    mapply( FUN = function( i, j ) 
      left_join(x = i, y = j, by = c("Age_1", "Emp_Stat_1", "Emp_Stat_2", "Health_Stat_1"))
      , .
      , Emp_Trans.02
      , SIMPLIFY = FALSE )
  
  # * * Finalize matrix of transition probabilities and save results ####
  
  # Create two new columns - one that combines probabilities for those years in which
  # the individual can modify their health-related behaviors (i.e., Years 0 and 9), and
  # another that combines probabilities relating to transitions in health status for
  # people who are currently healthy.
  TP <- TP %>% 
    lapply(., FUN = function(x) 
      # "HRD" refers to "health-related decision".
    {x <- x %>% mutate(., "p_HRD_change" = coalesce(p_Yr_0_Choice, p_Yr_9_Choice)) %>%
      mutate(., "p_HS_change" = coalesce(p_healthy_T2DM, p_healthy_healthy)) %>%
      # Eliminate the original columns - we don't want to carry around the dead 
      # weight.
      subset(., select = -c(23:26)); x})
  
  # Set NA rows equal to 1, so that even when the column is not pertinent to a
  # particular transition, we don't get an NA when we multiply through.
  # https://stackoverflow.com/questions/19379081/how-to-replace-na-values-in-a-table-for-selected-columns-data-frame-data-tab
  TP <- TP %>% 
    lapply(., FUN = function(x) 
    {x[, 23:25][is.na(x[, 23:25])] <- 1 ; x})
  
  TP <- TP %>% 
    lapply(., FUN = function(x) 
    {x <- x %>% 
      mutate("Probability" = format(p_DS_change*p_HS_change*p_HRD_change*p_ES_change, scientific = FALSE)) %>%
      subset(., select = c(1:21, 26)) %>%
      left_join(., State_List_TP.02a, 
            by = c("Age_1", "Duration_1", "Health_Stat_1", "Emp_Stat_1",
                   "Per_1_Adherence_LTPA_1", "Per_2_Adherence_LTPA_1",
                   "Per_1_Adherence_Diet_1", "Per_2_Adherence_Diet_1",
                   "Per_1_Adherence_Medication_1", "Per_2_Adherence_Medication_1")) %>% 
      left_join(., State_List_TP.02b, 
            by = c("Age_2", "Duration_2", "Health_Stat_2", "Emp_Stat_2",
                   "Per_1_Adherence_LTPA_2", "Per_2_Adherence_LTPA_2",
                   "Per_1_Adherence_Diet_2", "Per_2_Adherence_Diet_2",
                   "Per_1_Adherence_Medication_2", "Per_2_Adherence_Medication_2")); x})
  
  TP.03 <- TP %>% 
    lapply(., FUN = function(x) 
    {x <- x %>% 
      subset(., select = c(23:24, 22)); x})
  
  for (i in seq(TP.03))
  {TP.03[[i]][1:3] <- sapply(TP.03[[i]][1:3], as.numeric)}
 
  TP.03 <- lapply(TP.03, FUN = function(x) arrange(x, ID.x, ID.y))
   
}

# Reward Matrix -----------------------------------------------------------

# * Set Up Data Structures for Reward Matrix ====

# This list will ultimately hold the rows of the reward matrix.
RWD <- vector("list", 24)
# https://www.datamentor.io/r-programming/examples/odd-even

Factors <- 
  list(Gender = GEN_K[[Gender.U]], Ethnicity = ETH_K[[Ethnicity.U]], 
       Age_at_Diagnosis = NA, Age_1 = (seq(1:70)+Start_Age_K), Duration_1 = NA, 
       Health_Stat_1 = "Healthy", Per_1_Adherence_LTPA_1 = NA, Per_2_Adherence_LTPA_1 = NA, 
       Per_1_Adherence_Diet_1 = NA, Per_2_Adherence_Diet_1 = NA, 
       Per_1_Adherence_Medication_1 = NA, Per_2_Adherence_Medication_1 = NA, 
       Emp_Stat_1 = c("Employed", "Unemployed"))

RWD_Healthy.01 <- expand.grid(Factors, stringsAsFactors = FALSE)
RWD_Healthy.01$Adjusted_Utility <- 1
RWD_Healthy.01$Annual_Absences <- 0

# DOUBLE-CHECK TO MAKE SURE THE PROGRAM CORRECTLY FILTERED ON GENDER AND AGE.
RWD_T2DM <- Output6.02 %>%
  subset(., select = -c(11:14)) %>%
  add_column(., "Health_Stat_1" = "T2DM", .before = "Age_at_Diagnosis") %>%
  rename("Age_1" = "Age", "Duration_1" = "Duration",
         "Per_1_Adherence_LTPA_1" = "Per_1_Adherence_LTPA", 
         "Per_2_Adherence_LTPA_1" = "Per_2_Adherence_LTPA", 
         "Per_1_Adherence_Diet_1" = "Per_1_Adherence_Diet", 
         "Per_2_Adherence_Diet_1" = "Per_2_Adherence_Diet", 
         "Per_1_Adherence_Medication_1" = "Per_1_Adherence_Medication", 
         "Per_2_Adherence_Medication_1" = "Per_2_Adherence_Medication")

RWD_T2DM <- RWD_T2DM[,c(1:4, 14, 5:11, 15, 12:13)]
RWD_T2DM[RWD_T2DM == "Not Applicable"] <- NA
RWD_T2DM[,"Duration_1"] <- RWD_T2DM$Duration_1 - 1 
RWD_T2DM[,"Age_1"] <- RWD_T2DM$Age_1 - 1
RWD_T2DM[RWD_T2DM$Duration_1 == 0, c("Per_1_Adherence_LTPA_1", "Per_1_Adherence_Diet_1",
                                     "Per_1_Adherence_Medication_1")] <- NA
RWD_T2DM <- RWD_T2DM[!duplicated(RWD_T2DM[1:13]),]

RWD_Death <- 
  c(Gender = GEN_K[[Gender.U]], Ethnicity = ETH_K[[Ethnicity.U]], Age_at_Diagnosis = NA,
    Age_1 = NA, Duration_1 = NA, Health_Stat_1 = "Dead", Per_1_Adherence_LTPA_1 = NA, 
    Per_2_Adherence_LTPA_1 = NA, Per_1_Adherence_Diet_1 = NA, Per_2_Adherence_Diet_1 = NA, 
    Per_1_Adherence_Medication_1 = NA, Per_2_Adherence_Medication_1 = NA, 
    Emp_Stat_1 = NA, Adjusted_Utility = 0, Annual_Absences = 0)

# This is the list of possible states the individual can occupy
State_List_RWD <- rbind(RWD_Healthy.01, RWD_T2DM, RWD_Death)

rownames(State_List_RWD) <- seq(1, nrow(State_List_RWD))
State_List_RWD <- rownames_to_column(State_List_RWD, var = "ID_States")
State_List_RWD.02 <- State_List_RWD[, c(1, 5:14)]
Ref_States <- data.frame(lapply(State_List_RWD.02, function(x)
{x <- x %>%
  gsub("Non-Adherent", "DA",.) %>% 
  gsub("Adherent", "A",.) %>%
  gsub("Unemployed", "U", .) %>% 
  gsub("Employed", "E", .)}), stringsAsFactors = FALSE)

write.table(Ref_States,
            file = paste0(RDA, "Ref_States", ".dat"), quote = FALSE,
            sep = ",", row.names = FALSE)

for (i in 1:length(RWD)) {
  RWD[[i]] = State_List_RWD
}

RWD <- RWD %>%
  # For whatever reason, I get an error when I try to apply AS.INTEGER or AS.NUMERIC
  # to several columns at once, so for the time being I apply them to each column
  # individually (I will modify this later, if I can).
  lapply(., function(x) {x[,1] <- as.integer(x[,1]); x}) %>%
  lapply(., function(x) {x[,15] <- as.numeric(x[,15]); x}) %>%
  lapply(., function(x) {x[,16] <- as.numeric(x[,16]); x}) %>%
  # The warnings appear traceable to this column (i.e., Column 3).
  lapply(., function(x) {x[,4] <- as.integer(x[,4]); x}) %>%
  lapply(., function(x) {x[,5] <- as.integer(x[,5]); x}) %>%
  lapply(., function(x) {x[,6] <- as.integer(x[,6]); x}) 

# * Linking Income Data ====

Emp_Income.01 <- Emp_Income %>%
  filter(Sex == GEN_K[[Gender.U]], 
         Highest_Degree == EDU_K[Education.U,2]) %>%
  subset(., select = c("Age_Cat", "Med_Emp_Inc_2017.FT", "Med_Emp_Inc_2017.PT")) %>%
  rename("Age_Group" = "Age_Cat")

Emp_Income.02 <- data.frame("Age_1" = as.integer(seq(25:94)+24))
Emp_Income.02 <- Emp_Income.02 %>%
  mutate("Age_Group" =
           case_when(Age_1 < 35 ~ "25_34",
                     Age_1 < 45 ~ "35_44",
                     Age_1 < 55 ~ "45_54",
                     Age_1 < 65 ~ "55_64",
                     Age_1 >= 65 ~ "65_")) %>%
  left_join(., Emp_Income.01, by = "Age_Group") %>%
  rename("Gross_Income.FT" = "Med_Emp_Inc_2017.FT", "Gross_Income.PT" = 
           "Med_Emp_Inc_2017.PT") %>%
  mutate("Emp_Stat_1" = "Employed")

Emp_Income.02a <- Emp_Income.02 %>% subset(., select = -c(Gross_Income.PT))
Emp_Income.02a <- rename(Emp_Income.02a, "Gross_Income" = "Gross_Income.FT")
Emp_Income.02a$Age_1 <- as.character(Emp_Income.02a$Age_1)
Emp_Income.02b <- Emp_Income.02 %>% subset(., select = -c(Gross_Income.FT))
Emp_Income.02b <- rename(Emp_Income.02b, "Gross_Income" = "Gross_Income.PT")
Emp_Income.02b$Age_1 <- as.character(Emp_Income.02b$Age_1)

# The next step is to calculate CPP and EI contributions.
Emp_Income.02a$Age_1 <- as.integer(Emp_Income.02a$Age_1)
Emp_Income.02b$Age_1 <- as.integer(Emp_Income.02b$Age_1)

# 052518: We assume the individual participates in some scheme whereby a fixed
# proportion of their employment income is siphoned out into an income-earning
# retirement fund from which they are required to withdraw, beginning at age 65.
# A fixed proportion of income (established by the parameter Savings_Per_K) goes into
# the fund each year.
# 
# There's one really important implication of all this to bear in mind, which is 
# that annual savings depend ONLY on the individual's productivity/employment
# income, which is itself determined by the individual's gender, age (including
# age-related productivity decline), and level of educational attainment. It does
# NOT depend on the individual's state of health, or even on employment history
# (individuals who are employed part-time or who are unemployed at various points
# throughout their lives save the same amount [in absolute terms] throughout their 
# lives as those who remain employed full-time over their entire careers). This is
# in turn necessitated by the model's inability to keep track of what was done
# when. Ideally, we'd have a complete record of what the agent at did at each point
# in time, which would in turn allow us to determine how much was contributed to
# the fund and, consequently, how much will have been accumulated when the person
# retires.
Emp_Income.02a <- Emp_Income.02a %>%
  # The order of these calculations is important:
  # 1) ARPD
  # 2) Savings
  # 3) TO COME
  # 4) TO COME
  # We apply age-related productivity decline before anything else. Originally, this
  # was done because it was reasoned that ARPD could influence private savings
  # and so should be calculated first. However, ARPD CANNOT influence savings, because 
  # it only begins to affect the individual in the year after the fund begins paying out.
  # Now it is done largely out of convenience.
  mutate(Net_Emp_Inc.FT.00 = 
           ifelse(Age_1<65, 
                  Gross_Income, 
                  Gross_Income*(1-(Age_1-65)*ARPD_K))) %>%
  # We then calculate household savings, which, as noted above, are assumed to represent
  # a fixed proportion of household income. Importantly, this is done before applying
  # income tax. Why?  The rationale here is that private saving in this context involves 
  # an RRSP-like mechanism where earnings aren't taxed until they are withdrawn from the 
  # individual's retirement fund.
  mutate(Savings.New = 
           ifelse(Age_1<65,
                  Net_Emp_Inc.FT.00*Savings_Per_K,
                  0))
# When it comes to tracking accumulation of wealth in the fund, timing matters. When
# should we assume interest is applied to the principal? We shall assume this occurs
# at the end of each period (i.e., following each decision epoch). This means that
# when it is time to begin withdrawing from the fund (i.e., during the decision epoch
# corresponding to age 65), the fund will already have achieved its maximum value.
# The other implication (which is actually related) is that from the second year onward,
# accumulated savings at the time of the decision epoch equal the sum of this year's 
# savings (which have not yet had time to accrue interest) and savings accumulated up
# up the present (which accrued interest over the interval between the last decision
# epoch, and this one).

# Nothing has been accumulated in Year 1, so set this value manually.
Emp_Income.02a$Savings.Accum <- 0
for (i in 2:40) {
  Emp_Income.02a$Savings.Accum[i] <-
    Emp_Income.02a$Savings.New[i] + Emp_Income.02a$Savings.Accum[i-1]*(1+Rate_K)
}

# The next step is to calculate the size of annual withdrawals. This formula comes from:
# Tallarida, R. J. (2015). Pocket Book of Integrals and Mathematical Formulas (5th ed.). 
# Philadelphia, PA: CRC Press. See, in particular, "Amounts to Withdraw for a Specified
# Number of Withdrawals II: Payments at the Beginning of Each Year" (pp. 211-213).
Withdrawals <- 
  (Rate_K*Emp_Income.02a$Savings.Accum[Emp_Income.02a$Age_1==(Year_Ret_K-1)]/
     ((1+Rate_K)-(1/((1+Rate_K)^(Life_Exp[Gender.U,2]-Year_Ret_K-1)))))

# Create a new column that characterizes the impact of the private retirement fund on
# the stream of income available to employed individuals.
Emp_Income.02a <- Emp_Income.02a %>%
  mutate(Income_Stream = 
           ifelse(Age_1<65,
                  # If the individual is younger than age 65, their employment income 
                  # is reduced (in the short-term) by the amount allocated to the fund,
                  # hence the negative sign.
                  -Savings.New,
                  # If the individual is 65 or older, however, their employment income
                  # is SUPPLEMENTED by withdrawals from the fund.
                  ifelse(Age_1<ceiling(Life_Exp[Gender.U,2]),
                         Withdrawals,
                         0))) %>%
  mutate(Net_Emp_Inc.FT.01 = Net_Emp_Inc.FT.00 + Income_Stream) %>%
  mutate(Adj_Pens_Earnings = 
           (pmin(Net_Emp_Inc.FT.00,Max_Pen_Earn_2017_K)/Max_Pen_Earn_2017_K)*
           Max_Pen_Earn_K) %>%
  select(-c(Savings.New, Savings.Accum))

# Now do something similar for part-time workers, EXCEPT that we don't need any of the
# intermediate calculations - we can just borrow them from the first dataframe (since
# we assume the same absolute amount of income is siphoned from both FT and PT workers).
Emp_Income.02b <- Emp_Income.02b %>%
  mutate(Net_Emp_Inc.PT.00 = 
           ifelse(Age_1<65, 
                  Gross_Income, 
                  Gross_Income*(1-(Age_1-65)*ARPD_K))) %>%
  mutate(Income_Stream = 
           Emp_Income.02a$Income_Stream[Emp_Income.02b$Age_1 == 
                                          Emp_Income.02a$Age_1]) %>%
  mutate(Net_Emp_Inc.PT.01 = Net_Emp_Inc.PT.00 + Income_Stream) %>%
  mutate(Adj_Pens_Earnings = 
           (pmin(Net_Emp_Inc.PT.00,Max_Pen_Earn_2017_K)/Max_Pen_Earn_2017_K)*
           Max_Pen_Earn_K)

# In calculating CPP benefits, we note that according to CRA itself (as indicated by
# a CRA representative, since this is nowhere explicitly stated on the department's
# website), RRSP contributions affect neither CPP contributions nor CPP benefits.
# Note that we effectively assume that all individuals earn CPP benefits consistent 
# with a career spent working entirely full-time.
CPP.00 = sum(Emp_Income.02a$Adj_Pens_Earnings[Emp_Income.02a$Age_Group!="65_" &
                                                Emp_Income.02a$Age_1!=25])*(0.25/39)

# The reason we needed to calculate private savings now is that eligibility for OAS/
# GIS depends on the size of other sources of household income, which, in this model,
# includes both private savings and CPP.

# For retirement income, we'll begin by copying the template we used for employment
# income. Note that we don't carry out the same calculations for individuals who work
# part time. Why? This is because the model cannot track employment history. We
# therefore cannot distinguish between people who worked PT for 1 period or 10, or who
# were even unemployed for significant stretches of time. We therefore impose the
# simplying assumption that the same level of retirement income (i.e., the level
# associated with FT employment over the entirety of an individual's career) applies to
# all decision-makers.
Ret_Income <- Emp_Income.02a %>%
  # The CPP is the simplest to provide, since it is just a fixed amount, and depends
  # only on how much was contributed over the person's career (i.e., there is no
  # clawback).
  mutate(CPP = ifelse(Age_1 <65, 0, CPP.00)) %>%
  # Now we move to the OAS. The OAS has a relatively simple structure - every dollar
  # of income in excess of the threshold of $72,809 (for 2017) is clawed back at a
  # rate of 15%. It's not entirely clear whether the clawback affects all OAS income, or
  # only amounts exceeding the threshold. We'll assume the latter.
  mutate(OAS = ifelse(Age_1 < 65, 
                      0,
                      # Noting that the individual may not make RRSP contributions
                      # after age 65 (and hence cannot intentionally decrease their 
                      # income to qualify for more OAS benefits), for the purposes
                      # of calculating OAS amounts, we compare the unsaved proportion
                      # income from private savings and CPP and maximum OAS against the
                      #  OAS threshold.
                      ifelse(Income_Stream + CPP + OAS_Amt_K <= OAS_Threshold_K,
                             # If the sum is less than the threshold, the individual
                             # receives the maximum OAS amount.
                             12*OAS_Amt_K,
                             # If the sum is more than the threshold, there are two
                             # possibilities. One is that the amount withdrawn from the
                             # retirement fund plus CPP is greater than the threshold
                             # (i.e., even before considering OAS). In this case, the
                             # clawback applies to every dollar of OAS received (but
                             # note the PMAX function, which prevents too much from
                             # being deducted).
                             ifelse(Income_Stream + CPP > OAS_Threshold_K,
                                    12*pmax(OAS_Amt_K - 0.15*
                                              (Income_Stream + CPP + OAS_Amt_K - 
                                                 OAS_Threshold_K),0),
                                    # The second possibility is that it is the OAS
                                    # itself which pushes the individual over the
                                    # threshold. In this case, the clawback applies
                                    # only to OAS dollars in excess of the threshold.
                                    # This means the full difference between the
                                    # threshold and CPP/retirement fund income is
                                    # retained.
                                    12*OAS_Threshold_K - (Income_Stream + CPP) + 
                                      # The clawback then applies only to the 
                                      # remainder, which is the difference between the
                                      # maximum OAS amount and the expression in the
                                      # previous line.
                                      12*pmax(OAS_Amt_K - (OAS_Threshold_K - 
                                                             (Income_Stream + CPP)) -
                                                0.15*(Income_Stream + CPP + OAS_Amt_K - 
                                                        OAS_Threshold_K),0))))) %>%
  # Now we do something similar with the GIS. We note that the OAS does not appear
  # in these calculations; this is because neither OAS pensions nor GIS benefits are 
  # not included as income for OAS income-tested benefits (this is assumed to also be
  # true of the OAS, which is why the GIS does not appear in the above calculations).
  # https://www.theglobeandmail.com/globe-investor/personal-finance/how-rrsp-payments-can-help-seniors-with-benefits/article548651/
  # https://www.canada.ca/en/services/benefits/publicpensions/cpp/old-age-security/guaranteed-income-supplement/benefit-amount.html
  # http://client.advisor.ca/tax/understanding-oas-and-gis-clawbacks-11037
  # http://www.advisor.ca/tax/tax-news/how-to-minimize-oas-clawback-234417
  # http://www.moneysense.ca/save/retirement/pensions/oas-clawback-retirement/
  mutate(GIS = ifelse(Age_1 < 65, 
                      0,
                      # The GIS threshold is extremely low (i.e., $3,500), so the 
                      # following situation (i.e., where at least some clawback does 
                      # not take place) is extremely unlikely to occur.
                      ifelse(Income_Stream + CPP <= GIS_Threshold_K,
                             12*GIS_Amt_K,
                             # NOTE 090818: In earlier versions of this file, the next
                             # line read "12*pmax(GIS_Amt_K ...". This was an error. The
                             # GIS parameter is expressed in terms of a monthly amount,
                             # whereas Income_Stream, CPP and GIS_Threshold_K are all 
                             # expressed in terms of annual amounts. Because the first
                             # argument in pmax would have been much smaller than it
                             # should be (whereas the second argument would remain the
                             # correct magnitude), the implication of the error would 
                             # have been that many individuals who are in fact eligible
                             # for the GIS would not have received it.
                             pmax(12*GIS_Amt_K - 0.50*(Income_Stream + CPP - 
                                                         GIS_Threshold_K),0)))) %>%
  # Private savings, CPP, and OAS are all considered taxable income; GIS, however, is
  # not.
  mutate(Total_Income = Income_Stream + CPP + OAS + GIS, 
         Taxable_Income = Total_Income - GIS, Non_Taxable_Income = GIS)

# This is the first reference to the MAPPLY function, which appears frequently throughout
# this script. The essential idea here is that the position of an element within RWD
# tells us which action the individual has taken and, consequently, what reward they will
# obtain - conditional, of course, upon their state.

RWD.02 <-
  mapply(FUN = function(i,j)
    # This first case refers to individuals who have chosen to retire.
    if(ceiling(i/8) %% 3 == 0){
      # 060118: It is important to recognize that the reward matrix calculations for
      # people who have chosen to retire or work are slightly different here. This is
      # because an added step (i.e., adjustment for absenteeism and presenteeism) must
      # be applied to employment, but not retirement income.
      left_join( x = j
             # 060218: Removed "Emp_Stat_1"
             , y = Ret_Income[,c("Age_1", "Total_Income",
                                 "Taxable_Income", "Non_Taxable_Income")]
             , by = c("Age_1"))
      # 051018 - IMPORTANT: MERGING DATA.FRAMES REARRANGES THE COLUMNS. This segment of
      # code caused issues because some list elements (those involving retirement) were
      # sorted or otherwise manipulated differently than the remaining elements. This
      # created a situation where subsequent code dependent on the position of particular
      # columns did not operate the same way for part of the list.
      # ORIGINAL CODE: merge( x = j, y = Gross_Ret_Income, all.x = TRUE )
    } else{
      if(ceiling(i/8) %% 3 == 2){
        left_join( x = j
               , y = Emp_Income.02b[,c("Age_1", "Emp_Stat_1", "Net_Emp_Inc.PT.00",
                                       "Income_Stream")]
               , by = c("Age_1", "Emp_Stat_1"))
      } else {
        left_join( x = j
               , y = Emp_Income.02a[,c("Age_1", "Emp_Stat_1", "Net_Emp_Inc.FT.00",
                                       "Income_Stream")]
               , by = c("Age_1", "Emp_Stat_1"))
      }
    }
    , seq_along(along.with = RWD)
    , RWD
    , SIMPLIFY = FALSE)

RWD.02 <-
  lapply(RWD.02, FUN = function(x) arrange(x, ID_States))

# * Linking Presenteeism Data ====

# Transpose data so we can filter on a particular subset of rows.
CT_2.02 <- data.frame(t(CT_2))
# Convert rownames to first column.
CT_2.02 <- rownames_to_column(CT_2.02, var = "Category")

# https://stackoverflow.com/questions/11121385/repeat-rows-of-a-data-frame
CT_2.02 <- CT_2.02 %>%
  .[c(rep(seq_len(10), each = 2), rep(seq_len(10) + 10, each = 3), 
      rep(seq_len(10) + 20, each = 2), rep(seq_len(10) + 30, each = 3)),]

rownames(CT_2.02) <- seq(1:nrow(CT_2.02))

CT_2.02 <- CT_2.02 %>%
  .[c(1:50, rep(seq_len(50) + 50, each = 2)),]

rownames(CT_2.02) <- seq(1:nrow(CT_2.02))
CT_2.02[, 1] <- gsub("F", "Female.", CT_2.02[, 1])
CT_2.02[, 1] <- gsub("M", "Male.", CT_2.02[, 1])
CT_2.02[c(seq(1, 100, by = 2) + 50), 1] <- 
  gsub("NC", "Afro-Caribbean", CT_2.02[c(seq(1, 100, by = 2) + 50), 1])
CT_2.02[c(seq(2, 100, by = 2) + 50), 1] <- 
  gsub("NC", "Asian-Indian", CT_2.02[c(seq(2, 100, by = 2) + 50), 1])
CT_2.02[, 1] <- gsub("C.", "Caucasian.", CT_2.02[, 1], fixed = TRUE)

CT_2.02 <- CT_2.02 %>% arrange(., Category)

CT_2.02[c(seq(31, 50, by = 2), seq(81, 100, by = 2), seq(131, 150, by = 2)), 1] <- 
  gsub("NoDeg", "Educ_1", CT_2.02[c(seq(31, 50, by = 2), seq(81, 100, by = 2), seq(131, 150, by = 2)), 1])
CT_2.02[c(seq(32, 50, by = 2), seq(82, 100, by = 2), seq(132, 150, by = 2)), 1] <- 
  gsub("NoDeg", "Educ_2", CT_2.02[c(seq(32, 50, by = 2), seq(82, 100, by = 2), seq(132, 150, by = 2)), 1])
CT_2.02[c(seq(1, 30, by = 3), seq(51, 80, by = 3), seq(101, 130, by = 3)), 1] <- 
  gsub("Deg", "Educ_3", CT_2.02[c(seq(1, 30, by = 3), seq(51, 80, by = 3), seq(101, 130, by = 3)), 1])
CT_2.02[c(seq(2, 30, by = 3), seq(52, 80, by = 3), seq(102, 130, by = 3)), 1] <- 
  gsub("Deg", "Educ_4", CT_2.02[c(seq(2, 30, by = 3), seq(52, 80, by = 3), seq(102, 130, by = 3)), 1])
CT_2.02[c(seq(3, 30, by = 3), seq(53, 80, by = 3), seq(103, 130, by = 3)), 1] <- 
  gsub("Deg", "Educ_5", CT_2.02[c(seq(3, 30, by = 3), seq(53, 80, by = 3), seq(103, 130, by = 3)), 1])

CT_2.02 <- CT_2.02 %>%
  filter(., grepl(GEN_K[[Gender.U]], Category, fixed = FALSE) & 
           grepl(ETH_K[[Ethnicity.U]], Category) & 
           grepl(EDU_K[Education.U, 3], Category)) %>%
  .[c(2:nrow(.)),]

CT_2.02 <- data.frame(t(CT_2.02))
rownames(CT_2.02) <- seq(1:nrow(CT_2.02))
# https://stackoverflow.com/questions/32566763/convert-row-values-into-column-names-in-r
colnames(CT_2.02) <- unlist(CT_2.02[row.names(CT_2.02)=='1',])
CT_2.02 <- CT_2.02 %>%
  .[!row.names(.)=='1',] %>%
  mutate(Duration = seq(1:nrow(.)))
CT_2.02[, 5] = CT_2.02[, 5] - 1 
CT_2.02 <- CT_2.02 %>% filter(., Duration < 68)

# https://stackoverflow.com/questions/7303322/apply-function-to-each-column-in-a-data-frame-observing-each-columns-existing-da
CT_2.03 <- as.data.frame(sapply(CT_2.02, function(x) as.character(x)), stringsAsFactors = FALSE)
# https://stackoverflow.com/questions/18503177/r-apply-function-on-specific-dataframe-columns
CT_2.03[1:4] <- sapply(CT_2.03[1:4], function(x) as.numeric(x))
CT_2.03[5] <- sapply(CT_2.03[5], function(x) as.integer(x))

CT_2.04 <- melt(CT_2.03, id.vars = "Duration", 
                variable.names = c(names(CT_2.03)[1:4]), value.name = "Pres")
CT_2.04$variable <- as.character(CT_2.04$variable)

CT_2.04 <- CT_2.04 %>%
  mutate("WC_P1.02" = 
           ifelse(variable == names(CT_2.03)[1] | variable == names(CT_2.03)[3], 
                  "Yes", "No")) %>%
  mutate("WC_P2.02" =
           ifelse(Duration < 10, NA,
                  ifelse(variable == names(CT_2.03)[1] | variable == names(CT_2.03)[4],
                         "Yes", "No")))

CT_2.04 <- CT_2.04[!duplicated(CT_2.04[,c(1,4:5)]),]

Levels <- c("Adherent", "Non-Adherent")

Factors <- 
  list(Duration <- as.character(seq(1:68)-1),
       AD_LTPA_P1 <- Levels, AD_LTPA_P2 <- Levels, 
       AD_DIET_P1 <- Levels, AD_DIET_P2 <- Levels,
       AD_MEDS_P1 <- Levels, AD_MEDS_P2 <- Levels)
Headers <- c("Duration", "Per_1_Adherence_LTPA_1", "Per_2_Adherence_LTPA_1",
             "Per_1_Adherence_Diet_1", "Per_2_Adherence_Diet_1", 
             "Per_1_Adherence_Medication_1", "Per_2_Adherence_Medication_1")

Presenteeism.01 <- expand.grid(Factors)
colnames(Presenteeism.01) <- Headers
# https://stackoverflow.com/questions/13482532/sum-cells-of-certain-columns-for-each-row
Presenteeism.01$Duration <- as.integer(as.character(Presenteeism.01$Duration))

Presenteeism.01 <- Presenteeism.01 %>%
  mutate("WC_P1.01" = rowSums(.[,c(2,4,6)]=="Adherent")) %>%
  mutate("WC_P2.01" = 
           ifelse(Duration < 10, NA, rowSums(.[,c(3,5,7)]=="Adherent"))) %>%
  mutate("WC_P1.02" = ifelse(WC_P1.01 >= 2, "Yes", "No")) %>%
  mutate("WC_P2.02" = ifelse(is.na(WC_P2.01), NA, 
                             ifelse(WC_P2.01 >= 2, "Yes", "No")))
Presenteeism.01[2:7] <- sapply(Presenteeism.01[2:7], function(x) as.character(x))

# 101518: Added "_1" to the end of each of these three variables, which was not included
# in earlier versions of this script. I suspect this may have caused issues during the
# merge between RWD.02 and Presenteeism.02 (i.e., because the values in these three
# fields would not have been modified as intended, and thus would not correctly match 
# with the corresponding variables in RWD.02, which all contain values of NA).
Presenteeism.01$Per_2_Adherence_LTPA_1[Presenteeism.01$Duration <10] <- NA
Presenteeism.01$Per_2_Adherence_Diet_1[Presenteeism.01$Duration <10] <- NA
Presenteeism.01$Per_2_Adherence_Medication_1[Presenteeism.01$Duration <10] <- NA

Presenteeism.02 <- Presenteeism.01[!duplicated(Presenteeism.01[,c(1:7)]),]
Presenteeism.02 <- 
  merge(x = Presenteeism.02, y = CT_2.04[,c("Duration", "WC_P1.02", "WC_P2.02", "Pres")], 
        by = c("Duration", "WC_P1.02", "WC_P2.02"), all.x = TRUE)
Presenteeism.02 <- rename(Presenteeism.02, "Duration_1" = "Duration")

RWD.02 <-
  lapply(RWD.02, FUN = function(x) 
    left_join(x, Presenteeism.02[,c("Duration_1", "Per_1_Adherence_LTPA_1", 
                                "Per_2_Adherence_LTPA_1", "Per_1_Adherence_Diet_1", 
                                "Per_2_Adherence_Diet_1", "Per_1_Adherence_Medication_1",
                                "Per_2_Adherence_Medication_1", "Pres",
                                "WC_P1.02","WC_P2.02")]
          , by = c("Duration_1","Per_1_Adherence_LTPA_1", "Per_2_Adherence_LTPA_1",
                   "Per_1_Adherence_Diet_1", "Per_2_Adherence_Diet_1",
                   "Per_1_Adherence_Medication_1", "Per_2_Adherence_Medication_1")))

RWD.02 <-
  lapply(RWD.02, FUN = function(x) arrange(x, ID_States))

# We assume that people who are healthy or dead do not experience presenteeism. This 
# snippet of code accurately converts all NAs in the column entitled "Pres" into zero
# values.
RWD.03 <- RWD.02 %>% 
  lapply(., FUN = function(x) 
  {x[, "Pres"][is.na(x[, "Pres"])] <- 0 ; x})

# NOTE 100218: Need to augment Ref_States dataframe with projected absences and
# presenteeism data. You can accomplish this by merging Ref_States with data found
# in Output6.02 (which already filters on individual characteristics) and
# Presenteeism.02. Importantly, both absenteeism and presenteeism depend ONLY on the
# individual's state (i.e., not on their actions), which means we don't need the
# optimal policy to estimate productivity losses attributable to T2DM. We may,
# however, need to adjust for the units of time in which absenteeism and
# presenteeism are expressed.
Absenteeism <- Output6.02 %>% mutate(Health_Stat_1 = "T2DM")
Absenteeism <- Absenteeism[,c(17, 4, 19, 5:10, 18, 15)]
colnames(Absenteeism) <- c(colnames(Ref_States[2:11]), "Annual_Absences")
Absenteeism$Emp_Stat_1[Absenteeism$Emp_Stat_1=="Employed"] <- "E"
Absenteeism$Emp_Stat_1[Absenteeism$Emp_Stat_1=="Unemployed"] <- "U"
Absenteeism[,"Duration_1"] <- Absenteeism$Duration_1 - 1 
Absenteeism[,"Age_1"] <- Absenteeism$Age_1 - 1
Absenteeism[Absenteeism$Duration_1 == 0, 
            c("Per_1_Adherence_LTPA_1", "Per_1_Adherence_Diet_1",
              "Per_1_Adherence_Medication_1")] <- NA
Absenteeism$Duration_1 <- as.character(Absenteeism$Duration_1)
Absenteeism$Age_1 <- as.character(Absenteeism$Age_1)

Absenteeism[4:9] %<>%
  sapply(., function(x) gsub("Non-Adherent", "DA", x)) %>%
  sapply(., function(x) gsub("Adherent", "A", x))

Absenteeism <- Absenteeism[!duplicated(Absenteeism[1:10]),]

Presenteeism <- Presenteeism.02 %>% mutate(Health_Stat_1 = "T2DM")
Presenteeism <- Presenteeism[,c(1, 13, 4:9, 12)]
colnames(Presenteeism) <- c(colnames(Ref_States[3:10]), "Pres")

Presenteeism$Duration_1 <- as.character(Presenteeism$Duration_1)
Presenteeism[3:8] %<>%
  sapply(., function(x) gsub("Non-Adherent", "DA", x)) %>%
  sapply(., function(x) gsub("Adherent", "A", x))

Merge_Vars <- c("Age_1", "Duration_1", "Health_Stat_1", "Per_1_Adherence_LTPA_1",
                "Per_2_Adherence_LTPA_1", "Per_1_Adherence_Diet_1",
                "Per_2_Adherence_Diet_1", "Per_1_Adherence_Medication_1",
                "Per_2_Adherence_Medication_1", "Emp_Stat_1")

Ref_States.02 <- left_join(Ref_States, Absenteeism, 
                       by = Merge_Vars)

Ref_States.02 <- left_join(Ref_States.02, Presenteeism, 
                       by = c("Duration_1", "Health_Stat_1", "Per_1_Adherence_LTPA_1",
                              "Per_2_Adherence_LTPA_1", "Per_1_Adherence_Diet_1",
                              "Per_2_Adherence_Diet_1", "Per_1_Adherence_Medication_1",
                              "Per_2_Adherence_Medication_1"))

# 111518: Previously, the ANNUAL_ABSENCES and PRES variables were multiplied by
# the term (Time_Available_K/Annual_Work_Hrs.FT_K) to reflect the fact that illness
# can and does affect productivity outside of work. However, TIME_AVAILABLE_K is 
# 5,840, which is nearly three times the value of the parameter associated with hours
# allocated to full time work (i.e., 2,000 hours); this implies dramatically higher
# estimates of productivity losses than we would obtain by considering work hours
# alone.
# 
# NOTE 112818: IMPORTANT - Eliminating (Time_Available_K/Annual_Work_Hrs.FT_K) when
# constructing Ref_States.02 does NOT affect decision-making, because we continue to 
# apply this term when calculating losses to the individual, in the "Allocation of
# Time" section.
# 
# Absenteeism in the original study by Sorensen and Ploug was presented in days
# per year, so we multiply by the length of one working day (i.e., 8 hours).
Ref_States.02$Annual_Absences <- Ref_States.02$Annual_Absences*8
# Presenteeism in the original study by Lavigne, Phelps, Mushlin and Lednar was
# presented in hours per month, so we multiply by the number of months in a year.
# 111518: The line of code below contained an error, whereby PRES was multiplied
# by 8 rather than 12, as should have been the case. The error has now been corrected.
Ref_States.02$Pres <- Ref_States.02$Pres*12
# Rename the columns, and add a new column that reflects total efficiency losses.
Ref_States.02 %<>% rename(Abs_Hours = Annual_Absences, Pres_Hours = Pres)
Ref_States.02[12:13] <- 
  sapply(Ref_States.02[12:13], FUN = function(x) {x[is.na(x)] <- 0; x})
Ref_States.02 %<>% mutate(Tot_Eff_Losses = Abs_Hours + Pres_Hours)

Time_Value.00 <- Emp_Income.02a[, c("Age_1", "Net_Emp_Inc.FT.00")]
Time_Value.00$Age_1 <- as.character(Time_Value.00$Age_1)
Ref_States.02 <- left_join(Ref_States.02, Time_Value.00,by = c("Age_1"))

# 101518: The RHS of this line was originally Ref_States.02$Net_Emp_Inc.FT.00*
# Time_Available.K. It's unclear what this was intended to accomplish - only a
# single row did NOT contain a value (i.e., the death state). Furthermore, the
# variable in question (i.e., Net_Emp_Inc.FT.00) is expressed as an annual
# figure, not an hourly one - so multiplying by hours worked would have led to
# erroneous results even if Net_Emp_Inc.FT.00 had contained a value. Finally,
# Net_Emp_Inc.FT.00 is empty, so the code effectively involved replacing a 
# blank cell by itself, which generated an error.
Ref_States.02$Net_Emp_Inc.FT.00[is.na(Ref_States.02$Net_Emp_Inc.FT.00)==T] <- 0
Ref_States.02 %<>%
  mutate(Val_Abs_Eff_Losses =
           (Net_Emp_Inc.FT.00/Annual_Work_Hrs.FT_K)*Abs_Hours) %>%
  mutate(Val_Pres_Eff_Losses = 
           (Net_Emp_Inc.FT.00/Annual_Work_Hrs.FT_K)*Pres_Hours) %>%
  mutate(Val_Tot_Eff_Losses =
           Val_Abs_Eff_Losses + Val_Pres_Eff_Losses)

Ref_States.02 <- Ref_States.02[, c(11, 1:10, 16:18)]
Ref_States.02$ID_States <- as.integer(Ref_States.02$ID_States)
Ref_States.02 <- arrange(Ref_States.02, ID_States)

# * Adjustment for Absenteeism/Presenteeism ====

# I would like to be able to use lapply to indicate that people who opt for retirement
# in a given year do not pay taxes on retirement earnings in that year. However, this
# turned out to be more complicated than I expected, so I will put a pin in it for the
# moment.

RWD.03 <- RWD.03 %>%
  # Note that these columns are also being added to the matrices associated with the
  # decision to retire. While not relevant in this specific context, not discriminating
  # simplifies the coding somewhat. We can always drop these columns later.
  lapply(., function(x) {x <- x %>%
    mutate("Gross_Hourly_Wage" = NA) %>%
    mutate("Adj_Absenteeism" = NA) %>%
    mutate("Adj_Presenteeism" = NA); x}) 

# NOTE 042019: Modified calculation of gross hourly wage to allow for adjustment (i.e.,
# in line with the conceptualization of the wage as the "price" of leisure).
RWD.03 <- RWD.03 %>%
  mapply(FUN = function(i, j)
    if (ceiling(i/8) %% 3 == 1)
    {j$Gross_Hourly_Wage <- Mult_w_K*(j$Net_Emp_Inc.FT.00/Annual_Work_Hrs.FT_K); j}
    else
    {if (ceiling(i/8) %% 3 == 2)
    {j$Gross_Hourly_Wage <- Mult_w_K*(j$Net_Emp_Inc.PT.00/Annual_Work_Hrs.PT_K); j} 
      else
        # NOTE 042019: If the person is retired, they do not earn income in the labour
        # market, so we leave Gross_Hourly_Wage as is (i.e., equal to NA).
      {j$Gross_Hourly_Wage <- j$Gross_Hourly_Wage; j}
    }
    , seq_along(along.with = .)
    , .
    , SIMPLIFY = FALSE) %>%
  mapply(FUN = function(i, j)
    # If the individual is not part time, then ...
    if (ceiling(i/8) %% 3 == 1) 
    {j$Adj_Absenteeism <- 
      # If paid sick days exceed annual absences, then there are no adjusted absences 
      # (i.e., the person's sick days entirely compensate the person for illness).
      ifelse(Paid_SDs.FT_K > j$Annual_Absences, 
             0,
             # However, if annual absences exceed available sick days, then the person 
             # cannot make up the difference, and works less than the intended number 
             # of hours.
             j$Annual_Absences - Paid_SDs.FT_K); j}
    else 
    {if (ceiling(i/8) %% 3 == 2)
    {j$Adj_Absenteeism <- 
      ifelse(Paid_SDs.PT_K > j$Annual_Absences/2, 
             0,
             # This follows essentially the same logic as above, with two essential 
             # differences: 
             # 1) Annual absences are divided by two; this is required because the
             # individual is assumed to work half as many hours (note that the original 
             # values upon which the estimated figures are based are assumed to reflect 
             # absenteeism in individuals who are working FT).
             # 2) We apply adifferent sick day allowance (it is assumed that part-
             # time individuals have correspondly less generous allowances).
             j$Annual_Absences/2 - Paid_SDs.PT_K); j} 
      else
      {j$Adj_Absenteeism <- j$Adj_Absenteeism; j}
    }
    , seq_along(along.with = .)
    , .
    , SIMPLIFY = FALSE) %>%
  mapply(FUN = function(i, j)
    # If the individual is not part time, then ...
    if (ceiling(i/8) %% 3 == 1)
      # In the first case (FT), losses associated with presenteeism are divided between 
      # workers and employers as dictated by the parameter "Pres_Share_K".
    {j$Adj_Presenteeism <- Pres_Share_K*j$Pres; j}
    else
      # The second case (PT) essentially follows the first, but we assume that the 
      # burden of presenteeism is less, because the person does not work as frequently 
      # (the implicit assumption being that individuals are not more likely to 
      # experience the symptoms that manifest in presenteeism at work than elsewhere).
      # Again, it is also assumed that the figures used to generate the presenteeism
      # estimates were derived primarily from individuals who work FT (review of the
      # estimates of mean individual annual income in the sample would appear to
      # corroborate this - although the income distribution could also conceivably be
      # skewed).
    {if (ceiling(i/8) %% 3 == 2)
    {j$Adj_Presenteeism <- Pres_Share_K*j$Pres/2; j}
      else
      {j$Adj_Presenteeism <- j$Adj_Presenteeism; j}
    }
    , seq_along(along.with = .)
    , .
    , SIMPLIFY = FALSE)

RWD.03 <- RWD.03 %>%
  lapply(., function(x) {x <- x %>%
    # Absenteeism in the original study by Sorensen and Ploug was presented in days
    # per year, so we multiply by the length of one working day (i.e., 8 hours).
    mutate("Adj_Abs_Hours" = Adj_Absenteeism*8) %>%
    # Presenteeism in the original study by Lavigne, Phelps, Mushlin and Lednar was
    # presented in hours per month, so we multiply by the number of months in a year.
    mutate("Adj_Pres_Hours" = Adj_Presenteeism*12) %>%
    mutate("Work_Eff_Losses" = Adj_Abs_Hours + Adj_Pres_Hours); x})

# Initialize adjusted income column.
RWD.03 <- lapply(RWD.03, FUN = function(x) {x <- x %>% 
  mutate("Adjusted_Income" = NA); x})

# NOTE 091318: Fixed error relating to calculation of adjusted income for
# individuals on EIA. Previously, the lines of code under notes marked by
# asterices read "j$Adjusted_Income", which, as per the line above, is
# simply equal to NA. This lead to a situation whereby unemployed
# individuals were incapable of generating utility.
# 
# This code says the following. If the individual does not decide to retire
# this year, then if they are employed, their wages are adjusted according
# to their levels of absenteeism and presenteeism, and otherwise (i.e., if
# they are unemployed) they simply receive EIA (EIA payments being unaffected
# by the individual's work productivity). However, if the individual DOES
# choose to retire, then if they are younger than 65 they receive a penalty
# of $1,000,000 (to discourage them from making precisely this choice), and
# otherwise they simply receive CPP/OAS/GIS (where again, payments do not
# depend on workplace productivity).
RWD.03 <-
  mapply( FUN = function( i, j )
    if(ceiling(i/8) %% 3 == 0)
      # If the individual retires or is retired, adjusted income has no real relevance
      # to them (as neither absenteeism nor presenteeism influences household 
      # [financial] resources).
    {j$Adjusted_Income <- j$Adjusted_Income; j}
    else
    {if(ceiling(i/8) %% 3 != 2)
    {j$Adjusted_Income <- 
      ifelse(j$Emp_Stat_1 == "Employed",
             (Annual_Work_Hrs.FT_K - j$Work_Eff_Losses)*j$Gross_Hourly_Wage,
             # There are two very important sub-cases in which Adjusted Income is
             # not meaningful even for people who choose to pursue/sustain
             # employment, both of which involve individuals who are initially
             # unemployed. If an unemployed individual is younger than 65, they
             # still have a chance of recovering employment in this or some other
             # subsequent period, but earn EIA in the interim. If the individual is
             # 65 or older, we have assumed they transition automatically to 
             # retirement if they were initially unemployed, and so begin 
             # enjoying public and/or private pensions. In both cases, no 
             # adjustment is required or appropriate.
             # *
             j$Adjusted_Income); j}
      else
      {j$Adjusted_Income <- 
        ifelse(j$Emp_Stat_1 == "Employed",
               (Annual_Work_Hrs.PT_K - j$Work_Eff_Losses)*j$Gross_Hourly_Wage, 
               # *
               j$Adjusted_Income); j}}
    , seq_along( along.with = RWD.03 )
    , RWD.03
    , SIMPLIFY = FALSE )

RWD.03 <- RWD.03 %>%
  mapply( FUN = function( i, j )
    if(ceiling(i/8) %% 3 == 0)
      # Taxable and non-taxable income have already been calculated for people deciding
      # to retire, so there is no need to add such a column.
    {j$Taxable_Income <- j$Taxable_Income; j}
    else
      # For individuals who are employed (either PT or FT), all income is taxable, BUT
      # "all income" includes the negative income streams attributable to the
      # allocation of a share of employed households' resources to the private pension.
      # That is to say, "all income" refers to income net of funds that are saved,
      # since we are implicitly assuming such funds are saved in RRSPs which reduce the
      # overall amount of taxable income.
    {j$Taxable_Income <- ifelse(j$Emp_Stat_1 == "Employed", 
                                j$Adjusted_Income + j$Income_Stream,
                                ifelse(j$Age_1 < 65,
                                       # If the individual is younger than 65, being
                                       # unemployed, we assume they will be receiving
                                       # EIA, which is assumed not to constitute
                                       # taxable income.
                                       0,
                                       # If the individual is age 65 or older and
                                       # unemployed, we assume they transition
                                       # automatically to retirement, beginning in the
                                       # current period, and so are eligible to begin
                                       # earning retirement income. We use the MATCH
                                       # function to extract the appropriate row from
                                       # the retirement income table, noting that the
                                       # order of the arguments in the function
                                       # significantly affects the results received.
                                       Ret_Income$Total_Income[match(j$Age_1,Ret_Income$Age_1)])); j}
    , seq_along( along.with = . )
    , .
    , SIMPLIFY = FALSE ) %>%
  mapply( FUN = function( i, j )
    if(ceiling(i/8) %% 3 == 0)
      # This follows the same logic as above.
    {j$Non_Taxable_Income <- j$Non_Taxable_Income; j}
    else
      # Following the rationale provided above, there is no non-taxable income,
      # strictly speaking.
    {j$Non_Taxable_Income <- ifelse(j$Emp_Stat_1 == "Employed", 
                                    0,
                                    ifelse(j$Age_1 < 65,
                                           # If the individual is younger than 65, being
                                           # unemployed, we assume they will be receiving
                                           # EIA, which is assumed not to constitute
                                           # taxable income.
                                           EIA_General_K,
                                           Ret_Income$GIS[match(j$Age_1,Ret_Income$Age_1)])); j}
    , seq_along( along.with = . )
    , .
    , SIMPLIFY = FALSE )

# * Adjustment for Taxes, and CPP/EI Contributions ====

# This code says the following: 
# 1) If the individual does NOT decide to retire this year, then 
#   a) If they are employed, they pay taxes on their labour market earnings
#   b) Otherwise (i.e., if they are unemployed) they do not pay taxes.
# 2) If the individual DOES retire, then 
#   a) If they are younger than 65, they do not pay taxes (the large
#   penalty in the adjusted income column should ensure they do not make 
#   this choice).
#   b) Otherwise (i.e., if they are age 65+), they pay income on retirement
#   income, which they are assumed to begin earning immediately. Note that
#   the original code had such individuals paying taxes on all of their
#   income (which would not have been appropriate, because the GIS is not
#   considered taxable income).

# 060218: Ran into major issues calculating tax paid. The source of the issue was
# ultimately determined to be tied to the naming of a single variable, LOWER.LIMIT. In
# some parts of the code, the variable was instead referred to as LOWER_LIMIT, which
# was not the variable's name. This, in turn, led to a state of affairs in which a
# relatively large number was subtracted from cumulative taxes paid, resulting in
# negative tax bills. The issue appeared to affect only individuals who chose to retire,
# suggesting the misspelling occurred in only a single place.
RWD.03 <- lapply(RWD.03, function(x)
{x <- x %>%
  mutate(Cum_Tax_Amt = Tax_Rates[
    findInterval(x$Taxable_Income,Tax_Rates$Lower.Limit,all.inside = TRUE),9]) %>%
  mutate(MR = Tax_Rates[
    findInterval(x$Taxable_Income,Tax_Rates$Lower.Limit,all.inside = TRUE),4]) %>%
  mutate(Lower_Limit = Tax_Rates[
    findInterval(x$Taxable_Income,Tax_Rates$Lower.Limit,all.inside = TRUE),1]); x})

RWD.03 <-
  mapply( FUN = function( i, j )
    if( ceiling(i/8) %% 3 == 0 )
    {j$Tax_Paid <- 
      ifelse(j$Age_1 < 65, 
             # If the individual is younger than 65, they don't have any retirement 
             # income, so they shouldn't pay taxes on income they don't have.
             0,
             # Find the tax bracket into which the individual falls, and 
             # determine the cumulative amount of tax they would have paid 
             # across all earlier brackets. We still need to calculate the 
             # amount of tax paid in this bracket. We do this by calculating 
             # the difference between the individual's income and the lower 
             # range of the current bracket, and then multiplying this amount
             # by the corresponding marginal tax rate.
             j$Cum_Tax_Amt + j$MR*(j$Taxable_Income - j$Lower_Limit)); j}
    else
    {j$Tax_Paid <- 
      ifelse(j$Emp_Stat_1 == "Employed",
             # If the individual is employed, the calculation of taxes is as above.
             j$Cum_Tax_Amt + j$MR*(j$Taxable_Income - j$Lower_Limit),
             ifelse(j$Age_1 < 65, 
                    # If the individual is unemployed and below the standard age of 
                    # retirement they earn EIA, which is assumed not to be taxable
                    # income.
                    0,
                    # Again, if the individual is unemployed and age 65 or older,
                    # they transition to retirement, effective immediately. They
                    # therefore pay taxes on their retirement income (excluding the
                    # GIS, of course).
                    j$Cum_Tax_Amt + j$MR*(j$Taxable_Income - j$Lower_Limit))); j} 
    , seq_along( along.with = RWD.03 )
    , RWD.03
    , SIMPLIFY = FALSE )

# RWD.03 <- lapply(RWD.03, function(x) {x$Tax_Paid <- unlist(x$Tax_Paid); x})

# This next line of code, which was part of an earlier version of this script,
# was originally required to remove certain attributes that accompanied the
# "Tax_Paid" columns. However, these attributes do not appear to stick around
# when the mapply function is applied (the reason for this is unclear), so
# this code is not necessary at this time.
# Output6.05$Tax_Paid.FT <- as.vector(Output6.05$Tax_Paid.FT)

# Initialize CPP & EI contributions
RWD.03 <- lapply(RWD.03, FUN = function(x) {x <- x %>% 
  mutate("CPP_EI_Contribs" = NA); x})

RWD.03 <-
  mapply( FUN = function( i, j )
    if( ceiling(i/8) %% 3 == 0 )
      # People who are retired do not pay CPP/EI contributions. These apply only to
      # employment income.
    {j$CPP_EI_Contribs <- 0; j}
    else
    {j$CPP_EI_Contribs <- 
      ifelse(j$Emp_Stat_1 == "Employed",
             # The model assumes that people who are employed pay CPP and EI
             # contributions. Importantly, the amount of the contributions depend on 
             # their EMPLOYMENT income, and not their TAXABLE income (this was
             # confirmed by CRA representatives). The appropriate base for the purposes
             # of these calculations is not gross employment income, but rather
             # employment income net of ARPD and losses attributable to absenteeism and
             # presenteeism - that is, adjusted income.
             pmin(Max_CPP_Contrib_K, CPP_Contrib_Rate_K*(j$Adjusted_Income-Basic_Exemption_K)) +
               pmin(Max_EI_Contrib_K, EI_Contrib_Rate_K*j$Adjusted_Income), 0); j}
    # If the person is not employed, there are again two possibilities. Either
    # the person is younger than 65 (in which case they are an EIA recipient),
    # or they are age 65 or older (in which case they transition to retirement).
    # In either case they do not pay either CPP or EI contributions.
    , seq_along( along.with = RWD.03 )
    , RWD.03
    , SIMPLIFY = FALSE )

# * Pharmacare Program Expenditures ====

RWD.03 <- lapply(RWD.03, FUN = function(x) {x <- x %>% 
  mutate("Rx_Deductible" = NA); x})

# Pharmacare deductibles are based on line 150 of an individual's income tax return,
# which refers to Total Income (see: 
# http://www.gov.mb.ca/health/pharmacare/estimator.html). It is not altogether clear
# whether line 150 includes or excludes CPP/QPP/EI contributions; for simplicity, we
# will assume that if such contributions have tax implications, these are applied
# after the calculation of Total Income. In other words, we assume line 150 refers
# simply to everything the individual earned.
RWD.03 <-
  mapply( FUN = function( i, j )
    if( ceiling(i/8) %% 3 == 0 )
    {j$Rx_Deductible <- 
      ifelse(j$Age_1 < 65, 
             # People who are younger than 65 cannot retire, and so earn no retirement
             # income. Therefore, there isn't much point dedicating computational
             # effort to the calculation of Pharmacare deductibles.
             0,
             # People 65 or older pay a Pharmacare deductible corresponding to their
             # total retirement income, which, as discussed above, is not reduced by
             # any amount of GIS they may earn.
             pmax(Min_Ded_K,
                  data.matrix(j$Total_Income*
                                Pharmacare_K[findInterval(j$Total_Income,
                                                          Pharmacare_K$Lower.Limit, 
                                                          all.inside = TRUE),3]))); j}
    else
      # for people who are not retiring, there are two possibilities, depending on
      # whether or not they are employed.
    {j$Rx_Deductible <- 
      # If the individual is employed and younger than 65, their total income would
      # presumably include any income they will subsequently save, so Adjusted Income,
      # which does not deduct amounts allocated to RRSPs, would be the appropriate
      # base.
      ifelse(j$Emp_Stat_1 == "Employed",
             ifelse(j$Age_1 < 65,
                    pmax(Min_Ded_K,
                         data.matrix(j$Adjusted_Income*
                                       Pharmacare_K[findInterval(j$Adjusted_Income,
                                                                 Pharmacare_K$Lower.Limit, 
                                                                 all.inside = TRUE),3])),
                    # If, however, the individual is aged 65 or older, their total
                    # income would presumably include any amount they receive from
                    # their retirement fund (in this respect there is something of an
                    # asymmetry, in the sense that amounts placed into RRSPs and
                    # amounts withdrawn from the retirement fund are BOTH counted as
                    # income for the purposes of calculating Pharmacare deductibles).
                    # Therefore, the appropriate base for the purposes of the
                    # calculations is not Adjusted Income, which does not include
                    # these amounts, but Taxable Income.
                    pmax(Min_Ded_K,
                         data.matrix(j$Taxable_Income*
                                       Pharmacare_K[findInterval(j$Taxable_Income,
                                                                 Pharmacare_K$Lower.Limit, 
                                                                 all.inside = TRUE),3]))),
             # Now what if the individual is unemployed? Again, there are two
             # possibilities. If they are younger than 65, they still have a chance of
             # recovering employment in this or a future period, but are earning EIA
             # in the interim. Again, people with EIA do not pay Pharmacare
             # deductibles, so this amount should be zero.
             ifelse(j$Age_1 < 65,
                    0,
                    # If the individual is 65 or older, however, being unemployed they are
                    # assumed to transition immediately to retirement, receiving annual
                    # transfers from their retirement fund, as well as from CPP, OAS & GIS,
                    # as applicable. Therefore, the appropriate base for the purposes of
                    # these calculations is the sum of Taxable Income and Non-Taxable
                    # Income, since GIS is considered part of Total Income.
                    pmax(Min_Ded_K,
                         data.matrix((j$Taxable_Income+j$Non_Taxable_Income)*
                                       Pharmacare_K[findInterval((j$Taxable_Income+j$Non_Taxable_Income),
                                                                 Pharmacare_K$Lower.Limit, 
                                                                 all.inside = TRUE),3])))); j}
    , seq_along( along.with = RWD.03 )
    , RWD.03
    , SIMPLIFY = FALSE )

# * Household Expenditures ====

# Start with expenditures associated with LTPA. First initialize the rows 
# that will hold the data.
RWD.03 <- lapply(RWD.03, FUN = function(x) {x <- x %>% 
  mutate("LTPA_Cost" = NA); x})

# We allow LTPA costs to vary by gender, recognizing that these could in principle
# vary according to a range of personal characteristics. It may be, for example, that
# other things equal, older individuals need to spend more to adhere to LTPA than do
# younger individuals (e.g., they may need to investment in special clothing and/or
# equipment to accommodate special needs associated with health conditions). It is
# unclear, however, whether it can be credibly argued that differences in LTPA
# expenditures are not largely driven by personal preferences, which is something we
# have elected not to model - instead adopting the pretense that it is meaningful to
# regard there as being a well-defined "cost of admission" associated with LTPA
# expenditures.
LTPA_Annual_Cost = LTPA_Cost_List_K[Gender.U]

# The rather convoluted structure of the code in this section is intended to force
# people who have previously selected particular patterns of health-related behaviour
# to adopt them consistently. In particular, if a person has chosen to adhere to a
# given behaviour between Year 0 and Year 1, then in every subsequent (until the 9th) year,
# if they affirm their commitment to that behaviour, they incur the associated costs;
# if, however, they break their commitment, they pay an infinitely large penalty.
# Conversely, if they previously chose not to commit to a behaviour, they will incur
# different (and possibly lower) costs so long as they follow this course of action
# consistently, but face an infinitely large penalty for changing their mind. Under
# these circumstances, obviously, the individuals would be expected to stick with 
# whatever choice they originally made.
RWD.03 <-
  mapply( FUN = function( i, j )
    # Recall that with ceiling(i/x), "x" refers to "block size" - that is, the number
    # of consecutive rows in "Possible_Decisions" that involve the same level of
    # adherence to a particular class of health-related behaviour. In this case,
    # every other row of the "Possible_Decisions" DF specifies the same level of
    # adherence to LTPA, so the block size is 1. Also note (recalling that adherence
    # to health-related behaviors is modeled as a dichotomous choice) that if 
    # ceiling(i/x) %% 2 != 0 (i.e., ceiling(i/1) is an odd number), it indicates
    # adherence to a particular behaviour, while ceiling(i/x) %% 2 == 0 denotes non-
    # adherence.
    if( ceiling(i/1) %% 2 != 0 ){
      j$LTPA_Cost <- 
        # If the individual is newly-diagnosed (i.e., duration = 0) OR they are healthy,
        # the cost of participating in LTPA is as defined above.
        ifelse((j$Duration_1 == 0 | j$Health_Stat_1 == "Healthy"), LTPA_Annual_Cost,
               # Otherwise (i.e., if the person has T2DM) then determine whether they
               # have had it for less than 9 years.
               ifelse(j$Duration_1 < 9
                      # If so, they only adhere now if they chose to adhere when they
                      # were originally diagnosed.
                      , ifelse(j$Per_1_Adherence_LTPA == "Adherent"
                               # If they did in fact choose to adhere, they pay the
                               # price of admission.
                               , LTPA_Annual_Cost
                               # Otherwise, as we discussed above, they cannot now
                               # choose to change their mind (i.e., cannot decide to
                               # adhere if they were not previously interested in doing
                               # so), so we penalize them to prevent this from 
                               # occurring.
                               , Inf)
                      # The other possibility is that the individual has had T2DM for
                      # at least nine years. If they have had T2DM for a full 9 years,
                      # they now have the option of revisiting their health-related
                      # behaviours.
                      , ifelse(j$Duration_1 == 9, 
                               LTPA_Annual_Cost,
                               # If not, they must have had T2DM for a decade or more. 
                               # The same logic applies as in the years following the
                               # initial diagnosis, except the choice to which they 
                               # are required to adhere is not the one they made upon 
                               # diagnosis, but the one they made between Years 9 
                               # and 10.
                               ifelse(j$Per_2_Adherence_LTPA == "Adherent"
                                      , LTPA_Annual_Cost
                                      , Inf)))); j} 
    # This follows the same logic as above, with one critical exception: here the
    # individual has elected NOT to adhere to LTPA, so they shouldn't pay any costs
    # for engaging in regular physical activity. As above, if they opt to make this
    # decision after having originally committed to remain physically active, they pay
    # an arbitrarily large penalty.
    else
    {j$LTPA_Cost <- 
      ifelse((j$Duration_1 == 0 | j$Health_Stat_1 == "Healthy"), 
             0,
             ifelse(j$Duration_1 < 9,
                    ifelse(j$Per_1_Adherence_LTPA == "Adherent",
                           Inf,
                           0),
                    ifelse(j$Duration_1 == 9,
                           0,
                           ifelse(j$Per_2_Adherence_LTPA == "Adherent",
                                  Inf,
                                  0)))); j}
    , seq_along( along.with = RWD.03 )
    , RWD.03
    , SIMPLIFY = FALSE )

# Now consider dietary expenditures. Begin by initializing the rows that will hold the 
# data.
RWD.03 <- lapply(RWD.03, FUN = function(x) {x <- x %>% 
  mutate("Diet_Cost" = NA); x})

# 090818: In this new version of the script, dietary costs depend not only on gender,
# but also upon ethnicity and age. Here, we draw upon the tables compiled as part of
# the dietary costing exercise, extracting specific rows relevant to the decision-
# maker - specifically, those relating to their gender and ethnicity. Because caloric
# needs depend also upon age, costs are presented as a table, the values of which we
# will draw upon below.
Diet_Annual_Cost.H <- 
  Diet_Cost.H.01_K[Diet_Cost.H.01_K$Gender==GEN_K[[Gender.U]] & 
                     Diet_Cost.H.01_K$Ethnicity == ETH_K[[Ethnicity.U]] & 
                     Diet_Cost.H.01_K$Age < 95,]
Diet_Annual_Cost.UH <- 
  Diet_Cost.UH.01_K[Diet_Cost.UH.01_K$Gender==GEN_K[[Gender.U]] & 
                      Diet_Cost.UH.01_K$Ethnicity == ETH_K[[Ethnicity.U]] & 
                      Diet_Cost.UH.01_K$Age < 95,]

RWD.03 <-
  mapply( FUN = function( i, j )
    # As above, an odd number (i.e., x %% 2 != 0) denotes adherence, while an even
    # number denotes non-adherence. Thus, we begin with the case where the individual
    # chooses to adhere to healthy eating habits.
    if( ceiling(i/2) %% 2 != 0 ){
      j$Diet_Cost <- 
        ifelse((j$Duration_1 == 0 | j$Health_Stat_1 == "Healthy"), 
               Diet_Annual_Cost.H[findInterval(j$Age_1, Diet_Annual_Cost.H$Age, 
                                               all.inside = T),5],
               ifelse(j$Duration_1 < 9,
                      # As above, if the person originally committed to healthy eating,
                      # then they are charged the price of adhering to this health-
                      # related behaviour. If, however, they did NOT previously commit
                      # to this, they are forced to stick with this (which is enforced
                      # through the application of an arbitrarily large penalty).
                      ifelse(j$Per_1_Adherence_Diet == "Adherent",
                             Diet_Annual_Cost.H[findInterval(j$Age_1, 
                                                             Diet_Annual_Cost.H$Age, 
                                                             all.inside = T),5], 
                             Inf),
                      ifelse(j$Duration_1 == 9,
                             Diet_Annual_Cost.H[findInterval(j$Age_1, 
                                                             Diet_Annual_Cost.H$Age, 
                                                             all.inside = T),5],
                             ifelse(j$Per_2_Adherence_Diet == "Adherent",
                                    Diet_Annual_Cost.H[findInterval(j$Age_1, 
                                                                    Diet_Annual_Cost.H$Age,
                                                                    all.inside = T),5],
                                    Inf)))); j} 
    else
    {j$Diet_Cost <- 
      # 060318: The dietary cost included here originally was Diet_Annual_Cost.U.FAH,
      # which I don't believe was correct (since if ceiling(i/2) %% 2 == 0), this
      # denotes non-adherence.
      ifelse((j$Duration_1 == 0 | j$Health_Stat_1 == "Healthy"), 
             Diet_Annual_Cost.UH[findInterval(j$Age_1, Diet_Annual_Cost.UH$Age, 
                                              all.inside = T),5],
             ifelse(j$Duration_1 < 9,
                    ifelse(j$Per_1_Adherence_Diet == "Adherent",
                           Inf,
                           Diet_Annual_Cost.UH[findInterval(j$Age_1, 
                                                            Diet_Annual_Cost.UH$Age, 
                                                            all.inside = T),5]),
                    ifelse(j$Duration_1 == 9,
                           Diet_Annual_Cost.UH[findInterval(j$Age_1, 
                                                            Diet_Annual_Cost.UH$Age, 
                                                            all.inside = T),5],
                           ifelse(j$Per_2_Adherence_Diet == "Adherent",
                                  Inf,
                                  Diet_Annual_Cost.UH[findInterval(j$Age_1,
                                                                   Diet_Annual_Cost.UH$Age, 
                                                                   all.inside = T),5])))); j}
    , seq_along( along.with = RWD.03 )
    , RWD.03
    , SIMPLIFY = FALSE )

# Finally, we consider pharmaceutical expenditures. As before, initialize the rows
# that will hold the data.
RWD.03 <- lapply(RWD.03, FUN = function(x) {x <- x %>% 
  mutate("Rx_Price" = NA); x})

RWD.03 <-
  mapply( FUN = function( i, j )
    # Observe that the number to which the CEILING function is applied keeps increasing.
    # This is deliberate, and reflects the fact that the block size (see above) itself
    # keeps growing, from 1, to 2, to 4.
    # Now, if the person chooses to adhere ...
    if( ceiling(i/4) %% 2 != 0 ){
      j$Rx_Price <- 
        # If the individual is healthy, there are no drug costs.
        ifelse(j$Health_Stat_1 == "Healthy", 
               0,
               # If, however, they have T2DM, then: 
               # 1) If just diagnosed, they pay costs reflecting expenses for people
               # at a relatively early stage of T2DM progression (RxExp.Early.Annual_K).
               ifelse(j$Duration_1 == 0, 
                      RxExp.Early.Annual_K,
                      # 2) In the years immediately following diagnosis but prior to
                      # Year 10, having committed to adhering with pharmacotherapy,
                      # they must remain faithful to this choice. However, if they
                      # did not previously commit to medication adherence, they cannot
                      # change their mind now, and, indeed, are penalized if they try.
                      ifelse(j$Duration_1 < 9,
                             ifelse(j$Per_1_Adherence_Medication_1 == "Adherent",
                                    RxExp.Early.Annual_K,
                                    Inf),
                             # 3) In Year 9, the individual can revisit their choice.
                             ifelse(j$Duration_1 == 9, 
                                    ifelse(j$WC_P1.02 == "Yes",
                                           RxExp.Adv.WC.Annual_K,
                                           RxExp.Adv.PC.Annual_K),
                                    # 4) In Year 10 and every year thereafter, the health-
                                    # related decisions the individual made during Year 9
                                    # become part of the individual's history of T2DM
                                    # management, in addition to the choices made upon
                                    # diagnosis with the condition. If the individual decided
                                    # to adhere to pharmacotherapy in Year 10, then they are,
                                    # as before, required to adhere to this choice. However,
                                    # the costs of adherence now depend on whether their
                                    # condition was "well-controlled" during the first decade.
                                    # If so, they pay one (lower) cost, whereas they otherwise
                                    # pay another (higher) cost.
                                    ifelse(j$Per_2_Adherence_Medication_1 == "Adherent",
                                           ifelse(j$WC_P1.02 == "Yes", 
                                                  RxExp.Adv.WC.Annual_K,
                                                  RxExp.Adv.PC.Annual_K),
                                           Inf))))); j} 
    # If the person chooses not to adhere ...
    else
    {j$Rx_Price <- 
      ifelse(j$Health_Stat_1 == "Healthy", 
             Inf,
             ifelse(j$Duration_1 == 0, 
                    0,
                    ifelse(j$Duration_1 < 9,
                           ifelse(j$Per_1_Adherence_Medication_1 == "Adherent",
                                  Inf,
                                  0),
                           ifelse(j$Duration_1 == 9, 
                                  ifelse(j$WC_P1.02 == "Yes",
                                         0,
                                         0),
                                  ifelse(j$Per_2_Adherence_Medication_1 == "Adherent",
                                         ifelse(j$WC_P1.02 == "Yes", 
                                                Inf,
                                                Inf),
                                         # 102318: Previously, this next line read
                                         # "Inf". It should, however, be 0, because
                                         # people who chose not to adhere with meds
                                         # in Year 9 should not be penalized for
                                         # doing so now. The effect of this error
                                         # appears to have been that utility was set
                                         # to Penalty_K for all states in which the
                                         # individual decided not to adhere, from
                                         # Year 10 on. As such, there was no clear
                                         # incentive to choose the course of action
                                         # to which they had committed.
                                         0))))); j}
    , seq_along( along.with = RWD.03 )
    , RWD.03
    , SIMPLIFY = FALSE )

# Now we calculate out-of-pocket expenditures on pharmaceuticals.
RWD.03 <- lapply(RWD.03, FUN = function(x) {x <- x %>% 
  mutate("Rx_Cost" = NA); x})

RWD.03 <- lapply(RWD.03, FUN = function(x) 
  # If the price of pharmaceuticals is less than Inf (implying that the individual has
  # made choices that are "consistent", in the sense that they do not deviate from prior
  # commitments), then out-of-pocket costs are simply the lesser of their Pharmacare
  # deductible and their actual medication costs.
{x$Rx_Cost <- ifelse(x$Rx_Price < Inf
                     , pmin(x$Rx_Deductible, x$Rx_Price)
                     # If the individual's choices were, however, "inconsistent" 
                     # (i.e., they violated earlier commitments), then their OOP drug
                     # costs are set to Inf (carrying along the information that
                     # such choices are not permissible).
                     , Inf); x})

# This just fills in rows in this section (i.e., relating to household expenditures)
# that consist of NAs - which, in principle - should affect exactly one row in each
# list element (i.e., the death state).
RWD.03 <- RWD.03 %>% 
  lapply(., FUN = function(x) 
  {x[, c("Diet_Cost", "LTPA_Cost", "Rx_Price", "Rx_Cost")][is.na(x[, c("Diet_Cost", "LTPA_Cost", "Rx_Price", "Rx_Cost")])] <- 0 ; x})

# * Consumption ====

RWD.03 <- lapply(RWD.03, FUN = function(x) {x <- x %>% 
  mutate("Consumption" = NA); x})

# NOTE 091418: The income people actually have available to allocate to consumption consists of
# taxable and non-taxable income. Previously, the expression below calculated consumption as
# the difference between total income (if retired) or adjusted income (if still in the labour
# force), and total household expenditures. This, however, caused issues in tabulating
# consumption among labour force participants, because
# 1) For individuals who are employed, adjusted income does not include amounts stocked away in
# their retirement fund (if younger than 65), or withdrawn from said fund (if 65+). As such, it
# overstates and understates the resources available for consumption by the former and latter
# groups, respectively.
# 2) For individuals who are unemployed but younger than 65, there is no adjusted income,
# because such income ultimately originates from labour market earnings, which by definition
# they do not earn; rather, their income is derived from EIA, which is classified as non-
# taxable income.
# 3) Individuals who are unemployed but 65 years or older are assumed to transition
# automatically to retirement. They, too, receive no adjusted income. They can, however (and
# must) withdraw money from their retirement fund, which is classified as taxable income -
# as are CPP and OAS benefits. They may also be eligible for GIS benefits, which are classified
# as non-taxable income. As such, resources available for consumption consist of both income
# categories.
# NOTE 042019: Divided consumption by its price (which in this model was assumed).
RWD.03 <-
  mapply( FUN = function( i, j )
  {j$Consumption <- (j$Taxable_Income + j$Non_Taxable_Income - j$Tax_Paid - j$CPP_EI_Contribs - 
                       j$LTPA_Cost - j$Diet_Cost - j$Rx_Cost)/P_c_K; j}
  , seq_along( along.with = RWD.03 )
  , RWD.03
  , SIMPLIFY = FALSE )

# * Allocation of Time ====

# Set up columns for time dedicated to various activities.
RWD.03 <- lapply(RWD.03, FUN = function(x) {x <- x %>% 
  mutate("Time_Work" = NA) %>%
  mutate("Time_LTPA" = NA) %>%
  mutate("Time_Diet" = NA) %>%
  mutate("Time_Sick" = NA); x})

# * * Allocate time for work and illness ####
RWD.03 <- RWD.03 %>%
  mapply( FUN = function( i, j )
    # If the individual retires, by definition they are not working, and time
    # allocated to work is therefore assumed to be 0. Employment state is not really
    # pertinent to this calculation, but I apply an ifelse statement just to keep my
    # options open.
    if( ceiling(i/8) %% 3 == 0 )
    {j$Time_Work <- ifelse(j$Emp_Stat_1 == "Employed"
                           , 0
                           , 0); j} 
    else
      # If the individual chooses to work FT ...
    {if( ceiling(i/8) %% 2 != 0 )
      # Note that we do NOT subtract presenteeism from these calculations,
      # as we did when we were endeavouring to calculate compensation. The
      # reason for this is that people are still, strictly speaking, working - 
      # albeit at less than full capacity.
    {j$Time_Work <- ifelse(j$Emp_Stat_1 == "Employed"
                           , Annual_Work_Hrs.FT_K - j$Adj_Abs_Hours
                           # If the individual is NOT employed, but still seeking
                           # FT hours, they allocate a pre-determined amount of
                           # time to search.
                           , Time_Req_Search_K); j}
      else
        # "Adj_Abs_Hours" should already reflect the individual's PT/FT status, so
        # we can use the value in the same column, without any need for 
        # modification.
      {j$Time_Work <- ifelse(j$Emp_Stat_1 == "Employed"
                             , Annual_Work_Hrs.PT_K - j$Adj_Abs_Hours
                             , Time_Req_Search_K); j}
    }
    , seq_along( along.with = . )
    , .
    , SIMPLIFY = FALSE ) %>%
  mapply( FUN = function( i, j )
    if( ceiling(i/8) %% 3 == 0 )
      # If a person is retiring, they bear the entire burden associated with 
      # absenteeism and presenteeism, beginning in the same year the decision is made.
    {j$Time_Sick <- 
      ifelse(j$Emp_Stat_1 == "Employed"
             , (Time_Available_K/Annual_Work_Hrs.FT_K)*(j$Annual_Absences*8 + j$Pres*12)
             # Again, there is no obvious need to account for employment status if
             # the individual is retiring, but we include this additional line just
             # in case we require it for some reason in the future.
             , (Time_Available_K/Annual_Work_Hrs.FT_K)*(j$Annual_Absences*8 + j$Pres*12)); j} 
    else
    {j$Time_Sick <- 
      ifelse(j$Emp_Stat_1 == "Employed"
             # This ugly statement says the following. We estimate that the
             # individual experiences a certain total number of days during which
             # they are "bed-ridden" (the generalization of absenteeism) or
             # functioning at reduced capacity (the generalization of presenteeism).
             # Some fall on work days, while others occur at other times (e.g.,
             # weekends, overnight, etc.). Sick time for which the individual is
             # compensated does not count as sick time, but rather as work time. In
             # addition, workplace presenteeism is never counted as sick time,
             # because the intervals over which the individual is experiencing
             # symptoms are actually allocated to work. Therefore, we calculate
             # sick time by applying the following methodology:
             # 1) Calculate total annual sick time by multiplying the inverse of the
             # proportion of available time a FT employee spends at work by the sum
             # of the hours associated with absenteeism and presenteeism (as derived
             # from the two studies cited in the thesis).
             # 2) Subtract all workplace presenteeism from this total. We are only
             # interested in impaired productivity experienced outside of work.
             # 3) Subtract all absences for which the individual was remunerated.
             # This is assumed to be the minimum of actual absences (if absences <
             # PSDs) and paid sick days (if absences > PSDs).
             # 
             # Note that the total amount of time the individual spends ill is not
             # influenced by how many hours they work (i.e., PT vs. FT); rather,
             # their employment status affects only how much of this time occurs when
             # they are (or should be) working.
             , (Time_Available_K/Annual_Work_Hrs.FT_K)*(j$Annual_Absences*8 + j$Pres*12) - 
               j$Pres*12/ceiling(i/8) - 
               # This expression preserves consistency with earlier calculations by
               # indicating that the number of workplace absences that are classified
               # as work time is the lesser of actual absences and PSDs. If the
               # individual does not require all their sick days, only what they used
               # is considered work. If, however, they exceeded their allotment of
               # PSDs, the entire allotment - but only that - is considered work;
               # everything in addition to this would be regarded as sick time, for
               # which the individual is not compensated.
               pmin(j$Annual_Absences/ceiling(i/8),Paid_SDs.FT_K/ceiling(i/8))*8
             # If the individual is unemployed, then they bear the entire burden of
             # absenteeism/presenteeism themselves, much the same as if they had
             # retired.
             , (Time_Available_K/Annual_Work_Hrs.FT_K)*(j$Annual_Absences*8 + j$Pres*12)); j}
    , seq_along( along.with = . )
    , .
    , SIMPLIFY = FALSE )

# * * Allocate time for LTPA and diet ####

# Note that we do not bother assigning penalties for individuals who fail to behave in
# a manner consistent with earlier commitments regarding health-related behaviour (i.e.,
# for individuals with T2DM). We do not due this largely because it is not necessary -
# imposing income penalties should provide sufficient incentive to ensure consistency is
# maintained.
RWD.03 <- RWD.03 %>%
  mapply( FUN = function( i, j )
    if( ceiling(i/1) %% 2 != 0 ){
      j$Time_LTPA <- Time_Req_LTPA_K; j
    } else{j$Time_LTPA <- 0; j}
    , seq_along( along.with = . )
    , .
    , SIMPLIFY = FALSE ) %>%
  mapply( FUN = function( i, j )
    if( ceiling(i/2) %% 2 != 0 ){
      j$Time_Diet <- Time_Req_Diet_H_K; j
    } else{j$Time_Diet <- Time_Req_Diet_UH_K; j}
    , seq_along( along.with = . )
    , .
    , SIMPLIFY = FALSE )

# * Calculate leisure

RWD.03 <- lapply(RWD.03, FUN = function(x) {x <- x %>% 
  mutate("Leisure" = NA); x})

RWD.03 <- lapply(RWD.03, FUN = function(x)
{x$Leisure <- Time_Available_K - x$Time_Work - x$Time_LTPA - x$Time_Diet - 
  x$Time_Sick; x})

# * Calculate Utility ====

RWD.03 <- lapply(RWD.03, FUN = function(x) {x <- x %>% 
  mutate("Utility" = NA); x})

# NOTE 042019: Modified utility function. We now use a more general functional form
# (constant elasticity of substitution utility function) than in previous versions
# of the model (Cobb-Douglas utility function). The essential parameter
# underpinning the CES utility function is Rho, which is related to the elasticity
# of substitution. If Rho is equal to one, the CES utility function collapses to
# the CD utility function, as before; otherwise, however, the resultant equation is
# more complicated, but is also able to reflect consumer preferences more flexibly.
# NOTE 091418: The condition marked by the asterix below is new. This is intended
# to address the issue that individuals younger than 65 who opt to retire have
# negative consumption levels, because they do not earn income but are required to
# save for retirement. Such individuals are assigned a utility level of Penalty_K, 
# which should always be lower than the available alternative (i.e., remaining in
# the labour force).

# NOTE 042319: Caught and corrected error concerning calculation of utility. Originally,
# I used an IFELSE statement, where the choice of utility function (i.e., Cobb-Douglas or
# CES) hinged on the value of Rho_K. I found that this resulted in a single value being
# assigned to all rows, excepting those assigned 0 or -1. I concluded that this was likely
# attributable to the vectorized nature of the IFELSE function. In particular, I theorize
# that for some reason, the first value generated through the application of the 
# conditional logic was applied to all rows, some of which were subsequently overwritten
# by 0 and -1. The customized function below seems to circumvent this issue by applying
# the if-else logic once (noting that the utility function should never vary within a 
# loop), and then applying the resultant function consistently.

# NOTE 052519: Caught and corrected another error relating to calculation of
# utility, revolving around the expression below. The nested IF statements
# did not function as intended, with the effect that (possibly among other
# things) penalties were applied to actions which should have been permissible.
# I ultimately addressed this by applying the penalty associated with negative
# levels of consumption and/or leisure right at the end (of this sub-section).
# if(Rho_K==0) 
# {Utility_Fxn <- function(x, y, z) 
#   {if(y <= 0 | z <= 0)
#     {Penalty_K
#   } else {
#       x*(y^Alpha_K)*(z^(1-Alpha_K))}}
# } else {
#   Utility_Fxn <- function(x, y, z) 
#   {if(y <= 0 | z <= 0)
#     {Penalty_K
#   } else {
#     x*(Alpha_K*y^Rho_K + (1-Alpha_K)*z^Rho_K)^(1/Rho_K)}}}

if(Rho_K==0) 
{Utility_Fxn <- function(x, y, z) x*(y^Alpha_K)*(z^(1-Alpha_K))
} else {
  Utility_Fxn <- function(x, y, z) 
    x*(Alpha_K*y^Rho_K + (1-Alpha_K)*z^Rho_K)^(1/Rho_K)}

RWD.03 <- RWD.03 %>%
  mapply( FUN = function(i, j)
    # If the person has chosen to retire ...
    if( ceiling(i/8) %% 3 == 0 ){
      j$Utility <- ifelse(j$Health_Stat_1 == "Dead"
                          , 0
                          # *
                          , ifelse(j$Age_1 < 65
                                   , Penalty_K
                                   , Utility_Fxn(j$Adjusted_Utility, j$Consumption, j$Leisure))); j
    } else {
      j$Utility <- ifelse(j$Health_Stat_1 == "Dead"
                          # Why the penalty? The intent here is to corral the dead
                          # into "choosing" retirement, primarily because it makes
                          # little sense to categorize them as being employed. The
                          # difference in rewards, however modest, should be
                          # sufficient to nudge people into the decisions we want
                          # them to make.
                          , Penalty_K
                          , Utility_Fxn(j$Adjusted_Utility, j$Consumption, j$Leisure)); j}
    , seq_along(along.with = .)
    , .
    , SIMPLIFY = FALSE) %>%
  mapply( FUN = function(i, j)
    # If the individual has chosen to work part time ...
    if( ceiling(i/8) %% 3 == 2 )
    {j$Utility <- ifelse(j$Health_Stat_1 == "Dead",
                         j$Utility,          
                         ifelse(j$Emp_Stat_1 == "Unemployed" 
                                # What is the meaning of this? Again, the intent is to force
                                # decision-makers to engage in actions that are not
                                # unreasonable, given their current state, and that display a
                                # degree of consistency. Here, we want all job seekers to 
                                # elect to choose "full-time hours" despite the fact that the
                                # distinction has virtually no significance to those who are
                                # unemployed (in part, because people who had been working FT
                                # or PT are required to invest the same amount of effort in
                                # search).
                                , Penalty_K
                                , j$Utility)); j}
    # NOTE 091518: Another possibility to consider is that the individual has elected to work 
    # full-time. However, not everyone can make this choice, since we assume individuals 65 
    # or older cannot re-enter the labour force after having previously left it. Interestingly, 
    # in prior versions of this code, the reward for working FT actually surpassed the reward
    # for remaining retired, likely reflecting the fact that private retirement savings are 
    # assumed to have been depleted by the time individuals reach their average life expectancy
    # (which, for men, is approximately 80).
    else
    {if( ceiling(i/8) %% 3 == 1 )
    {j$Utility <- ifelse(j$Health_Stat_1 == "Dead",
                         # If the individual is dead, they have already been assigned a 
                         # utility value above - there's no reason to overwrite it.
                         j$Utility,
                         ifelse(j$Emp_Stat_1 == "Unemployed",
                                ifelse(j$Age_1 > 64,
                                       # If the individual is unemployed and over the age
                                       # of retirement, they cannot re-enter the labour
                                       # force; if they try, we assign them a penalty.
                                       Penalty_K,
                                       Utility_Fxn(j$Adjusted_Utility, j$Consumption, j$Leisure)),
                                # If the individual is employed and opts to retire, they
                                # begin earning retirement income immediately, and 
                                # transition to unemployment the instant before the next
                                # decision epoch.
                                Utility_Fxn(j$Adjusted_Utility, j$Consumption, j$Leisure))); j}
      # NOTE 091618: The last possibility, of course, is that the individual opts to retire.
      else
      {j$Utility <- ifelse(j$Health_Stat_1 == "Dead",
                           j$Utility,
                           ifelse(j$Emp_Stat_1 == "Unemployed",
                                  ifelse(j$Age_1 > 64,
                                         # Here things are a little different. People above
                                         # the standard age of retirement who are 
                                         # unemployed or who choose to retire begin earning
                                         # retirement income immediately (and transition to
                                         # unemployment if they did not already occupy said
                                         # state).
                                         Utility_Fxn(j$Adjusted_Utility, j$Consumption, j$Leisure),
                                         # Meanwhile, people below the standard age of 
                                         # retirement cannot (yet) retire, so if they
                                         # attempt to do so, they receive a penalty.
                                         Penalty_K),
                                  # The individual might also be employed.
                                  ifelse(j$Age_1 > 64,
                                         # If they are employed and at or above the
                                         # standard age of retirement, they are eligible
                                         # to retire during any decision epoch, enabling
                                         # them to begin earning their public pension.
                                         Utility_Fxn(j$Adjusted_Utility, j$Consumption, j$Leisure),
                                         # If they are below the standard age of
                                         # retirement, they cannot yet retire, and
                                         # receive a penalty if they attempt to do so.
                                         Penalty_K))); j}}
    , seq_along(along.with = .)
    , .
    , SIMPLIFY = FALSE)

RWD.03 <- RWD.03 %>%
  mapply( FUN = function(i, j)
  {j$Utility <- ifelse((is.nan(j$Utility) == T | j$Consumption <= 0 | 
                         j$Leisure <= 0) & j$Health_Stat_1 != "Dead"
                       , Penalty_K
                       # *
                       , j$Utility); j}
  , seq_along(along.with = .)
  , .
  , SIMPLIFY = FALSE)

# * Worked Examples ====

# NOTE 052519: The function of this section is to carry out the calculations
# required to produce three worked examples, which are as follows (the
# content in parentheses indicate the specific combinations of gender and 
# ethnicity incorporated into the body of the thesis itself):
# 1) A retired 93 year-old recently diagnosed with T2DM (Asian-Indian Female
# with a university certificate, diploma or degree at bachelor level or above)
# 2) An employed 60 year-old who has had T2DM for 9 years (Caucasian Male
# with a high school diploma or equivalent)
# 3) An unemployed 30 year-old who is now healthy (Afro-Caribbean Female with
# a postsecondary certificate or diploma below bachelor level)

# * * Setup ####

# RWD.04 <- RWD.03 %>%
#   mapply( FUN = function(i, j)
#     if( ceiling(i/8) %% 3 == 0 )
#     {j <- j %>% subset(., select = c(1:3, 5, 14, 7, 6, 8:13, 15, 
#                                      26:29, 18:19, 33:34, 36:37, 39:46)); j}
#     else 
#     {j <- j %>% subset(., select = c(1:3, 5, 14, 7, 6, 8:13, 15,
#                                      25:30, 34:35, 37:38, 40:47)); j}
#     , seq_along(along.with = .)
#     , .
#     , SIMPLIFY = FALSE)
# 
# # * * Year 0 ####
# 
# RWD_WE_T0 <- vector("list", 3); TP_WE_T0 <- vector("list", 3)
# RWD_WE_T0[[1]] <- RWD.04 %>%
#   lapply(., FUN = function(x) {x <- x %>% 
#     subset(., Age_1 == 93 & Duration_1 == 0 & 
#              Emp_Stat_1 == "Unemployed"); x})
# TP_WE_T0[[1]] <- TP.03 %>% lapply(., FUN = function(x){x <- x %>% 
#   subset(., ID.x == 228347); x})
# RWD_WE_T0[[2]] <- RWD.04 %>%
#   lapply(., FUN = function(x) {x <- x %>% 
#     subset(., Age_1 == 60 & Duration_1 == 9 & 
#              Emp_Stat_1 == "Employed" & 
#              Per_1_Adherence_LTPA_1 == "Adherent" &
#              Per_1_Adherence_Diet_1 == "Non-Adherent" & 
#              Per_1_Adherence_Medication_1 == "Non-Adherent"); x})
# TP_WE_T0[[2]] <- TP.03 %>% lapply(., FUN = function(x){x <- x %>% 
#   subset(., ID.x == 73321); x})
# RWD_WE_T0[[3]] <- RWD.04 %>%
#   lapply(., FUN = function(x) {x <- x %>% 
#     subset(., Age_1 == 30 & Health_Stat_1 == "Healthy" & 
#              Emp_Stat_1 == "Unemployed"); x})
# TP_WE_T0[[3]] <- TP.03 %>% lapply(., FUN = function(x){x <- x %>% 
#   subset(., ID.x == 76); x})
# 
# for (i in 1:3) 
# {RWD_WE_T0[i] %<>% lapply(., FUN = function(x) do.call("rbind", x))}
# 
# # 061219: Illustrating the consequences of poor health decisions
# RWD_WE_T0_JW <- vector("list", 2); TP_WE_T0_JW <- vector("list", 2)
# 
# RWD_WE_T0_JW[[1]] <- RWD.04 %>%
#   lapply(., FUN = function(x) {x <- x %>% 
#     subset(., Age_1 == 70 & Duration_1 == 19 & 
#              Emp_Stat_1 == "Unemployed" & 
#              Per_1_Adherence_LTPA_1 == "Adherent" &
#              Per_1_Adherence_Diet_1 == "Non-Adherent" & 
#              Per_1_Adherence_Medication_1 == "Non-Adherent" & 
#              Per_2_Adherence_LTPA_1 == "Non-Adherent" &
#              Per_2_Adherence_Diet_1 == "Adherent" & 
#              Per_2_Adherence_Medication_1 == "Adherent"); x})
# TP_WE_T0_JW[[1]] <- TP.03 %>% lapply(., FUN = function(x){x <- x %>% 
#   subset(., ID.x == 188044); x})
# 
# RWD_WE_T0_JW[[2]] <- RWD.04 %>%
#   lapply(., FUN = function(x) {x <- x %>% 
#     subset(., Age_1 == 70 & Duration_1 == 19 & 
#              Emp_Stat_1 == "Unemployed" & 
#              Per_1_Adherence_LTPA_1 == "Adherent" &
#              Per_1_Adherence_Diet_1 == "Non-Adherent" & 
#              Per_1_Adherence_Medication_1 == "Non-Adherent" & 
#              Per_2_Adherence_LTPA_1 == "Non-Adherent" &
#              Per_2_Adherence_Diet_1 == "Non-Adherent" & 
#              Per_2_Adherence_Medication_1 == "Non-Adherent"); x})
# TP_WE_T0_JW[[2]] <- TP.03 %>% lapply(., FUN = function(x){x <- x %>% 
#   subset(., ID.x == 188049); x})
# 
# for (i in 1:2) 
# {RWD_WE_T0_JW[i] %<>% lapply(., FUN = function(x) do.call("rbind", x))}

# * * Year 1 ####

# RWD_WE_T1 <- vector("list", 3)
# RWD_WE_T1[[1]] <- RWD.04 %>%
#   lapply(., FUN = function(x) {x <- x %>% 
#     subset(., (Age_1 == 94 & Duration_1 == 1) | ID_States == 228357); x})
# RWD_WE_T1[[2]] <- RWD.04 %>%
#   lapply(., FUN = function(x) {x <- x %>% 
#     subset(., (Age_1 == 61 & Duration_1 == 10 & 
#                  Per_1_Adherence_LTPA_1 == "Adherent" &
#                  Per_1_Adherence_Diet_1 == "Non-Adherent" & 
#                  Per_1_Adherence_Medication_1 == "Non-Adherent") | 
#              ID_States == 228357); x})
# RWD_WE_T1[[3]] <- RWD.04 %>%
#   lapply(., FUN = function(x) {x <- x %>% 
#     subset(., (Age_1 == 31 & (is.na(Duration_1) == T | Duration_1 == 0)) | 
#              ID_States == 228357); x})
# 
# WE_T1 <- vector("list", 3)
# for (a in 1:3){
#   WE_T1[[a]] <- mapply( FUN = function( i, j ) 
#     merge(x = i, y = j, by.x = "ID_States", by.y = "ID.y", all = TRUE)
#     , RWD_WE_T1[[a]]
#     , TP_WE_T0[[a]]
#     , SIMPLIFY = FALSE )
# }
# 
# save(RWD_WE_T0, file = paste(RDA, "/Worked_Examples/","WE_T0", LM, "_",
#                              Gender.U, Ethnicity.U, Education.U, ".Rda", sep = ""))
# save(WE_T1, file = paste(RDA, "/Worked_Examples/", "WE_T1", LM, "_",
#                          Gender.U, Ethnicity.U, Education.U,".Rda", sep = ""))
# 
# # 061219: Illustrating the consequences of poor health decisions
# RWD_WE_T1_JW <- vector("list", 2)
# RWD_WE_T1_JW[[1]] <- RWD.04 %>%
#   lapply(., FUN = function(x) {x <- x %>% 
#     subset(., (Age_1 == 71 & Duration_1 == 20 & 
#                  Emp_Stat_1 == "Unemployed" & 
#                  Per_1_Adherence_LTPA_1 == "Adherent" &
#                  Per_1_Adherence_Diet_1 == "Non-Adherent" & 
#                  Per_1_Adherence_Medication_1 == "Non-Adherent" & 
#                  Per_2_Adherence_LTPA_1 == "Non-Adherent" &
#                  Per_2_Adherence_Diet_1 == "Adherent" & 
#                  Per_2_Adherence_Medication_1 == "Adherent") | 
#              ID_States == 228357); x})
# RWD_WE_T1_JW[[2]] <- RWD.04 %>%
#   lapply(., FUN = function(x) {x <- x %>% 
#     subset(., (Age_1 == 71 & Duration_1 == 20 & 
#                  Emp_Stat_1 == "Unemployed" & 
#                  Per_1_Adherence_LTPA_1 == "Adherent" &
#                  Per_1_Adherence_Diet_1 == "Non-Adherent" & 
#                  Per_1_Adherence_Medication_1 == "Non-Adherent" & 
#                  Per_2_Adherence_LTPA_1 == "Non-Adherent" &
#                  Per_2_Adherence_Diet_1 == "Non-Adherent" & 
#                  Per_2_Adherence_Medication_1 == "Non-Adherent") | 
#              ID_States == 228357); x})
# 
# WE_T1_JW <- vector("list", 2)
# for (a in 1:2){
#   WE_T1_JW[[a]] <- mapply( FUN = function( i, j ) 
#     merge(x = i, y = j, by.x = "ID_States", by.y = "ID.y", all = TRUE)
#     , RWD_WE_T1_JW[[a]]
#     , TP_WE_T0_JW[[a]]
#     , SIMPLIFY = FALSE )
# }
# 
# save(RWD_WE_T0_JW, file = paste(RDA, "/Worked_Examples/","WE_T0_JW", LM, "_",
#                                 Gender.U, Ethnicity.U, Education.U, ".Rda", sep = ""))
# save(WE_T1_JW, file = paste(RDA, "/Worked_Examples/", "WE_T1_JW", LM, "_",
#                             Gender.U, Ethnicity.U, Education.U,".Rda", sep = ""))

# * Finalize reward matrix and save results ====

RWD.03 <- RWD.03 %>% 
  lapply(., FUN = function(x) 
  {x <- x %>% subset(., select = c("ID_States", "Utility")); x})

# SUBSET doesn't seem to function here, even though it has been applied successfully
# elsewhere in the code.
# UPDATE: It doesn't work properly because RWD.03 is a list, not a dataframe.
# RWD.03 <- subset(RWD.03, subset = RWD.03[,], select = c(ID, Utility))