# Control Module for the Patient Decision Model
# Created: 09/10/18; Last Modified: 06/24/19

# Objective/Rationale -----------------------------------------------------

# NOTE 062419: Attempted to accelerate execution of code by employing the
# FWRITE function (part of the DATA.TABLE package) to accelerate generation
# of reward and transition probability matrices. The original approach,
# which utilized the functionality of the OPENXLSX package, took more than
# 10 minutes to accomplish this - the FWRITE function does it in seconds.

# 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 script passes parameters to RWD&TP_Matrices_Construction_090618.R,
# before executing said script. This is done primarily to automate the
# generation of TP and RWD matrix components for individuals with differing
# personal characteristics. It also serves as a template for the script
# operationalizing the sensitivity analysis.

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

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

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

# NOTE 042019: SensAnal is a switch, the intent of which is to indicate whether it is
# necessary to calculate the TP matrices. The motivation for this is the observation
# that none of the sensitivity analyses actually vary the transition probabilities,
# so there is no real advantage in effectively carrying out the same calculations
# repeatedly - and significant added cost, both in terms of processing time and CPU
# memory.
SensAnal <- 0

Model_Path <- here()
# String for date last modified
LM <- "_042019"
# Filename for output
Filename <- "RWD&TP_Matrices_Construction"
# Paths for output
RDA <- paste0(Model_Path, "/R_Objects/")
MATLAB <- paste0(Model_Path, "/MATLAB/")
# Paths for components of the transition probability and reward matrices
path_TP <- paste0(MATLAB, "Baseline_Results/TP_Matrix_Components/Case_")
path_RWD <- paste0(MATLAB, "Baseline_Results/RWD_Matrix_Components/Case_")

# Load RDA files derived from other scripts, beginning with the UKPDS-OM2 simulation
# results
filenames <- c(paste0(RDA, "UKPDS_OM2_Results_Processing", LM, ".Rda"), 
               # Load life tables for probabilities relating to risk of mortality.
               paste0(RDA, "Life_Tables", LM, ".Rda"), 
               # Load vector of life expectancies (which are derived from the life tables).
               paste0(RDA, "Life_Exp", LM, ".Rda"),
               # Load health state transition probabilities derived from application of the DPoRT
               # risk equations.
               paste0(RDA, "DPoRT_Calculations", LM, ".Rda"),
               # Loading our calculations from the analysis of the 2011 SLID data.
               paste0(RDA, "SLID_2011_Calculations", LM, ".Rda"),
               # Load estimates of gross individual employment income derived from analysis of 2011
               # NHS data
               paste0(RDA, "NHS_2011_Calculations", LM, "_1", ".Rda"),
               # Load tax rates used in analysis of 2011 NHS data
               paste0(RDA, "NHS_2011_Calculations", LM, "_2", ".Rda"),
               # Load presenteeism calculations
               paste0(RDA, "Presenteeism_Calculations", LM, ".Rda"))

# https://stackoverflow.com/questions/3764292/loading-many-files-at-once
lapply(filenames, load, environment())

# * 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_K <- c(53480)
# We also generate a separate variable for 2017 MPE.
Max_Pen_Earn_2017_K <- c(55300)
# The sources for the remainder of these parameters are cited in the thesis. The
# OAS amounts are effective October 2017.
OAS_Amt_K <- 585.49; OAS_Threshold_K <- 72809.00
GIS_Amt_K <- 874.48; GIS_Threshold_K <- 3500.00
CPP_Contrib_Rate_K <- 0.0495; Max_CPP_Contrib_K <- 2564.10; Basic_Exemption_K <- 3500.00
EI_Contrib_Rate_K <- 0.0163; Max_EI_Contrib_K <- 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_K <- 0.03; Savings_Per_K <- 0.10; Year_Ret_K <- 65

# * Parameters relating to household expenditures ====

# Annual pharmaceutical expenditures, based on results of the pharmacist consultation.
# https://stackoverflow.com/questions/7070173/r-friendly-way-to-convert-r-data-frame-column-to-a-vector
RxExp.Early.Annual_K <- read_xlsx(path = paste0(Model_Path, "/Raw_Data/",
                                                "Pharmaceutical_Costs_V2.xlsx"), 
                                  sheet = "ND", range = "I7", col_names = FALSE)[[1]]
RxExp.Adv.WC.Annual_K <- read_xlsx(path = paste0(Model_Path, "/Raw_Data/",
                                                 "Pharmaceutical_Costs_V2.xlsx"), 
                                   sheet = "WM", range = "I12", col_names = FALSE)[[1]] 
RxExp.Adv.PC.Annual_K <- read_xlsx(path = paste0(Model_Path, "/Raw_Data/",
                                                 "Pharmaceutical_Costs_V2.xlsx"), 
                                   sheet = "PM", range = "I15", col_names = FALSE)[[1]] 
Min_Ded_K <- 100

# Now we load the Pharmacare rates for the 2018/19 benefit year, which were collected
# from the province's website.
url <- 'https://www.gov.mb.ca/health/pharmacare/estimator.html'
webpage <- read_html(url)

Pharm.00 <- html_nodes(webpage, css = "p+ div tr+ tr td")
Pharm.01 <- as.matrix(html_text(Pharm.00))

Pharm.01_Income <- Pharm.01[1]
Pharm.01_Income <- unlist(strsplit(Pharm.01_Income, split = "\r\n"))
# https://stackoverflow.com/questions/5992082/how-to-remove-all-whitespace-from-a-string
# https://www.rforexcelusers.com/remove-currency-dollar-sign-r/
# https://stackoverflow.com/questions/33949945/replace-multiple-strings-in-one-gsub-or-chartr-statement-in-r
# https://stackoverflow.com/questions/11129107/removing-non-alpanumeric-characters-from-an-ordered-collection-of-objects-list

# Remove all non-numeric characters, excepting dashes
Pharm.01_Income %<>%
  gsub(pattern = "[^0-9-]", replacement = "", .) %>%
  data.frame(., stringsAsFactors = FALSE) %>%
  rename("Income" = ".") %>%
  separate(., col = "Income", into = c("Lower.Limit", "Upper.Limit"), sep = "-")

# https://stackoverflow.com/questions/15629885/replace-na-in-column-with-value-in-adjacent-column
Pharm.01_Income[is.na(Pharm.01_Income)] <- 1000000

Pharm.01_DedPer <- Pharm.01[2]
Pharm.01_DedPer <- unlist(strsplit(Pharm.01_DedPer, split = "%\r\n"))
Pharm.01_DedPer %<>%
  gsub(pattern = "[^0-9.]", replacement = " ", .) %>%
  data.frame(., stringsAsFactors = FALSE) %>%
  rename("Rate" = ".")

Pharmacare_K <- cbind(Pharm.01_Income, Pharm.01_DedPer)
# https://stackoverflow.com/questions/18503177/r-apply-function-on-specific-dataframe-columns
Pharmacare_K[1:3] <- sapply(Pharmacare_K[1:3], as.numeric)
Pharmacare_K$Rate <- Pharmacare_K$Rate/100

# Dietary expenditures

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

Diet_Cost <- data.table(read_xlsx(path = paste0(Model_Path, "/Raw_Data/",
                                                "Meal_Plan_MS_v9.xlsm"), 
                                  sheet = "EER", range = "A21:H33", col_names = TRUE))
Diet_Cost <- Diet_Cost %>% rename(Diet_Type = "Diet Type")
# https://stackoverflow.com/questions/9202413/how-do-you-delete-a-column-by-name-in-data-table
# https://stackoverflow.com/questions/16638484/remove-multiple-columns-from-data-table
Diet_Cost.H <- Diet_Cost[grep("AHEI Score: 102.5", Diet_Type),]
Diet_Cost.UH <- Diet_Cost[grep("AHEI Score: 2.5", Diet_Type),]
Diet_Cost.H <- Diet_Cost.H[, c("Diet_Type") := NULL]
Diet_Cost.UH <- Diet_Cost.UH[, c("Diet_Type") := NULL]

Diet_Cost.H <- melt(Diet_Cost.H, id.vars = c("Age_Group"), variable.name = c("Chars"),
                    value.name = "Cost")
Diet_Cost.H <- separate(data = Diet_Cost.H, col = Chars, into = c("Ethnicity", "Gender")
                        , sep = "_")
Diet_Cost.UH <- melt(Diet_Cost.UH, id.vars = c("Age_Group"), variable.name = c("Chars"),
                     value.name = "Cost")
Diet_Cost.UH <- separate(data = Diet_Cost.UH, col = Chars, into = c("Ethnicity", "Gender")
                         , sep = "_")

Diet_Cost.H.01_K <- Diet_Cost.UH.01_K <- 
  data.frame(expand.grid("Age" = seq(from = 25, to = 95), "Gender" = c("Male", "Female"),
                         Ethnicity = c("Caucasian", "Afro-Caribbean", "Asian-Indian"), 
                         stringsAsFactors = F))
Diet_Cost.H.01_K <- Diet_Cost.H.01_K %>%
  mutate(Age_Group = Age_Lookup_K[findInterval(Age,Age_Lookup_K$Age),2])
Diet_Cost.H.01_K <- merge(Diet_Cost.H.01_K, Diet_Cost.H, 
                          by = c("Age_Group", "Gender", "Ethnicity"), all.x = T)
Diet_Cost.UH.01_K <- Diet_Cost.UH.01_K %>%
  mutate(Age_Group = Age_Lookup_K[findInterval(Age,Age_Lookup_K$Age),2])
Diet_Cost.UH.01_K <- merge(Diet_Cost.UH.01_K, Diet_Cost.UH, 
                           by = c("Age_Group", "Gender", "Ethnicity"), all.x = T)

# LTPA expenditures
LTPA_Cost_List_K <- c(0,0)

# * Parameters relating to time allocation ====

# Time allocated to sleep each year
Time_Req_Sleep_K <- 8*365
# Time available for any other activity
Time_Available_K <- 24*365 - Time_Req_Sleep_K
# Time allocated to LTPA (for those adherent to it)
Time_Req_LTPA_K <- (150/60)*52
# Time allocated to food preparation
Time_Req_Diet_H_K <- 90*365/60; Time_Req_Diet_UH_K <- 30*365/60
# Time allocated to job search (for individuals younger than 65 who are
# currently unemployed)
Time_Req_Search_K <- 9.5*52

# * Other miscellaneous parameters ====

# 052519: Penalty for non-admissible actions
Penalty_K <- -1000000
# 042019: Price of consumption, the multiplier applied to the wage, and Rho, which is 
# a key parameter of the CES utility function, and which relates to the elasticity of
# substitution in the following way: sigma = (1/(1-Rho)).
P_c_K <- 1; Mult_w_K <- 1; Rho_K <- 0;
# Impact of healthy, prediabetic, and T2DM states on probability of sustaining and/or
# recovering employment
Emp_Diff_Healthy_K = 0; Emp_Diff_Pred_K = 0; Emp_Diff_T2DM_K = 0
# Starting age for the decision problems
Start_Age_K = 24
# Monthly EIA for individuals with no dependents in the General Assistance category: 
# https://www.gov.mb.ca/fs/eia/eia_general.html
EIA_General_K = 12*806
# Monthly EIA for individuals with no dependents who is considered to have a disability:
# https://www.gov.mb.ca/fs/eia/EIA_Disability_K_K.html
EIA_Disability_K = 12*1027
Annual_Work_Hrs.FT_K <- 2000; Annual_Work_Hrs.PT_K <- 1000
Pres_Share_K <- 0.5
Paid_SDs.FT_K <- 7; Paid_SDs.PT_K <- Paid_SDs.FT_K/2
# "ARPD" = "age-related productivity decline"
ARPD_K <- 0.025
# Parameters for the utility function
Alpha_K <- 0.5

# * Personal characteristics ====
GEN_K = list("Male", "Female")
ETH_K = list("Caucasian", "Afro-Caribbean", "Asian-Indian")

EDU_K <- 
  data.frame(Index = as.character(seq(1:5)), 
             Levels = 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",
                        "5. Postsecondary certificate, diploma or degree"), 
             Pres_Code = c("Educ_1", "Educ_2", "Educ_3", "Educ_4", "Educ_5"),
             # 081718: The next line effectively defines PSE for the purposes of our
             # calculations. The paper in which the DPoRT 2.0 Risk Equations were published
             # included a variable that denoted any PS attendance. Reviewing the data dictionary
             # and derived variable guide for the CCHS 2012 confirms the authors' interpretation
             # was correct.
             # 
             # However, the NHS data indicates only the individual's level of educational
             # attainment (i.e., whether they actually progressed far enough in their education 
             # to receive some kind of credential). We bridge the gap by assuming a person must
             # have obtained a PS certificate or diploma below bachelor level to satisfy the 
             # definition. This is conservative in that it potentially excludes many individuals
             # who attended PS educational institutions but did not progress far enough to 
             # receive a credential. However, it avoids the risk of including a potentially 
             # large segment of individuals who did not attend PSE after completing high school 
             # (NOTE: the next line of code been revised - an earlier version of the code 
             # included everyone with at least a HS education in the definition of "some PSE").
             Attended_PS = c(rep("NoPS", 2), rep("PS", 3)))

Age_Align_K <- 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 = FALSE)

rm(Diet_Cost, Diet_Cost.H, Diet_Cost.UH, Pharm.00, Pharm.01, Pharm.01_DedPer, Pharm.01_Income,
   webpage, url)

# Generation of Reward and Transition Probability Matrices ----------------

for (delta in seq(4)){
  for (gamma in seq(3)){
    for (theta in seq(2)){

      # These integer values index positions in the three lists constructed above. BY
      # MODIFYING THESE THREE VALUES, IT IS POSSIBLE TO GENERATE TRANSITION PROBABILITY AND
      # REWARD MATRICES FOR ALL COMBINATIONS OF FIXED PERSONAL CHARACTERISTICS.
      Gender.U = theta; Ethnicity.U = gamma; Education.U = delta
      # This next variable is included specifically to faciliate integration of 2011 SLID
      # data. Anyone who has at minimum a high school diploma or equivalent is assigned a
      # value of 1 (note that this does not align precisely with the 2011 SLID codebook,
      # since many individuals completing secondary school will not go on to participate in
      # ANY post-secondary education).
      Education_SLID = as.character(ifelse(Education.U > 1, 1, 0))
      
      IsTopQuin <- Emp_Income[Emp_Income$Sex == GEN_K[[Gender.U]] & Emp_Income$Highest_Degree ==
                                EDU_K[Education.U, 2],c(1:3, 5)]
      IsTopQuin <- merge(IsTopQuin, Age_Align_K, by = "Age_Cat", all.x = TRUE)
      IsTopQuin$IsTopQuin <- ifelse(IsTopQuin$IsTopQuin==FALSE, "NoTQ", "TQ")

      # NOTE 092318: Now we call the appropriate script to initiate the generation of the reward
      # and transition probability matrices (note that at present, this is slighly different
      # than the script used for the baseline analysis - for reasons we discuss below).      
      source(paste0(Model_Path, "/RWD&TP_Matrices_Construction", LM, ".R"))
      
      # NOTE 092318: For each model run, we save the reward and transition probability matrix
      # components in the appropriate place. This code was actually migrated here from the
      # RWD&TP_Matrices_Construction script.
      # https://stackoverflow.com/questions/48063019/formatting-percentages-in-r-package-openxlsx
      # https://stackoverflow.com/questions/50413253/openxlsx-currency-style-is-there-a-rounding-option?rq=1
      # https://stackoverflow.com/questions/20648149/what-are-numberformat-options-in-excel-vba
      # 
      # 091018: I tried several approaches before finding one that worked. Originally,
      # each component of the TP matrix was saved to CSV; for reasons I cannot now
      # recall, I determined that MATLAB doesn't work well with CSV files, and had
      # developed an Excel macro to convert these to XLSX files. I subsequently learned 
      # (as the note below suggests), that there are faster ways to save to CSV.
      # However, this added an unnecessary extra step. Why not just save to XLSX
      # directly from R? I found that this could be done relatively easily using the
      # OPENXLSX package, but encountered some new issues. In particular, it became
      # clear that MATLAB categorized matrices imported from these XLSX files as cell
      # arrays rather than as double-precision. Initially, the only way I could
      # successfully implement the conversion was by manually setting data types in
      # Excel, which is not feasible for large numbers of workbooks. The issue
      # persisted even after I programmed OPENXLSX to format the content of each
      # component of the TP matrix as numeric. I ultimately addressed the issue by
      # applying the AS.NUMERIC function to each element of TP.03. Additional details
      # about the issue and its resolution are presented in the MATLAB M-File.
      # 
      # 090818: This has been revised. The earlier version of the code below used the
      # WRITE.CSV function rather than the FWRITE function, which I believe to be a
      # less efficient way of writing to CSV. In addition, the file name now includes
      # a series of integers indicating the gender, ethnicity, and education of the
      # individual for whom the decision problem is being set up.
      
      for (i in seq(TP.03))
      {fwrite(x = TP.03[[i]], sep = ",", col.names = FALSE, file = 
                paste0(path_TP, Gender.U, Ethnicity.U, Education.U, "/TP_", i, "_",
                       Gender.U, Ethnicity.U, Education.U, ".csv"))}
      
      for (i in seq(RWD.03))
      {fwrite(x = RWD.03[[i]], sep = ",", col.names = FALSE, file = 
                paste0(path_RWD, Gender.U, Ethnicity.U, Education.U, "/RWD_", i, "_",
                       Gender.U, Ethnicity.U, Education.U, ".csv"))}
      
      # NOTE 112318: Save efficiency loss calculations for individuals who have developed
      # T2DM. Whereas the individual's accounts for losses outside the workplace in their
      # own decision-making, these figures include only losses incurred on the job (since
      # calculated values are otherwise significantly higher than estimates from other
      # sources, such as the Canadian Diabetes Association).
      write.table(Ref_States.02,
                  file = paste0(MATLAB, "Simulation_Results/Efficiency_Loss_Calculations/",
                                "Eff_Losses_Alive_", Gender.U, Ethnicity.U, Education.U, 
                                ".dat"), quote = FALSE,
                  sep = ",", row.names = FALSE)
      
      # NOTE 112318: Save efficiency loss calculations for individuals who are deceased.
      # These losses include everything the individual WOULD HAVE earned if they had
      # survived to a given point in time and remained working.
      write.table(Emp_Income.02a,
                  file = paste0(MATLAB, "Simulation_Results/Efficiency_Loss_Calculations/",
                                "Eff_Losses_Dead_", Gender.U, Ethnicity.U, Education.U, 
                                ".dat"), quote = FALSE,
                  sep = ",", row.names = FALSE)
      
      # Clear variables used in preparation of TP matrix for next iteration. Note that in the
      # original code, the construct Factors was deleted twice, since it was redefined for
      # the purposes of RWD matrix component generation after having previously been used to
      # generate TP matrix components. When the code was migrated here from the script in
      # which it previously appeared (i.e., RWD&TP_Matrices_Construction), it initially
      # returned an error, because the same construct was effectively being deleted twice
      # (i.e., in the original code, it was defined, deleted, and then defined again - here
      # it is defined, redefined [overwritten], and then deleted). 
      
      if (SensAnal == 0) {
        
        rm(Decision_Combos, Decision_List, Emp_Diff, Emp_List.01, Emp_List.02, Emp_Stat_1, Emp_Trans.01,
           Emp.02, Factors, Healthy_Healthy_List, Healthy_Healthy.02, Healthy_T2DM_List, Healthy_T2DM.01,
           Healthy_T2DM.02, Healthy_Trans_A.01, Healthy_Trans_A.02, Healthy_Trans_B.01, Healthy_Trans_B.02,
           Healthy_Trans_C.01, Healthy_Trans_C.02, Healthy_Trans.01, Mortality.01, Mortality.02,
           Mortality.03, Mortality.03a, Mortality.03b, Mortality.03c, Output6_TEST, Possible_Decisions,
           Pr_Risk.01, Pr_Risk.02, Ref_Actions, SLID_2011_S6.01, State_List_TP, State_List_TP.02a,
           State_List_TP.02b, T2DM_Trans_0_9.01, T2DM_Trans_0_9.02, T2DM_Trans_10_94.01, T2DM_Trans_10_94.02,
           T2DM_Trans_9_10.01, T2DM_Trans_9_10.02, T2DM_Trans_Dead.01, T2DM_Trans_Dead.02, TP, TP_Dead.02,
           TP_Healthy_Dead, TP_Healthy_Live, TP_Healthy.01, TP_Healthy.02, TP_T2DM, TP_T2DM_0_9.01,
           TP_T2DM_0_9.02, TP_T2DM_10_94.01, TP_T2DM_10_94.02, TP_T2DM_9_10.01, TP_T2DM_9_10.02,
           TP_T2DM_Dead_0_1, TP_T2DM_Dead_1_94, TP_T2DM_Dead.02, TP_T2DM_Live_0_1, TP_T2DM_Live_1_94,
           TP_T2DM.02, TP.03, Transition_List.01, Transition_List.02, Yr_0_Choices.01, Yr_0_Choices.02, 
           Yr_9_Choices.02)
        
        rm(Activity, AD_Diet, AD_LTPA, AD_Medication, Diet, Emp_List_Names, LMP, TP_Death)
        
      } else {
        
      }
      
      # Clear variables used in preparation of RWD matrix for next iteration.
      # NOTE 042019: I now also clear the three new variables I added (i.e., the price of a unit of
      # consumption, the wage multiplier, and the parameter Rho)
      rm(CT_2.02, CT_2.03, CT_2.04, Diet_Annual_Cost.H, Diet_Annual_Cost.UH, Emp_Income.01, 
         Emp_Income.02, Emp_Income.02a, Emp_Income.02b, Presenteeism.01, Presenteeism.02, Ref_States, 
         Ret_Income, RWD, RWD_Healthy.01, RWD_T2DM, RWD.02, RWD.03, State_List_RWD, State_List_RWD.02)
      
      rm(AD_DIET_P1, AD_DIET_P2, AD_LTPA_P1, AD_LTPA_P2, AD_MEDS_P1, AD_MEDS_P2, CPP.00, Duration,
         Headers, Levels, LTPA_Annual_Cost, RWD_Death, Withdrawals)
      
      # Clear variables associated with preparation of dataframe illustrating the calculation
      # of the value function.
      # rm(TP.INT, INT, RWD.INT)
      
    } 
  }
}

# Scrapyard ---------------------------------------------------------------

# NOTE 062419: This code reflects an earlier approach to generating the
# components of the reward and transition probability matrix, using the
# functionality of the OPENXLSX package. I subsequently replaced this with
# simpler code that employs the FWRITE function from the DATA.TABLE package,
# as this reduces the time required for execution of the code substantially
# (> 10 minutes).
# for (i in seq(TP.03))
# {# Create an Excel workbook object and add a worksheet
#   wb = createWorkbook()
#   sht = addWorksheet(wb, "Sheet_1")
#   dat = TP.03[[i]]
#   # Create a percent style
#   numFmt = createStyle(numFmt="#,##0.00")
#   # Add data to the worksheet we just created
#   writeData(wb, sht, dat, colNames = F)
#   # Add the percent style to the desired cells
#   addStyle(wb, sht, style=numFmt, cols=c(1:3), rows=1:(nrow(dat)), gridExpand=TRUE)
#   saveWorkbook(wb, paste0(path_TP, Gender.U, Ethnicity.U, Education.U, "/TP_", i, "_",
#                           Gender.U, Ethnicity.U, Education.U, ".xlsx"), overwrite = T)}
# 
# for (i in seq(RWD.03))
# {# Create an Excel workbook object and add a worksheet
#   wb = createWorkbook()
#   sht = addWorksheet(wb, "Sheet_1")
#   dat = RWD.03[[i]]
#   # Create a percent style
#   numFmt = createStyle(numFmt="#,##0.00")
#   # Add data to the worksheet we just created
#   writeData(wb, sht, dat, colNames = F)
#   # Add the percent style to the desired cells
#   addStyle(wb, sht, style=numFmt, cols=c(1:2), rows=1:(nrow(dat)), gridExpand=TRUE)
#   saveWorkbook(wb, paste0(path_RWD, Gender.U, Ethnicity.U, Education.U, "/RWD_", i, "_",
#                           Gender.U, Ethnicity.U, Education.U, ".xlsx"), overwrite = T)}