library(heRomod) library(ggplot2)
This example is an implementation of the assessment of a new total hip replacement (THR) technology described in chapter 3.5 of Decision Modelling for Health Economic Evaluation. A more detailed report is available at this location. This reports goes a bit further in the analysis. For the sake of simplicity we will not reproduce exactly the analysis from the book. See vignette vignette("i-reproduction", "heRomod")
for an exact reproduction.
This model has 5 states:
omrPTHR
);omrRTHR
);Two transition probabilities are time-varying in this model:
Other-cause death probabilities (mortality rate mr
) for the United Kingdom is taken from WHO databases using the get_who_mr()
function. The variable sex
, taking values 0
and 1
, must be recoded in "FMLE"
and "MLE"
before being passed to this function.
standardRR
increases with time with the following formula (a Weibull distribution):$$ P_{revision} = 1 - \exp(\lambda \times ((t-1)^\gamma-t^\gamma)) $$
Where $t$ is the time since revision, $\gamma = 1.45367786$ and:
$$ \lambda = exp(cons + ageC \times age + maleC \times sex) $$
Where $age$ and $sex$ (female = 0, male = 1) are individual characteristics, $cons = -5.49094$, $ageC = -0.0367$ and $maleC = 0.768536$.
standardRR
is modified by the relative risk $rrNP1 = 0.260677$.$$ P_{revision} = 1 - \exp(\lambda \times rr \times NP1 \times ((t-1)^\gamma-t^\gamma)) $$
rrr
) probability is set to be constant at 0.04 per year.The key element to specify time-varying elements in heRomod
is through the use of the package-defined variables markov_cycle
and state_cycle
. See vignette vignette("b-time-dependency", "heRomod")
for more details.
In order to build this more complex Markov model, parameters need to be defined through define_parameters()
(for 2 reasons: to keep the transition matrix readable, and to avoid repetition by re-using parameters between strategies).
The equations decribed in the previous section can be written easily, here for a female population (sex
= 0) starting at 60 years old (age_init
= 60).
param <- define_parameters( age_init = 60, sex = 0, # age increases with cycles age = age_init + markov_cycle, # operative mortality rates omrPTHR = .02, omrRTHR = .02, # re-revision mortality rate rrr = .04, # parameters for calculating primary revision rate cons = -5.49094, ageC = -.0367, maleC = .768536, lambda = exp(cons + ageC * age_init + maleC * sex), gamma = 1.45367786, rrNP1 = .260677, # revision probability of primary procedure standardRR = 1 - exp(lambda * ((markov_cycle - 1) ^ gamma - markov_cycle ^ gamma)), np1RR = 1 - exp(lambda * rrNP1 * ((markov_cycle - 1) ^ gamma - markov_cycle ^ gamma)), # age-related mortality rate sex_cat = ifelse(sex == 0, "FMLE", "MLE"), mr = get_who_mr(age, sex_cat, country = "GBR", local = TRUE), # state values u_SuccessP = .85, u_RevisionTHR = .30, u_SuccessR = .75, c_RevisionTHR = 5294 ) param
Now that parameters are defined, the probability transitions can be easily written:
mat_standard <- define_transition( state_names = c( "PrimaryTHR", "SuccessP", "RevisionTHR", "SuccessR", "Death" ), 0, C, 0, 0, omrPTHR, 0, C, standardRR, 0, mr, 0, 0, 0, C, omrRTHR+mr, 0, 0, rrr, C, mr, 0, 0, 0, 0, 1 ) mat_standard mat_np1 <- define_transition( state_names = c( "PrimaryTHR", "SuccessP", "RevisionTHR", "SuccessR", "Death" ), 0, C, 0, 0, omrPTHR, 0, C, np1RR, 0, mr, 0, 0, 0, C, omrRTHR+mr, 0, 0, rrr, C, mr, 0, 0, 0, 0, 1 ) mat_np1
While it is possible to plot the matrix thanks to the diagram
package, the results may not always be easy to read.
plot(mat_standard)
Utilities and costs are then associated to states. In this model costs are discounted at a rate of 6% and utilities at a rate of 1.5%.
Now that parameters, transition matrix and states are defined we can define the strategies for the control group and the NP1 treatment.
We use define_starting_values()
to take into account the cost of surgery.
strat_standard <- define_strategy( transition = mat_standard, PrimaryTHR = define_state( utility = 0, cost = 0 ), SuccessP = define_state( utility = discount(u_SuccessP, .015), cost = 0 ), RevisionTHR = define_state( utility = discount(u_RevisionTHR, .015), cost = discount(c_RevisionTHR, .06) ), SuccessR = define_state( utility = discount(u_SuccessR, .015), cost = 0 ), Death = define_state( utility = 0, cost = 0 ), starting_values = define_starting_values( cost = 394 ) ) strat_standard strat_np1 <- define_strategy( transition = mat_np1, PrimaryTHR = define_state( utility = 0, cost = 0 ), SuccessP = define_state( utility = discount(u_SuccessP, .015), cost = 0 ), RevisionTHR = define_state( utility = discount(u_RevisionTHR, .015), cost = discount(c_RevisionTHR, .06) ), SuccessR = define_state( utility = discount(u_SuccessR, .015), cost = 0 ), Death = define_state( utility = 0, cost = 0 ), starting_values = define_starting_values( cost = 579 ) ) strat_np1
Both strategies can now be run for 60 years. By default models are computed for 1000 person starting in PrimaryTHR
.
res_mod <- run_model( standard = strat_standard, np1 = strat_np1, parameters = param, cycles = 60, cost = cost, effect = utility )
A comparison of both strategies can be done with summary()
. The incremental cost and effect are displayed in columns Cost
and Effect
.
summary(res_mod)
The new treatment costs £r round(summary(res_mod)$res_comp$.icer[2])
more per QALY gained.
It should be noted that this result differs from the original study. This difference is explained by higher population-level all-causes mortality rates in the original study than in the WHO database (used here). See vignette vignette("i-reproduction", "heRomod")
for an exact reproduction of the analysis.
We can plot the counts per state:
plot(res_mod, type = "counts", panel = "by_state", free_y = TRUE) + theme_bw() + scale_color_brewer( name = "Strategy", palette = "Set1" )
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.