# Script for Generation of Tables Incorporated Into the Thesis
# Created: 08/10/18; Last Modified: 05/06/19

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

# NOTE 042019: Changes in the RSTUDIOAPI package appear to have affected the
# getSourceEditorContext()$path operation, causing issues when referencing
# project components. Ultimately, I determined that by setting the default
# working directory to the DIABETES_MODEL subfolder, it is possible to
# eliminate a lot of unnecessary complexity, since the path to all project
# components can now defined relative to the output generated by the 
# HERE function. I will update the User Guide to communicate the need for
# this added step to prospective users. 
# 
# 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 081818: Much of this syntax is borrowed directly from the RWD&TP Matrices
# Construction script.
# 
# This script generates 3 key tables incorporated in Chapter 3 of the 
# dissertation, namely:
# 1) Estimates of employment income for healthy individuals
# 2) Estimates of retirement income
# 3) Estimates of the probabilities of transitions from normoglycemia/
# prediabetes to T2DM

# It should be noted that these tables do not replace any of the tabulations
# in the RWD&TP Matrices Construction script, since they are intended purely
# to facilitate reporting, not for use in the actual model.

# LEGEND
# A <- Net_Emp_Inc.FT.00
# A_B <- Savings.New
# A_B_C <- Savings.Accum
# A_B_C_D <- Withdrawals
# A_B_C_D_E <- Income_Stream
# A_B_C_D_E_F <- Net_Emp_Inc.FT.01
# A_B_C_D_E_F_G <- Adj_Pens_Earnings

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

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

# Call required libraries
library(dplyr); library(stringr); library(tidyr); library(httr); library(magrittr)
library(purrr); library(reshape2); library(readxl); library(tibble); library(here)
library(data.table); library(stringr); library(openxlsx)

Model_Path <- here()

# String for date last modified
LM <- "_042019"
# This is where the model is located, but we'll need to concatenate this with the
# date the model was last modified to complete the path.
# Filename for output
Filename <- "RPT_Tables"
# Paths for RDA file outputs
RDA <- paste0(Model_Path, "/R_Objects/")
# Paths for Excel table outputs
Excel <- paste0(Model_Path, "/Excel_Tables/")

# We need life expectancy information, because this determines how long
# peoples' private savings will last.
load(paste0(RDA, "Life_Exp", LM, ".Rda"))
load(paste0(RDA, "NHS_2011_Calculations", LM, "_1", ".Rda"))
load(paste0(RDA, "NHS_2011_Calculations", LM, "_2", ".Rda"))
load(paste0(RDA, "DPoRT_Calculations", LM, ".Rda"))

# The following is a critical function we call upon when compiling the table
# using the DPoRT risk equations to estimate the one-year probabilities of
# developing T2DM. We discuss its application in further detail below.
# https://stackoverflow.com/questions/11693599/alternative-to-expand-grid-for-data-frames
expand.grid.df <- function(...) Reduce(function(...) merge(..., by=NULL), list(...))

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

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

# * Other miscellaneous parameters ====

# Monthly EIA for individuals with no dependents in the General Assistance category: 
# https://www.gov.mb.ca/fs/eia/eia_general.html
EIA_General = 12*806
# Monthly EIA for individuals with no dependents who is considered to have a disability
# (not used in the current version of the model):
# https://www.gov.mb.ca/fs/eia/eia_disability.html
EIA_Disability = 12*1027
Annual_Work_Hrs.FT <- 2000; Annual_Work_Hrs.PT <- 1000
# "ARPD" = "age-related productivity decline"
ARPD <- 0.025

# Table Generation: Employment Income --------------------------------------------------------

# * Full-Time Employment ====

Emp_Inc_RPT.FT.00 <- melt(Emp_Income, id.vars = c("Age_Cat", "Sex", "Highest_Degree"), 
                              measure.vars = c("Med_Emp_Inc_2017.FT"))

# The descriptors used to denote educational attainment are far too long. Here we
# shorten them significantly.
Emp_Inc_RPT.FT.00 <- Emp_Inc_RPT.FT.00 %>%
  mutate("Highest_Degree_1" = 
           case_when(Highest_Degree == "1. No certificate, diploma or degree" ~ "E1_None",
                     Highest_Degree == "2. High school diploma or equivalent" ~ "E2_HS",
                     Highest_Degree == "3. Postsecondary certificate or diploma below bachelor level" ~ "E3_No_Bach",
                     Highest_Degree == "4. University certificate, diploma or degree at bachelor level or above" ~ "E4_Bach_Plus"))

Emp_Inc_RPT.FT.01 <- dcast(Emp_Inc_RPT.FT.00, Age_Cat + variable ~ Highest_Degree_1 + Sex)
Emp_Inc_RPT.FT.01 <- rename(Emp_Inc_RPT.FT.01, "Age_Group" = "Age_Cat")
Emp_Inc_RPT.FT.02 <- data.frame("Age_1" = as.integer(seq(25:94)+24))

Emp_Inc_RPT.FT.02 <- Emp_Inc_RPT.FT.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_Inc_RPT.FT.01, by = "Age_Group")

# We use this set later to refer to columns referring to calculations for
# individuals of different sexes with different levels of education.
Emp_Inc_RPT.FT_Set.00 = 
  c("E1_None_Female", "E1_None_Male", "E2_HS_Female", "E2_HS_Male", 
    "E3_No_Bach_Female", "E3_No_Bach_Male", "E4_Bach_Plus_Female",
    "E4_Bach_Plus_Male")

Emp_Inc_RPT.FT_Set.01 = 
  paste(Emp_Inc_RPT.FT_Set.00, "_A", sep = "")
Emp_Inc_RPT.FT_Set.02 = 
  paste(Emp_Inc_RPT.FT_Set.01, "_B", sep = "")
Emp_Inc_RPT.FT_Set.02.F = Emp_Inc_RPT.FT_Set.02[grepl("Female",Emp_Inc_RPT.FT_Set.02)]
Emp_Inc_RPT.FT_Set.02.M = setdiff(Emp_Inc_RPT.FT_Set.02, Emp_Inc_RPT.FT_Set.02.F)
Emp_Inc_RPT.FT_Set.03 = 
  paste(Emp_Inc_RPT.FT_Set.02, "_C", sep = "")
# https://stackoverflow.com/questions/7187442/how-to-filter-a-vector-of-strings-in-r-based-on-string-matching
Emp_Inc_RPT.FT_Set.03.F = Emp_Inc_RPT.FT_Set.03[grepl("Female",Emp_Inc_RPT.FT_Set.03)]
Emp_Inc_RPT.FT_Set.03.M = setdiff(Emp_Inc_RPT.FT_Set.03, Emp_Inc_RPT.FT_Set.03.F)
Emp_Inc_RPT.FT_Set.04 = 
  paste(Emp_Inc_RPT.FT_Set.03, "_D", sep = "")
Emp_Inc_RPT.FT_Set.04.F = Emp_Inc_RPT.FT_Set.04[grepl("Female",Emp_Inc_RPT.FT_Set.04)]
Emp_Inc_RPT.FT_Set.04.M = setdiff(Emp_Inc_RPT.FT_Set.04, Emp_Inc_RPT.FT_Set.04.F)

# We confront a new challenge here. In the RWD&TD Matrices Construction script,
# we only had to carry out calculations for one person "type" (i.e., each 
# combination of sex and level of educational attainment) at a time. We could
# simply carry out these calculations independently, but that is extremely
# inefficient. For simple calculations, one alternative approach is to use the
# MUTATE_AT function (part of the DPLYR package), which creates a series of new
# columns based on the contents of another predefined set of columns.
# https://stackoverflow.com/questions/27027347/mutate-each-summarise-each-in-dplyr-how-do-i-select-certain-columns-and-give
# 
# There is a catch, however, which relates to variable naming. In particular,
# the naming process effectively tacks on the given variable name to the label
# of the columns upon which the mutation was based. This means, in effect, that
# column names will gradually grow lengthier as more calculations are performed,
# which, in turn, implies that it is probably advisable to keep names as short
# as possible. What we've done above is assign a new letter to each step of the
# calculations. Thus the variables resulting from Step #1 have the suffix "_A",
# those resulting from Step #2 have the suffix "
Emp_Inc_RPT.FT.03 <- Emp_Inc_RPT.FT.02 %>%
  mutate_at(.vars = vars(Emp_Inc_RPT.FT_Set.00),
            .funs = funs(A = ifelse(Age_1<65, 
                                    ., 
                                    .*(1-(Age_1-65)*ARPD)))) %>%
  mutate_at(.vars = vars(Emp_Inc_RPT.FT_Set.01),
            .funs = funs(B = ifelse(Age_1<65, 
                                    .*Savings_Per, 
                                    0))) %>%
  mutate_at(.vars = vars(Emp_Inc_RPT.FT_Set.02),
            .funs = funs(C = ifelse(Age_1 > 64 | Age_1 == 25, 
                                    0, 
                                    NA)))

for (k in which(colnames(Emp_Inc_RPT.FT.03)=="E1_None_Female_A_B_C"):ncol(Emp_Inc_RPT.FT.03)){
  for (i in 2:40){
    Emp_Inc_RPT.FT.03[i,k] <- Emp_Inc_RPT.FT.03[i,k-8] + Emp_Inc_RPT.FT.03[i-1,k]*(1+Rate)
  }
}

Emp_Inc_RPT.FT.04 <- Emp_Inc_RPT.FT.03 %>%
  mutate_at(.vars = vars(Emp_Inc_RPT.FT_Set.03.F),
            .funs = funs(D = 
                           ifelse(Age_1<65, 
                                  0, 
                                  Rate*.[Age_1==(Year_Ret-1)]/
                                    ((1+Rate)-(1/((1+Rate)^
                                                    (Life_Exp[2,2]-Year_Ret-1))))))) %>%
  mutate_at(.vars = vars(Emp_Inc_RPT.FT_Set.03.M),
            .funs = funs(D = 
                           ifelse(Age_1<65, 
                                  0, 
                                  Rate*.[Age_1==(Year_Ret-1)]/
                                    ((1+Rate)-(1/((1+Rate)^
                                                    (Life_Exp[1,2]-Year_Ret-1)))))))

Emp_Inc_RPT.FT.04 <- Emp_Inc_RPT.FT.04[,c(1:36,40,37,41,38,42,39,43)]

Emp_Inc_RPT.FT.05 <- Emp_Inc_RPT.FT.04 %>%
  mutate(E1_None_Female_A_B_C_D_E = 
           ifelse(Age_1<65, -E1_None_Female_A_B,
                  ifelse(Age_1<ceiling(Life_Exp[2,2]), E1_None_Female_A_B_C_D, 0))) %>%
  mutate(E1_None_Male_A_B_C_D_E = 
           ifelse(Age_1<65, -E1_None_Male_A_B,
                  ifelse(Age_1<ceiling(Life_Exp[1,2]), E1_None_Male_A_B_C_D, 0))) %>%
  mutate(E2_HS_Female_A_B_C_D_E = 
           ifelse(Age_1<65, -E2_HS_Female_A_B,
                  ifelse(Age_1<ceiling(Life_Exp[2,2]), E2_HS_Female_A_B_C_D, 0))) %>%
  mutate(E2_HS_Male_A_B_C_D_E = 
           ifelse(Age_1<65, -E2_HS_Male_A_B,
                  ifelse(Age_1<ceiling(Life_Exp[1,2]), E2_HS_Male_A_B_C_D, 0))) %>%
  mutate(E3_No_Bach_Female_A_B_C_D_E = 
           ifelse(Age_1<65, -E3_No_Bach_Female_A_B,
                  ifelse(Age_1<ceiling(Life_Exp[2,2]), E3_No_Bach_Female_A_B_C_D, 0))) %>%
  mutate(E3_No_Bach_Male_A_B_C_D_E = 
           ifelse(Age_1<65, -E3_No_Bach_Male_A_B,
                  ifelse(Age_1<ceiling(Life_Exp[1,2]), E3_No_Bach_Male_A_B_C_D, 0))) %>%
  mutate(E4_Bach_Plus_Female_A_B_C_D_E = 
           ifelse(Age_1<65, -E4_Bach_Plus_Female_A_B,
                  ifelse(Age_1<ceiling(Life_Exp[2,2]), E4_Bach_Plus_Female_A_B_C_D, 0))) %>%
  mutate(E4_Bach_Plus_Male_A_B_C_D_E = 
           ifelse(Age_1<65, -E4_Bach_Plus_Male_A_B,
                  ifelse(Age_1<ceiling(Life_Exp[1,2]), E4_Bach_Plus_Male_A_B_C_D, 0)))

CS_FT_Emp_Inc <- Emp_Inc_RPT.FT.05[,grepl("Age_1|Age_Group|A$",colnames(Emp_Inc_RPT.FT.05))]
CS_FT_Income_Stream <- Emp_Inc_RPT.FT.05[,grepl("Age_1|Age_Group|A_B_C_D_E$",colnames(Emp_Inc_RPT.FT.05))]
names(CS_FT_Income_Stream)[3:10] <- names(CS_FT_Emp_Inc)[3:10] <- Emp_Inc_RPT.FT_Set.00

for (n in names(CS_FT_Emp_Inc[3:10])){
  Emp_Inc_RPT.FT.05[[paste0(n, "_Pre_Tax")]] <- 
    CS_FT_Emp_Inc[[n]] + CS_FT_Income_Stream[[n]]}

for (n in names(CS_FT_Emp_Inc[3:10])){
  Emp_Inc_RPT.FT.05[[paste0(n, "_Adj_Pens_Earnings")]] <-
    (pmin(CS_FT_Emp_Inc[[n]],Max_Pen_Earn_2017)/Max_Pen_Earn_2017)*Max_Pen_Earn
}

CS_FT_Pre_Tax <- 
  Emp_Inc_RPT.FT.05[, grepl("Age_1|Age_Group|Pre_Tax", colnames(Emp_Inc_RPT.FT.05))]
CS_FT_Adj_Pens_Earnings <- 
  Emp_Inc_RPT.FT.05[, grepl("Age_1|Age_Group|Adj_Pens_Earnings", colnames(Emp_Inc_RPT.FT.05))]
names(CS_FT_Pre_Tax)[3:10] <- names(CS_FT_Adj_Pens_Earnings)[3:10] <- Emp_Inc_RPT.FT_Set.00

# 081918: Originally, we used CS_FT_Emp_Inc in these three calculations.
# This approach was incorrect, since taxes are assumed to be calculated based
# upon all income that is not saved (i.e., because since we are implicitly 
# assuming savings are invested in RRSPs, reducing the overall amount of 
# taxable income).

# 082418: Note that we need to calculate pre-tax income BEFORE we carry out the tabulations
# below (previously, for some reason this was done afterwards, resulting in an error - I'm
# uncertain how this eluded detection before now).
for (n in names(CS_FT_Emp_Inc[3:10])){
  Emp_Inc_RPT.FT.05[[paste0(n, "_Cum_Tax_Amt")]] <- 
    Tax_Rates[findInterval(CS_FT_Pre_Tax[[n]],Tax_Rates$Lower.Limit,all.inside = TRUE),9]}

for (n in names(CS_FT_Emp_Inc[3:10])){
  Emp_Inc_RPT.FT.05[[paste0(n, "_MR")]] <- 
    Tax_Rates[findInterval(CS_FT_Pre_Tax[[n]],Tax_Rates$Lower.Limit,all.inside = TRUE),4]}

for (n in names(CS_FT_Emp_Inc[3:10])){
  Emp_Inc_RPT.FT.05[[paste0(n, "_LowerLim")]] <- 
    Tax_Rates[findInterval(CS_FT_Pre_Tax[[n]],Tax_Rates$Lower.Limit,all.inside = TRUE),1]}

CS_FT_Cum_Tax_Amt <- 
  Emp_Inc_RPT.FT.05[, grepl("Age_1|Age_Group|Cum_Tax_Amt", colnames(Emp_Inc_RPT.FT.05))]
CS_FT_MR <- 
  Emp_Inc_RPT.FT.05[, grepl("Age_1|Age_Group|MR", colnames(Emp_Inc_RPT.FT.05))]
CS_FT_LowerLim <- 
  Emp_Inc_RPT.FT.05[, grepl("Age_1|Age_Group|LowerLim", colnames(Emp_Inc_RPT.FT.05))]
names(CS_FT_Cum_Tax_Amt)[3:10] <- names(CS_FT_MR)[3:10] <- Emp_Inc_RPT.FT_Set.00
names(CS_FT_LowerLim)[3:10] <- Emp_Inc_RPT.FT_Set.00

# 084218: Taxes paid should be based upon PRE-TAX income, not EMPLOYMENT income
# (as was previously the case).
for (n in names(CS_FT_Emp_Inc[3:10])){
  Emp_Inc_RPT.FT.05[[paste0(n, "_TaxPaid")]] <- 
    CS_FT_Cum_Tax_Amt[[n]] + CS_FT_MR[[n]]*(CS_FT_Pre_Tax[[n]] - 
                                                  CS_FT_LowerLim[[n]])}

CS_FT_TaxPaid <- 
  Emp_Inc_RPT.FT.05[, grepl("Age_1|Age_Group|TaxPaid", colnames(Emp_Inc_RPT.FT.05))]
names(CS_FT_TaxPaid)[3:10] <- Emp_Inc_RPT.FT_Set.00

for (n in names(CS_FT_Emp_Inc[3:10])){
  Emp_Inc_RPT.FT.05[[paste0(n, "_CPP_Contrib")]] <- 
    pmin(Max_CPP_Contrib, CPP_Contrib_Rate*(CS_FT_Emp_Inc[[n]]-Basic_Exemption))}

for (n in names(CS_FT_Emp_Inc[3:10])){
  Emp_Inc_RPT.FT.05[[paste0(n, "_EI_Contrib")]] <- 
    pmin(Max_EI_Contrib, EI_Contrib_Rate*CS_FT_Emp_Inc[[n]])}

CS_FT_CPP_Contrib <- 
  Emp_Inc_RPT.FT.05[, grepl("Age_1|Age_Group|CPP_Contrib", colnames(Emp_Inc_RPT.FT.05))]
CS_FT_EI_Contrib <- 
  Emp_Inc_RPT.FT.05[, grepl("Age_1|Age_Group|EI_Contrib", colnames(Emp_Inc_RPT.FT.05))]
names(CS_FT_CPP_Contrib)[3:10] <- names(CS_FT_EI_Contrib)[3:10] <- Emp_Inc_RPT.FT_Set.00

for (n in names(CS_FT_Emp_Inc[3:10])){
  Emp_Inc_RPT.FT.05[[paste0(n, "_Net_Inc")]] <- 
    CS_FT_Pre_Tax[[n]] - CS_FT_TaxPaid[[n]] - CS_FT_CPP_Contrib[[n]] - 
    CS_FT_EI_Contrib[[n]]}

# For validation step
TEST_CASE_E1_M <- 
  Emp_Inc_RPT.FT.05[, grepl("Age_1|Age_Group|E1_None_Male", colnames(Emp_Inc_RPT.FT.05))]

# NOTE 050719: I want to include in the thesis income net of ARPD. To do this, I
# filter on all columns whose labels end in "_A" (and also age, as denoted by 
# "Age_1").
RPT_TBL_FT_GrossEmp <- 
  Emp_Inc_RPT.FT.05[, grepl(
    "Age_1|_A$", colnames(
      Emp_Inc_RPT.FT.05))][RPT_TBL_FT$Age_1 %in% c(seq(25, 85, 10), 94), 
                           c(1,3,2,5,4,7,6,9,8)]
RPT_TBL_FT_Savings <- 
  Emp_Inc_RPT.FT.05[, grepl(
    "Age_1|A_B_C_D_E$", colnames(
      Emp_Inc_RPT.FT.05))][RPT_TBL_FT$Age_1 %in% c(seq(25, 85, 10), 94), 
                           c(1,3,2,5,4,7,6,9,8)]
wb = createWorkbook()
sht1 = addWorksheet(wb, "GrossEmp"); dat1 = RPT_TBL_FT_GrossEmp
sht2 = addWorksheet(wb, "Savings"); dat2 = RPT_TBL_FT_Savings
writeData(wb, sht1, dat1, colNames = T); writeData(wb, sht2, dat2, colNames = T)
saveWorkbook(wb, paste0(Model_Path, "/Report_Tables/", "Employment_Income",
                        LM, ".xlsx"), overwrite = T)

# * Part-Time Employment ====

Emp_Inc_RPT.PT.00 <- melt(Emp_Income, id.vars = c("Age_Cat", "Sex", "Highest_Degree"), 
                              measure.vars = c("Med_Emp_Inc_2017.PT"))
Emp_Inc_RPT.PT.00 <- Emp_Inc_RPT.PT.00 %>%
  mutate("Highest_Degree_1" = 
           case_when(Highest_Degree == "1. No certificate, diploma or degree" ~ "E1_None",
                     Highest_Degree == "2. High school diploma or equivalent" ~ "E2_HS",
                     Highest_Degree == "3. Postsecondary certificate or diploma below bachelor level" ~ "E3_No_Bach",
                     Highest_Degree == "4. University certificate, diploma or degree at bachelor level or above" ~ "E4_Bach_Plus"))

Emp_Inc_RPT.PT.01 <- dcast(Emp_Inc_RPT.PT.00, Age_Cat + variable ~ Highest_Degree_1 + Sex)
Emp_Inc_RPT.PT.01 <- rename(Emp_Inc_RPT.PT.01, "Age_Group" = "Age_Cat")
Emp_Inc_RPT.PT.02 <- data.frame("Age_1" = as.integer(seq(25:94)+24))
# https://stackoverflow.com/questions/27027347/mutate-each-summarise-each-in-dplyr-how-do-i-select-certain-columns-and-give
Emp_Inc_RPT.PT.02 <- Emp_Inc_RPT.PT.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_Inc_RPT.PT.01, by = "Age_Group") 
Emp_Inc_RPT.PT_Set.00 <- Emp_Inc_RPT.FT_Set.00
Emp_Inc_RPT.PT_Set.01 = 
  paste(Emp_Inc_RPT.PT_Set.00, "_A", sep = "")

Emp_Inc_RPT.PT.03 <- Emp_Inc_RPT.PT.02 %>%
  mutate_at(.vars = vars(Emp_Inc_RPT.PT_Set.00),
            .funs = funs(A = ifelse(Age_1<65, 
                                    ., 
                                    .*(1-(Age_1-65)*ARPD))))

# Income streams (i.e., savings and withdrawals) are the same for people working FT and PT,
# so we simply extract these values from calculations we have already undertaken (i.e., for
# FT people).
Emp_Inc_RPT.PT.05 <- Emp_Inc_RPT.PT.03 %>%
  mutate(E1_None_Female_A_B_C_D_E = 
           Emp_Inc_RPT.FT.05$E1_None_Female_A_B_C_D_E[Emp_Inc_RPT.PT.03$Age_1==Emp_Inc_RPT.FT.05$Age_1]) %>%
  mutate(E1_None_Male_A_B_C_D_E = 
           Emp_Inc_RPT.FT.05$E1_None_Male_A_B_C_D_E[Emp_Inc_RPT.PT.03$Age_1==Emp_Inc_RPT.FT.05$Age_1]) %>%
  mutate(E2_HS_Female_A_B_C_D_E = 
           Emp_Inc_RPT.FT.05$E2_HS_Female_A_B_C_D_E[Emp_Inc_RPT.PT.03$Age_1==Emp_Inc_RPT.FT.05$Age_1]) %>%
  mutate(E2_HS_Male_A_B_C_D_E = 
           Emp_Inc_RPT.FT.05$E2_HS_Male_A_B_C_D_E[Emp_Inc_RPT.PT.03$Age_1==Emp_Inc_RPT.FT.05$Age_1]) %>%
  mutate(E3_No_Bach_Female_A_B_C_D_E = 
           Emp_Inc_RPT.FT.05$E3_No_Bach_Female_A_B_C_D_E[Emp_Inc_RPT.PT.03$Age_1==Emp_Inc_RPT.FT.05$Age_1]) %>%
  mutate(E3_No_Bach_Male_A_B_C_D_E = 
           Emp_Inc_RPT.FT.05$E3_No_Bach_Male_A_B_C_D_E[Emp_Inc_RPT.PT.03$Age_1==Emp_Inc_RPT.FT.05$Age_1]) %>%
  mutate(E4_Bach_Plus_Female_A_B_C_D_E = 
           Emp_Inc_RPT.FT.05$E4_Bach_Plus_Female_A_B_C_D_E[Emp_Inc_RPT.PT.03$Age_1==Emp_Inc_RPT.FT.05$Age_1]) %>%
  mutate(E4_Bach_Plus_Male_A_B_C_D_E = 
           Emp_Inc_RPT.FT.05$E4_Bach_Plus_Male_A_B_C_D_E[Emp_Inc_RPT.PT.03$Age_1==Emp_Inc_RPT.FT.05$Age_1])

# By "employment income", we refer to income after accounting for ARPD.
CS_PT_Emp_Inc <- Emp_Inc_RPT.PT.05[,grepl("Age_1|Age_Group|A$",colnames(Emp_Inc_RPT.PT.05))]
CS_PT_Income_Stream <- Emp_Inc_RPT.PT.05[,grepl("Age_1|Age_Group|A_B_C_D_E$",colnames(Emp_Inc_RPT.PT.05))]
names(CS_PT_Income_Stream)[3:10] <- names(CS_PT_Emp_Inc)[3:10] <- Emp_Inc_RPT.PT_Set.00

# "Pre-tax income" refers to income after accounting for income flows attributable to private
# savings.
# Note that we use new for loops to carry out these calculations, when in principle we could
# accomplish the same with much less code using a single statement. The reason we haven't done
# this is that it influences the order of the columns; in particular, columns associated with 
# each formula would be separated rather than clustered together, as is currently the case.
for (n in names(CS_PT_Emp_Inc[3:10])){
  Emp_Inc_RPT.PT.05[[paste0(n, "_Pre_Tax")]] <- 
    CS_PT_Emp_Inc[[n]] + CS_PT_Income_Stream[[n]]}

for (n in names(CS_PT_Emp_Inc[3:10])){
  Emp_Inc_RPT.PT.05[[paste0(n, "_Adj_Pens_Earnings")]] <-
    (pmin(CS_PT_Emp_Inc[[n]],Max_Pen_Earn_2017)/Max_Pen_Earn_2017)*Max_Pen_Earn}

CS_PT_Pre_Tax <- 
  Emp_Inc_RPT.PT.05[, grepl("Age_1|Age_Group|Pre_Tax", colnames(Emp_Inc_RPT.PT.05))]
CS_PT_Adj_Pens_Earnings <- 
  Emp_Inc_RPT.PT.05[, grepl("Age_1|Age_Group|Adj_Pens_Earnings", colnames(Emp_Inc_RPT.PT.05))]
names(CS_PT_Pre_Tax)[3:10] <- names(CS_PT_Adj_Pens_Earnings)[3:10] <- Emp_Inc_RPT.PT_Set.00

for (n in names(CS_PT_Emp_Inc[3:10])){
  Emp_Inc_RPT.PT.05[[paste0(n, "_Cum_Tax_Amt")]] <- 
    Tax_Rates[findInterval(CS_PT_Pre_Tax[[n]],Tax_Rates$Lower.Limit,all.inside = TRUE),9]}

for (n in names(CS_PT_Emp_Inc[3:10])){
  Emp_Inc_RPT.PT.05[[paste0(n, "_MR")]] <- 
    Tax_Rates[findInterval(CS_PT_Pre_Tax[[n]],Tax_Rates$Lower.Limit,all.inside = TRUE),4]}

for (n in names(CS_PT_Emp_Inc[3:10])){
  Emp_Inc_RPT.PT.05[[paste0(n, "_LowerLim")]] <- 
    Tax_Rates[findInterval(CS_PT_Pre_Tax[[n]],Tax_Rates$Lower.Limit,all.inside = TRUE),1]}

# When we create new tables that include results from application of one or more formulas, it
# indicates that none of the formulas in that group were interdependent, in the sense that the
# results of one depended on the results of any of the others. The break therefore signifies
# that the next set of formulas will depend upon recent tabulations.
CS_PT_Cum_Tax_Amt <- 
  Emp_Inc_RPT.PT.05[, grepl("Age_1|Age_Group|Cum_Tax_Amt", colnames(Emp_Inc_RPT.PT.05))]
CS_PT_MR <- 
  Emp_Inc_RPT.PT.05[, grepl("Age_1|Age_Group|MR", colnames(Emp_Inc_RPT.PT.05))]
CS_PT_LowerLim <- 
  Emp_Inc_RPT.PT.05[, grepl("Age_1|Age_Group|LowerLim", colnames(Emp_Inc_RPT.PT.05))]
names(CS_PT_Cum_Tax_Amt)[3:10] <- names(CS_PT_MR)[3:10] <- Emp_Inc_RPT.FT_Set.00
names(CS_PT_LowerLim)[3:10] <- Emp_Inc_RPT.PT_Set.00

for (n in names(CS_PT_Emp_Inc[3:10])){
  Emp_Inc_RPT.PT.05[[paste0(n, "_TaxPaid")]] <- 
    CS_PT_Cum_Tax_Amt[[n]] + CS_PT_MR[[n]]*(CS_PT_Pre_Tax[[n]] - 
                                                  CS_PT_LowerLim[[n]])}

CS_PT_TaxPaid <- 
  Emp_Inc_RPT.PT.05[, grepl("Age_1|Age_Group|TaxPaid", colnames(Emp_Inc_RPT.PT.05))]
names(CS_PT_TaxPaid)[3:10] <- Emp_Inc_RPT.PT_Set.00

for (n in names(CS_PT_Emp_Inc[3:10])){
  Emp_Inc_RPT.PT.05[[paste0(n, "_CPP_Contrib")]] <- 
    pmin(Max_CPP_Contrib, CPP_Contrib_Rate*(CS_PT_Emp_Inc[[n]]-Basic_Exemption))}

for (n in names(CS_PT_Emp_Inc[3:10])){
  Emp_Inc_RPT.PT.05[[paste0(n, "_EI_Contrib")]] <- 
    pmin(Max_EI_Contrib, EI_Contrib_Rate*CS_PT_Emp_Inc[[n]])}

CS_PT_CPP_Contrib <- 
  Emp_Inc_RPT.PT.05[, grepl("Age_1|Age_Group|CPP_Contrib", colnames(Emp_Inc_RPT.PT.05))]
CS_PT_EI_Contrib <- 
  Emp_Inc_RPT.PT.05[, grepl("Age_1|Age_Group|EI_Contrib", colnames(Emp_Inc_RPT.PT.05))]
names(CS_PT_CPP_Contrib)[3:10] <- names(CS_PT_EI_Contrib)[3:10] <- Emp_Inc_RPT.PT_Set.00

for (n in names(CS_PT_Emp_Inc[3:10])){
  Emp_Inc_RPT.PT.05[[paste0(n, "_EI_Contrib")]] <- 
    pmin(Max_EI_Contrib, EI_Contrib_Rate*CS_PT_Emp_Inc[[n]])}

for (n in names(CS_PT_Emp_Inc[3:10])){
  Emp_Inc_RPT.PT.05[[paste0(n, "_Net_Inc")]] <- 
    CS_PT_Pre_Tax[[n]] - CS_PT_TaxPaid[[n]] - CS_PT_CPP_Contrib[[n]] - 
    CS_PT_EI_Contrib[[n]]}

RPT_TBL_PT <- Emp_Inc_RPT.PT.05[, grepl("Age_1|Net_Inc", colnames(Emp_Inc_RPT.PT.05))]
RPT_TBL_PT <- RPT_TBL_PT[RPT_TBL_PT$Age_1 %in% c(25, 35, 45, 55, 65:94),c(1,3,2,5,4,7,6,9,8)]
fwrite(RPT_TBL_PT, paste(Excel, "RPT_TBL_PT", LM, ".csv", sep = ""))

# Table Generation: Retirement Income ------------------------------------------------

Emp_Inc_RPT.RT_Set.00 <- Emp_Inc_RPT.FT_Set.00

# https://stackoverflow.com/questions/25391975/grepl-in-r-to-find-matches-to-any-of-a-list-of-character-strings
Emp_Inc_RPT.RT <- Emp_Inc_RPT.FT.05[,grepl("Age_1|Age_Group|A_B_C_D_E|Adj_Pens_Earnings", colnames(Emp_Inc_RPT.FT.05))]

# https://stackoverflow.com/questions/5487164/r-how-to-replace-parts-of-variable-strings-within-data-frame
colnames(Emp_Inc_RPT.RT) %<>% gsub(., pattern = "A_B_C_D_E", replacement = "Income_Stream")

# https://stackoverflow.com/questions/14324562/creating-new-named-variable-in-dataframe-using-loop-and-naming-convention
RT_COLS <- c("E1_None_Female", "E1_None_Male", "E2_HS_Female", "E2_HS_Male", "E3_No_Bach_Female", 
             "E3_No_Bach_Male", "E4_Bach_Plus_Female", "E4_Bach_Plus_Male")

# CPP benefits are assumed to correspond to those that would be available to an individual
# who worked full-time over the entire length of their careers.
for (n in names(CS_FT_Adj_Pens_Earnings[3:10])){
  Emp_Inc_RPT.RT[[paste0(n, "_CPP")]] <- 
    ifelse(Emp_Inc_RPT.RT$Age_1 < 65, 
           0,
           sum(CS_FT_Adj_Pens_Earnings[[n]][CS_FT_Adj_Pens_Earnings$Age_Group!="65_" & CS_FT_Adj_Pens_Earnings$Age_1!=25])*(0.25/39))
}

CS_RT_Income <- Emp_Inc_RPT.RT[, grepl("Age_1|Age_Group|Income_Stream", colnames(Emp_Inc_RPT.RT))]
CS_RT_CPP <- Emp_Inc_RPT.RT[, grepl("Age_1|Age_Group|CPP", colnames(Emp_Inc_RPT.RT))]
names(CS_RT_Income)[3:10] <- names(CS_RT_CPP)[3:10] <- Emp_Inc_RPT.RT_Set.00

for (n in names(CS_RT_Income[3:10])){
  Emp_Inc_RPT.RT[[paste0(n, "_OAS")]] <- 
    ifelse(Emp_Inc_RPT.RT$Age_1 < 65, 
           0,
           ifelse(CS_RT_Income[[n]] + CS_RT_CPP[[n]] + OAS_Amt <= OAS_Threshold,
                  12*OAS_Amt,
                  ifelse(CS_RT_Income[[n]] + CS_RT_CPP[[n]] > OAS_Threshold,
                         12*pmax(OAS_Amt - 0.15*
                                   (CS_RT_Income[[n]] + CS_RT_CPP[[n]] + OAS_Amt - 
                                      OAS_Threshold),0),
                         12*OAS_Threshold - (CS_RT_Income[[n]] + CS_RT_CPP[[n]]) + 
                           12*pmax(OAS_Amt - (OAS_Threshold - 
                                                (CS_RT_Income[[n]] + CS_RT_CPP[[n]])) -
                                     0.15*(CS_RT_Income[[n]] + CS_RT_CPP[[n]] + OAS_Amt - 
                                             OAS_Threshold),0))))
}

for (n in names(CS_RT_Income[3:10])){
  Emp_Inc_RPT.RT[[paste0(n, "_GIS")]] <- 
    ifelse(Emp_Inc_RPT.RT$Age_1 < 65, 
           0,
           ifelse(CS_RT_Income[[n]] + CS_RT_CPP[[n]] <= GIS_Threshold,
                  12*GIS_Amt,
                  pmax(12*GIS_Amt - 0.50*(CS_RT_Income[[n]] + CS_RT_CPP[[n]] - 
                                            GIS_Threshold),0)))
}

CS_RT_OAS <- Emp_Inc_RPT.RT[, grepl("Age_1|Age_Group|OAS", colnames(Emp_Inc_RPT.RT))]
CS_RT_GIS <- Emp_Inc_RPT.RT[, grepl("Age_1|Age_Group|GIS", colnames(Emp_Inc_RPT.RT))]
names(CS_RT_OAS)[3:10] <- names(CS_RT_GIS)[3:10] <- Emp_Inc_RPT.RT_Set.00

for (n in names(CS_RT_Income[3:10])){
  Emp_Inc_RPT.RT[[paste0(n, "_Pre_Tax")]] <- CS_RT_Income[[n]] + CS_RT_CPP[[n]] +
    CS_RT_OAS[[n]] + CS_RT_GIS[[n]]}

CS_RT_Pre_Tax <- Emp_Inc_RPT.RT[, grepl("Age_1|Age_Group|Pre_Tax", colnames(Emp_Inc_RPT.RT))]
names(CS_RT_Pre_Tax)[3:10] <- Emp_Inc_RPT.RT_Set.00

# These expressions isolate taxable and non-taxable income, respectively.
for (n in names(CS_RT_Income[3:10])){
  Emp_Inc_RPT.RT[[paste0(n, "_TaxY_Inc")]] <- CS_RT_Pre_Tax[[n]] - CS_RT_GIS[[n]]}
for (n in names(CS_RT_Income[3:10])){
  Emp_Inc_RPT.RT[[paste0(n, "_TaxN_Inc")]] <- CS_RT_GIS[[n]]}

CS_RT_TaxY <- Emp_Inc_RPT.RT[, grepl("Age_1|Age_Group|TaxY_Inc", colnames(Emp_Inc_RPT.RT))]
CS_RT_TaxN <- Emp_Inc_RPT.RT[, grepl("Age_1|Age_Group|TaxN_Inc", colnames(Emp_Inc_RPT.RT))]
names(CS_RT_TaxY)[3:10] <- names(CS_RT_TaxN)[3:10] <- Emp_Inc_RPT.RT_Set.00

for (n in names(CS_RT_Income[3:10])){
  Emp_Inc_RPT.RT[[paste0(n, "_Cum_Tax_Amt")]] <- 
    Tax_Rates[findInterval(CS_RT_TaxY[[n]],Tax_Rates$Lower.Limit,all.inside = TRUE),9]}

for (n in names(CS_RT_Income[3:10])){
  Emp_Inc_RPT.RT[[paste0(n, "_MR")]] <- 
    Tax_Rates[findInterval(CS_RT_TaxY[[n]],Tax_Rates$Lower.Limit,all.inside = TRUE),4]}

for (n in names(CS_RT_Income[3:10])){
  Emp_Inc_RPT.RT[[paste0(n, "_LowerLim")]] <- 
    Tax_Rates[findInterval(CS_RT_TaxY[[n]],Tax_Rates$Lower.Limit,all.inside = TRUE),1]}

CS_RT_Cum_Tax_Amt <- 
  Emp_Inc_RPT.RT[, grepl("Age_1|Age_Group|Cum_Tax_Amt", colnames(Emp_Inc_RPT.RT))]
CS_RT_MR <- 
  Emp_Inc_RPT.RT[, grepl("Age_1|Age_Group|MR", colnames(Emp_Inc_RPT.RT))]
CS_RT_LowerLim <- 
  Emp_Inc_RPT.RT[, grepl("Age_1|Age_Group|LowerLim", colnames(Emp_Inc_RPT.RT))]
names(CS_RT_Cum_Tax_Amt)[3:10] <- names(CS_RT_MR)[3:10] <- Emp_Inc_RPT.RT_Set.00
names(CS_RT_LowerLim)[3:10] <- Emp_Inc_RPT.RT_Set.00

for (n in names(CS_RT_Income[3:10])){
  Emp_Inc_RPT.RT[[paste0(n, "_TaxPaid")]] <- 
    ifelse(Emp_Inc_RPT.RT$Age_1 < 65,
           0,
           CS_RT_Cum_Tax_Amt[[n]] + CS_RT_MR[[n]]*(CS_RT_TaxY[[n]] - 
                                                         CS_RT_LowerLim[[n]]))}

CS_RT_TaxPaid <- 
  Emp_Inc_RPT.RT[, grepl("Age_1|Age_Group|TaxPaid", colnames(Emp_Inc_RPT.RT))]
names(CS_RT_TaxPaid)[3:10] <- Emp_Inc_RPT.RT_Set.00

for (n in names(CS_RT_Income[3:10])){
  Emp_Inc_RPT.RT[[paste0(n, "_Net_Inc")]] <- 
    ifelse(Emp_Inc_RPT.RT$Age_1 < 65,
           0,
           CS_RT_Pre_Tax[[n]] - CS_RT_TaxPaid[[n]])}

for (n in names(CS_RT_Income[3:10])){
  Emp_Inc_RPT.RT[[paste0(n, "_LowerLim")]] <- 
    Tax_Rates[findInterval(CS_RT_TaxY[[n]],Tax_Rates$Lower.Limit,all.inside = TRUE),1]}

for (n in names(CS_RT_Income[3:10])){
  Emp_Inc_RPT.RT[[paste0(n, "_Pub_Pen")]] <- CS_RT_CPP[[n]] + CS_RT_OAS[[n]] + 
    CS_RT_GIS[[n]]}

# For validation step
TEST_CASE_E3_F <- 
  Emp_Inc_RPT.RT[, grepl("Age_1|Age_Group|E3_No_Bach_Female", colnames(Emp_Inc_RPT.RT))]

# 050719 - I want to include in the report "before" and "afters" based upon 
# retirement income. The "Pre_Tax" and "Net_Tax" suffixes mark columns that relate
# to the former and latter, respectively. Let's take this information, and cram it
# into a single Excel workbook, so that we can included them in the report.
RPT_TBL_RT_Pub_Pen <-  
  Emp_Inc_RPT.RT[Emp_Inc_RPT.RT$Age_1 %in% c(65, 94), 
                 grepl("Age_1|Pub_Pen", 
                       colnames(Emp_Inc_RPT.RT))][, c(3,2,5,4,7,6,9,8)]
RPT_TBL_RT_Pub_Pen %<>% mutate(Status = c("Savings", "No_Savings")) %>%
  melt(., id.vars = "Status", variable = "Education", value.name = "Pub_Pen") %>% 
  mutate(Gender = ifelse(grepl("Female", Education), "F", "M"))
RPT_TBL_RT_Pub_Pen$Education %<>% 
  gsub(., "_Male_Pub_Pen", "") %>% gsub(., "_Female_Pub_Pen", "")
RPT_TBL_RT_Pub_Pen %<>% 
  dcast(., Education ~ Gender + Status, value.var = "Pub_Pen")
RPT_TBL_RT_Pub_Pen <- RPT_TBL_RT_Pub_Pen[, c(1,5,4,3,2)]

RPT_TBL_RT_Assets <-  
  Emp_Inc_RPT.RT[Emp_Inc_RPT.RT$Age_1 %in% c(65, 94), 
                 grepl("Age_1|Income_Stream", 
                       colnames(Emp_Inc_RPT.RT))][, c(3,2,5,4,7,6,9,8)]
RPT_TBL_RT_Assets %<>% mutate(Status = c("Savings", "No_Savings")) %>%
  melt(., id.vars = "Status", variable = "Education", value.name = "Assets")
RPT_TBL_RT_Assets %<>% mutate(Gender = ifelse(grepl("Female", Education), "F", "M"))
RPT_TBL_RT_Assets$Education %<>% 
  gsub(., pattern = "_Male_Income_Stream", replacement = "") %>%
  gsub(., pattern = "_Female_Income_Stream", replacement = "")
RPT_TBL_RT_Assets %<>% 
  dcast(., Education ~ Gender + Status, value.var = "Assets")
RPT_TBL_RT_Assets <- RPT_TBL_RT_Assets[, c(1,5,4,3,2)]

RPT_TBL_RT_Pre_Tax <-  
  Emp_Inc_RPT.RT[Emp_Inc_RPT.RT$Age_1 %in% c(65, 94), 
                 grepl("Age_1|Pre_Tax", 
                       colnames(Emp_Inc_RPT.RT))][, c(3,2,5,4,7,6,9,8)]
RPT_TBL_RT_Pre_Tax %<>% mutate(Status = c("Savings", "No_Savings")) %>%
  melt(., id.vars = "Status", variable = "Education", value.name = "Pre_Tax")
RPT_TBL_RT_Pre_Tax %<>% mutate(Gender = ifelse(grepl("Female", Education), "F", "M"))
RPT_TBL_RT_Pre_Tax$Education %<>% 
  gsub(., pattern = "_Male_Pre_Tax", replacement = "") %>%
  gsub(., pattern = "_Female_Pre_Tax", replacement = "")
RPT_TBL_RT_Pre_Tax %<>% 
  dcast(., Education ~ Gender + Status, value.var = "Pre_Tax")
RPT_TBL_RT_Pre_Tax <- RPT_TBL_RT_Pre_Tax[, c(1,5,4,3,2)]
  
RPT_TBL_RT_Net_Inc <-  
  Emp_Inc_RPT.RT[Emp_Inc_RPT.RT$Age_1 %in% c(65, 94), 
                 grepl("Age_1|Net_Inc", 
                       colnames(Emp_Inc_RPT.RT))][, c(3,2,5,4,7,6,9,8)]
RPT_TBL_RT_Net_Inc %<>% mutate(Status = c("Savings", "No_Savings")) %>%
  melt(., id.vars = "Status", variable = "Education", value.name = "Net_Inc")
RPT_TBL_RT_Net_Inc %<>% mutate(Gender = ifelse(grepl("Female", Education), "F", "M"))
RPT_TBL_RT_Net_Inc$Education %<>% 
  gsub(., pattern = "_Male_Net_Inc", replacement = "") %>%
  gsub(., pattern = "_Female_Net_Inc", replacement = "")
RPT_TBL_RT_Net_Inc %<>% 
  dcast(., Education ~ Gender + Status, value.var = "Net_Inc")
RPT_TBL_RT_Net_Inc <- RPT_TBL_RT_Net_Inc[, c(1,5,4,3,2)]

wb = createWorkbook()
sht = list(addWorksheet(wb, "Pub_Pen"), addWorksheet(wb, "Assets"), 
           addWorksheet(wb, "Pre_Tax"), addWorksheet(wb, "Net_Inc"))
dat = list(RPT_TBL_RT_Pub_Pen, RPT_TBL_RT_Assets, RPT_TBL_RT_Pre_Tax, 
           RPT_TBL_RT_Net_Inc)
for (i in 1:4){writeData(wb, sht[[i]], dat[[i]], colNames = T)}
saveWorkbook(wb, paste0(Model_Path, "/Report_Tables/", "Retirement_Income",
                        LM, ".xlsx"), overwrite = T)

# DPoRT 2.0 Table ---------------------------------------------------------

# This table reports the probabilities of transitioning from normoglycemia or
# prediabetes to T2DM, using estimates obtained using the DPoRT 2.0 Risk
# Equations (the table generated through these calculations corresponds to 
# Table 11 in V10 of the report [dated 060919], entitled, "One-year 
# probabilities (%) of developing T2DM, by gender, ethnicity, level of 
# education, age, adherence to LTPA, and dietary quality").

# The development of the table illustrates a more general approach to 
# transferring information from one data frame from to another. This approach 
# involves carrying out the following three steps:
# 1) "Flatten" both tables, using the MELT function from the RESHAPE2 package.
# This enables us to see where there is misalignment in the definition of 
# particular variables.
# 2) Develop "bridging tables". This involves compiling tables that show how
# variables that are defined in different ways across tables are related. It
# should be noted that problems can emerge if the categories used to define
# one variable do not wholly encompass or are not wholly encompassed by the
# categories representing the corresponding variable in the second table.
# 3) Using the "bridging tables" to redefine variables in one or both tables
# to ensure alignment. In R, this can readily be accomplished using one of
# DPLYR's join functions, although I believe the MERGE function in the BASE
# package should work equally well.
# 4) "Transfer" the data by merging the tables. This is made possible by the
# fact that variables are defined similarly in each. Practically, this is
# achieved by performing a left join from the table containing the data of
# interest to the one into which we'd like this information to be integrated.
# With reference to the syntax for DPLYR's join functions, we join by all of
# the variables that are now defined in the same way across tables.
load(paste0(RDA, "DPoRT_Calculations", LM, "_One_Year_Probs", ".Rda"))
DPoRT.00 <- data.table(RPT_TBL)
Emp_Income.00 <- data.table(Emp_Income)
  
# STEP 1: Flatten data
DPoRT.01 <- melt(DPoRT.00, id = 1:2, measure = 3:14)
# Eliminate the "Very Active" and "Low Active" activity levels, since the
# calculation of equilibrium BMI/probabilities does not utilize these.
DPoRT.02 <- DPoRT.01[!grepl("Very_Active|Low_Active",DPoRT.01$variable),,]
# In two of the columns, the variables of interest are concatenated into a
# string. Fortunately, we can use the SEPARATE function (from the TIDYR)
# package to split each string everywhere there is an underscore.
DPoRT.03 <- DPoRT.02 %>%
  # https://stackoverflow.com/questions/4350440/split-data-frame-string-column-into-multiple-columns
  separate(variable, c("LTPA","PostSec","TopQuin","Risk"),"_") %>%
  separate(variable.1, c("Ethnicity", "Sex", "Diet"),"_")
DPoRT.03 <- DPoRT.03[,!c("Risk"),]
# We now classify DPoRT.03 as a data frame, rather than a data table. This
# was done because some functions appear not to work well with data tables.
DPoRT.03 <- data.frame(DPoRT.03)
# For some reason, the variable Age_Group was defined as a factor. We need
# to rectify this because it is defined as a character in Emp_Income.00. 
DPoRT.03$Age_Group <- as.character(DPoRT.03$Age_Group)

# Get rid of columns we don't need. The actual income level isn't relevant
# for the purposes of the calculations, only the determination that a
# person is, or is not in the top income quintile. One implication of this
# is that we can quickly update the table if we change how median income is
# defined (i.e., because this determination has already taken place).
Emp_Income.01 <- Emp_Income.00[,!c("Med_Emp_Inc_2017.FT", "Med_Emp_Inc_2017.PT"),]

# Here, we test out a new function, EXPAND.GRID.DF (defined above), which is 
# drawn directly from the following link:
# https://stackoverflow.com/questions/11693599/alternative-to-expand-grid-for-data-frames
# The function enables us to compile a new data frame consisting of all
# possible combinations of variables across several more primitive data
# frames. It should be noted that the expand.grid.df function will not work 
# if Emp_Income.01 is defined as a data.table, which is why we converted it
# to a data frame above.
LTPA <- data.frame(LTPA = c("Sedentary", "Active"), stringsAsFactors = F); 
Diet <- data.frame(Diet = c("Recommended","Excessive"), stringsAsFactors = F)
Ethnicity <- data.frame(Ethnicity = c("Afro-Caribbean", "Asian-Indian", "Caucasian"), 
                        stringsAsFactors = F)
Emp_Income.01 <- data.frame(Emp_Income.01)
Emp_Income.02 <- expand.grid.df(Emp_Income.01, LTPA, Diet, Ethnicity)
Emp_Income.02$IsTopQuin <- as.character(Emp_Income.02$IsTopQuin)

# STEP 3: Construct Bridging Tables

# We need three tables to faciliate variable redefinition. The first aligns
# the highest degree ("Highest_Degree") variable (which originates from the
# 2011 NHS data) with the post-secondary attendance ("PostSec") variable 
# (which originates from the DPoRT 2.0 Risk Equations).

# The first table relates to educational attainment. We have elected to
# define post-secondary attendance in terms of attainment of any PS degree,
# certificate or diploma, acknowledging that at least some HS graduates will
# have attended PS studies without obtaining a credential (this is discussed
# in greater depth in the RWD&TP_Matrices Construction script).
Educ_Att_Levels <- 
  data.frame(Highest_Degree = c("1. No certificate, diploma or degree", 
                        "2. High school diploma or equivalent", 
                        "3. Postsecondary certificate or diploma below bachelor level", 
                        "4. University certificate, diploma or degree at bachelor level or above"),
             PostSec = c(rep("NoPS",2), rep("PS", 2)), stringsAsFactors = F)

# The second table aligns age categories. the DPoRT 2.0 Risk Equations
# employ fewer categories than the NHS data (although, now that I think of it,
# the categories in the NHS data may be artificial). Fortunately, however, the
# categories in the latter source do not overlap those in the former, making
# the process of alignment much simpler.
Age_Align <- data.frame("Age_Cat" = c("25_34", "35_44", "45_54", "55_64", "65_"),
                        "Age_Group" = c("Age_19_44", "Age_19_44", "Age_45_64", "Age_45_64",
                                        "Age_65_95"), stringsAsFactors = F)

# The third and last table aligns the definitions of top quintile occupation
# across data frames. Since all that differentiates the two is the column
# label and the actual phrase used to denote occupation/non-occupation, this
# is the simplest of all.
IncQuin <- data.frame("IsTopQuin" = c("FALSE", "TRUE"), "TopQuin" = c("NoTQ", "TQ"),
                      stringsAsFactors = F)

# Now we augment the data frame to which the information will be transferred
# by carrying out a series of left joins. I originally tried to skip over this
# step using complicated referencing techniques, but this did not work as
# intended (we can try again later).
Emp_Income.03 <- Emp_Income.02 %>% left_join(Age_Align, by = "Age_Cat")
Emp_Income.03 <- Emp_Income.03 %>% left_join(Educ_Att_Levels, by = "Highest_Degree")
Emp_Income.03 <- Emp_Income.03 %>% left_join(IncQuin, by = "IsTopQuin")

# STEP 4: Merge Tables
RPT_TBL_DPoRT.00 <- Emp_Income.03 %>% 
  left_join(DPoRT.03, by = c("Age_Group", "LTPA", "PostSec", "TopQuin", "Ethnicity", "Sex", "Diet"))

RPT_TBL_DPoRT.00 <- data.table(RPT_TBL_DPoRT.00)
RPT_TBL_DPoRT.01 <- RPT_TBL_DPoRT.00[,!c("IsTopQuin", "Age_Group","PostSec","TopQuin"),]
RPT_TBL_DPoRT.01 <- RPT_TBL_DPoRT.01[,c(1,4,5,3,6,2,7),]
RPT_TBL_DPoRT.01 <- RPT_TBL_DPoRT.01 %>% rename("Prob" = "value")
# https://stackoverflow.com/questions/1296646/how-to-sort-a-dataframe-by-multiple-columns
# https://stackoverflow.com/questions/7793295/how-to-order-a-data-frame-by-one-descending-and-one-ascending-column
# https://stackoverflow.com/questions/13685295/sort-a-data-table-fast-by-ascending-descending-order
setorder(RPT_TBL_DPoRT.01, Age_Cat, LTPA, -Diet, Highest_Degree, Ethnicity, -Sex)

# For some reason, the descending sorts performed in the previous step are undone
# when I apply the DCAST function, necessating some manual manipulation in Excel
# after the table is exported.
RPT_TBL_DPoRT.01 <- dcast.data.table(RPT_TBL_DPoRT.01, Age_Cat + LTPA + Diet ~ 
                                       Highest_Degree + Ethnicity + Sex, value.var = "Prob")
setorder(RPT_TBL_DPoRT.01, Age_Cat, LTPA, -Diet)
  
RPT_TBL_DPoRT.01 <- 
  RPT_TBL_DPoRT.01[, c(1:3, 5, 4, 7, 6, 9, 8, 11, 10, 13, 12, 15, 14, 17, 16, 19, 18, 21, 20, 
                     23, 22, 25, 24, 27, 26)]
RPT_TBL_DPoRT.02 <- RPT_TBL_DPoRT.01 %>% 
  transmute_if(is.numeric, .funs = function(x) {x <- x*100})
RPT_TBL_DPoRT.02 <- cbind(RPT_TBL_DPoRT.01[, c(1:3)], RPT_TBL_DPoRT.02)

fwrite(RPT_TBL_DPoRT.02, paste(Excel, Filename, LM, ".csv", sep = ""))

# Equilibrium EER and BMI Tables ------------------------------------------------------

# Now we want to compile this information into a single table that we can
# use for the purposes of reporting. In addition, the table may be referenced
# directly in the decision model (although additional transformations may be
# required before it is fit for this latter purpose).

Pred_EER_1 <- Pr_Risk[Pr_Risk$Caloric_Intake=="Recommended",grepl("Gender|Age_Group|Ethnicity|Caloric_Intake|^EER.Sedentary$|^EER.Active$",colnames(Pr_Risk))]
Pred_EER_2 <- melt(Pred_EER_1, id=c(1:3), measure=c(5:6))
Pred_EER_3 <- dcast(Pred_EER_2, variable + Age_Group ~ Ethnicity + Gender
                    , value.var = "value")
Pred_EER_4 <- Pred_EER_3 %>% rename("Activity_Level" = "variable")

fwrite(Pred_EER_4, paste(Excel, Filename, "_EER", LM, ".csv", sep = ""))

Pred_BMI_1 <- Pr_Risk[,grepl("Gender|Age_Group|Ethnicity|Caloric_Intake|^Sedentary.BMI$|^Active.BMI$",colnames(Pr_Risk))]
Pred_BMI_2 <- melt(Pred_BMI_1, id=c(1:4), measure=c(5:6))
Pred_BMI_3 <- dcast(Pred_BMI_2, variable + Caloric_Intake + Age_Group ~ Ethnicity + Gender
                    , value.var = "value")
Pred_BMI_4 <- Pred_BMI_3 %>%
  rename("LTPA" = "variable", "Healthy_Eating" = "Caloric_Intake")

fwrite(Pred_BMI_4, paste(Excel, Filename, "_BMI", LM, ".csv", sep = ""))