library(heemod)
The purpose of this vignette is to present how to reproduce exactly the results from Decision Modelling for Health Economic Evaluation. While other vignettes such as vignette("c-homogeneous", "heemod")
, vignette("d-non-homogeneous", "heemod")
or vignette("e-probabilistic", "heemod")
are greatly inspired from this book, key elements differ (mostly for the sake of clarity) and thus results differ, sometimes significantly, from the book. Here we show how to exactly reproduce the results with the heemod
package.
Key differences in DMHEE:
It is possible to reproduce 1. and 2. by making transition happen at the end of each year with method = "end". Since with this method the transition occur 1 year after the beginning, costs should be discounted from the first cycle with the argument first = TRUE in discount().
Point 3. is reproduced by making rr
and cost_lami
a time changing variable like this rr = ifelse(model_time <= 2, .509, 1.00)
.
The last point is reproduced by writing transition probabilities as fractions.
par_mod <- define_parameters( rr = ifelse(model_time <= 2, .509, 1), cost_lami = ifelse(model_time <= 2, 2086.5, 0), cost_zido = 2278 ) mat_mono <- define_transition( 1251/1734, 350/1734, 116/1734, 17/1734, 0, 731/1258, 512/1258, 15/1258, 0, 0, 1312/1749, 437/1749, 0, 0, 0, 1.00 ) mat_comb <- define_transition( C, 350/1734*rr, 116/1734*rr, 17/1734*rr, 0, C, 512/1258*rr, 15/1258*rr, 0, 0, C, 437/1749*rr, 0, 0, 0, 1.00 ) A_mono <- define_state( cost_health = 2756, cost_drugs = cost_zido, cost_total = discount( cost_health + cost_drugs, .06, first = T), life_year = 1 ) B_mono <- define_state( cost_health = 3052, cost_drugs = cost_zido, cost_total = discount( cost_health + cost_drugs, .06, first = T), life_year = 1 ) C_mono <- define_state( cost_health = 9007, cost_drugs = cost_zido, cost_total = discount( cost_health + cost_drugs, .06, first = T), life_year = 1 ) D_mono <- define_state( cost_health = 0, cost_drugs = 0, cost_total = discount( cost_health + cost_drugs, .06, first = T), life_year = 0 ) A_comb <- define_state( cost_health = 2756, cost_drugs = cost_zido + cost_lami, cost_total = discount( cost_health + cost_drugs, .06, first = T), life_year = 1 ) B_comb <- define_state( cost_health = 3052, cost_drugs = cost_zido + cost_lami, cost_total = discount( cost_health + cost_drugs, .06, first = T), life_year = 1 ) C_comb <- define_state( cost_health = 9007, cost_drugs = cost_zido + cost_lami, cost_total = discount( cost_health + cost_drugs, .06, first = T), life_year = 1 ) D_comb <- define_state( cost_health = 0, cost_drugs = 0, cost_total = discount( cost_health + cost_drugs, .06, first = T), life_year = 0 ) mod_mono <- define_strategy( transition = mat_mono, A_mono, B_mono, C_mono, D_mono ) mod_comb <- define_strategy( transition = mat_comb, A_comb, B_comb, C_comb, D_comb ) res_mod <- run_model( mono = mod_mono, comb = mod_comb, parameters = par_mod, cycles = 20, cost = cost_total, effect = life_year, method = "end", init = c(1, 0, 0, 0) ) summary(res_mod)
Key difference in DMHEE:
This can be corrected by using a user-specified mortality table and then fetch the values with:
look_up(death_prob, age = age, sex = sex, bin = "age")
# a function to return age-related mortality rate # given age and sex death_prob <- data.frame( age = rep(seq(35, 85, 10), each = 2), sex = rep(1:0, 6), value = c( 1.51e-3, .99e-3, 3.93e-3, 2.6e-3, 10.9e-3, 6.7e-3, 31.6e-3, 19.3e-3, 80.1e-3, 53.5e-3, 187.9e-3, 154.8e-3 ) ) death_prob 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.490935, ageC = -.0367022, maleC = .768536, lambda = exp(cons + ageC * age_init + maleC * sex), lngamma = 0.3740968, gamma = exp(lngamma), lnrrNP1 = -1.344474, rrNP1 = exp(lnrrNP1), # 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 = look_up(death_prob, age = age, sex = sex, bin = "age"), u_successP = .85, u_revisionTHR = .30, u_successR = .75, c_revisionTHR = 5294 ) 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_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 ) mod_standard <- define_strategy( transition = mat_standard, PrimaryTHR = define_state( utility = 0, cost = 394 ), 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 ) ) mod_np1 <- define_strategy( transition = mat_np1, PrimaryTHR = define_state( utility = 0, cost = 579 ), 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 ) ) res_mod <- run_model( standard = mod_standard, np1 = mod_np1, parameters = param, cycles = 60, cost = cost, effect = utility, method = "beginning", init = c(1, 0, 0, 0, 0) ) summary(res_mod)
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.