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(markov_cycle <= 2, .509, 1.00)
.
The last point is reproduced by writing transition probabilities as fractions.
<- define_parameters(
par_mod rr = ifelse(markov_cycle <= 2, .509, 1),
cost_lami = ifelse(markov_cycle <= 2, 2086.5, 0),
cost_zido = 2278
)
<- define_transition(
mat_mono 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
)
## No named state -> generating names.
<- define_transition(
mat_comb 350/1734*rr, 116/1734*rr, 17/1734*rr,
C, 0, C, 512/1258*rr, 15/1258*rr,
0, 0, C, 437/1749*rr,
0, 0, 0, 1.00
)
## No named state -> generating names.
<- define_state(
A_mono cost_health = 2756,
cost_drugs = cost_zido,
cost_total = discount(
+ cost_drugs, .06, first = T),
cost_health life_year = 1
)<- define_state(
B_mono cost_health = 3052,
cost_drugs = cost_zido,
cost_total = discount(
+ cost_drugs, .06, first = T),
cost_health life_year = 1
)<- define_state(
C_mono cost_health = 9007,
cost_drugs = cost_zido,
cost_total = discount(
+ cost_drugs, .06, first = T),
cost_health life_year = 1
)<- define_state(
D_mono cost_health = 0,
cost_drugs = 0,
cost_total = discount(
+ cost_drugs, .06, first = T),
cost_health life_year = 0
)
<- define_state(
A_comb cost_health = 2756,
cost_drugs = cost_zido + cost_lami,
cost_total = discount(
+ cost_drugs, .06, first = T),
cost_health life_year = 1
)<- define_state(
B_comb cost_health = 3052,
cost_drugs = cost_zido + cost_lami,
cost_total = discount(
+ cost_drugs, .06, first = T),
cost_health life_year = 1
)<- define_state(
C_comb cost_health = 9007,
cost_drugs = cost_zido + cost_lami,
cost_total = discount(
+ cost_drugs, .06, first = T),
cost_health life_year = 1
)<- define_state(
D_comb cost_health = 0,
cost_drugs = 0,
cost_total = discount(
+ cost_drugs, .06, first = T),
cost_health life_year = 0
)
<- define_strategy(
mod_mono transition = mat_mono,
A_mono,
B_mono,
C_mono,
D_mono )
## No named state -> generating names.
<- define_strategy(
mod_comb transition = mat_comb,
A_comb,
B_comb,
C_comb,
D_comb )
## No named state -> generating names.
<- run_model(
res_mod 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)
## 2 strategies run for 20 cycles.
##
## Initial state counts:
##
## A = 1
## B = 0
## C = 0
## D = 0
##
## Counting method: 'end'.
##
## Values:
##
## cost_health cost_drugs cost_total life_year
## mono 45541.24 18203.97 44663.45 7.991207
## comb 48082.83 24492.28 50601.65 8.937389
##
## Efficiency frontier:
##
## mono -> comb
##
## Differences:
##
## Cost Diff. Effect Diff. ICER Ref.
## comb 5938.198 0.9461822 6275.956 mono
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
<- data.frame(
death_prob 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
## age sex value
## 1 35 1 0.00151
## 2 35 0 0.00099
## 3 45 1 0.00393
## 4 45 0 0.00260
## 5 55 1 0.01090
## 6 55 0 0.00670
## 7 65 1 0.03160
## 8 65 0 0.01930
## 9 75 1 0.08010
## 10 75 0 0.05350
## 11 85 1 0.18790
## 12 85 0 0.15480
<- define_parameters(
param 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.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 * ((markov_cycle - 1) ^ gamma -
^ gamma)),
markov_cycle np1RR = 1 - exp(lambda * rrNP1 * ((markov_cycle - 1) ^ gamma -
^ gamma)),
markov_cycle
# 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
)
<- define_transition(
mat_standard 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
)
<- define_transition(
mat_np1 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
)
<- define_strategy(
mod_standard 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
)
)
<- define_strategy(
mod_np1 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
)
)
<- run_model(
res_mod 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)
## 2 strategies run for 60 cycles.
##
## Initial state counts:
##
## PrimaryTHR = 1
## SuccessP = 0
## RevisionTHR = 0
## SuccessR = 0
## Death = 0
##
## Counting method: 'beginning'.
##
## Values:
##
## utility cost
## standard 14.65283 512.4327
## np1 14.69734 610.3112
##
## Efficiency frontier:
##
## standard -> np1
##
## Differences:
##
## Cost Diff. Effect Diff. ICER Ref.
## np1 97.87858 0.04451095 2198.977 standard