Time-varying Markov Models (Non-Homogeneous)

library(heemod)
library(ggplot2)

Model description

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", "heemod") for an exact reproduction.

This model has 5 states:

Two transition probabilities are time-varying in this model:

Other-Cause death

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.

THR revision

$$ 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$.

$$ P_{revision} = 1 - \exp(\lambda \times rr \times NP1 \times ((t-1)^\gamma-t^\gamma)) $$

Parameter definition

The key element to specify time-varying elements in heemod is through the use of the package-defined variables model_time and state_time. See vignette vignette("b-time-dependency", "heemod") 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 + model_time,

    # 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 * ((model_time - 1) ^ gamma -
                                     model_time ^ gamma)),
    np1RR = 1 - exp(lambda * rrNP1 * ((model_time - 1) ^ gamma - 
                                        model_time ^ gamma)),

    # age-related mortality rate
    sex_cat = ifelse(sex == 0, "FMLE", "MLE"),
    mr = get_who_mr(age, sex_cat, local = TRUE),

    # state values
    u_SuccessP = .85,
    u_RevisionTHR = .30,
    u_SuccessR = .75,
    c_RevisionTHR = 5294
)
param

Transition matrix definition

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)

State & strategy definition

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

Model analysis

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", "heemod") 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"
  )


Try the heemod package in your browser

Any scripts or data that you put into this service are public.

heemod documentation built on July 26, 2023, 5:45 p.m.