1 Overview

An N-state partitioned survival model (PSM) simulates the probabilitythat a patient is in each of N distinct health states at a given pointof time when treated by a particular therapy. State membership isestimated from a set of non-mutually exclusive survival curves; for anN-state model, N-1 survival curves are required.

The cumulative survival function,$S_n(t)$, represents the probability that apatient survives to health state$n$or to a lower indexed state beyond time$t$. The probability that a patient is inhealth state 1 is always just$S_1(t)$. State membership in health states$2,\ldots, n-1$ is calculated as$S_{n}(t) - S_{n-1}(t)$. Finally, theprobability of being in the final health state$N$ (i.e., the dead state) is$1-S_{N-1}(t)$, or one minus the overallsurvival curve.

Inhesim, an N-state PSM consists of three types ofseparate statistical models: a set of N-1 survival models, a model ofthe utility values to assign to each health state, and a set of modelsfor costs to assign to each health state. Multiple cost models can beused to categorize costs into different categories (e.g., costs ofmedical care, drug costs). All models can either be fit using R orconstructed with external data.

Separate survival curves as a function of time,$t$, are predicted for each treatmentstrategy$k$ and patient$i$ as a function of$L$ parameters,$\alpha_l$, which can depend on covariates,$x_l$ (seeparams_surv()for more details),

\[\begin{aligned}S_{n}(t|\alpha_{1n}(x_{1n, ik}), \ldots, \alpha_{Ln}(x_{Ln, ik})).\end{aligned}\] Quality-adjusted life-years (QALYs) and total costs associatedwith a given health state for a given treatment strategy and patient arecalculated by integrating the “weighted” probability of survival, whereweights are a function of the discount factor and predicted statevalues. Letting$h$ index a healthstate, QALYs and total costs for the$m$th cost category are calculated as,

\[\begin{aligned}\rm{QALYs}_{hik} &= \int_{0}^{T} q_{hik}(t) e^{-rt} p_{hik}(t)dt, \\\rm{Costs}_{m,hik} &= \int_{0}^{T} c_{m,hik}(t) e^{-rt}p_{hik}(t)dt,\end{aligned}\] where$q_{hik}$ is aquality-of-life weight,$c_{m,hik}$is (assuming that time is measured in years) annualized costs,$r$ is the discount rate, and$p_{hik}(t)$ is the probability of being ina given health state at a given time, and QALYs and costs are calculatedover$T$ time periods.

2 An example 4-statePSM

To illustrate a partitioned survival analysis, consider theevaluation of a two-line sequential treatment strategy in oncology usinga 4-state PSM. The 4 health states might be:

Progression free on first line treatment
Progression free on second line treatment
Progressed on second line treatment
Death

We begin by defining the population, treatment strategies, and modelstructure. In this example, we model three treatment strategies, threedistinct patients (that differ in age and gender), and four healthstates (death and three non-death states). This information is combinedusinghesim_data(), which creates a general object of classhesim_data.

library("data.table")library("hesim")set.seed(101)strategies<-data.table(strategy_id =c(1,2,3),strategy_name =paste0("Strategy ",1:3))patients<-data.table(patient_id =seq(1,3),age =c(45,50,60),female =c(0,0,1))states<-data.frame(state_id =seq(1,3),state_name =c("Progression free (1st line)","Progression free (2nd line)","Progressed (2nd line)"),stringsAsFactors =FALSE)hesim_dat<-hesim_data(strategies = strategies,patients = patients,states = states)print(hesim_dat)

## $strategies##    strategy_id strategy_name##          <num>        <char>## 1:           1    Strategy 1## 2:           2    Strategy 2## 3:           3    Strategy 3## ## $patients##    patient_id   age female##         <int> <num>  <num>## 1:          1    45      0## 2:          2    50      0## 3:          3    60      1## ## $states##   state_id                  state_name## 1        1 Progression free (1st line)## 2        2 Progression free (2nd line)## 3        3       Progressed (2nd line)## ## attr(,"class")## [1] "hesim_data"

We useget_labels() to obtain nice labels for plots andsummary tables.

labs<-get_labels(hesim_dat)print(labs)

## $strategy_id## Strategy 1 Strategy 2 Strategy 3 ##          1          2          3 ## ## $state_id## Progression free (1st line) Progression free (2nd line) ##                           1                           2 ##       Progressed (2nd line)                       Death ##                           3                           4

3 Parameterestimation

3.1 Survival models

The parameters of the survival models would typically be estimated inone of two ways. First, if access to patient level data from a clinicaltrial were available, survival models would be fit using available Rpackages. Second, a formal evidence synthesis—such as a networkmeta-analysis—might be conducted. Here, we provide an example of ananalysis with trial level data.

We consider an estimation dataset (surv_est_data) withthree survival endpoints with endpoint 1 denoting the lowest indexedstate, endpoint 2 the next indexed state, and endpoint 3 representingoverall survival. In our example, endpoint 1 represents progression on1st line treatment, endpoint 2 represents progression on 2nd linetreatment, and endpoint 3 represents death.

Survival models inhesim can be fit using eitherflexsurvreg orflexsurvspline from theflexsurv package.hesim currently supportsparametric (exponential, Weibull, Gompertz, gamma, log-logistic,lognormal, and generalized gamma), splines, and fractional polynomialsurvival models (seeparams_surv).In this example, we fit a Weibull model to each of the three survivalendpoints.

library("flexsurv")surv_est_data<- psm4_exdata$survivalfit1<-flexsurvreg(Surv(endpoint1_time, endpoint1_status)~ age+ female+factor(strategy_id),data = surv_est_data,dist ="weibull")fit2<-flexsurvreg(Surv(endpoint2_time, endpoint2_status)~ age+ female+factor(strategy_id),data = surv_est_data,dist ="weibull")fit3<-flexsurvreg(Surv(endpoint3_time, endpoint3_status)~ age+ female+factor(strategy_id),data = surv_est_data,dist ="weibull")psfit_wei<-flexsurvreg_list(fit1, fit2, fit3)

3.2 Utility and costmodels

Although utility and cost models can be fit using regression models,in many cases, we won’t fit a model directly, but will use estimatedmean values from published sources. In this case, we can use astateval_tbl() to store estimates. For instance, consideran example where we are only given a range of values for utility and wetherefore approximate the posterior distribution of mean utility valuesby health state using a uniform distribution.

utility_tbl<-stateval_tbl(tbl =data.frame(state_id = states$state_id,min = psm4_exdata$utility$lower,max = psm4_exdata$utility$upper),dist ="unif")

Similar to utility, we will not simulate drug costs based on astatistical model but will instead use published estimates. Drug costswill only vary by treatment strategy (i.e., they do not vary by healthstate or across samples in the PSA).

drugcost_tbl<-stateval_tbl(tbl =data.frame(strategy_id = strategies$strategy_id,est = psm4_exdata$costs$drugs$costs),dist ="fixed")

For cases in which we want to use a regression model to simulatecosts,hesim supports parameterizing a cost or utilitymodel using linear models. Here, we fit a model for medical costs as afunction of the three non-death health states using thelm() function in R.

medcosts_fit<-lm(costs~ female+ state_name,data = psm4_exdata$costs$medical)

4 Simulation

4.1 Constructing theeconomic model

A PSM inhesim is anR6 object of classPsm and comprises of three separate submodels: (1) a set ofsurvival models for generating survival curves (of classPsmCurves), (2) a utility model (of classStateVals), and (3) a set of cost models for each costcomponent (a list ofStateVals objects). Since analyses inhesim are, by default, based on sampled values of theparameters from their joint posterior distribution in order tofacilitate probabilistic sensitivity analysis (PSA), we will specify thenumber of sampled parameter sets that are desired.

n_samples<-100

4.1.1 Survivalmodels

Survival curves are simulated usingPsmCurves objects.We can create aPsmCurves object usingcreate_PsmCurves() as a function of our fitted Weibullmodels and input data describing the target population and treatmentstrategies.

To conduct the PSA, we must sample random draws of the regressioncoefficients from suitable probability distribution, which can be donein one of two ways. First, the parameters of each survival model can besampled separately using a multivariate normal distribution. However,this option does not preserve the correlation in the survival endpoints.By default,create_PsmCurves() consequently samples theparameters via bootstrapping, whereby the survival models are refitrepeatedly to resamples of the estimation dataset.

surv_input_data<-expand(hesim_dat,by =c("strategies","patients"))survmods<-create_PsmCurves(psfit_wei,input_data = surv_input_data,n = n_samples,uncertainty ="bootstrap",est_data = surv_est_data)

4.1.2 Utility and costmodels

Utility and cost values are simulated using models of classStateVals. We can instantiate aStateValsobject from a ’stateval_tbl` object, which creates a “mean model” thatpredicts mean state values for each treatment strategy, patient, andhealth state combination. In the utility model, utility only variesacross health states and is constant within patients and treatmentstrategies conditional on health state.

utilitymod<-create_StateVals(utility_tbl,n = n_samples,hesim_data = hesim_dat)

We take a similar approach for drug costs, in which costs only varyacross treatment strategies.

drugcostmod<-create_StateVals(drugcost_tbl,n = n_samples,hesim_data = hesim_dat)

Finally, we can also instantiate aStateVals() from afitted linear model and a corresponding dataset of classexpanded_hesim_data.

medcost_data<-expand(hesim_dat,by =c("strategies","patients","states"))medcostmod<-create_StateVals(medcosts_fit,input_data = medcost_data,n = n_samples)

4.1.3 Combining thestatistical models

The PSM model is instantiated using$new().

psm<- Psm$new(survival_models = survmods,utility_model = utilitymod,cost_models =list(medical = medcostmod,drug = drugcostmod))

4.2 Simulatingoutcomes

4.2.1 Simulating survivalcurves

Survival curves are predicted by treatment strategy and patient withthe function$sim_survival(). We can quickly generate aplot of survival curves (averaged across patients) by treatment strategyusingautoplot.survival().

# Simulatetimes<-seq(0,12,by = .1)psm$sim_survival(t = times)head(psm$survival_)

##    sample strategy_id patient_id grp_id curve     t  survival##     <num>       <int>      <int>  <int> <num> <num>     <num>## 1:      1           1          1      1     1   0.0 1.0000000## 2:      1           1          1      1     1   0.1 0.9189061## 3:      1           1          1      1     1   0.2 0.8411512## 4:      1           1          1      1     1   0.3 0.7688137## 5:      1           1          1      1     1   0.4 0.7019982## 6:      1           1          1      1     1   0.5 0.6405107

# Plotlibrary("ggplot2")labs$curve<-c("Curve 1"=1,"Curve 2"=2,"Curve 3"=3)autoplot(psm$survival_,labels = labs,ci =TRUE,ci_style ="ribbon")+scale_x_continuous(breaks =seq(0,max(times),2))

4.2.2 Simulating healthstate probabilities

After simulating survival, we can calculate the probability of beingin each of the four health states using$sim_stateprobs().At a given point in time, the probability of being in state 1 is thesurvival probability from the first survival curve, the probability ofbeing in state 2 is the difference in survival probabilities between the2nd and 1st curves, the probability of being in state 3 is thedifference in survival probabilities between the 3rd and 2nd curves, andthe probability of death is 1 minus the survival probability from curve3.

While plots can be quickly generated withautoplot.stateprobs(), the output is returned in a tidyformat (and stored as a field in the class namedstateprobs_) that facilitates creation of custom plots.We’ll plot state probabilities for the 1st patient with the 30thrandomly sampled parameter set.

# Simulatepsm$sim_stateprobs()head(psm$stateprobs_)

##    sample strategy_id patient_id grp_id state_id     t      prob##     <num>       <int>      <int>  <int>    <int> <num>     <num>## 1:      1           1          1      1        1   0.0 1.0000000## 2:      1           1          1      1        1   0.1 0.9189061## 3:      1           1          1      1        1   0.2 0.8411512## 4:      1           1          1      1        1   0.3 0.7688137## 5:      1           1          1      1        1   0.4 0.7019982## 6:      1           1          1      1        1   0.5 0.6405107

# Plotstateprobs<- psm$stateprobs_[sample==30& patient_id==1]stateprobs[, state:=factor(state_id,levels =rev(unique(state_id)))]stateprobs[, strategy:=factor(strategy_id,labels =names(labs$strategy_id))]ggplot(stateprobs[strategy_id%in%c(1,2,3)],aes(x = t,y = prob,fill = state,group = state))+geom_area(alpha = .65)+facet_wrap(~strategy)+xlab("Years")+ylab("Proportion in state")+scale_fill_manual(name ="Health state",values =c("gray92","green4","orange","purple"),labels =c("Death",rev(hesim_dat$states$state_name)))+scale_x_continuous(breaks =seq(0,max(times),2))+guides(fill =guide_legend(reverse =TRUE,nrow =2,byrow =TRUE))+theme(legend.position ="bottom")

4.2.3 Simulating costsand QALYs

Finally, we can simulate discounted costs and QALYs by numericallyintegrating the “weighted” probabilities in$stateprobs_ asdescribed above. Costs and QALYs are produced for each treatmentstrategy, patient, health state, and sampled parameter set.

# Costspsm$sim_costs(dr = .03)head(psm$costs_)

##    sample strategy_id patient_id grp_id state_id    dr category    costs##     <num>       <int>      <int>  <int>    <int> <num>   <char>    <num>## 1:      1           1          1      1        1  0.03  medical 28707.06## 2:      1           1          1      1        2  0.03  medical 24729.36## 3:      1           1          1      1        3  0.03  medical 24914.92## 4:      1           1          2      1        1  0.03  medical 27384.46## 5:      1           1          2      1        2  0.03  medical 23700.48## 6:      1           1          2      1        3  0.03  medical 23230.48

# QALYspsm$sim_qalys(dr = .03)head(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 0.8913423 1.0479377## 2:      1           1          1      1        2  0.03 0.5448242 0.7391230## 3:      1           1          1      1        3  0.03 0.5277850 0.7778092## 4:      1           1          2      1        1  0.03 0.8502762 0.9996569## 5:      1           1          2      1        2  0.03 0.5221563 0.7083712## 6:      1           1          2      1        3  0.03 0.4921027 0.7252234

Movatterモバイル変換

Partitioned survival models

2025-10-10