# R Code for Development of the Simulation Cohort
# Created: 10/05/18; Last Modified: 11/18/18

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

# NOTE 111818: The code now detects the path where the script is located, eliminating 
# the need for the user to change this manually upon downloading the script for the 
# first time.
# 
# This code develops the cohort for the simulation exercise. The information
# necessary to achieve this is drawn from two main sources:
# 1) The 2011 NHS PUMF Individuals File
# 2) Demographic statistics (Population Report - June 1, 2017) produced by
# Manitoba Health, Seniors and Active Living (see page 11: 
# https://www.gov.mb.ca/health/population/pr2017.pdf)

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

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

# Call requisite libraries
library(dplyr); library(tidyr); library(readxl); library(xlsx); library(CANSIM2R)
library(reshape2); library(tibble); library(rvest); library(pdftools)
library(magrittr); library(data.table); library(here)

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

# Preparing the Data ------------------------------------------------------

# We begin with the 2011 NHS PUMF data (99-014-X2011041). Using the FREAD
# command ensures the information will be imported as a data table, enabling
# us to take advantage of the data manipulation functionality of the
# DATA.TABLE package.
# http://www12.statcan.gc.ca/nhs-enm/2011/dp-pd/prof/index.cfm?Lang=E
# http://www.dummies.com/programming/r/how-to-use-read-csv-to-import-data-in-r/
Demogr_NHS.00 <- 
  fread(paste0(Data_Path, "NHS-99M001X-E-2011-pumf-individuals_F1.csv"))

Demogr_NHS.01 <- Demogr_NHS.00 %>% 
  # Focus on individuals in the age 25-29 age range, who reported their level
  # of educational attainment.
  filter(AGEGRP == 9, HDGREE %in% c(1:13), PR == 46) %>%
  # Select three variables for the remainder of the analysis:
  # 1) Sex/Gender
  # 2) Highest certificate, diploma or degree
  # 3) Population group (i.e., ethnicity)
  select(SEX, HDGREE, DPGRSUM) %>%
  # Aggregate responses relating to level of educational attainment (this
  # preserves consistency with the approach applied elsewhere in the model).
  mutate(HDGREE.01 = case_when(
    HDGREE %in% 1 ~ "1. No certificate, diploma or degree",
    HDGREE %in% 2 ~ "2. High school diploma or equivalent",
    HDGREE %in% 3:8 ~ "3. Postsecondary certificate or diploma below bachelor level",
    HDGREE > 8 ~ "4. University certificate, diploma or degree at bachelor level or above")) %>%
  mutate(SEX.01 = case_when(SEX == 1 ~ "Female", SEX == 2 ~ "Male")) %>%
  # This next piece is important, because of the shortcomings it reveals in 
  # the methodology underlying the simulations. Because the UKPDS-OM2
  # includes only three ethnic groups (i.e., Caucasian, Afro-Caribbean, and
  # Asian-Indian), it is necessary to somehow align these with the wide
  # range of ethnic groups recorded in the 2011 NHS. Here's how I've chosen
  # to do this.
  mutate(DPGRSUM.01 = case_when(
    # A person is classified as Caucasian if they are recorded in the PUMF
    # as being "white" OR "white and visible minority(ies)" OR if this
    # information is not available. 
    DPGRSUM %in% c(1, 14, 88) ~ "Caucasian",
    # A person is classified as Afro-Caribbean if they are recorded in the
    # PUMF as being "black".
    DPGRSUM == 4 ~ "Afro-Caribbean",
    # Everyone else is classified as being Asian-Indian. This includes
    # individuals who are categorized as "Arab" OR "visible minority, 
    # n.i.e." OR "multiple visible minorities, n.i.e.". Critically, it also 
    # includes those who self-identify as Aboriginal. This is defensible 
    # only in that Asian-Indian and Indigenous Peoples are both at very high 
    # risk for T2DM, and that this is tied to non-modifiable risk factors 
    # (i.e., genetic predisposition).
    TRUE ~ "Asian-Indian")) %>%
  # Having labelled the contents of this excerpt of the dataset in a
  # manner than facilitates straightforward interpretation, we now filter
  # out the original variables.
  select(SEX.01, HDGREE.01, DPGRSUM.01)

# For some reason, the manipulations described above converted the 
# data table to a dataframe. We now reverse this, simultaneously
# filtering on gender so as to divide the dataset in half.
Demogr_NHS.02_M <- data.table(Demogr_NHS.01)[SEX.01=="Male"]
Demogr_NHS.02_F <- data.table(Demogr_NHS.01)[SEX.01=="Female"]

# These next lines of code, which involve a series of manipulations using the
# subset of the sub-sample of PUMF records for males, illustrate the power of 
# the DATA.TABLE package.
Demogr_NHS.02_M %<>% 
  # We first generate record counts for every combination of ethnicity,
  # gender, and level of educational attainment. Note that the order in
  # which the variables appear in the statement below determine the order
  # in which they will then appear in the data table.
  .[, .N, by = .(DPGRSUM.01, SEX.01, HDGREE.01)] %>%
  # Now sort the variables
  .[order(DPGRSUM.01, SEX.01, HDGREE.01)] 
# In the DATA.TABLE equivalent to MUTATE (which is part of the DPLYR package), 
# we create a new variable that converts these counts to proportions.
# http://brooksandrew.github.io/simpleblog/articles/advanced-data-table/
Demogr_NHS.02_M %<>% .[, Prop := N/sum(N)]

# Now we carry out the same calculations using the subset of the sub-
# sample of PUMF records for females.
Demogr_NHS.02_F %<>% 
  .[, .N, by = .(DPGRSUM.01, SEX.01, HDGREE.01)] %>%
  .[order(DPGRSUM.01, SEX.01, HDGREE.01)] 
# http://brooksandrew.github.io/simpleblog/articles/advanced-data-table/
Demogr_NHS.02_F %<>% .[, Prop := N/sum(N)]

# These counts come directly from demographic statistics produced by
# Manitoba Health, Seniors and Active Living (refer to the "Objective
# section above for a link to the original source), and represent the
# number of males and females who were aged 25 as of June 1, 2017.
# Ordinarily, my preference would have been to scrape this data directly
# from the PDF file; however, the page from which the data is taken is 
# rotated. Given how little information we actually require from this
# document, investing the time and effort necessary to automate its
# collection is not justified.
Num_M <- 10208 ; Num_F <- 9927
# Now that we know how many males and females comprise the cohort, we
# multiply these with the proportions derived using the sub-sample of 
# the 2011 NHS PUMF. Implicitly, the assumption is that the ethnic and
# educational characteristics of the individuals in the sub-sample -
# which, to be clear, represent Canada as a whole - are reasonably 
# similar to those of Manitoba youth.
Demogr_NHS.02_M %<>% .[, Cohort_Size.00 := round(Prop*Num_M)]
Demogr_NHS.02_F %<>% .[, Cohort_Size.00 := round(Prop*Num_F)]

Demogr_NHS.03_M <- Demogr_NHS.02_M[, CumProp := cumsum(Prop)]
Rand_M <- runif(1); Index <- 1:nrow(Demogr_NHS.03_M)
Demogr_NHS.03_M <- cbind(Index, Demogr_NHS.03_M)
Demogr_NHS.03_M %<>% .[, Lottery := Index == 
                       findInterval(Rand_M, CumProp)]
Demogr_NHS.03_M %<>% 
  .[, Cohort_Size.01 := 
      ifelse(Lottery==T, Cohort_Size.00 - 1, Cohort_Size.00)]

Demogr_NHS.03_F <- Demogr_NHS.02_F[, CumProp := cumsum(Prop)]
Rand_F <- runif(2); Index <- 1:nrow(Demogr_NHS.03_F)
Demogr_NHS.03_F <- cbind(Index, Demogr_NHS.03_F)
Demogr_NHS.03_F %<>% .[, c("Lottery_1", "Lottery_2") := 
                         list(Index == findInterval(Rand_F[1], CumProp),
                              Index == findInterval(Rand_F[2], CumProp))]
Demogr_NHS.03_F %<>% 
  .[, Lottery := Lottery_1 + Lottery_2] %>%
  .[, Cohort_Size.01 := Cohort_Size.00 - Lottery]

Demogr_NHS.04_M <- Demogr_NHS.03_M[, c("DPGRSUM.01", "SEX.01", "HDGREE.01", 
                                       "Cohort_Size.01")]
Demogr_NHS.04_M %<>% .[, .(Ethnicity = DPGRSUM.01, Gender = SEX.01,
                           Education = HDGREE.01, Number = Cohort_Size.01)]
Demogr_NHS.04_F <- Demogr_NHS.03_F[, c("DPGRSUM.01", "SEX.01", "HDGREE.01", 
                                       "Cohort_Size.01")]
Demogr_NHS.04_F %<>% .[, .(Ethnicity = DPGRSUM.01, Gender = SEX.01,
                           Education = HDGREE.01, Number = Cohort_Size.01)]

Demogr_NHS <- rbind(Demogr_NHS.04_M, Demogr_NHS.04_F)
write.table(Demogr_NHS,
            file = paste0(MATLAB, "Simulation_Results/", "Starting_Cohort", ".dat"), 
            quote = FALSE, sep = ";", row.names = FALSE)

# Prepare Table for Report ------------------------------------------------

Demogr_NHS_RPT.00 <- Demogr_NHS
Demogr_NHS_RPT.01 <- 
  dcast.data.table(Demogr_NHS_RPT.00, Ethnicity + Education ~ Gender, 
                   value.var = "Number")

Demogr_NHS_RPT.01 %<>% 
  .[, c(1:2, 4, 3)] %<>%
  .[, Total := Male + Female]

Demogr_NHS_RPT.01.T <- 
  Demogr_NHS_RPT.01[, .(Ethnicity = "Total", Education = "Total", Male = sum(Male), 
                        Female = sum(Female), Total = sum(Total))]

Demogr_NHS_RPT.02 <- rbind(Demogr_NHS_RPT.01, Demogr_NHS_RPT.01.T)
write.table(Demogr_NHS_RPT.02,
            file = paste0(MATLAB, "Simulation_Results/", "Starting_Cohort_RPT", 
                          ".dat"), quote = FALSE, sep = ";", row.names = FALSE)