Simulate health state probabilities from a survival object using partitioned survival analysis.

# S3 method for survival
sim_stateprobs(x, ...)

## Arguments

x An object of class survival. Further arguments passed to or from other methods.

## Value

A stateprobs object.

## Details

In an $$N$$-state partitioned survival model there are $$N-1$$ survival curves and $$S_n(t)$$ is the cumulative survival function denoting the probability of survival to health state $$n$$ or a lower indexed state beyond time $$t$$. The probability that a patient is in health state 1 at time $$t$$ is simply $$S_1(t)$$. State membership in health states $$2,\ldots, N -1$$ is calculated as $$S_n(t) - S_{n-1}(t)$$. Finally, the probability of being in the final health state $$N$$ (i.e., the death state) is $$1-S_{N-1}(t)$$, or one minus the overall survival curve.

In some cases, the survival curves may cross. hesim will issue a warning but the function will still run. Probabilities will be set to 0 in a health state if the prior survival curve lies above the curve for state $$n$$; that is, if $$S_n(t) < S_{n-1}(t)$$, then the probability of being in state $$n$$ is set to 0 and $$S_n(t)$$ is adjusted to equal $$S_{n-1}(t)$$. The probability of being in the final health state is also adjusted if necessary to ensure that probabilities sum to 1.

library("data.table") library("survival") # This example shows how to simulate a partitioned survival model by # manually constructing a "survival" object. We will consider a case in which # Cox proportional hazards models from the survival package---which are not # integrated with hesim---are used for parameter estimation. We will use # point estimates in the example, but bootstrapping, Bayesian modeling, # or other techniques could be used to draw samples for a probabilistic # sensitivity analysis. # (0) We first setup our model per usual by defining the treatment strategies, # target population, and health states hesim_dat <- hesim_data( strategies = data.table(strategy_id = 1:3, strategy_name = c("SOC", "New 1", "New 2")), patients = data.table(patient_id = 1:2, female = c(0, 1), grp_id = 1), states = data.table(state_id = 1:2, state_name = c("Stable", "Progression")) ) # (1) Next we will estimate Cox models with survival::coxph(). We illustrate # by predicting progression free survival (PFS) and overall survival (OS) ## Fit models onc3_pfs_os <- as_pfs_os(onc3, patient_vars = c("patient_id", "female", "strategy_name")) fit_pfs <- coxph(Surv(pfs_time, pfs_status) ~ strategy_name + female, data = onc3_pfs_os) fit_os <- coxph(Surv(os_time, pfs_status) ~ strategy_name + female, data = onc3_pfs_os) ## Predict survival on input data surv_input_data <- expand(hesim_dat) times <- seq(0, 14, 1/12) predict_survival <- function(object, newdata, times) { surv <- summary(survfit(object, newdata = newdata, se.fit = FALSE), t = times) pred <- newdata[rep(seq_len(nrow(newdata)), each = length(times)), ] pred[, sample := 1] # Point estimates only in this example pred[, time := rep(surv$time, times = nrow(newdata))] pred[, survival := c(surv$surv)] return(pred[, ]) } pfs <- predict_survival(fit_pfs, newdata = surv_input_data, times = times) os <- predict_survival(fit_os, newdata = surv_input_data, times = times) surv <- rbind( as.data.table(pfs)[, curve := 1L], as.data.table(os)[, curve := 2L] ) ## Convert predictions to a survival object surv <- survival(surv, t = "time") if (FALSE) autoplot(surv) # (2) We can then compute state probabilities from the survival object stprobs <- sim_stateprobs(surv) # (3) Finally, we can use the state probabilities to compute QALYs and costs ## A dummy utility model to illustrate utility_tbl <- stateval_tbl( data.table(state_id = 1:2, est = c(1, 1) ), dist = "fixed" ) utilitymod <- create_StateVals(utility_tbl, hesim_data = hesim_dat, n = 1) ## Instantiate Psm class and compute QALYs psm <- Psm$new(utility_model = utilitymod) psm$stateprobs_ <- stprobs psm$sim_qalys() psm$qalys_