# Control Module for the Patient Decision Model - Sensitivity Analysis
# Created: 09/22/18; Last Modified: 06/24/19

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

# NOTE 042019: Extended the sensitivity analysis to allow for examination of the impact
# of changes to the price of consumption and leisure (i.e., wages), as well as to the
# parameters of the CES utility function. In addition, by setting the default (RSTUDIO) 
# working directory to the DIABETES_MODEL subfolder, I have simplified the process of 
# detecting the location of the active script. 
# 
# 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 <- 1

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/")

# 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 sensitivity analysis ====

# NOTE 092018: This next line is VERY important, because it determines the characteristics of
# the decision-maker for the purposes of the decision analysis, AS WELL AS the the parameter
# being varied.
Person_Type <- c(2, 3, 4); Sens_Analysis_ID <- 5
# 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 = Person_Type[1]; Ethnicity.U = Person_Type[2]; Education.U = Person_Type[3]

# NOTE 092018: This line indicates the parameter we are currently considering, and will need
# to be modified every time we examine a new one by changing the value of Sens_Analysis_ID
# below (note that said variable simply indexes a location in this vector).
Sens_Analysis_Types <- 
  c("1_Functional_Form", "2_Alpha", "3_Discount_Factor", "4_Price_Consump", "5_Wage", 
    "6_ARPD", "7_Savings_Percentage", "8_Dietary_Costs")
Sens_Analysis_Vars <- 
  c("Rho_K", "Alpha_K", "Disc_Factor_K", "P_c_K", "Mult_w_K", "ARPD_K", "Savings_Per_K", 
    "Diet_Cost_K")
# NOTE 092318: Sens_Analysis_Table indicates the levels of all of the parameters 
# included in the sensitivity analysis. The levels for Diet_Cost_K will need to be
# revisited when we determine how this parameter will be worked into the model.
Sens_Analysis_Table <- 
  data.table(SA_Type = Sens_Analysis_Types,
             SA_Vars = Sens_Analysis_Vars,
             Level_1 = c(-10.00, 0.00, 0.00, 0.50, 0.50, 0.00, 0.00, 0.00),
             Level_2 = c(-1.00, 0.25, 0.50, 0.75, 0.75, 0.01, 0.05, 0.50), 
             Level_3 = c(+0.00, 0.50, 0.75, 1.00, 1.00, 0.025, 0.10, 1.00),
             Level_4 = c(+0.50, 0.75, 0.95, 1.50, 1.50, 0.05, 0.15, 1.50), 
             Level_5 = c(+1.00, 1.00, 1.00, 2.00, 2.00, 0.10, 0.20, 2.00))

# "ARPD" = "age-related productivity decline";
Rho_K <- as.numeric(Sens_Analysis_Table[SA_Vars == "Rho_K", "Level_3"])
Alpha_K <- as.numeric(Sens_Analysis_Table[SA_Vars == "Alpha_K", "Level_3"])
Disc_Factor_K <- as.numeric(Sens_Analysis_Table[SA_Vars == "Disc_Factor_K", "Level_4"])
P_c_K <- as.numeric(Sens_Analysis_Table[SA_Vars == "P_c_K", "Level_3"])
Mult_w_K <- as.numeric(Sens_Analysis_Table[SA_Vars == "Mult_w_K", "Level_3"]) 
ARPD_K <- as.numeric(Sens_Analysis_Table[SA_Vars == "ARPD_K", "Level_3"])
Savings_Per_K <- as.numeric(Sens_Analysis_Table[SA_Vars == "Savings_Per_K", "Level_3"])
Diet_Cost_K <- as.numeric(Sens_Analysis_Table[SA_Vars == "Diet_Cost_K", "Level_3"])

Baseline_Values <- 
  data.frame(Baseline = c(Rho_K, Alpha_K, Disc_Factor_K, P_c_K,
                          Mult_w_K, ARPD_K, Savings_Per_K, Diet_Cost_K))

Sens_Analysis_Table <- data.frame(cbind(Sens_Analysis_Table, Baseline_Values))

# * 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; 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: http://www.gov.mb.ca/health/pharmacare/estimator.html

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)
# This is unique to the sensitivity analysis. We need to store the values in the dataframe
# containing our original estimates of the cost of sustaining an unhealthy diet, because
# the values in Diet_Cost.UH.01_K can change in every loop.
Diet_Cost.UH.01_BASE <- Diet_Cost.UH.01_K

# 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
# 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

# * 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 ----------------

# NOTE 042219: We no longer require the path to the transition probability matrices, as
# we are now drawing upon the matrices generated through the baseline analysis rather
# than needlessly recreating them on every loop (i.e., because no part of the sensitivity
# analysis actually modifies any element of the matrices).
# 
# Paths for components of the reward matrices (i.e., for the output generated by this script).
path_RWD <- paste0(MATLAB, "Sensitivity_Analysis/", Sens_Analysis_Types[Sens_Analysis_ID],
                   "/RWD_Matrix_Components/Case_")

for (delta in 1:5){
  
  # NOTE 042419: Custom function for determination of parameter values.
  # https://stackoverflow.com/questions/4108577/cast-function-argument-as-a-character-string
  Parameters <- function(x) {
    ifelse(Sens_Analysis_ID==which(Sens_Analysis_Table$SA_Vars==deparse(substitute(x))), 
           Sens_Analysis_Table[Sens_Analysis_ID, delta + 2], 
           Sens_Analysis_Table[which(Sens_Analysis_Table$SA_Vars==deparse(substitute(x))), 
                               "Baseline"])}
  
  # NOTE 042219: I identified a much cleaner way of defining the values of these
  # parameters during each run, by defining a custom function (see above).
  Rho_K <- Parameters(Rho_K); Alpha_K <- Parameters(Alpha_K)
  Disc_Factor_K <- Parameters(Disc_Factor_K); P_c_K <- Parameters(P_c_K)
  Mult_w_K <- Parameters(Mult_w_K); ARPD_K <- Parameters(ARPD_K)
  Savings_Per_K <- Parameters(Savings_Per_K); Diet_Cost_K <- Parameters(Diet_Cost_K)
  
  # NOTE 113018: The value of the unhealthy diet is uniformly multiplied by the value
  # of Diet_Cost_K.
  Diet_Cost.UH.01_K <- Diet_Cost.UH.01_BASE 
  Diet_Cost.UH.01_K$Cost <- Diet_Cost.UH.01_K$Cost*Diet_Cost_K
  
  # 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")
  
  if (Sens_Analysis_ID %in% 4:5 == FALSE){
    
    # 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.
    for (i in seq(RWD.03))
    {fwrite(x = RWD.03[[i]], sep = ",", col.names = FALSE, file = 
              paste0(path_RWD, delta, "/RWD_", i, "_",
                     Gender.U, Ethnicity.U, Education.U, ".csv"))}
    
    # NOTE 042119: Eliminate code clearing variables used in preparation of TP matrix 
    # for next iteration (this code still appears in the control module for the baseline analysis).
    
    # Clear variables used in preparation of RWD matrix for next iteration.
    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)
    
  } else {
    
    for (phi in seq(from = 1, to = 5, by = 2)){
      
      Rho_K <- Sens_Analysis_Table[1, phi + 2]
      
      # 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.
      for (i in seq(RWD.03))
      {fwrite(x = RWD.03[[i]], sep = ",", col.names = FALSE, file = 
                paste0(path_RWD, delta, "/Rho_", Sens_Analysis_Table[1, phi + 2], "/RWD_", i, "_",
                       Gender.U, Ethnicity.U, Education.U, ".csv"))}
      
      # NOTE 042119: Eliminated code clearing variables used in preparation of TP matrix 
      # for next iteration (this code still appears in the control module for the baseline analysis).
      
      # Clear variables used in preparation of RWD matrix for next iteration.
      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)
      
    }
  }
}