vignettes/markov-inhomogeneous.Rmd
markov-inhomogeneous.Rmd
We will illustrate a time inhomogeneous Markov cohort model by replicating the total hip replacement (THR) example from the Decision Modeling for Health Economic Evaluation textbook. The analysis compares two treatment strategies, a “standard” prosthesis and a “new” prosthesis. A probabilistic sensitivity analysis (PSA) will be used to quantify parameter uncertainty.
The model consists 5 health states: (i) primary THR, (ii) successful primary, (iii) revision THR, (iv) successful revision, and (v) death.
The model is time-inhomogeneous because (i) the prosthesis survival time is modeled using a Weibull survival model so that the probability of failure is time-varying and (ii) background mortality rates increase as patients age.
We set up the model for two treatment strategies and 10 subgroups defined by age and sex.
library("hesim") library("data.table") # Treatment strategies strategies <- data.table(strategy_id = 1:2, strategy_name = c("Standard prosthesis", "New prosthesis")) # Patients ages <- seq(55, 75, 5) age_weights <- c(.05, .1, .4, .25, .20) gender_weights <- c(.65, .35) weights <- rep(age_weights, times = 2) * rep(gender_weights, each = length(ages)) patients <- data.table(patient_id = 1:10, grp_id = 1:10, gender = rep(c("Female", "Male"), each = length(ages)), age = rep(ages, times = 2), patient_wt = weights) hesim_dat <- hesim_data(strategies = strategies, patients = patients) print(hesim_dat)
## $strategies
## strategy_id strategy_name
## 1: 1 Standard prosthesis
## 2: 2 New prosthesis
##
## $patients
## patient_id grp_id gender age patient_wt
## 1: 1 1 Female 55 0.0325
## 2: 2 2 Female 60 0.0650
## 3: 3 3 Female 65 0.2600
## 4: 4 4 Female 70 0.1625
## 5: 5 5 Female 75 0.1300
## 6: 6 6 Male 55 0.0175
## 7: 7 7 Male 60 0.0350
## 8: 8 8 Male 65 0.1400
## 9: 9 9 Male 70 0.0875
## 10: 10 10 Male 75 0.0700
##
## attr(,"class")
## [1] "hesim_data"
The model will be simulated for each treatment strategy and patient combination. Data containing treatment strategy and patient characteristics can be generated using expand()
.
## strategy_id patient_id strategy_name grp_id gender age patient_wt
## 1: 1 1 Standard prosthesis 1 Female 55 0.0325
## 2: 1 2 Standard prosthesis 2 Female 60 0.0650
## 3: 1 3 Standard prosthesis 3 Female 65 0.2600
## 4: 1 4 Standard prosthesis 4 Female 70 0.1625
## 5: 1 5 Standard prosthesis 5 Female 75 0.1300
## 6: 1 6 Standard prosthesis 6 Male 55 0.0175
Transition probabilities from state \(i\) (rows) to state \(j\) (columns) are determined using the parameters shown in the table below:
Primary THR | Successful primary | Revision THR | Successful revision | Death | |
---|---|---|---|---|---|
Primary THR | 0 | C | 0 | 0 | omrPTHR |
Successful primary | 0 | C | rr | 0 | mr |
Revision THR | 0 | 0 | 0 | C | omrRTHR + mr |
Successful revision | 0 | 0 | rrr | C | mr |
Death | 0 | 0 | 0 | 0 | 1 |
Three parameters of the transition probability matrix are assumed to be constant in the model: (i) the operative mortality rate following primary THR, omrPTHR, (ii) the operative mortality rate following revision THR, omrRTHR, and the re-revision rate, rrr. In a sample of 100 patients receiving primary THR 2 died implying that omrPTHR can be characterized by a beta distribution with \(\alpha = 2\) and \(\beta = 98\). Similarly, in a sample of 100 patients experiencing a revision procedure, four patients had another procedure within one year suggesting that rrr can be characterized by a beta distribution with \(\alpha = 4\) and \(\beta = 96\). Finally. there was information available for omrRTHR, so it assumed to follow the same beta distribution (\(\alpha = 2\) and \(\beta = 98\)) as omrPTHR.
The remaining parameters, rr and mr, are time-varying. The yearly mortality rates, mr, are stratified by age and sex.
print(mort_tbl)
## age_lower age_upper male female
## 1 35 45 0.00151 0.00099
## 2 45 55 0.00393 0.00260
## 3 55 65 0.01090 0.00670
## 4 65 75 0.03160 0.01930
## 5 75 85 0.08010 0.05350
## 6 85 Inf 0.18790 0.15480
Revision rates, rr, were modeled using a proportional hazards Weibull model. The scale parameter was modeled as a function of age and indicators for male sex and and whether a new prosthesis was used. That is, given a scale parameter \(\lambda\) and a shape parameter \(\gamma\), the survival function is,
\[ S(t) = \exp(-\lambda t^\gamma) \\ \lambda = \exp(\beta_0 + \beta_1 \cdot age + \beta_2 \cdot male + \beta_3 \cdot new) \\ \]
Note that this Weibull distribution corresponds to the distribution implemented in flexsurv::pweibullPH()
with \(a = \gamma\) and \(m = \lambda\). It can be shown (see Appendix below) that the transition probabilities (with a one-year model cycle) are given by,
\[
1 - \exp \left\{\lambda\left[(t-1)^\gamma - t^\gamma \right] \right\}
\] The coefficients from the regression model and the variance-covariance matrix used for the PSA are stored in rr_coef
and rr_vcov
, respectively.
print(rr_coef)
## lngamma cons age male np1
## 0.3740968 -5.4909350 -0.0367022 0.7685360 -1.3444740
print(rr_vcov)
## lngamma cons age male np1
## lngamma 0.002251512 -0.00569100 2.800000e-08 0.00000510 0.0002590
## cons -0.005691000 0.04321908 -7.830000e-04 -0.00724700 -0.0006420
## age 0.000000028 -0.00078300 2.715661e-05 0.00003300 -0.0001110
## male 0.000005100 -0.00724700 3.300000e-05 0.01189539 0.0001840
## np1 0.000259000 -0.00064200 -1.110000e-04 0.00018400 0.1463686
The mean (standard error) of utility was estimated to be 0.85 (0.03), 0.30 (0.03), and 0.75 (0.04) in the successful primary, revision, and successful revision health states, respectively.
The largest costs are the cost of the prostheses themselves. The standard prosthesis costs \(£394\) while the new prosthesis costs \(£579\). Both are assumed to be known with certainty.
The model assumes that there are no ongoing medical costs. The only remaining cost is therefore the cost of the revision procedure, which was estimated to have a mean of \(£5,294\) and standard error of \(£1,487\).
All underlying parameter estimates are stored in a list.
params <- list( # Transition probabilities ## Operative mortality following primary THR omrPTHR_shape1 = 2, omrPTHR_shape2 = 98, ## Revision rate for prosthesis rr_coef = rr_coef, rr_vcov = rr_vcov, ## Mortality_rates mr = mort_tbl, ## Operative mortality following revision THR omrRTHR_shape1 = 2, omrRTHR_shape2 = 98, ## re-revision rate rrr_shape1 = 4, rrr_shape2 = 96, # Utility u_mean = c(PrimaryTHR = 0, SuccessP = .85, Revision = .30, SuccessR = .75), u_se = c(PrimaryTHR = 0, SuccessP = .03, Revision = .03, SuccessR = .04), # Costs c_med_mean = c(PrimaryTHR = 0, SuccessP = 0, Revision = 5294, SuccessR = 0), c_med_se = c(PrimaryTHR = 0, SuccessP = 0, Revision = 1487, SuccessR = 0), c_Standard = 394, c_NP1 = 579 )
As noted above, omrPTHR, omrRTHR, and rrr are drawn from beta distributions with \(\alpha\) and \(\beta\) (i.e. shape1
and shape2
) specified. Similarly, utility is drawn from a beta distribution, but shape1
and shape2
are derived from the mean and standard error using the method of moments. The mortality rate and the cost of the prostheses are assumed to be known with certainty. The medical costs associated with health states are drawn from a gamma distribution, for which, like utility, the underlying parameters are derived from the mean and standard error using the method of moments. Finally, the parameters of the Weibull survival model are drawn from a multivariate normal distribution using the point estimates and the variance-covariance matrix. 500 samples will be drawn for the PSA.
rng_def <- define_rng({ list( omrPTHR = beta_rng(shape1 = omrPTHR_shape1, shape2 = omrPTHR_shape2), rr_coef = multi_normal_rng(mu = rr_coef, Sigma = rr_vcov), mr_male = fixed(mr$male, names = mr$age_lower), mr_female = fixed(mr$female, names = mr$age_lower), omrRTHR = beta_rng(shape1 = omrRTHR_shape1, shape2 = omrRTHR_shape2), rrr = beta_rng(shape1 = rrr_shape1, shape2 = rrr_shape2), u = beta_rng(mean = u_mean, sd = u_se), c_med = gamma_rng(mean = c_med_mean, sd = c_med_se), c_Standard = c_Standard, c_NP1 = c_NP1 ) }, n = 500)
The sampled parameter values are “transformed” as a function of the input data. When evaluating the define_tparams()
expression, operations are vectorized and performed for each combination of the sampled parameters, input data, and time periods. The model will be simulated for 60 model cycles (i.e., 60 years) so we will create 60 time intervals of length 1.
Separate transformation functions are used for the transition model and the cost/utility models. This is done for computational efficiency since only the transition model depends on cycle time. In the transition model, the revision rate (rr) depends on the scale and the shape parameters, which, in turn, depend on the sampled values of the parameters from the Weibull model. The mortality rate depends on gender and a patient’s age during each model cycle.
transitions_def <- define_tparams({ # Regression for revision risk male <- ifelse(gender == "Female", 0, 1) np1 <- ifelse(strategy_name == "Standard prosthesis", 0, 1) scale <- exp(rr_coef$cons + rr_coef$age * age + rr_coef$male * male + rr_coef$np1 * np1) shape <- exp(rr_coef$lngamma) rr <- 1 - exp(scale * ((time - 1)^shape - time^shape)) # Mortality rate age_new <- age + time mr <- mr_female[["35"]] * (gender == "Female" & age_new >= 35 & age_new < 45) + mr_female[["45"]] * (gender == "Female" & age_new >= 45 & age_new < 55) + mr_female[["55"]] * (gender == "Female" & age_new >= 55 & age_new < 65) + mr_female[["65"]] * (gender == "Female" & age_new >= 65 & age_new < 75) + mr_female[["75"]] * (gender == "Female" & age_new >= 75 & age_new < 85) + mr_female[["85"]] * (gender == "Female" & age_new >= 85) + mr_male[["35"]] * (gender == "Male" & age_new >= 35 & age_new < 45) + mr_male[["45"]] * (gender == "Male" & age_new >= 45 & age_new < 55) + mr_male[["55"]] * (gender == "Male" & age_new >= 55 & age_new < 65) + mr_male[["65"]] * (gender == "Male" & age_new >= 65 & age_new < 75) + mr_male[["75"]] * (gender == "Male" & age_new >= 75 & age_new < 85) + mr_male[["85"]] * (gender == "Male" & age_new >= 85) list( tpmatrix = tpmatrix( 0, C, 0, 0, omrPTHR, 0, C, rr, 0, mr, 0, 0, 0, C, omrRTHR + mr, 0, 0, rrr, C, mr, 0, 0, 0, 0, 1) ) }, times = 1:60) statevals_def <- define_tparams({ c_prosthesis <- ifelse(strategy_name == "Standard prosthesis", c_Standard, c_NP1) list( utility = u, costs = list( prosthesis = c_prosthesis, medical = c_med ) ) })
The model is defined by specifying parameter estimates, a random number generation expression for the PSA, and the transformed parameter expressions.
mod_def <- define_model(tparams_def = list(transitions_def, statevals_def), rng_def = rng_def, params = params)
The economic model is created from the input data and the defined model. Arguments for instantiating the cost models (a list of objects of class StateVals
) can be specified with a named list containing the values of arguments to pass to StateVals$new()
. In this case, we note that the model for the cost of the prostheses is based on one-time costs accrued at the start of the model (with no discounting); conversely, medical costs accrue by weighting state values by time spent in the health states (in the current model simply the cost of revision over 1 year).
cost_args <- list( prosthesis = list(method = "starting"), medical = list(method = "wlos") ) econmod <- create_CohortDtstm(mod_def, data, cost_args = cost_args)
Health state probabilities are simulated for 60 model cycles.
econmod$sim_stateprobs(n_cycles = 60)
The expected number of patients per 1,000 in each health state is plotted below. The counts in the plot are for patient_id = 2
, a 60 year old female as in the original textbook example. Since a PSA was performed both mean values and 95 percent credible intervals can be displayed. Patients are more likely to have a revision with a standard prosthesis than with the new prosthesis.
Subgroup analyses can be performed with $sim_summarize()
using the by_grp = TRUE
option. Cost-effectiveness analyses are then performed for each subgroup (grp_id
) from input data. We will begin by reporting results for the 60-year old female, which are consistent with those from the textbook.
ce_sim <- econmod$summarize(by_grp = TRUE) wtp <- seq(0, 25000, 500) icea_pw_out <- icea_pw(ce_sim, comparator = 1, dr_qalys = 0.015, dr_costs = .06, k = wtp) ce_sim$qalys[grp_id == 2, .(mean = mean(qalys)), by = c("strategy_id", "grp_id")]
## strategy_id grp_id mean
## 1: 1 2 14.61403
## 2: 2 2 14.65681
ce_sim$costs[grp_id == 2 & category == "total", .(mean = mean(costs)), by = c("strategy_id", "grp_id")]
## strategy_id grp_id mean
## 1: 1 2 508.9896
## 2: 2 2 612.7832
We can then compute incremental cost-effectiveness ratios for each subgroup. There is considerable variation which is not unexpected since the the revision rate for prostheses (rr) depends on covariates.
icer_tbl(icea_pw_out, output = "data.table")
## strategy_id grp_id iqalys icosts inmb
## 1: 2 1 0.07 (0.02, 0.13) 72 (-30, 133) 3,386 (1,074, 6,664)
## 2: 2 2 0.04 (0.02, 0.08) 104 (28, 148) 2,035 (626, 4,114)
## 3: 2 3 0.03 (0.01, 0.05) 129 (74, 160) 1,149 (312, 2,388)
## 4: 2 4 0.01 (0.01, 0.03) 147 (109, 169) 590 (107, 1,347)
## 5: 2 5 0.01 (0.00, 0.02) 161 (135, 175) 255 (-19, 694)
## 6: 2 6 0.11 (0.04, 0.21) -21 (-213, 93) 5,541 (1,937, 10,292)
## 7: 2 7 0.07 (0.02, 0.12) 39 (-104, 120) 3,315 (1,132, 6,221)
## 8: 2 8 0.04 (0.01, 0.07) 85 (-16, 141) 1,883 (599, 3,651)
## 9: 2 9 0.02 (0.01, 0.04) 119 (51, 156) 1,003 (294, 2,038)
## 10: 2 10 0.01 (0.00, 0.02) 143 (98, 168) 484 (82, 1,091)
## icer conclusion
## 1: 1,042 Cost-effective
## 2: 2,426 Cost-effective
## 3: 5,028 Cost-effective
## 4: 9,979 Cost-effective
## 5: 19,325 Cost-effective
## 6: Dominates Cost-effective
## 7: 583 Cost-effective
## 8: 2,167 Cost-effective
## 9: 5,298 Cost-effective
## 10: 11,379 Cost-effective
Estimates aggregated across all patients can also be computed with $sim_summarize()
by using by_grp = FALSE
. In this case, overall costs and QALYs are a weighted average (using the column patient_wt
) of costs and QALYs for each patient in the input data.
ce_sim <- econmod$summarize(by_grp = FALSE) icea_pw_out <- icea_pw(ce_sim, comparator = 1, dr_qalys = .015, dr_costs = .06, k = wtp) icer_tbl(icea_pw_out)
## 1 2
## Incremental QALYs "-" "0.03 (0.01, 0.05)"
## Incremental costs "-" "121 (58, 157)"
## Incremental NMB "-" "1,268 (381, 2,580)"
## ICER "-" "4,341"
## Conclusion "-" "Cost-effective"
In general, the transition probability of an event given a model cycle is the probability a patient has an event during the cycle conditional on not having an event up until that cycle. That is, given a survivor function, \(S(t)\) and a cycle of length \(u\), the transition probability, \(tp(t_u)\) of an event is,
\[ \begin{align} tp(t_u) &= \frac{1 - S(t)}{S(t-u)} \\ &= \frac{1 - exp\{-H(t)\}}{exp\{-H(t-u)\}} \\ &= 1 - exp\{H(t-u) - H(t)\}, \end{align} \]
where \(H(t)\) is the cumulative hazard function.
The proportional hazard Weibull distribution with scale parameter \(\lambda\) and shape parameter, \(\gamma\) has the cumulative hazard function,
\[ H(t) = \lambda t^\gamma \] The transition probability of an event for the Weibull model is thus,
\[ \begin{align} tp(t_u) &= 1 - exp\{\lambda (t-u)^\gamma -\lambda t^\gamma\} \\ &= 1 - exp\{\lambda[(t-u)^\gamma - t^\gamma]\} \end{align} \]