# R Code for Application of DPoRT Risk Equations 
# Created: 02/10/18; Last Modified: 12/22/18

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

# NOTE 122218: The code now enables the user to directly modify the 
# definition of recommended and excessive caloric consumption through the
# new parameters Rec_CI and Exc_CI (where "CI" refers to "caloric intake").

# 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 060818: The code has been now been revised to draw upon the most
# recent DPoRT risk equations (i.e., DPoRT 2.0) presented in: Rosella, L. 
# C., Lebenbaum, M., Li, Y., Wang, J., & Manuel, D. G. (2014). Risk
# distribution and its influence on the population targets for diabetes 
# prevention. Preventive Medicine, 58(Supplement C), 17-21. 
# https://doi.org/10.1016/j.ypmed.2013.10.007
# 
# NOTE 021018: Some of the calculations in this code are based on results
# presented in: Rosella, L. C., Manuel, D. G., Burchill, C., & Stukel, T. A. 
# (2011). A population-based risk algorithm for the development of diabetes: 
# development and validation of the Diabetes Population Risk Tool (DPoRT). 
# Journal of Epidemiology and Community Health, 65(7), 613. 
# https://doi.org/10.1136/jech.2009.102244. Additionally, SAS code used as 
# the basis for this R script may be found here: 
# http://www.rrasp-phirn.ca/index.php?option=com_content&view=article&id=305%3Adport-diabetes-population-risk-tool-&catid=41%3Aother-tools&Itemid=75&lang=en

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

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

start.time <- Sys.time()
# Call DPLYR, TIDYR, READXL, XLSX, RESHAPE2, and TIBBLE
library(dplyr); library(tidyr); library(readxl); library(reshape2); 
library(tibble); library(ggplot2); library(data.table); library(openxlsx)
library(magrittr); library(here)

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

# Load estimates of gross individual employment income derived from analysis of 2011
# NHS data:
load(paste(RDA, "NHS_2011_Calculations", LM, "_1", ".Rda", sep = ""))

# Assiging Key Parameter Values #1: DPoRT Parameters-------------------------------------------

# NOTE 122218: These new parameters Rec_CI and Exc_CI enable the user to
# directly manipulate the definition of recommended and excessive caloric 
# consumption, which, by default, correspond to the level of caloric intake 
# required to sustain a stable body weight when a person is classified as "low 
# active" or "very active", respectively. Although the user could previously
# have modified this on their own, it was buried deep in the code and 
# difficult to find.
EER_List <- c("EER.Sedentary", "EER.Low_Active", "EER.Active", "EER.Very_Active")
Rec_CI = EER_List[1]; Exc_CI = EER_List[3]

# These parameters are drawn directly from the original paper. The first set
# are specific to males. It is important to note that our analysis will focus
# primarily on the influence of gender, ethnicity, education, and income; as 
# such, although we include the estimated values of the parameters here, they 
# do not really come into play. More concretely, we assume that individuals do
# not have hypertension or heart disease (although they may develop it later),
# and are non-smokers. This is consistent with the approach taken in the
# remainder of the thesis.
Const.M = 10.3062; Hypert.M = -0.3629; NC.M = -0.5697; HeartD.M = -0.3483
Smoke.M = -0.0585; PostSec.M = 0.1884; TopQuin.M = 0.1173; Scale.M = 0.7994
# We insert the parameters associated with BMI/age into a vector. This was 
# done because it was perceived that it would simplify referencing later in 
# the code (choosing seperate names for each variable was found to be 
# confusing). Note that the first entry is zero, reflecting the fact that
# BMI<23 & Age<45 is the base category in the regressions formulated by the
# authors of the original paper.
BMI_Age.M = c(0.0000, -0.5520, -0.9521, -1.7162, -2.3310, 
        -1.3602, -1.6537, -2.0563, -2.5513, -2.9353)

# The next set of parameters are specific to females. It should be noted
# that some of them are different than the ones for males (i.e., not just
# the values, but the parameters themselves), implying that different risk
# factors may be important in the development of T2DM for men and women.
Const.F = 10.5777; Hypert.F = -0.4098; NC.F = -0.4528; Immigr.F = -0.1477;
PostSec.F = 0.1939; Scale.F = 0.8419

# We note that the original paper also estimated parameters for individuals
# from each age category for individuals with missing BMI observations. These
# are not pertinent to our analysis, and so are not included here.
BMI_Age.F = c(0.0000, -0.7432, -1.1521, -1.8479, -2.0562,
             -0.7100, -1.2338, -1.8357, -2.3742, -2.6631,
             -1.5956, -1.6144, -1.9830, -2.2148, -2.6448)

# We set out the row labels ahead of time. These are intended primarily to
# improve the readability of the tables we will produce.
Row_Labels = c("BMI<23 & Age<45", "BMI 23-25 & Age<45", "BMI 25-30 & Age<45",
               "BMI 30-35 & Age<45", "BMI>=35 & Age<45", "BMI<23 & Age 45-64", 
               "BMI 23-25 & Age 45-64", "BMI 25-30 & Age 45-64",
               "BMI 30-35 & Age 45-64", "BMI>=35 & Age 45-64", 
               "BMI<23 & Age>=65", "BMI 23-25 & Age>=65", 
               "BMI 25-30 & Age>=65", "BMI 30-35 & Age>=65", 
               "BMI>=35 & Age>=65")

# The third set of parameters are used during the final stage of the
# calculations, when we estimate the likelihood of T2DM. 365.25 is evidently
# a reference to the average number of days in a year. The original paper
# estimated 9-year predicted risk of developing T2DM; we, however, are
# interested in the annual risk of being diagnosed with the condition, and
# modify the parameter on YRS accordingly.
DAYS <- 365.25; YRS <- 1

# Assigning Key Parameter Values #2: Other Parameters ---------------------

# DRI levels originate from: Otten, J. J., Hellwig, J. P., & Meyers, L. D. 
# (Eds.). (2006). Dietary reference intakes: DRI; the essential guide to 
# nutrient requirements. Washington, DC: National Academies Press.
DRI.m = c(1.00, 1.11, 1.25, 1.48)
DRI.f = c(1.00, 1.12, 1.27, 1.45)
PA_Levels = c("Sedentary", "Low_Active", "Active", "Very_Active")
Ethnicities = c("Afro-Caribbean","Asian-Indian","Caucasian")
Genders = c("Male","Female")

# These parameters are again derived directly from Otten et al. (2006), and
# are used to estimate EERs for (non-pregnant or lactating) adults aged 19
# or older.
Cons.m = 662; Age.m = -9.53; Weight.m = 15.91; Height.m = 539.6;
Cons.f = 354; Age.f = -6.91; Weight.f = 9.36; Height.f = 726;
# The next parameter is starting BMI, the level of which (i.e., 23 kg/m^2)
# falls within the normal range.
Starting_BMI = 22.9999
# The next three parameters are estimates of height for males and females of
# varying ethnicities. Since anthropometric data like this appears to be
# difficult to find for Canada (although the Canadian Health Measures Survey
# may merit some additional investigation), we used data from the National 
# Health and Nutrition Examination Survey (NHANES):
# Fryar, C. D., Gu, Q., Ogden, C. L., & Flegal, K. M. (2016, August 3). 
# Anthropometric Reference Data for Children and Adults: United States, 
# 2011-2014. Vital and Health Statistics, 3(39). Retrieved from 
# https://www.cdc.gov/nchs/data/series/sr_03/sr03_039.pdf
Heights.Caucasian = c(1.771,1.629); Heights.Afro_Caribbean = c(1.764,1.630)
Heights.Asian_Indian = c(1.703,1.570)
# Now we stack these on top of each other using rbind, and package the result
# into a matrix.
Heights <- data.frame(rbind(Heights.Afro_Caribbean,Heights.Asian_Indian,
                            Heights.Caucasian))
colnames(Heights) <- Genders

# The equations published by Otten et al. (2006) incorporate age as a 
# continuous variable. However, we are working with age ranges. Therefore, for
# each age range we derive a single value - the simple average of the two end-
# points of the range.
Ages <- c("Age_19_44","Age_45_64","Age_65_95")
Ages_DPoRT <- c("Age<45","Age 45-64","Age>=65")
Point_Age <- as.numeric(c(((44-19)/2)+19,((64-45)/2)+45,((95-65)/2)+65))

# This bit of code generates a data frame, the row entries of which are based
# on all possible combinations and permutations of vectors defined by the
# user. This is MUCH faster than the alternative.
Pr_Risk <- expand.grid(Gender = Genders, Age_Group = Ages,
                      Ethnicity = Ethnicities,
                      Caloric_Intake = c("Recommended","Excessive"))

# We now generate three lookup tables we will need later on. For the first
# and third of these, the intent here is to generate "cross-walks" between 
# the categories used in the EER tables and the CCHS 2004, and those used in 
# the DPoRT. Eventually, we will use the table to generate strings we will 
# subsequently employ to reference predictions generated by the DPoRT (but 
# more on this later).
Age_Lookup <- data.frame(cbind(Ages,Point_Age,Ages_DPoRT))
Age_Lookup$Point_Age <- as.numeric(as.character(Point_Age))
Age_Lookup$Ages_DPoRT <- as.character(Age_Lookup$Ages_DPoRT)
Height_Lookup <- data.frame(cbind(Ethnicities,Heights))
rownames(Height_Lookup) <- c()

# The lookup table we construct next will enable us to link the BMI
# predictions obtained using the various information sources outlined above
# to string fragments which, when combined with similar fragments associated
# with age category, will enable us to reference DPoRT-generated predictions
# of the annual risk of developing T2DM.
BMI_Lookup_Calc <- c(0,23,25,30,35)
BMI_Lookup_DPoRT <- c("BMI<23 &","BMI 23-25 &","BMI 25-30 &",
                      "BMI 30-35 &", "BMI>=35 &")
BMI_Lookup <- data.frame(cbind(BMI_Lookup_Calc,BMI_Lookup_DPoRT))
# The "as.numeric(as.character(BMI_Lookup$BMI_Lookup_Calc))" line of code is
# actually critically important. For some reason, R originally regarded these
# as factors, such that when I attempted to reference this column later on,
# the functions used returned the factor level, rather than the value with
# which it was associated. This code converts the column into a vector of 
# numeric values, so that this doesn't happen. Note that (again, for reasons
# which are not altogether clear), it is not sufficient to simply specify 
# "as.numeric" - the operation in parentheses (i.e., "as.character") is also
# necessary. It may in the future be worthwhile researching why this is so.
BMI_Lookup$BMI_Lookup_Calc <- 
  as.numeric(as.character(BMI_Lookup$BMI_Lookup_Calc))
BMI_Lookup$BMI_Lookup_DPoRT <- as.character(BMI_Lookup$BMI_Lookup_DPoRT)

# Generating the Tables: DPoRT ---------------------------------------------------

# We now generate a series of vectors that will ultimately be used to
# produce the columns of our data tables. Each vector represents a different
# column, and indicates the combination of personal characteristics being
# considered in each one. For example, under "Headings.M", "c(0,1,0,0,1)"
# denotes estimated probabilities for non-Caucasians with post-secondary
# education. It should be noted that several values in each vector do not
# change, reflecting the fact that this analysis focuses exclusively on 
# gender, ethnicity, and educational attainment. In principle, however, it
# would be relatively straightforward to extend the tables to include all
# possible combinations of individual characteristics.
# 
# Not coincidentally, it is straightforward to index an individual element of
# a vector contained in a list thusly: Headings.M[[x]][y]. This will come in
# very handy later.

Headings.M <- list(c(0,0,0,0,0,0), c(0,1,0,0,0,0), c(0,0,0,0,1,0), c(0,1,0,0,1,0),
                   c(0,0,0,0,0,1), c(0,1,0,0,0,1), c(0,0,0,0,1,1), c(0,1,0,0,1,1))
Headings.F <- list(c(0,0,0,0), c(0,1,0,0), c(0,0,0,1), c(0,1,0,1),
                   c(0,0,0,0), c(0,1,0,0), c(0,0,0,1), c(0,1,0,1))

# Now we create an empty data frame to hold our calculations, adapting the
# approach described here:
# https://stackoverflow.com/questions/10689055/create-an-empty-data-frame
MU.M <- data.frame(C.NoPS.NoTQ.M = numeric(), NC.NoPS.NoTQ.M = numeric(),
                   C.PS.NoTQ.M = numeric(), NC.PS.NoTQ.M = numeric(),
                   C.NoPS.TQ.M = numeric(), NC.NoPS.TQ.M = numeric(),
                   C.PS.TQ.M = numeric(), NC.PS.TQ.M = numeric())
# The transitivity of the "<-" operator enables us to clone the DF we just
# constructed. This is made possible by the fact that each set of three
# tables has exactly the same dimensions.
p.M <- m.M <- MU.M
MU.F <- data.frame(C.NoPS.NoTQ.F = numeric(), NC.NoPS.NoTQ.F = numeric(),
                   C.PS.NoTQ.F = numeric(), NC.PS.NoTQ.F = numeric(),
                   C.NoPS.TQ.F = numeric(), NC.NoPS.TQ.F = numeric(),
                   C.PS.TQ.F = numeric(), NC.PS.TQ.F = numeric())
p.F <- m.F <- MU.F
# Generate a diagonal matrix. This will act as a switch that enables us to
# select the BMI/age parameter corresponding to each row of the table.
DM <- diag(1,15,15)

# This is the core of the calculations, and closely mirrors the published SAS
# code. There is a substantial amount of indexing going on in the parentheses.
# The first entry in the parentheses for the first set of parameters (i.e.,
# those unrelated to BMI/age) draws upon the binary values within a series of
# vectors that are themselves nested in a list (i.e., Headings.M). For the
# second set of parameters, the first entry reflects a value in the diagonal
# matrix we created earlier. The rationale for this approach is that each
# BMI/age parameter applies (i.e., is equal to 1) in exactly one row of the 
# table, and not the others (i.e., it is otherwise equal to 0). The second
# value in each pair of parenthesis references parameters estimated from the 
# CCHS 4.1 PUMF, the interpretation of which is not clearly laid out in the
# original article or in the publicly-available SAS code.
for (j in 1:length(Headings.M)){
  # There are 15 combinations of BMI/age, so why do we only populate the first
  # 10 rows? The answer is that for males, the original analysis estimated
  # parameters for two age categories (i.e., < 45 years, >= 45 years), whereas
  # for females there are three age categories (i.e., < 45 years, 45-64 years,
  # >=65 years). To simplify future coding, we split the ">= 45 years"
  # category into two, paralleling the breakdown for females, and simply apply
  # the same parameter values to each sub-category.
  
  for (i in 1:10){
    MU.M[i,j] = Const.M + Hypert.M*(Headings.M[[j]][1]) + 
      NC.M*(Headings.M[[j]][2]) + HeartD.M*(Headings.M[[j]][3]) + 
      Smoke.M*(Headings.M[[j]][4]) +PostSec.M*(Headings.M[[j]][5]) + 
      TopQuin.M*(Headings.M[[j]][6]) + BMI_Age.M[1]*(DM[i,1]) + 
      BMI_Age.M[2]*(DM[i,2]) + BMI_Age.M[3]*(DM[i,3]) +
      BMI_Age.M[4]*(DM[i,4]) + BMI_Age.M[5]*(DM[i,5]) + 
      BMI_Age.M[6]*(DM[i,6]) + BMI_Age.M[7]*(DM[i,7]) +
      BMI_Age.M[8]*(DM[i,8]) + BMI_Age.M[9]*(DM[i,9]) +
      BMI_Age.M[10]*(DM[i,10]) 
  }
  # Rather than repeating the calculations for the last five rows, we simply
  # duplicate the values appearing in the five preceding rows (the rationale
  # for this is described above).
  for (i in 11:15){
    MU.M[i,j] = MU.M[i-5,j]
  }
}

# To faciliate readability, we now insert the row labels we developed earlier.
rownames(MU.M) <- Row_Labels

for (j in 1:length(Headings.F)){
  for (i in 1:length(BMI_Age.F)){
    MU.F[i,j] = Const.F + Hypert.F*(Headings.F[[j]][1]) + 
      NC.F*(Headings.F[[j]][2]) + Immigr.F*(Headings.F[[j]][3]) + 
      PostSec.F*(Headings.F[[j]][4]) + BMI_Age.F[1]*(DM[i,1]) + 
      BMI_Age.F[2]*(DM[i,2]) + BMI_Age.F[3]*(DM[i,3]) + 
      BMI_Age.F[4]*(DM[i,4]) + BMI_Age.F[5]*(DM[i,5]) + 
      BMI_Age.F[6]*(DM[i,6]) + BMI_Age.F[7]*(DM[i,7]) + 
      BMI_Age.F[8]*(DM[i,8]) + BMI_Age.F[9]*(DM[i,9]) + 
      BMI_Age.F[10]*(DM[i,10]) + BMI_Age.F[11]*(DM[i,11]) + 
      BMI_Age.F[12]*(DM[i,12]) + BMI_Age.F[13]*(DM[i,13]) + 
      BMI_Age.F[14]*(DM[i,14]) + BMI_Age.F[15]*(DM[i,15]) 
  }
}
rownames(MU.F) <- Row_Labels

# The next step in the calculations is to derive values for "m", using the 
# following formula: m = (log(follow - uptime) - MU)/scale
for (i in 1:nrow(MU.M)){
  for (j in 1:ncol(MU.M)){
    m.M[i,j] = (log(DAYS*YRS) - MU.M[i,j])/Scale.M
    m.F[i,j] = (log(DAYS*YRS) - MU.F[i,j])/Scale.F
  }
}
rownames(m.M) <- rownames(m.F) <- Row_Labels

# Finally, we calculate the actual probabilities, using the following formula:
# Predicted risk (P) = 1 - exp(-exp^m)
for (i in 1:nrow(m.M)){
  for (j in 1:ncol(m.M)){
    p.M[i,j] = 1-exp(-exp(m.M[i,j]))
    p.F[i,j] = 1-exp(-exp(m.F[i,j]))
  }
}
rownames(p.M) <- rownames(p.F) <- Row_Labels

# Generating the Tables: Other --------------------------------------------

# We proceed to calculate equilibrium body weights, which are then converted 
# into BMIs, and then finally combined with other individual characteristics
# to reference information from the DPoRT tables generated in the previous
# section of the code. We begin by recalling point ages, mean heights and
# starting weights from the various reference tables constructed above.
Pr_Risk <- Pr_Risk %>%
  mutate(Age_Value =
           Age_Lookup[findInterval(Age_Group,Age_Lookup$Ages),2]) %>%
  mutate(Height_Value =
           ifelse(Gender=="Male",
                  Height_Lookup[findInterval(Ethnicity,
                                             Height_Lookup$Ethnicities),2],
                  Height_Lookup[findInterval(Ethnicity,
                                             Height_Lookup$Ethnicities),3])) %>%
  mutate(Starting_Weight = Starting_BMI*Height_Value^2)

Pr_Risk$EER.Low_Active <- Pr_Risk$EER.Sedentary <- 0
Pr_Risk$EER.Very_Active <- Pr_Risk$EER.Active <- 0

for (j in 1:4){
  Pr_Risk[,7+j] = 
    ifelse(Pr_Risk$Gender=="Male",
           Cons.m + Age.m*Pr_Risk$Age_Value + DRI.m[j]*
             (Weight.m*Pr_Risk$Starting_Weight + Height.m*Pr_Risk$Height_Value),
           Cons.f + Age.f*Pr_Risk$Age_Value + DRI.f[j]*
             (Weight.f*Pr_Risk$Starting_Weight + Height.f*Pr_Risk$Height_Value))
}

# This fairly intimidating bit of code calculates equilibrium body weight in
# kilograms according to gender, caloric intake and level of physical activity
# by combining the EER estimates published by Health Canada, figures on the
# 75th percentile of caloric intake among Manitobans from the 2004 CCHS, EER
# formulas and Physical Activity Coefficients from the National Academy of
# Sciences (i.e., Otten et al. (2006)), and the Variable Calorie Expenditure
# model outlined in: Barton, J. T. (2016). Models for life: an introduction to  
# discrete mathematical modeling with Microsoft Office Excel. Hoboken, New
# Jersey: Wiley.
# 
# The equations themselves are described in greater detail in the body of the
# thesis. The most important thing to acknowledge is that they are not exactly
# the same as the ones Barton (2016) uses, because he relies upon the Mifflin-
# St. Jeor equations for resting energy expenditure (REE).
WT_List <- data.frame("Sedentary.WT" = 1, "Low_Active.WT" = 2, "Active.WT" = 3, 
                      "Very_Active.WT" = 4)

DRI <- data.frame(rbind(DRI.m, DRI.f))
names(DRI) <- c("Sedentary.WT", "Low_Active.WT", "Active.WT", "Very_Active.WT")

for (i in names(WT_List)){Pr_Risk[[i]] <- 1}

for (i in names(WT_List)){
  for (j in 1:nrow(Pr_Risk)){
    Pr_Risk[[i]][j] = 
      with(Pr_Risk,
           ifelse(Gender[j] == "Male",
                  ifelse(Caloric_Intake[j] == "Recommended", 
                         Pr_Risk[[Rec_CI]][j]/(DRI[[i]][1]*Weight.m) -
                           (Height.m*Height_Value[j])/Weight.m - 
                           Cons.m/(DRI[[i]][1]*Weight.m) -
                           (Age.m*Age_Value[j])/(DRI[[i]][1]*Weight.m), 
                         Pr_Risk[[Exc_CI]][j]/(DRI[[i]][1]*Weight.m) - 
                           (Height.m*Height_Value[j])/Weight.m -
                           Cons.m/(DRI[[i]][1]*Weight.m) - 
                           (Age.m*Age_Value[j])/(DRI[[i]][1]*Weight.m)),
                  ifelse(Caloric_Intake[j] == "Recommended", 
                         Pr_Risk[[Rec_CI]][j]/(DRI[[i]][2]*Weight.f) -
                           (Height.f*Height_Value[j])/Weight.f - 
                           Cons.f/(DRI[[i]][2]*Weight.f) -
                           (Age.f*Age_Value[j])/(DRI[[i]][2]*Weight.f),
                         Pr_Risk[[Exc_CI]][j]/(DRI[[i]][2]*Weight.f) - 
                           (Height.f*Height_Value[j])/Weight.f -
                           Cons.f/(DRI[[i]][2]*Weight.f) - 
                           (Age.f*Age_Value[j])/(DRI[[i]][2]*Weight.f))))
  }
}

Pr_Risk <- Pr_Risk %>%
  # Now we just need to convert the values we just derived back into BMIs.
  mutate(Sedentary.BMI = Sedentary.WT/Height_Value^2) %>%
  mutate(Low_Active.BMI = Low_Active.WT/Height_Value^2) %>%
  mutate(Active.BMI = Active.WT/Height_Value^2) %>%
  mutate(Very_Active.BMI = Very_Active.WT/Height_Value^2)

# Now we sort the columns, purely to facilitate presentation.
Pr_Risk <- Pr_Risk %>%
  arrange(Caloric_Intake,Gender,Age_Group,Ethnicity)

# We would like to be able to draw values from the DPoRT tables directly into
# this table. I opted to do this by concatenating text strings containing
# information about BMI and age category. These strings then serve as the
# basis for the lookup by matching on the first row of the DPoRT tables. Note
# that the concatenation is undertaken by applying the "paste" function. Also
# note that when this function is used, R by default inserts a space between
# the strings that are being merged together. Originally, I inserted a space
# myself at the end of the first string, which caused the matching process to
# fail (i.e., because there were now two spaces between the strings that
# were being concatenated, and only one in the source being referenced).
Pr_Risk <- Pr_Risk %>% 
  mutate(Sedentary.BMI_Cat = paste(BMI_Lookup[findInterval(
    Sedentary.BMI,BMI_Lookup$BMI_Lookup_Calc),2],
    Age_Lookup[findInterval(Age_Group,Age_Lookup$Ages),3])) %>%
  mutate(Low_Active.BMI_Cat = paste(BMI_Lookup[findInterval(
    Low_Active.BMI,BMI_Lookup$BMI_Lookup_Calc),2],
    Age_Lookup[findInterval(Age_Group,Age_Lookup$Ages),3])) %>%
  mutate(Active.BMI_Cat = paste(BMI_Lookup[findInterval(
    Active.BMI,BMI_Lookup$BMI_Lookup_Calc),2],
    Age_Lookup[findInterval(Age_Group,Age_Lookup$Ages),3])) %>%
  mutate(Very_Active.BMI_Cat = paste(BMI_Lookup[findInterval(
    Very_Active.BMI,BMI_Lookup$BMI_Lookup_Calc),2],
    Age_Lookup[findInterval(Age_Group,Age_Lookup$Ages),3])) 

# The next step is to duplicate the DPoRT tables, since we will be 
# manipulating them somewhat (and because the original tables will ultimately
# be exported into Excel).
p.F_Link <- p.F; p.M_Link <- p.M
# As we did previously, we now convert the rownames into column values using a
# special function in the DPLYR package. This will give us a field from the
# DPoRT data that we can match on.
p.F_Link <- rownames_to_column(p.F_Link, var="BMI_Age")
# p.F_Link <- p.F_Link[order(p.F_Link$BMI_Age),]
p.M_Link <- rownames_to_column(p.M_Link, var="BMI_Age")
# p.M_Link <- p.M_Link[order(p.M_Link$BMI_Age),]

# Finally, we link the equilibrium BMIs to the DPoRT transition probabilities.
# Since the DPoRT model recognizes an association between educational
# attainment (i.e., whether they do or do not have any post-secondary
# education) and the risk of T2DM, we generate four pairs of sets of TPs. Note
# that we use "match" rather than "findInterval" to perform the lookup. This
# was done because it appears that "findInterval" does not work with strings.
# In fact, it appears it was designed to work only with numeric values, so it
# is somewhat surprising it functioned properly when a lookup was carried out
# using the Age_Lookup table. One explanation for this is that R used the
# factor values as the basis for the lookup, rather than relying upon the
# strings with which they are associated; this theory may, however, merit
# further examination. Finally, it should be noted that the code below is
# also somewhat unwieldy, and that there may be opportunities to condense it.

Pr_Risk <- Pr_Risk %>%
  mutate(C_OR_NC = ifelse(Ethnicity=="Caucasian", 1, 2)) %>%
  mutate(Sedentary_NoPS_NoTQ_Risk = NA, Sedentary_PS_NoTQ_Risk = NA, 
         Sedentary_NoPS_TQ_Risk = NA, Sedentary_PS_TQ_Risk = NA, 
         Low_Active_NoPS_NoTQ_Risk = NA, Low_Active_PS_NoTQ_Risk = NA, 
         Low_Active_NoPS_TQ_Risk = NA, Low_Active_PS_TQ_Risk = NA,
         Active_NoPS_NoTQ_Risk = NA, Active_PS_NoTQ_Risk = NA,
         Active_NoPS_TQ_Risk = NA, Active_PS_TQ_Risk = NA,
         Very_Active_NoPS_NoTQ_Risk = NA, Very_Active_PS_NoTQ_Risk = NA,
         Very_Active_NoPS_TQ_Risk = NA, Very_Active_PS_TQ_Risk = NA)

BMI_Cat_List <- list("Sedentary.BMI_Cat", "Low_Active.BMI_Cat", "Active.BMI_Cat",
                     "Very_Active.BMI_Cat")
Link_List <- list("Male" = p.M_Link, "Female" = p.F_Link)
# Don't need to create a separate list for genders, as we already have a vector
# containing this information.
# 
# Order:
# 1) Level of physical activity
# 2) Income level (which is binary, since a person either is, or is not in the
# top income quintile).
# 3) PSE (which is also binary, since a person either has, or has not attended 
# PSE).
for (i in 1:length(BMI_Cat_List)){
  for (j in 1:2){
    for (k in 1:2){
      for (n in 1:nrow(Pr_Risk)){
        Pr_Risk[n, 24 + length(BMI_Cat_List)*(i-1) + 2*(j-1) + k] <-
          Link_List[[Pr_Risk$Gender[n]]][match(Pr_Risk[n, BMI_Cat_List[[i]]],
                                               # NOTE 060818: When I tried to reference
                                               # the BMI_Age column by number (i.e., 
                                               # [1]), the match uniformly failed. It
                                               # is evidently necessary to refer to
                                               # column name by label instead.
                                               Link_List[[Pr_Risk$Gender[n]]]$BMI_Age),
                                         1 + 4*(j-1) + 2*(k-1) + Pr_Risk$C_OR_NC[n]]
      }}}}

save(Pr_Risk, file = paste(RDA, Filename, LM, ".Rda", sep = ""))
write.xlsx(Pr_Risk, file = paste0(UKPDS_OM2, "Pr_Risk_NEW", LM, ".xlsx"))

# For the tables incorporated into the report.
RPT_TBL <- Pr_Risk[, c(1:4, 25:40)]
RPT_TBL %<>%
  mutate(Chars = paste(Ethnicity, Gender, Caloric_Intake, sep = "_")) %>%
  select(c(2, 21, 5:20))

RPT_TBL %<>% 
  melt(., id = c("Age_Group", "Chars"), variable.names = c(3:18),
       value.name = "Prob") %>%
  dcast(., Age_Group + variable ~ Chars, value.var = "Prob")

RPT_TBL %<>% .[, c(1:2, 14:11, 6:3, 10:7)]
RPT_TBL$variable %<>%
  gsub("Low_Active", "LowActive", .) %>%
  gsub("Very_Active", "VeryActive", .)
  
RPT_TBL %<>% separate(., variable, into = c("LTPA", "PS", "TQ", "Risk"), sep = "_")
RPT_TBL$LTPA <- as.factor(RPT_TBL$LTPA)
RPT_TBL$PS <- as.factor(RPT_TBL$PS)
RPT_TBL$TQ <- as.factor(RPT_TBL$TQ)
# https://www.r-bloggers.com/reorder-factor-levels/
RPT_TBL$LTPA = factor(RPT_TBL$LTPA, levels(RPT_TBL$LTPA)[c(3,2,1,4)])
RPT_TBL %<>% arrange(., TQ, PS, Age_Group, LTPA)
RPT_TBL[2:5] %<>% sapply(., as.character)
RPT_TBL %<>% mutate(variable = paste(LTPA, PS, TQ, Risk, sep = "_"))
RPT_TBL %<>% .[, c(1, 18, 6:17)]

RPT_TBL$variable %<>%
  gsub("LowActive", "Low_Active", .) %>%
  gsub("VeryActive", "Very_Active", .)

save(RPT_TBL, file = paste0(RDA, Filename, LM, "_One_Year_Probs", ".Rda"))

# For the tables incorporated into the Meal Plan template.
RPT_TBL2 <- Pr_Risk[Pr_Risk$Caloric_Intake == "Recommended", c(1:4, 8:11)]
RPT_TBL2$Caloric_Intake <- NULL
RPT_TBL2 %<>%
  mutate(Chars = paste(Ethnicity, Gender, sep = "_")) %>%
  select(c(2, 8, 4:7))

RPT_TBL2 %<>% 
  melt(., id = c("Age_Group", "Chars"), variable.names = c(3:6),
       value.name = "Prob") %>%
  dcast(., Age_Group + variable ~ Chars, value.var = "Prob")

RPT_TBL2 %<>% .[, c(1:2, 4:3, 6:5, 8:7)]

write.xlsx(RPT_TBL2, file = paste0(UKPDS_OM2, "EERs", LM, ".xlsx"))

# Foundation for UKPDS Data Template Builder -----------------------------------

# 032519: This code uses the calculations employed to generate the equilibrium
# BMIs to fill in the BMI tab of the UKPDS-OM2 template ("C:\Users\starg\Documents\ 
# PhD Dissertation\Model Calibration\Transition Probabilities\UKPDS\Dissertation\
# UKPDS_Data_Template_V5.xlsm"). This content is relatively new, and is intended to
# ensure a consistent approach to the calibration of the decision model across the
# modules drawing from the DPoRT 2.0 Risk Equations and the UKPDS-OM2 simulation
# results.

Ethnicity = c("White", "Afro-Caribbean", "Asian-Indian")
Gender = c("Male", "Female"); Age = c(27:95)
UKPDS.00 <- expand.grid(Ethnicity, Gender, Age)
names(UKPDS.00) <- c("Ethnicity", "Gender", "Age")
UKPDS.01 <- data.frame(matrix(NA, ncol = 70, nrow = nrow(UKPDS.00)))
names(UKPDS.01) <- paste("Year", seq(1:70), sep = " ")

UKPDS_Age_Lookup <- data.frame(Age = seq(from = 25, to = 170))
UKPDS_Age_Lookup <- UKPDS_Age_Lookup %>%
  mutate(Age_Group = case_when(
    Age < 45 ~ "Age_19_44",
    Age < 65 ~ "Age_45_64",
    TRUE ~ "Age_65_95"
  ))

# 081718: Note that UKPDS.02 is effectively an empty dataframe, serving as
# a template for UKPDS_Temp - which will actually contain the calculated
# values used to calibrate the UKPDS-OM2.
UKPDS.02 <- merge(UKPDS.00, UKPDS_Age_Lookup, by = "Age")
UKPDS.02 <- cbind(UKPDS.02, UKPDS.01)

levels(UKPDS.02$Ethnicity) <- c("White", "Afro-Caribbean", "Asian-Indian")
UKPDS.02 %<>% arrange(., Gender, Ethnicity, Age)
# https://stackoverflow.com/questions/19508256/how-to-add-new-column-to-an-dataframe-to-the-front-not-end
UKPDS.02 <- add_column(UKPDS.02, ID = rownames(UKPDS.02), .after = "Age_Group")

UKPDS_Temp <- vector(mode = "list", length = 64)
for (i in 1:length(UKPDS_Temp)){UKPDS_Temp[[i]] = UKPDS.02}

# 032519: This previously linked to a location on the author's C: drive, which is not
# inherently problematic, but would have prevented others from using the code without
# implementing modifications.
UKPDS_String <- 
  read_xlsx(paste0(UKPDS_OM2, "/Decisions.xlsx"), "Code", range = "K1:K65")

# 063018: This doesn't work properly. For some reason, the first part of the
# string is omitted. I cannot think of a logical explanation for this.
# for (i in 1:length(UKPDS_Temp)){names(UKPDS_Temp)[[i]] <- 
# assign(paste(c("UKPDS_", i), collapse = ""), UKPDS_Temp_Names.00$UKPDS_String[i])}

# 070118: This code DOES work properly. The AS.CHARACTER function converts the value
# from UKPDS_string into a string, whereupon it can be readily combined with the other
# string we specified.
for (i in 1:length(UKPDS_Temp)){names(UKPDS_Temp)[[i]] <- 
  paste("UKPDS_", as.character(UKPDS_String$CODE[i]), sep = "")}

Pr_Risk.01 <- Pr_Risk %>% select(., c(1:4, 12:15))
Pr_Risk.01 <- data.table(Pr_Risk.01)
Pr_Risk.01 <- melt(Pr_Risk.01, id = c(1:4), measure = c(5:ncol(Pr_Risk.01)), 
                   value.name = "Weight", variable.name = "LTPA")
# https://stackoverflow.com/questions/11121385/repeat-rows-of-a-data-frame
Pr_Risk.01 <- Pr_Risk.01[rep(seq_len(nrow(Pr_Risk.01)), each=2),]
rownames(Pr_Risk.01) <- seq(1:nrow(Pr_Risk.01))
levels(Pr_Risk.01$Ethnicity) <- c("White", "Afro-Caribbean", "Asian-Indian")

Pr_Risk.01 <- Pr_Risk.01 %>% 
  mutate(ID = seq(1:nrow(Pr_Risk.01))) %>%
  mutate(Med_Adherence = ifelse(ID %% 2 == 0, 
                                "Adhere", 
                                "Don't Adhere")) %>%
  arrange(Med_Adherence, Caloric_Intake, LTPA, Gender, Age_Group, Ethnicity)

Pr_Risk.01 <- Pr_Risk.01[, c(7, 1:3, 5, 4, 8, 6)]
Pr_Risk.01$Caloric_Intake <- as.character(Pr_Risk.01$Caloric_Intake)
Pr_Risk.01$LTPA <- as.character(Pr_Risk.01$LTPA)
Pr_Risk.01$Med_Adherence <- as.character(Pr_Risk.01$Med_Adherence)

Combinations <- expand.grid(LTPA = c("Active.WT", "Sedentary.WT"),
                       Caloric_Intake = c("Recommended","Excessive"),
                       Med_Adherence = c("Adhere", "Don't Adhere"), stringsAsFactors = F)

for (i in 1:nrow(Combinations)){
  for (j in 1:nrow(Combinations)){
    for (k in 1:nrow(UKPDS.02)){
      for (l in 1:8){
        UKPDS_Temp[[(i-1)*8+j]][k, l + 5] = 
          Pr_Risk.01[
            Pr_Risk.01$Ethnicity==UKPDS_Temp[[(i-1)*8+j]]$Ethnicity[k] & 
              Pr_Risk.01$Gender==UKPDS_Temp[[(i-1)*8+j]]$Gender[k] &
              Pr_Risk.01$Age_Group==UKPDS_Age_Lookup[UKPDS_Age_Lookup==(UKPDS_Temp[[(i-1)*8+j]]$Age[k]+l),2] &
              Pr_Risk.01$LTPA==Combinations$LTPA[j] & 
              Pr_Risk.01$Caloric_Intake==Combinations$Caloric_Intake[j] &
              Pr_Risk.01$Med_Adherence==Combinations$Med_Adherence[j], 8]
      }
      for (l in 9:70){
        UKPDS_Temp[[(i-1)*8+j]][k, l + 5] = 
          Pr_Risk.01[
            Pr_Risk.01$Ethnicity==UKPDS_Temp[[(i-1)*8+j]]$Ethnicity[k] & 
              Pr_Risk.01$Gender==UKPDS_Temp[[(i-1)*8+j]]$Gender[k] &
              Pr_Risk.01$Age_Group==UKPDS_Age_Lookup[UKPDS_Age_Lookup==(UKPDS_Temp[[(i-1)*8+j]]$Age[k]+l),2] &
              Pr_Risk.01$LTPA==Combinations$LTPA[i] & 
              Pr_Risk.01$Caloric_Intake==Combinations$Caloric_Intake[i] &
              Pr_Risk.01$Med_Adherence==Combinations$Med_Adherence[i], 8]
      }
    }
  }
}

# Order the list by name (we'll need this later).
UKPDS_Temp <- UKPDS_Temp[order(names(UKPDS_Temp))]

# https://stackoverflow.com/questions/27713310/r-easy-way-to-export-multiple-data-frame-to-multiple-excel-worksheets
# This is a great function (associated with the OPENXLSX package) - it enables the user
# to copy an arbitrarily large number of dataframes to Excel as separate worksheets. 
# NOTE: Changes to the file path presented here MUST be reflected in the macros in the
# data template (see file path at the beginning of this section) for the latter to work
# properly.
write.xlsx(UKPDS_Temp, file = paste0(UKPDS_OM2, "BMI_Calcs", LM, ".xlsx"))

end.time <- Sys.time()
duration <- end.time - start.time