Nothing
## ---- echo=FALSE, include=FALSE-----------------------------------------------
library(heemod)
## -----------------------------------------------------------------------------
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)
## -----------------------------------------------------------------------------
# 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.