R heemod package issue with reprex

[mmaughan0]

Hello all

Not sure if this was the right spot to post this question. I have quite a bit of code here, but my issue is only with the very last 5 lines of it.

Background: I am using R heemod package to run a cost effectiveness model

Issue: THe Probabilistic Sensitvity Analysis (Monte Carlo) at the end (under the "def_psa" command) will only vary the point estimate for effect related variables, and not cost related variables. In the reprex below I include one of each, but I tried many of each and it was a consistent issue - only the effect (x axis) would vary on the coste effectiveness plane. The cost (y-axis) stayed at the original point estimate before I ran the PSA

More context: the survival tables are dummy tables to get the reprex to run. The output values therefore are not meaningful. Only the lack of variance in the cost ouptus for the PSA/Monte Carlo is. Also my real version that is not a reprex has variance in effect. This doesnt because I didnt vary the survival tables for the sake of the reprex. The Resupply_Cost ETU should still show variance in the reprex as shown and it does not. I have also tried Tx_Cost, as it only affects the Mab114 and REGNeb3 arms and not the base case arms. Still no variance

[mmaughan0]

library(heemod)
library(flexsurv)
#> Loading required package: survival
library(rgho)
library(dplyr)
#> 
#> Attaching package: 'dplyr'
#> The following objects are masked from 'package:stats':
#> 
#>     filter, lag
#> The following objects are masked from 'package:base':
#> 
#>     intersect, setdiff, setequal, union
library(readxl)
library(reprex)

##setting up directory, as we have many tables to put in that are too large to practically code##

setwd("C:/Users/Matt/OneDrive/Documents/Portfolio Project")

##reading in necessary tables from directory##


objects
#> function (name, pos = -1L, envir = as.environment(pos), all.names = FALSE, 
#>     pattern, sorted = TRUE) 
#> {
#>     if (!missing(name)) {
#>         pos <- tryCatch(name, error = function(e) e)
#>         if (inherits(pos, "error")) {
#>             name <- substitute(name)
#>             if (!is.character(name)) 
#>                 name <- deparse(name)
#>             warning(gettextf("%s converted to character string", 
#>                 sQuote(name)), domain = NA)
#>             pos <- name
#>         }
#>     }
#>     all.names <- .Internal(ls(envir, all.names, sorted))
#>     if (!missing(pattern)) {
#>         if ((ll <- length(grep("[", pattern, fixed = TRUE))) && 
#>             ll != length(grep("]", pattern, fixed = TRUE))) {
#>             if (pattern == "[") {
#>                 pattern <- "\\["
#>                 warning("replaced regular expression pattern '[' by  '\\\\['")
#>             }
#>             else if (length(grep("[^\\\\]\\[<-", pattern))) {
#>                 pattern <- sub("\\[<-", "\\\\\\[<-", pattern)
#>                 warning("replaced '[<-' by '\\\\[<-' in regular expression pattern")
#>             }
#>         }
#>         grep(pattern, all.names, value = TRUE)
#>     }
#>     else all.names
#> }
#> <bytecode: 0x000001e3ede92290>
#> <environment: namespace:base>

##par_mod <- define_parameters(
##age_base = 30*365,##cycle in this model is 1 day
##age_cycle = model_time + age_base)

## I muted this code above because I don't think its necessary to import WHO life tables - as my cycle length is 1 day and the lifetables
#are based on years I do not want to muddle the two. I can realize the lost life years during transition to the death state

reprex
#> function (x = NULL, input = NULL, wd = NULL, venue = c("gh", 
#>     "r", "rtf", "html", "slack", "so", "ds"), render = TRUE, 
#>     advertise = NULL, session_info = opt(FALSE), style = opt(FALSE), 
#>     comment = opt("#>"), tidyverse_quiet = opt(TRUE), std_out_err = opt(FALSE), 
#>     html_preview = opt(TRUE), outfile = deprecated(), show = deprecated(), 
#>     si = deprecated()) 
#> {
#>     if (lifecycle::is_present(show)) {
#>         html_preview <- show
#>         lifecycle::deprecate_warn(when = "1.0.0", what = "reprex(show)", 
#>             with = "reprex(html_preview)")
#>     }
#>     if (lifecycle::is_present(si)) {
#>         session_info <- si
#>     }
#>     reprex_impl(x_expr = substitute(x), input = input, wd = wd, 
#>         venue = venue, render = render, new_session = TRUE, advertise = advertise, 
#>         session_info = session_info, style = style, html_preview = html_preview, 
#>         comment = comment, tidyverse_quiet = tidyverse_quiet, 
#>         std_out_err = std_out_err, outfile = outfile)
#> }
#> <bytecode: 0x000001e3fa02b410>
#> <environment: namespace:reprex>
  
par_mod <- define_parameters(
  
  p_death_disease_base = compute_surv(
    fit_death_disease_base, ##fit death disease base will be defined later, and will pull in TTE tables
    time = state_time, ##state time means time in that given state, in this case it will be in the "conf-sick" state
    km_limit = 28)) ## km_limit is 28 because the tables I have go to 28 days

##The above code is calculating probability based on a life curve that I will define in a later piece of code##


par_mod <- modify(
  par_mod,

  p_death_disease_mAb114 = compute_surv(
    fit_death_disease_mAb114,
    time = state_time,
    km_limit = 28))
  
##the above code does evertything the p_death_disease_base does, but for the REGNEB3 treatment


par_mod <- modify(
  par_mod,
  
  p_death_disease_REGN_EB3 = compute_surv(
    fit_death_disease_REGNEB3,
    time = state_time,
    km_limit = 28))

## the above code does everything p_death_disease_base does, but for the mAb114 treatment
##pulling in induct to treat table##



par_mod <- modify(
  par_mod,
  
  p_tx_screen = compute_surv(
    fit_screen_to_treatment,
    time = state_time,
    km_limit = 28))

## the above code is setting a probability of moving to the next stage for the time between induction(screen) and treatment. 
##fit_screen_treatment will be defined later

###pulling in discharge tables###

par_mod <- modify(
  par_mod,
  
  p_discharge_base = compute_surv(
    fit_discharge_base,
    time = state_time,
    km_limit = 28))

par_mod <- modify(
  par_mod,
  
  p_discharge_mAb114 = compute_surv(
    fit_discharge_mAb114,
    time = state_time,
    km_limit = 28))

par_mod <- modify(
  par_mod,
  
  p_discharge_REGNEB3 = compute_surv(
    fit_discharge_REGNEB3,
    time = state_time,
    km_limit = 28))

## the above probabilities function just like the death probabilities, but for discharge

##dummy table to allow the reprex to work. actual model has different tables for each##

tab_surv <- structure(list(time = c(0.4, 8.7, 7, 5.1, 9.2, 1, 0.5, 3.3, 1.8,
                                    3, 6.7, 3.7, 1.1, 5.9, 5.1, 10, 10, 10, 10, 10, 10, 10, 10, 10,10), 
                           status = c(1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
              1L, 1L, 1L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L)), .Names = c("time", "status"), row.names = c(NA, -25L), class = "data.frame")

##The below code is what fits the curves into probabilities

fit_death_disease_base <- flexsurv::flexsurvreg(
  survival::Surv(time, status) ~ 1,
  dist = "gompertz",
  data = tab_surv)


fit_death_disease_mAb114 <- flexsurv::flexsurvreg(
  survival::Surv(time, status) ~ 1,
  dist = "gompertz",
  data = tab_surv)

fit_death_disease_REGNEB3 <- flexsurv::flexsurvreg(
  survival::Surv(time, status) ~ 1,
  dist = "gompertz",
  data = tab_surv)


fit_screen_to_treatment <- flexsurv::flexsurvreg(
  survival::Surv(time, status) ~ 1,
  dist = "gompertz",
  data = tab_surv)

fit_discharge_base <- flexsurv::flexsurvreg(
  survival::Surv(time, status) ~ 1,
  dist = "gompertz",
  data = tab_surv)

fit_discharge_mAb114 <- flexsurv::flexsurvreg(
  survival::Surv(time, status) ~ 1,
  dist = "gompertz",
  data = tab_surv)

fit_discharge_REGNEB3 <- flexsurv::flexsurvreg(
  survival::Surv(time, status) ~ 1,
  dist = "gompertz",
  data = tab_surv)


###define matrices####

##putting in point estimates from Levine et. al. paper. All diagnostic rates can be derived from these 3 data points##
# prevalence may correlate with sens/spec in an epidemiological sense, but not in a mathematical sense. For our purposes they are considered independent

#Definitions:
#TP = True Positive
#TN = True Negative
#FP = False Positive
#FN = False Negative

#Sensitivity = TP/(TP+FN) --> Think "If I am positive, what are the chances I test positive"
#Specificity = TN/(TN+FP) --> Think "If I am negative, what are the chances I test negative"
# Prevalence = (TP+FN)/(TP+TN+FP+FN) --> What are the chances a person in the given population has the disease



##drew the below from Levine et. al. 2015 - which helpfully has confidence intervals, and thus can be coded probabilistically 
sensitivity_algo = .93
specificity_algo = .23
prevalence = .42

## the below just rearranges all the terms and definitions above algebraically. no new information is here, just rearranging

positive_algo = sensitivity_algo * prevalence + (1- specificity_algo) * (1- prevalence)

p_pos_low = ((1-sensitivity_algo)*prevalence)/(((1- sensitivity_algo)* prevalence) + (specificity_algo * (1- prevalence)))
p_pos_hi = (sensitivity_algo * prevalence)/((sensitivity_algo * prevalence)+((1- prevalence)*(1- specificity_algo)))


##define base transition matrix##



par_mod <- modify(
  par_mod,
  
  n_days = 11,
  death_transition_prob_base = ifelse(
    state_time < n_days, 
    p_death_disease_base, 
    0)
)

par_mod <- modify(
  par_mod,
  
  n_days = 11,
  death_transition_prob_mAb114 = ifelse(
    state_time < n_days, 
    p_death_disease_mAb114, 
    0)
)

par_mod <- modify(
  par_mod,
  
  n_days = 11,
  death_transition_prob_REGNEB3 = ifelse(
    state_time < n_days, 
    p_death_disease_REGN_EB3, 
    0)
)

par_mod <- modify(
  par_mod,
  
  n_days = 11,
  screen_transition_prob = ifelse(
    state_time < n_days - 5, 
    p_tx_screen, 
    1)
)


mat_base <- define_transition(
  state_names = c("induction", "triage_hi", "triage_lo", "conf_sick", "death", "terminal"),
  
  0, positive_algo, C, 0, 0, 0, ##you can only go from induction (where you start) to two states - hi or low triage. the positive algo gives a point estimate
  0, C, 0, screen_transition_prob*p_pos_hi, 0, screen_transition_prob*(1-p_pos_hi),  ##p_tx_screen*p_pos_hionce in triage_hi, the p_tx_screen gives the chances of moving to next stage by day, but multiplies this by chance of being positive. 
  0, 0, C, screen_transition_prob*p_pos_low, 0, screen_transition_prob*(1-p_pos_low), ##p_tx_screen*p_pos_low same as above but starting from low 
  0, 0, 0, C, death_transition_prob_base, p_discharge_base, ##there are two curves here from the confirmed sick state, you can either die or be discharged. the "plug" is remaining in that state. I  am honestly not sure if this will work.
  0, 0, 0, 0, 0, 1, ##death covers death through burial. it is not a terminal state - leaving it as a terminal state would count the death costs every remaining cycle, vastly overcounting QALY loss and burial cost
  0, 0, 0, 0, 0, 1) #terminal is the only terminal state.

plot(fit_death_disease_mAb114)

mat_mAb_114 <- define_transition(
  state_names = c("induction", "triage_hi", "triage_lo", "conf_sick", "death", "terminal"),
  
  0, positive_algo, C, 0, 0, 0, ##you can only go from induction (where you start) to two states - hi or low triage. the positive algo gives a point estimate
  0, C, 0, screen_transition_prob*p_pos_hi, 0, screen_transition_prob*(1-p_pos_hi),  ##p_tx_screen*p_pos_hionce in triage_hi, the p_tx_screen gives the chances of moving to next stage by day, but multiplies this by chance of being positive. 
  0, 0, C, screen_transition_prob*p_pos_low, 0, screen_transition_prob*(1-p_pos_low), ##p_tx_screen*p_pos_low same as above but starting from low 
  0, 0, 0, C, death_transition_prob_mAb114, p_discharge_mAb114, ##there are two curves here from the confirmed sick state, you can either die or be discharged. the "plug" is remaining in that state. I  am honestly not sure if this will work.
  0, 0, 0, 0, 0, 1, ##death covers death through burial. it is not a terminal state - leaving it as a terminal state would count the death costs every remaining cycle, vastly overcounting QALY loss and burial cost
  0, 0, 0, 0, 0, 1) #terminal is the only terminal state.

mat_REGN_EB3 <- define_transition(
  state_names = c("induction", "triage_hi", "triage_lo", "conf_sick", "death", "terminal"),
  
  0, positive_algo, C, 0, 0, 0, ##you can only go from induction (where you start) to two states - hi or low triage. the positive algo gives a point estimate
  0, C, 0, screen_transition_prob*p_pos_hi, 0, screen_transition_prob*(1-p_pos_hi),  ##p_tx_screen*p_pos_hionce in triage_hi, the p_tx_screen gives the chances of moving to next stage by day, but multiplies this by chance of being positive. 
  0, 0, C, screen_transition_prob*p_pos_low, 0, screen_transition_prob*(1-p_pos_low), ##p_tx_screen*p_pos_low same as above but starting from low 
  0, 0, 0, C, death_transition_prob_REGNEB3, p_discharge_REGNEB3, ##there are two curves here from the confirmed sick state, you can either die or be discharged. the "plug" is remaining in that state. I  am honestly not sure if this will work.
  0, 0, 0, 0, 0, 1, ##death covers death through burial. it is not a terminal state - leaving it as a terminal state would count the death costs every remaining cycle, vastly overcounting QALY loss and burial cost
  0, 0, 0, 0, 0, 1) #terminal is the only terminal state.

###define state values###


state_induction <- define_state(
  cost_total = 0, 
  qaly = 0
)
## costs are still dummy values, BUT the following should be included
##SoC costs (Bartsh et. al) (all strategies, all states within the ETU)
##incremental cold chain costs (Mvundura et. al. WHO/UNICEF) (experimental arms)
## Incremental Training costs for mAb (conf_sick state in experimental arms, first day, spread across all pos patients)
##burial costs (for death state, all arms)

## qalys are negligible in every state but death, where the lost life years over the course of a lifetime are considered
## can absolutely code probabilistically based on life tables, have found literature of average age at induction

state_triage_hi <- define_state(
  cost_total = Avg_2024_SoC_Pt_Cost,
  qaly = 0
)

state_triage_low <- define_state(
  cost_total = Avg_2024_SoC_Pt_Cost_low,
  qaly = 0
)

state_terminal <- define_state(
  cost_total = 0,
  qaly = 0
)






## step 1 - Bartsch et. al. 2014 defines ppe and drug cost per episode for both those who recover and those who die, and duration of treatment.
##step 2: divide this cost for those who recovered (830/11.2) = 74.1 and those who die 321/4.2 = 76.4
##step 3: average out to 75
##step 4: convert to 2024 dollars.https://www.bls.gov/data/inflation_calculator.htm Decided to use exchange rates rather than PPP due to Bartsch methodology and variety of possible payors and locations

Per_Survivor_Cost = 830
Per_Dead_Cost = 321
Avg_2014_Pt_Cost = (Per_Survivor_Cost/11.2 + Per_Dead_Cost/4.2)/2

CPI_Multiplier_2014_2024 = 1.32

Avg_2024_SoC_Pt_Cost = CPI_Multiplier_2014_2024 * Avg_2014_Pt_Cost

### we are going to do same calculation for the not as severe supportive care from Bartsch 2014 so we can estimate SoC low triage

Per_Survivor_Cost_low = 446
Per_Dead_Cost_low = 185
Avg_2014_Pt_Cost_low = (Per_Survivor_Cost_low/11.2 + Per_Dead_Cost_low/4.2)/2

Avg_2024_SoC_Pt_Cost_low = CPI_Multiplier_2014_2024 * Avg_2014_Pt_Cost_low

CI_Per_Surv_cost = 862-800
CI_Per_Dead_cost = 351-292
CI_Per_Surv_cost_low = 466-428
CI_Per_Dead_cost_low = 202-169

##step 1: Wilson 2014 suggests ~100 patients on a continuing basis (80 when she arrived, 120 at peak). She mentions 102 direct care personnel (78/24 registered vs. certified) total
##step 2: assuming incremental training on mAbs will take all of these nurses 1 day of wages
##step 3: average wage calc - McCoy et al. 2008 - 4 different countries - assuming registered nurse closer to the 78 nurses and certified closer to assistants
##3a: Burkina Faso 2454 RN 1871 CN (2006) Ghana RN 4896 (2005) Zambia 4128 (2005) Nigeria 5097 (2001). Use BF ratio for other countries 

## the below have already been run through the CPI calculator in Jan of their respective year to Jan 23##

BF_RN_ann_wage = 3702
BF_CN_ann_wage = 2823
Zambia_RN_ann_wage = 6476
Nigeria_RN_ann_wage = 8709
Ghana_RN_ann_wage = 7680

BF_ratio_RN_CN = BF_RN_ann_wage/BF_CN_ann_wage

Zambia_CN_ann_wage = Zambia_RN_ann_wage/BF_ratio_RN_CN
Nigeria_CN_ann_wage = Nigeria_RN_ann_wage/BF_ratio_RN_CN
Ghana_CN_ann_wage = Ghana_RN_ann_wage/BF_ratio_RN_CN

##Step 4: do weighted average of CN/RN wage based on wilson paper

RN_pct_DC = 78/102

Ghana_avg_wage = RN_pct_DC*Ghana_RN_ann_wage + (1-RN_pct_DC)*Ghana_CN_ann_wage
Zambia_avg_wage = RN_pct_DC*Zambia_RN_ann_wage + (1-RN_pct_DC)*Zambia_CN_ann_wage
Nigeria_avg_wage = RN_pct_DC*Nigeria_RN_ann_wage + (1-RN_pct_DC)*Nigeria_CN_ann_wage
BF_avg_wage = RN_pct_DC*BF_RN_ann_wage + (1-RN_pct_DC)*BF_CN_ann_wage



Regional_avg_wage = mean(Ghana_avg_wage, Zambia_avg_wage, Nigeria_avg_wage, BF_avg_wage)

working_days_per_year = 250
##5 days per week, 50 weeks per year, accounting for travel time sickness etc.##


Regional_daily_wage = Regional_avg_wage/working_days_per_year

##step 5: add wage back up for all health workers, allocate it at a patient level, using total patients in WIlson paper



Total_ETU_Patients = 600
Confirmed_cases = prevalence* Total_ETU_Patients

Training_cost_pp = Regional_daily_wage*102/Confirmed_cases


###DEFINE STATES####

#First define DALYs so that we can place into distribution for sensitivity analysis later#




LY_point_est = 43.2

class(LY_dist)
#> Error in eval(expr, envir, enclos): object 'LY_dist' not found
###All Prev Figures taken from Shantha et. al.###
###All DW Figures taken from Salomon et. al.###

Mild_Prev = 14/191
Mild_DW = .004

Moderate_Prev = 2/191
Moderate_DW = .033

Blindness_Prev = 10/191
Blindness_DW = .195



DW_Vision_Impair = Mild_Prev*Mild_DW + Moderate_Prev*Moderate_DW + Blindness_Prev*Blindness_DW



state_sick_base <- define_state(
  cost_total = Avg_2024_SoC_Pt_Cost, 
  qaly = 0
)




###before defining the sick state, we need to allow it to only charge for mAb on day 1###

##going to first define Soc per day cost##

Tx_cost = 500

##this is a discreet decision point less than a probabilistic sensitivity. the types of people using this model as a decision tool often are also negotiating the price, so it can be an output of the decsion making process, not an input



TLC_Cost= 6890
Days_Operation = 120
Resupply_Cost_ETU = 92


Amortization_Rate = Days_Operation/(365*5)
Refrigeration_Cost = TLC_Cost * Amortization_Rate
CCL_cost_ETU = Resupply_Cost_ETU + Refrigeration_Cost

CCL_Cost_pp = CCL_cost_ETU/Confirmed_cases



##McCoy et. al.

par_mod <- modify(
  par_mod,
  
  cost_tx_start = Avg_2024_SoC_Pt_Cost + Tx_cost + CCL_Cost_pp + Training_cost_pp,
  cost_tx_end = Avg_2024_SoC_Pt_Cost,
  n_days = 1,
  cost_tx_cycle_mAb114 = ifelse(
    state_time < n_days-9.9,
    cost_tx_start,
    cost_tx_end)
)


state_sick_mAb114 <- define_state(
  cost_total = cost_tx_cycle_mAb114,
  qaly = 0
)

par_mod <- modify(
  par_mod,
  
  cost_tx_start = Avg_2024_SoC_Pt_Cost + Tx_cost + CCL_Cost_pp + Training_cost_pp,
  cost_tx_end = Avg_2024_SoC_Pt_Cost,
  n_days = 11,
  cost_tx_cycle_REGN_EB3 = ifelse(
    state_time < n_days-9.9, ## we first defined in n days in the probability matrix to cut off the KM curve at 10 days. For that reason we have to use it as a reference point
    cost_tx_start,
    cost_tx_end)
)


state_sick_REGN_Eb3 <- define_state(
  cost_total = cost_tx_cycle_REGN_EB3,
  qaly = 0
)


state_death <- define_state(
  cost_total = 50,
  qaly = (LY_point_est*-1)*(1-DW_Vision_Impair)
)


##define strategies##

strat_base <- define_strategy(
  transition = mat_base,
  
  induction = state_induction,
  triage_hi = state_triage_hi,
  triage_lo = state_triage_low,
  conf_sick = state_sick_base,
  death = state_death,
  terminal = state_terminal
)


strat_mAb114 <- define_strategy(
  transition = mat_mAb_114,
  
  induction = state_induction,
  triage_hi = state_triage_hi,
  triage_lo = state_triage_low,
  conf_sick = state_sick_mAb114,
  death = state_death,
  terminal = state_terminal
)

strat_REGN_EB3 <- define_strategy(
  transition = mat_REGN_EB3,
  
  induction = state_induction,
  triage_hi = state_triage_hi,
  triage_lo = state_triage_low,
  conf_sick = state_sick_REGN_Eb3,
  death = state_death,
  terminal = state_terminal
)

res_mod <- run_model(
  parameters = par_mod,
  
  mAb114 = strat_mAb114,
  REGN_EB3 = strat_REGN_EB3,
  base = strat_base,
  
  cycles = 40,
  
  cost = cost_total,
  effect = qaly,
  
  method = "life-table"
)
#> mAb114: detected use of 'state_time', expanding states: conf_sick, triage_hi, triage_lo.
#> REGN_EB3: detected use of 'state_time', expanding states: conf_sick, triage_hi, triage_lo.
#> base: detected use of 'state_time', expanding states: conf_sick, triage_hi, triage_lo.



summary(res_mod, threshold = c(1000,5000,15000))
#> 3 strategies run for 40 cycles.
#> 
#> Initial state counts:
#> 
#> induction = 1000L
#> triage_hi = 0L
#> triage_lo = 0L
#> conf_sick = 0L
#> death = 0L
#> terminal = 0L
#> 
#> Counting method: 'life-table'.
#> 
#> Values:
#> 
#>          cost_total    qaly
#> mAb114     935086.6 -7664.1
#> REGN_EB3   935086.6 -7664.1
#> base       719247.6 -7664.1
#> 
#> Net monetary benefit difference:
#> 
#>             1000    5000   15000
#> mAb114     0.000   0.000   0.000
#> REGN_EB3   0.000   0.000   0.000
#> base     215.839 215.839 215.839
#> 
#> Efficiency frontier:
#> 
#> base
#> 
#> Differences:
#> 
#>          Cost Diff. Effect Diff. ICER   Ref.
#> REGN_EB3      0.000            0  NaN mAb114
#> base       -215.839            0 -Inf mAb114




def_psa <- define_psa(
  LY_point_est ~ triangle(lower = 12, 65, peak = LY_point_est),
  Resupply_Cost_ETU ~ triangle(lower = 0, upper = 1000000, peak = Resupply_Cost_ETU)
  ) 
#> Loading required namespace: triangle
  

res_psa <- run_psa(res_mod, psa = def_psa, N = 20)
#> Resampling strategy 'mAb114'...
#> Resampling strategy 'REGN_EB3'...
#> Resampling strategy 'base'...

plot(res_psa, type = "ce")

^{Created on 2024-03-03 with reprex v2.1.0}

[mmaughan0]

The below is what you get when I use my actual code with my real survival tables

[mmaughan0]

I figured it out - kind of. It seems that the psa only works on variables explicitly defined in the "total_cost" or "qaly" step when defining each health state. Anything used to derive that will not work as a sensitivity.