R/sim-general.R
sim_stateprobs.survival.RdSimulate health state probabilities from a survival object using partitioned
survival analysis.
# S3 method for class 'survival'
sim_stateprobs(x, ...)An object of class survival.
Further arguments passed to or from other methods.
A stateprobs object.
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) # \dontrun{}
# (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_
#> sample strategy_id patient_id grp_id state_id dr qalys lys
#> <num> <int> <int> <int> <int> <num> <num> <num>
#> 1: 1 1 1 1 1 0.03 4.874007 4.874007
#> 2: 1 1 1 1 2 0.03 3.070008 3.070008
#> 3: 1 1 2 1 1 0.03 4.267061 4.267061
#> 4: 1 1 2 1 2 0.03 3.166468 3.166468
#> 5: 1 2 1 1 1 0.03 5.640397 5.640397
#> 6: 1 2 1 1 2 0.03 2.605955 2.605955
#> 7: 1 2 2 1 1 0.03 4.961563 4.961563
#> 8: 1 2 2 1 2 0.03 2.789254 2.789254
#> 9: 1 3 1 1 1 0.03 6.141992 6.141992
#> 10: 1 3 1 1 2 0.03 2.399936 2.399936
#> 11: 1 3 2 1 1 0.03 5.426923 5.426923
#> 12: 1 3 2 1 2 0.03 2.639381 2.639381