# R Code for Prediction of Presenteeism Attributable to T2DM
# Created: 02/04/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 020418: The calculations in this code are based on results presented in:
# Lavigne, J. E., Phelps, C. E., Mushlin, A., & Lednar, W. M. (2003). 
# Reductions in individual work productivity associated with type 2 
# diabetes mellitus. Pharmacoeconomics, 21(15), 1123-1134.

# 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(tibble); 
library(reshape2); library(here); library(openxlsx)

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

# Assiging Key Parameter Values -------------------------------------------

# These parameters are drawn directly from the original paper. I have elected
# not to include parameters that were determined not to be statistically
# significant at p < 0.10
DUR_T2DM <- 1.08; DUR_T2DM_2 <- -0.04; Male <- -8.27; Male_T2DM <- 12.77
Minority <- -7.29; Minority_T2DM <- 18.76; PS_Degree_T2DM <- -11.85
# Introduce a series of parameters that reflect the impact of the patient's 
# T2DM management activities; the parameters will be stored in a vector.
MGMT = c(1.00, 0.50, 1.50)

# Generating the Tables ---------------------------------------------------

# 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

WEL_1.00 <- data.frame(C.NoDeg.M = numeric(), C.NoDeg.F = numeric(),
                 C.Deg.M = numeric(), C.Deg.F = numeric(),
                 NC.NoDeg.M = numeric(), NC.NoDeg.F = numeric(),
                 NC.Deg.M = numeric(), NC.Deg.F = numeric())  

# In the year of diagnosis, our assumption is that the individual has not
# yet begun to incur productivity losses; thus, all values are set to 0.
WEL_1.00[1,] = 0

# Interestingly, the "<-" operator is transitive. We can take advantage of
# this to "clone" the original dataframe.
WEL_1.50B <- WEL_1.50A <- WEL_0.50B <- WEL_0.50A <- WEL_1.00

# We then generate a list containing each of the dataframes. This is done
# because it enables us to automate the next set of calculations by indexing
# each element of the list in turn.
Scenarios <- list(WEL_1.00, WEL_0.50A, WEL_1.50A, WEL_0.50B, WEL_1.50B)

# Now we use the parameters reported in the paper and presented above to
# predict productivity losses for individuals with all possible combinations
# of gender, ethnicity, and educational attainment. We do this using the
# original values of the parameters, and again for the two scenarios that
# reflect "good" or "poor" management of T2DM.
# 
# NOTE 020718: The way we do this below is clearly not the most elegant 
# approach possible, but will do until a superior alternative is identified.

for (h in 1:3){
  # First, we index element h from the list we compiled earlier (each element
  # of which is a dataframe). Then we perform the calculations on each one. It
  # is noteworthy that it is surprisingly straightforward to apply
  # manipulations to the dataframes included in the list above. In particular,
  # R somehow knows that "Scenarios[[h]]" refers to a dataframe, and therefore
  # substitutes in directly for the name of that dataframe. This is immensely
  # useful when, as in this case, we are trying to apply code to a series of
  # dataframes contained in a list (or some other structure) in a loop.
    for (dur in 2:76){
    Scenarios[[h]][dur,"C.NoDeg.M"] = MGMT[h]*DUR_T2DM*dur + 
      DUR_T2DM_2*dur^2 + (Male + Male_T2DM)
    Scenarios[[h]][dur,"C.NoDeg.F"] = MGMT[h]*DUR_T2DM*dur + DUR_T2DM_2*dur^2
    Scenarios[[h]][dur,"C.Deg.M"] = MGMT[h]*DUR_T2DM*dur + DUR_T2DM_2*dur^2 + 
      (Male + Male_T2DM) + PS_Degree_T2DM
    Scenarios[[h]][dur,"C.Deg.F"] = MGMT[h]*DUR_T2DM*dur + DUR_T2DM_2*dur^2 + 
      PS_Degree_T2DM
    Scenarios[[h]][dur,"NC.NoDeg.M"] = MGMT[h]*DUR_T2DM*dur + 
      DUR_T2DM_2*dur^2 + (Male + Male_T2DM) + (Minority + Minority_T2DM)
    Scenarios[[h]][dur,"NC.NoDeg.F"] = MGMT[h]*DUR_T2DM*dur + 
      DUR_T2DM_2*dur^2 + (Minority + Minority_T2DM)
    Scenarios[[h]][dur,"NC.Deg.M"] = MGMT[h]*DUR_T2DM*dur + DUR_T2DM_2*dur^2 + 
      (Male + Male_T2DM) + PS_Degree_T2DM + (Minority + Minority_T2DM)
    Scenarios[[h]][dur,"NC.Deg.F"] = MGMT[h]*DUR_T2DM*dur + DUR_T2DM_2*dur^2 + 
      PS_Degree_T2DM + (Minority + Minority_T2DM)
  }
}

# Two things become immediately apparent when we run the code above: 
# 1) The value of the parameter associated with DUR_T2DM_2 is such that at 
# some point, productivity losses associated with T2DM themselves begin to 
# decline and are ultimately converted into productivity gains. There is no 
# reason to expect that this should occur in reality. 
# 2) Second, the parameter associated with having a college or graduate degree
# and T2DM (i.e., PS_Degree_T2DM) is sufficiently large that Caucasians of both
# genders only ever incur negative productivity losses (i.e., productivity 
# gains), which, again, is not a meaningful result.
# 
# Thus, we adapt the results in the following way:
# 1) For individuals who never experience productivity losses, set all values
# to zero. Intuitively, higher education in some sense insulates them from
# work efficiency losses.
# 2) For individuals who at some point do experience productivity losses, all
# data points following the largest recorded value are set equal to the 
# largest recorded value. Intuitively, we assume that work efficiency losses
# ultimately plateau. It is admittedly unclear to what extent this approach is
# appropriate, since it seems at least as plausible that productivity losses
# should continue to mount indefinitely. Nonetheless, it is supported to some
# extent by the data. It should be noted that this set of individuals fall
# into two groups:
#   a) People who initially (i.e., before adjustment) experience productivity
#   losses, and these losses accumulate over time before reaching a maximum.
#   The process outlined above applies here without any need for modification.
#   b) People who initially (i.e., before adjustment) experience productivity
#   gains, but subsquently incur losses. In this case we apply one change to 
#   the process described above. In particular, negative values in the rows 
#   preceding the maximum for each column are set to zero (negative values
#   in all subsequent rows are already addressed by employing the approach
#   described in 2)
# 
# This next code segment implements the algorithm outlined above.

for (h in 1:3) {
  # First, we index element h from the list we compiled earlier (each element
  # of which is a dataframe), and assign it to the temporary dataframe, "df".
  df <- Scenarios[[h]]
  # The next "layer" of for loops actually implements the algorithm.
    for (j in 1:ncol(df))
    {
      # This says that if the value in the second row of each column is negative and
      # the maximum value in the column is negative, set all values in that column
      # to zero.
      if (df[2,j]<0 & max(df[,j])<0)
      {df[,j]=0}
      for (i in 1:nrow(df))
      {
        # Now we deal with the more complicated cases. If the current row index is
        # larger than the one associated with the largest value in a given column,
        # set the value in each row equal to the maximum (i.e., if it is less than
        # this value, it means that we are past the maximum and that WELs have begun
        # to decline).
        if (i>which.max(df[,j]))
        {df[i,j]=max(df[,j])}
        else
          # If the current row index is smaller than the one associated with the largest
          # value in a given column, there are two possibilities to consider. If the
          # value in this row is positive or equal to 0, we leave it as is. If it is
          # negative, however, we set it equal to zero.
          if (i<which.max(df[,j]))
          {
            if (df[i,j]<0)
            {df[i,j]=0}
          }
      }
      # We note that we have disregarded a third possibility, namely one in which
      # the current row index is exactly equal to the one associated with the
      # largest value. In this case, no adjustment is required.
    }
  # Finally, we take the transformed dataframes, and insert them back into the
  # list generated earlier.
  Scenarios[[h]] <- df
}

# To prepare the tables for export into Excel, we withdraw each dataframe from
# the list, and assign it an appropriate name.
WEL_1.00A <- Scenarios[[1]]
WEL_0.50A <- Scenarios[[2]]
WEL_1.50A <- Scenarios[[3]]

# Compiling the next two tables (i.e., WEL_0.50B and WEL_1.50B) is a little 
# trickier. These tables examine the impact of modifying one's approach to 
# T2DM management after having lived with the behaviour for the first
# decade. We assume that the impact of the change in behaviour, should the 
# change occur, would be observed beginning in Year 10 (which is actually the
# interval between Years 10 and 11). There are two possibilities here:
#   1) The individual was previously managing T2DM well, but now chooses to
#   be less vigilant. In this case, we assume the value of the duration
#   parameter increases from 0.50 to 1.50, beginning in the 11th year, and
#   that WELs then accumulate until they reach the level they WOULD have
#   reached if the individual had always managed their condition poorly.
#   2) The individual was previously managing T2DM poorly, but now chooses to
#   manage it more aggressively. In this case, we assume the value of the
#   duration parameter decreases from 1.50 to 0.50. We further assume that
#   any WELs accumulated so far are permanent, which limits the impact of the
#   the change in behaviour.
WEL_0.50B <- WEL_0.50A
for (j in 1:ncol(WEL_0.50B))
{
  # Recall that the change in behaviour (if it happens) occurs in the 11th 
  # year (which, again, is classified as Year 10). There are two possible 
  # situations to deal with. The first is that maximum WELs are experienced
  # prior to when the behavioural shift takes place. In this event, we use 
  # that value as the baseline, and then apply the new, higher rate of WEL 
  # growth to every subsequent year. Note that we start the "counter" back at
  # one, rather than basing the calculations on the integer associated with
  # the current year. The rationale for this is that we assume WELs should
  # rise to the level they would otherwise have attained (i.e., if the
  # individual had always managed their condition poorly) - merely at a later
  # time. Also note that the function of the "min" operator is to avoid a
  # situation in which WELs rise above the levels seen in the counterfactual.
  if (which.max(WEL_0.50A[,j])<=10)
    for (i in 11:nrow(WEL_0.50B)){
      WEL_0.50B[i,j] = min((max(WEL_0.50A[,j]) + MGMT[3]*DUR_T2DM*(i-10) + 
        DUR_T2DM_2*(i-10)^2),max(WEL_1.50A[,j]))
    }
    else
      # In the event that maximum WELs had not yet been achieved by Year 10,
      # the Year 9 value serves the baseline.
      for (i in 11:nrow(WEL_0.50B)){
        WEL_0.50B[i,j] = 
          min((max(WEL_0.50A[1:10,j]) + MGMT[3]*DUR_T2DM*(i-10) + 
                                DUR_T2DM_2*(i-10)^2),max(WEL_1.50A[,j]))
      }
}
Scenarios[[4]] <- WEL_0.50B

# In the second of the two tables examining the implications upon WELs of Year
# 11 behavioural changes, we look at what happens if a person who did not
# manage their condition well now becomes more vigilant.
WEL_1.50B <- WEL_1.50A
for (j in 1:ncol(WEL_1.50B))
{
  # The first possibility is that maximum WELs are achieved prior to the end
  # of Year 10. We assume that the progressive nature of T2DM and the 
  # irreversible nature of some complications (i.e., organ and tissue damage)
  # manifests in permanent loss of WELs. As such, if maximum WELs were
  # already achieved by the time the behavioural change takes place, it has
  # no real impact on individual work productivity.
  if (which.max(WEL_1.50A[,j])<=10)
    for (i in 11:nrow(WEL_1.50B)){
      WEL_1.50B[i,j] = WEL_1.50A[i,j]
    }
  else
    # The second possibility, however, is that maximum WELs have not yet been
    # achieved prior to the end of Year 10 (which appears to be the case,
    # given the assumptions we applied). In this event there are two sub-
    # possibilities. First, the individual was already relatively close to
    # maximum WELs, so that applying the new, lower rate of WEL growth would
    # result in a value lower than that observed for the previous year. In
    # this case, we apply the value observed for the previous year to this and
    # all future years.
    if (abs(DUR_T2DM_2^2) > abs(MGMT[2]*DUR_T2DM))
      for (i in 11:nrow(WEL_1.50B)){
        WEL_1.50B[i,j] = WEL_1.50A[10,j]
        }
      else
        # If, on the other hand, we had not and were not close to reaching the
        # maximum value, we apply essentially the same process as we did for
        # the first of these two tables.
        for (i in 11:nrow(WEL_1.50B)){
          {WEL_1.50B[i,j] = 
            min(max(WEL_1.50A[1:10,j]) + DUR_T2DM_2*(i-10)^2 + 
                      MGMT[2]*DUR_T2DM*(i-10))}
        }
}
Scenarios[[5]] <- WEL_1.50B

# The last thing to do is to ensure WELs not begin to decline after 
# plateauing. As above, this is done by setting all observations after the
# one associated with the maximum value, to that value. Note that we prepared
# the tables corresponding to the behavioural change scenarios by applying a
# two-step process, rather than trying to cram everything into a single step.
# Once we complete the manipulations, we pack the new tables back into the
# list in preparation for reporting.
for (h in 4:5) {
  df <- Scenarios[[h]]
  for (j in 1:ncol(df)){
    for (i in 1:nrow(df)){
      if (i>which.max(df[,j]))
      {df[i,j]=max(df[,j])}
      }
    }
  Scenarios[[h]] <- df
}

# Saving the Results ------------------------------------------------------

# Now we want to compile this information into a series of tables. In 
# particular, we want to generate five separate tables, four of which
# will ultimately be referenced by the decision model, and a fourth that 
# synthesizes the results from the first three, and that can be incorporated 
# directly into the body of the thesis.

# To prepare the tables for export into Excel, we withdraw each dataframe from
# the list, and assign it an appropriate name.
WEL_1.00_TBL <- Scenarios[[1]]
WEL_0.50A_TBL <- Scenarios[[2]]
WEL_1.50A_TBL <- Scenarios[[3]]
WEL_0.50B_TBL <- Scenarios[[4]]
WEL_1.50B_TBL <- Scenarios[[5]]

# Now we use the function "bind_cols" (which is included in the DPLYR package)
# to stack the three dataframes next to each other.
CT_1 <- bind_cols(WEL_1.00_TBL, WEL_0.50A_TBL, WEL_1.50A_TBL, 
                  WEL_0.50B_TBL, WEL_1.50B_TBL)
# Let's rearrange the columns of the table to allow the values across the 3
# scenarios to be more easily compared.

# To produce the tables we need, we must effectively place every 8th column
# next to one another. This created a situation where it was necessary to 
# diretly reference all 40 columns individually, which is not only inelegant,
# but also seriously prone to error. Fortunately, there is a (slightly) better
# way to do this, using an approach that was inspired from the following link:
# https://stackoverflow.com/questions/20531125/how-to-replace-values-in-every-nth-column-with-value-from-adjoining-column
# 
# Generate a vector consisting of every 8th integer between 1 and the number
# rows in CT_1.
Seq_1 <- seq(1,ncol(CT_1),8)
# Clone Seq_1 8 times.
Seq_8 <- Seq_7 <- Seq_6 <- Seq_5 <- Seq_4 <- Seq_3 <- Seq_2 <- Seq_1
# Generate a list containing all 8 vectors.
Sequences <- list(Seq_1, Seq_2, Seq_3, Seq_4, Seq_5, Seq_6, Seq_7, Seq_8)

# Recursively run through each element of the list just constructed, and
# increment by one each time (note that we use the same trick we applied to
# simplify the construction of the original tables). This generates 40 unique 
# integers, corresponding to the 40 columns in CT_1.
for (s in 2:length(Sequences)){
  vec <- Sequences[[1]] + (s-1)
  Sequences[[s]] <- vec
}

# Now we rearrange the columns of CT_1 by referring to the 8 sets of integers
# we produced above, each of which consists of 5 numbers. Note that I
# originally tried to reference "Sequences[[1:8]]" to avoid referring to
# elements of the list individually; this, unfortunately, appears to generate
# an error.
CT_2 <- CT_1[,c(Sequences[[1]],Sequences[[2]],Sequences[[3]],Sequences[[4]],
                Sequences[[5]],Sequences[[6]],Sequences[[7]],Sequences[[8]])]
# Let's select a subset of rows that would serve to illustrate the differences
# between the scenarios (i.e., the values for years 0-5, and for every 5th 
# year thereafter).
CT_3 <- CT_2[c(1:6,11,16,21,26,31,36,41,46,51,56,61,66,71,76),]

save(CT_2, file = paste(RDA, Filename, LM, ".Rda", sep = ""))

# Finally, we write the tables to XLSX, so that it can be formatted for insertion
# into the thesis. We do this using functions embedded into the OPENXLSX package.

# Create the workbook
wb = createWorkbook()

# Create the (blank) worksheets that will hold our data 
sht1 = addWorksheet(wb, "WEL_1.00_TBL"); sht2 = addWorksheet(wb, "WEL_0.50A_TBL")
sht3 = addWorksheet(wb, "WEL_1.50A_TBL"); sht4 = addWorksheet(wb, "WEL_0.50B_TBL")
sht5 = addWorksheet(wb, "WEL_1.50B_TBL"); sht6 = addWorksheet(wb, "CT_3")

# Add data to the worksheet we just created
writeData(wb, sht1, WEL_1.00_TBL); writeData(wb, sht2, WEL_0.50A_TBL); 
writeData(wb, sht3, WEL_1.50A_TBL); writeData(wb, sht4, WEL_0.50B_TBL); 
writeData(wb, sht5, WEL_1.50B_TBL); writeData(wb, sht6, CT_3)

# Save the workbook
saveWorkbook(wb, paste0(RPT_TBL, Filename, LM, ".xlsx"), overwrite = T)