Markov Cost Effectiveness Model - Surival curve is using a (parametric) best fit line. I need it to use the (non-parametric) Kaplan Meier curve

Consult reprex below, but my issue is not that the code will not run - its that doing something I do not want it to do.

In the disease I am modeling, the risk of death greatly decreases after a certain amount of days, and thus the Kaplan_meier/survival curve flattens out. I want my cost effectiveness model to follow the probabilities related to that curve. What the model does is estimate a parametric best fit line that CONTINUES to decline after that 10 day period. I very much do not want it to do this, and I thought that is what I was specifying by the "km-limit" in the p_death_disease_xx objects. It is not. I have tested it out by altering the cycles and although it should stop reducing the risk of death after those 10 cycles, it continues to add risk of death, clearly following that parametric best fit curve.

Is there anything I can do to specify it NOT to do that, and to only follow the KM/non-parametric curve?

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

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: 0x000001b08c87fe80>
#> <environment: namespace:base>

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

par_mod <- define_parameters(

  p_death_disease_base = compute_surv(
    fit_death_disease_base,
    time = state_time,
    km_limit = 14))

par_mod <- modify(
  par_mod,

  p_death_disease_mAb = compute_surv(
    fit_death_disease_mAb,
    time = state_time,
    km_limit = 14))

#establishing dummy table only for purpose of reprex. In reality each strategy would have a different table

tab_surv_base <- 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, 10, 10, 10,10,10,12,14, 15, 15, 15), status = c(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")

##code below creates a variable in which mortality curves reside. This variable (p_death_disease_strat) 
##can then be used in transition probability matrices etc.##

##pulling in induct to treat table##

plot(fit_death_disease_mAb)
#> Error in eval(expr, envir, enclos): object 'fit_death_disease_mAb' not found

##This code is what imports WHO life tables##

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

tab_surv_mAb <- 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), status = c(1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
                                                                            0L, 0L, 0L )), .Names = c("time","status"), row.names = c(NA, -14L), class = "data.frame")

fit_death_disease_mAb <- flexsurv::flexsurvreg(
  survival::Surv(time, status) ~ 1,
  dist = "weibull",
  data = tab_surv_mAb)

###define matrices####

##putting in point estimates from Levine et. al. paper. All diagnostic rates can be derived from these 3 data points##

p_hi_triage <-.7
p_pos_hi <- .8
p_pos_low <-.3

##define base transition matrix##

plot(fit_death_disease_base)


mat_base <- define_transition(
  state_names = c("unknown", "triage_hi", "triage_lo", "conf_sick", "death", "terminal"),

  0, p_hi_triage, C, 0, 0, 0,
  0, 0, C, p_pos_hi, 0, 0, 
  0, 0, C, p_pos_low, 0, 0,
  0, 0, 0, C, p_death_disease_base, 0,
  0, 0, 0, 0, 0, 1,
  0, 0, 0, 0, 0, 1)

mat_mAb<- define_transition(
  state_names = c("unknown", "triage_hi", "triage_lo", "conf_sick", "death", "terminal"),

  0, p_hi_triage, C, 0, 0, 0,
  0, 0, C, p_pos_hi, 0, 0, 
  0, 0, C, p_pos_low, 0, 0,
  0, 0, 0, C, p_death_disease_mAb, 0,
  0, 0, 0, 0, 0, 1,
  0, 0, 0, 0, 0, 1)

###define state values###

##unknown represents induction - is a transient state that occurs between presentation and triage. 
##Patients do not spend a full cycle (day) there,so no qaly assigned
##some marginal ppe usage assumed as well as clinician time that always has opportunity cost in this setting

state_unknown <- define_state(
  cost_total = 50,
  qaly = 0)

state_hi <- define_state(
  cost_total = 150,
  qaly = 1/365)

state_low <- define_state(
  cost_total =100, 
  qaly = 1/365)

state_sick_base <- define_state(
  cost_total = 200,
  qaly = 1/365
)

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

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

par_mod <- modify(
  par_mod,

  cost_tx_start = 300,
  cost_tx_end = 200,
  n_days = 2,
  cost_tx_cycle = ifelse(
    state_time < n_days,
    cost_tx_start,
  cost_tx_end)
  )

state_sick_mAb <- define_state(
  cost_total = cost_tx_cycle,
  qaly = 1/365
)

state_death <- define_state(
  cost_total = 50,
  qaly = -30
)

strat_base <- define_strategy(
  transition = mat_base,

  unknown = state_unknown,
  triage_hi = state_hi,
  triage_lo = state_low,
  conf_sick = state_sick_base,
  death = state_death,
  terminal = state_terminal
)

strat_mAb <- define_strategy(
  transition = mat_mAb,

  unknown = state_unknown,
  triage_hi = state_hi,
  triage_lo = state_low,
  conf_sick = state_sick_mAb,
  death = state_death,
  terminal = state_terminal
)

res_mod <- run_model(
  parameters = par_mod,

  mAb = strat_mAb,
  base = strat_base,

  cycles = 14,

  cost = cost_total,
  effect = qaly,

  method = "life-table"
)
#> mAb: detected use of 'state_time', expanding state: conf_sick.
#> base: detected use of 'state_time', expanding state: conf_sick.

summary(res_mod, threshold = c(1000, 5000, 15000))
#> 2 strategies run for 14 cycles.
#> 
#> Initial state counts:
#> 
#> unknown = 1000L
#> triage_hi = 0L
#> triage_lo = 0L
#> conf_sick = 0L
#> death = 0L
#> terminal = 0L
#> 
#> Counting method: 'life-table'.
#> 
#> Values:
#> 
#>      cost_total      qaly
#> mAb     1431582 -26425.86
#> base    1810375 -16237.00
#> 
#> Net monetary benefit difference:
#> 
#>          1000     5000    15000
#> mAb     0.000     0.00      0.0
#> base 9810.064 50565.49 152454.1
#> 
#> Efficiency frontier:
#> 
#> mAb -> base
#> 
#> Differences:
#> 
#>      Cost Diff. Effect Diff.     ICER Ref.
#> base   378.7928     10.18886 37.17716  mAb

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

r-lib / vctrs

Markov Cost Effectiveness Model - Surival curve is using a (parametric) best fit line. I need it to use the (non-parametric) Kaplan Meier curve #1911