# Processing Results from Life Tables for Manitoba
# Created: 03/30/18; Last Modified: 04/20/19

# Introduction & Rationale ------------------------------------------------

# NOTE 042019: Modified detection of script location. For this to work, the user
# must set the default working directory under Global Options to the base
# directory of the folder in which the decision model is stored.
# 
# NOTE 111818: The code now detects the path where the script is located, eliminating 
# the need for the user to change this manually upon downloading the script for the 
# first time.
# 
# NOTE 090618: The code has been updated to draw upon the most recent life tables
# available from Statistics Canada (2014/16)
# 
# The purpose of this script is to apply the RVEST package to download, manipulate, and
# ultimately aggregate data from the life tables for Manitoba for 2014-2016. The ultimate
# goal is to develop the foundation for a vectorized approach to the construction of the
# transition probability and reward matrices that will ultimately be incorporated into the
# model. A vectorized approach is preferred because experience suggests such an approach
# can expedite operations by a factor of magnitude or more, as implied by the observation
# that an operation that took a half-hour or more (probably more) to perform using for-
# loops could be carried out in under a minute using alternative approaches.

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

# The next line of code clears the variable list, if required.
rm(list=ls())

library(rvest);library(dplyr); library(stringr); library(tidyr); library(httr)
library(gsubfn); library(magrittr); library(purrr); library(CANSIM2R); 
library(data.table); library(here)

Model_Path <- here()
# String for date last modified
LM <- "_042019"
# Filename for output
Filename <- "Life_Tables"
# Paths for output
RDA <- paste0(Model_Path, "/R_Objects/")

# Extract Life Tables -----------------------------------------------------

x.00 <- getCANSIM(13100114)
x.00 <- data.table(x.00)

# https://stackoverflow.com/questions/11940605/selecting-a-subset-of-columns-in-a-data-table
x.01 <- x.00[i=="Canada" & t=="2014/2016", ]

# https://stackoverflow.com/questions/44057930/access-hidden-attribute-in-r-data-frame
z.00 <- lapply(x.01, function(x.01) attributes(x.01)$label)
z.01 <- unlist(z.00)
colnames(x.01) <- z.01

x.02 <- x.01[, -grep("Both sexes|Margin of error", colnames(x.01)), with = FALSE]
x.03 <- x.02[, grep("qx|(ex)", colnames(x.02)), with = FALSE]

z.02 <- z.01[-grep("Both sexes|Margin of error", z.01)]
z.03 <- z.02[grep("qx|(ex)", z.02)]

# https://stackoverflow.com/questions/28653867/best-way-to-transpose-data-table
x.04 <- data.table::transpose(x.03)
x.04$V2 <- z.03
x.04 <- separate(data = x.04, col = V2, into = c("Age", "Gender", "Metric"), sep = "; ")
x.04$Age <- gsub("[^0-9]", "", x.04$Age)
x.04$Age <- as.numeric(x.04$Age)

x.05 <- data.table(x.04)
# https://stackoverflow.com/questions/14630335/how-to-select-r-data-table-rows-based-on-substring-match-a-la-sql-like
x.05.M.Pr <- x.05[Gender=="Males" & Metric %like% "Death probability", c(1:3)]
x.05.M.LE <- x.05[Gender=="Males" & Metric %like% "Life expectancy", c(1:3)]
x.05.F.Pr <- x.05[Gender=="Females" & Metric %like% "Death probability", c(1:3)]
x.05.F.LE <- x.05[Gender=="Females" & Metric %like% "Life expectancy", c(1:3)]

x.06.M <- left_join(x.05.M.Pr, x.05.M.LE, by = "Age")
x.06.F <- left_join(x.05.F.Pr, x.05.F.LE, by = "Age")

x.06.M <- x.06.M[, c("Gender.x", "Age", "V1.x", "V1.y")]
x.06.F <- x.06.F[, c("Gender.x", "Age", "V1.x", "V1.y")]
x.06.M <- x.06.M %>% rename("Gender" = "Gender.x", "Prob" = "V1.x", "LifeExp" = "V1.y")
x.06.F <- x.06.F %>% rename("Gender" = "Gender.x", "Prob" = "V1.x", "LifeExp" = "V1.y")

# https://www.statmethods.net/management/sorting.html
x.06.M <- x.06.M[order(x.06.M$Age),]; x.06.F <- x.06.F[order(x.06.F$Age),]
# This would also work: setorder(x.06, Age)

x.06.M$Gender <- "Male"; x.06.F$Gender <- "Female"

# 052518: Save life expectancies for each gender from the original life tables; we will
# require this information to calculate the duration over which the individual will
# withdraw money from their private pensions.
Life_Exp <- data.frame("Gender" = c("Male", "Female"), 
                       "Life_Exp" = c(x.06.M$LifeExp[1], x.06.F$LifeExp[1]))
save(Life_Exp, file = paste(RDA, "Life_Exp", LM, ".Rda", sep = ""))

# Set up table of transition probabilities --------------------------------

# Now we set up the table of transition probabilities. We start by explicitly listing all
# possible combinations of each individual factor in which we are interested.
Factors <- list(Age <- as.character(seq(1:111)-1), Gender <- c("Male","Female"), 
                Employment <- c("Employed", "Unemployed"), Ethnicity <- 
                  c("Caucasian","Afro-Caribbean","Asian-Indian"))
# We also produce labels for each column.
Headers2 <- c("Age", "Gender", "Employment", "Ethnicity")

# The remarkable function EXPAND.GRID now generates a dataframe, the rows of which include
# all possible combinations of the factors listed above.
TP <- expand.grid(Factors, stringsAsFactors = FALSE)
# Apply the column labels produced above.
colnames(TP) <- Headers2
# Insert the matrix of transition probabilities and the two dataframes with values from
# the life tables into a list. This positions us to apply the reduce function momentarily.
Tables <- list(TP, "Male" = x.06.M, "Female" = x.06.F)

# This is where the magic happens. We ultimately want a single column containing all
# probabilities of interest, which we can then reshape to produce the matrix of transition
# probabilities that will ultimately be used in the decision model. So we need some way of
# linking each unique set of factors to its associated probability. We do this by modifying
# an approach described here: https://stackoverflow.com/questions/15291613/merge-3-data-frames-by-column-names
# This approach relies on two central functions, namely MERGE, which actually carries out
# the linking process, and REDUCE, from the PURRR package, which applies the linking
# recursively so that we don't have to MERGE multiple times (at least, not explicitly - 
# note that the rationale for this approach is that we could in principle be linking an
# arbitrarily large number of tables). With respect to the MERGE, note that we need to
# include from the RHS table not only the field(s) we want to link, but also the fields we
# will be merging on. Also note that specifying "all.x" ensures that we include every row
# from the LHS table - EVEN those that do not have a value in the RHS table (which in the
# present context will be the case 99% of the time).
TP.01 <- Reduce(function(x,y) merge(x,y[,c("Age","Gender","Prob")], 
                                    by = c("Age", "Gender" = "Gender"), 
                                    all.x = TRUE), Tables, accumulate = FALSE)

# Unfortunately, there doesn't seem to be a straightforward way of merging from multiple
# tables into a single column, which is why we engineered a sequence of merges. This means,
# however, that each merged field ultimately ends up in its own column. Fortunately, by
# construction, there is only ever one probability value in any row of the table, with the
# remainder being NAs. We can take advantage of this to condense everything into a single
# column using the COALESCE function, which is part of the DPLYR package. Essentially,
# COALESCE merges multiple columns, retaining the first non-NA value in each row. This is
# perfect, since, as we noted, there is only ever one such value in each row of our table.
TP.01 <- TP.01 %>%
  mutate("Probability" = 1 - coalesce(Prob.x,Prob.y))
  
# Having merged the fields of interest into a single column, we can now eliminate the
  # former using the SUBSET function. We now have what we originally wanted, which is a
  # single column that contains all probabilities.
TP.01 <- subset(TP.01, select = -c(Prob.x, Prob.y))
TP.01$Age <- as.integer(TP.01$Age)
TP.01 %<>% arrange(., Age, Gender, Ethnicity, Employment)

# https://stackoverflow.com/questions/8345759/how-to-save-a-data-frame-in-r
save(TP.01, file = paste0(RDA, Filename, LM, ".Rda"))