## ----warning = FALSE, message = FALSE-----------------------------------------
library("hesim")
library("data.table")
library("kableExtra")
library("flexsurv")
library("ggplot2")
## ----out.width = "700px", echo = FALSE----------------------------------------
knitr::include_graphics("markov-inhomogeneous.png")
## ----warning = FALSE, message = FALSE-----------------------------------------
# Treatment strategies
strategies <- data.table(
strategy_id = 1:2,
strategy_name = c("Standard prosthesis", "New prosthesis")
)
n_strategies <- nrow(strategies)
# Patients
n_patients <- 1000
patients <- data.table(
patient_id = 1:n_patients,
gender = "Female",
age = 60
)
# States
states <- data.table( # Non-death health states
state_id = 1:4,
state_name = c("Primary THR", "Sucessful primary", "Revision THR",
"Successful revision")
)
n_states <- nrow(states)
# "hesim data"
hesim_dat <- hesim_data(strategies = strategies,
patients = patients,
states = states)
print(hesim_dat)
## -----------------------------------------------------------------------------
labs <- get_labels(hesim_dat)
print(labs)
## -----------------------------------------------------------------------------
tmat <- rbind(c(NA, 1, NA, NA, 2),
c(NA, NA, 3, NA, 4),
c(NA, NA, NA, 5, 6),
c(NA, NA, 7, NA, 8),
c(NA, NA, NA, NA, NA))
colnames(tmat) <- rownames(tmat) <- names(labs$state_id)
tmat %>%
kable() %>%
kable_styling()
## -----------------------------------------------------------------------------
ttrrPTHR <- 2 # Time to recovery rate implies mean time of 1/2 years
## -----------------------------------------------------------------------------
# 2 out of 100 patients receiving primary THR died
omrPTHR_shape1 <- 2
omrPTHR_shape2 <- 98
## -----------------------------------------------------------------------------
prob_to_rate <- function(p, t = 1){
(-log(1 - p))/t
}
## ----include = FALSE, results = "hide"----------------------------------------
rr_coef <- c(0.3740968, -5.490935, -0.0367022, 0.768536, -1.344474)
names(rr_coef) <- c("lngamma", "cons", "age", "male", "np1")
# Variance-covariance matrix
rr_vcov <- matrix(
c(0.0474501^2, -0.005691, 0.000000028, 0.0000051, 0.000259,
-0.005691, 0.207892^2, -0.000783, -0.007247, -0.000642,
0.000000028, -0.000783, 0.0052112^2, 0.000033, -0.000111,
0.0000051, -0.007247, 0.000033, 0.109066^2, 0.000184,
0.000259, -0.000642, -0.000111, 0.000184, 0.3825815^2),
ncol = 5, nrow = 5, byrow = TRUE
)
rownames(rr_vcov) <- colnames(rr_vcov) <- names(rr_coef)
## -----------------------------------------------------------------------------
print(rr_coef)
print(rr_vcov)
## ----echo = FALSE-------------------------------------------------------------
mort_tbl <- rbind(
c(35, 45, .00151, .00099),
c(45, 55, .00393, .0026),
c(55, 65, .0109, .0067),
c(65, 75, .0316, .0193),
c(75, 85, .0801, .0535),
c(85, Inf, .1879, .1548)
)
colnames(mort_tbl) <- c("age_lower", "age_upper", "male", "female")
mort_tbl <- data.frame(mort_tbl)
## -----------------------------------------------------------------------------
mr <- c(.0067, .0193, .0535, .1548)
mr_times <- c(0, 5, 15, 25)
## -----------------------------------------------------------------------------
ttrRTHR <- 1 # There is no rate, the time is fixed
## -----------------------------------------------------------------------------
# 4 out of 100 patients with a successful revision needed another procedure r
omrRTHR_shape1 <- 4
omrRTHR_shape2 <- 96
## -----------------------------------------------------------------------------
n_samples <- 500
## -----------------------------------------------------------------------------
vec_to_dt <- function(v, n = NULL){
if (length(v) == 1) v <- rep(v, n_samples)
dt <- data.table(v)
colnames(dt) <- "cons"
return(dt)
}
## -----------------------------------------------------------------------------
transmod_coef_def <- define_rng({
omrPTHR <- prob_to_rate(beta_rng(shape1 = omrPTHR_shape1,
shape2 = omrPTHR_shape2))
mr <- fixed(mr)
mr_omrPTHR <- omrPTHR + mr
colnames(mr) <- colnames(mr_omrPTHR) <- paste0("time_", mr_times)
rr <- multi_normal_rng(mu = rr_coef, Sigma = rr_vcov)
rrr <- prob_to_rate(beta_rng(shape1 = 4, shape2 = 96))
list(
log_omrPTHR = vec_to_dt(log(omrPTHR)),
log_mr = log(mr),
log_ttrrPTHR = vec_to_dt(log(ttrrPTHR)),
log_mr_omrPTHR = log(mr_omrPTHR),
rr_shape = vec_to_dt(rr$lngamma),
rr_scale = rr[, -1,],
log_rrr = vec_to_dt(log(rrr))
)
}, n = n_samples)
transmod_coef <- eval_rng(transmod_coef_def)
## -----------------------------------------------------------------------------
summary(transmod_coef)
## -----------------------------------------------------------------------------
head(transmod_coef$log_omrPTHR)
head(transmod_coef$rr_scale)
## -----------------------------------------------------------------------------
as_dt_list <- function(x) {
lapply(as.list(x), vec_to_dt)
}
## -----------------------------------------------------------------------------
transmod_params <- params_surv_list(
# 1. Primary THR:Successful primary (1:2)
params_surv(coefs = list(rate = transmod_coef$log_ttrrPTHR),
dist = "exp"),
# 2. Primary THR:Death (1:5)
params_surv(coefs = list(rate = transmod_coef$log_omrPTHR),
dist = "exp"),
# 3. Successful primary:Revision THR (2:3)
params_surv(coefs = list(shape = transmod_coef$rr_shape,
scale = transmod_coef$rr_scale),
dist = "weibullPH"),
# 4. Successful primary:Death (2:5)
params_surv(coefs = as_dt_list(transmod_coef$log_mr),
aux = list(time = c(0, 5, 15, 25)),
dist = "pwexp"),
# 5. Revision THR:Successful revision (3:4)
params_surv(coefs = list(est = vec_to_dt(ttrRTHR)),
dist = "fixed"),
# 6. Revision THR:Death (3:5)
params_surv(coefs = as_dt_list(transmod_coef$log_mr_omrPTHR),
aux = list(time = c(0, 5, 15, 25)),
dist = "pwexp"),
# 7. Successful revision:Revision THR (4:3)
params_surv(coefs = list(rate = transmod_coef$log_rrr),
dist = "exp"),
# 8. Successful revision:Death (4:5)
params_surv(coefs = as_dt_list(transmod_coef$log_mr),
aux = list(time = c(0, 5, 15, 25)),
dist = "pwexp")
)
## -----------------------------------------------------------------------------
utility_tbl <- stateval_tbl(
data.table(state_id = states$state_id,
mean = c(0, .85, .3, .75),
se = c(0, .03, .03, .04)),
dist = "beta"
)
head(utility_tbl)
## -----------------------------------------------------------------------------
drugcost_tbl <- stateval_tbl(
data.table(strategy_id = rep(strategies$strategy_id, each = n_states),
state_id = rep(states$state_id, times = n_strategies),
est = c(394, 0, 0, 0,
579, 0, 0, 0)),
dist = "fixed"
)
medcost_tbl <- stateval_tbl(
data.table(state_id = states$state_id,
mean = c(0, 0, 5294, 0),
se = c(0, 0, 1487, 0)),
dist = "gamma",
)
## -----------------------------------------------------------------------------
transmod_data <- expand(hesim_dat, by = c("strategies", "patients"))
head(transmod_data)
## -----------------------------------------------------------------------------
transmod_data[, cons := 1]
transmod_data[, male := ifelse(gender == "Male", 1, 0)]
transmod_data[, np1 := ifelse(strategy_name == "New prosthesis", 1, 0)]
## -----------------------------------------------------------------------------
# Transition model
transmod <- create_IndivCtstmTrans(transmod_params,
input_data = transmod_data,
trans_mat = tmat,
clock = "forward",
start_age = patients$age)
## -----------------------------------------------------------------------------
# Utility
utilitymod <- create_StateVals(utility_tbl, n = transmod_coef_def$n,
hesim_data = hesim_dat)
# Costs
drugcostmod <- create_StateVals(drugcost_tbl, n = transmod_coef_def$n,
method = "starting", hesim_data = hesim_dat)
medcostmod <- create_StateVals(medcost_tbl, n = transmod_coef_def$n,
hesim_data = hesim_dat)
costmods <- list(Drug = drugcostmod,
Medical = medcostmod)
## -----------------------------------------------------------------------------
econmod <- IndivCtstm$new(trans_model = transmod,
utility_model = utilitymod,
cost_models = costmods)
## -----------------------------------------------------------------------------
ptm <- proc.time()
econmod$sim_disease(max_t = 60, max_age = 120)
proc.time() - ptm
head(econmod$disprog_)
## -----------------------------------------------------------------------------
econmod$sim_stateprobs(t = 0:60)
## ----echo = FALSE, fig.width = 7, fig.height = 4------------------------------
mean_stateprobs <- econmod$stateprobs_[, .(prob_mean = mean(prob)),
by = c("strategy_id", "state_id", "t")]
mean_stateprobs[, state_name := factor(state_id,
levels = 1:nrow(tmat),
labels = colnames(tmat))]
ggplot(mean_stateprobs, aes(x = t, y = prob_mean, col = factor(strategy_id))) +
geom_line() + facet_wrap(~state_name, scales = "free_y") +
xlab("Years") + ylab("Probability in health state") +
scale_color_discrete(name = "Strategy") +
theme(legend.position = "bottom") +
theme_bw()
## -----------------------------------------------------------------------------
econmod$sim_qalys(dr = .015)
econmod$sim_costs(dr = .06)
## -----------------------------------------------------------------------------
ce_sim <- econmod$summarize()
summary(ce_sim, labels = labs) %>%
format()
## -----------------------------------------------------------------------------
cea_pw_out <- cea_pw(ce_sim, comparator = 1, dr_qalys = 0.015, dr_costs = .06,
k = seq(0, 25000, 500))
icer(cea_pw_out, labels = labs) %>%
format(digits_qalys = 3)
## ----eval = FALSE, echo = FALSE-----------------------------------------------
# # Rather than use a piecewise exponential distribution, approximate
# # mortality rates with a Weibull distribution
# fit_from_pwexp <- function(n = 10000, rate, time = mr_times){
# sim_pwexp <- rpwexp(10000, rate = rate, time = time)
#
# fit_wei <- flexsurvreg(formula = Surv(sim_pwexp) ~ 1, dist = "weibull")
# sim_wei <- rweibull(n,
# shape = exp(fit_wei$res.t["shape", "est"]),
# scale = exp(fit_wei$res.t["scale", "est"]))
#
# res <- fit_wei
# attr(res, "sim_dist") <- list(wei = summary(sim_wei),
# pwexp = summary(sim_pwexp))
# return(res)
# }
# fit_mort_wei <- fit_from_pwexp(rate = mr)
# fit_mort2_wei <- fit_from_pwexp(rate = mr + .02) # .02 = mean of omrPTHR
#
# # Sample the parameters
# transmod_coef2 <- define_rng({
# mr <- multi_normal_rng(mu = fit_mort_wei$res.t[, "est"],
# Sigma = vcov(fit_mort_wei))
# mr_omrPTHR <- multi_normal_rng(mu = fit_mort2_wei$res.t[, "est"],
# Sigma = vcov(fit_mort2_wei))
# list(
# mr_shape = vec_to_dt(mr$shape),
# mr_scale = vec_to_dt(mr$scale),
# mr_omrPTHR_shape = vec_to_dt(mr_omrPTHR$shape),
# mr_omrPTHR_scale = vec_to_dt(mr_omrPTHR$scale)
# )
# }, n = n_samples)
# transmod_coef2 <- eval_rng(transmod_coef2)
#
# # Set parameters of transitions
# transition4_params <- params_surv(
# coefs = list(shape = transmod_coef2$mr_shape,
# scale = transmod_coef2$mr_scale),
# dist = "weibull")
#
# transition6_params <- params_surv(
# coefs = list(shape = transmod_coef2$mr_omrPTHR_shape,
# scale = transmod_coef2$mr_omrPTHR_scale),
# dist = "weibull")
#
# # Simulation
# ## Construct
# econmod2 <- econmod$clone(deep = TRUE)
# econmod2$trans_model$params[[4]] <- transition4_params
# econmod2$trans_model$params[[6]] <- transition6_params
# econmod2$trans_model$params[[8]]<- transition4_params
#
# ## Simulate
# econmod2$sim_disease(max_t = 60, max_age = 120)
# econmod2$sim_qalys(dr = .015)
# econmod2$sim_costs(dr = .06)
# ce_sim2 <- econmod2$summarize()
# ce_sim2$qalys[, .(mean = mean(qalys)),
# by = c("strategy_id")]
# ce_sim2$costs[category == "total",
# .(mean = mean(costs)),
# by = c("strategy_id")]