Nothing
## ---- include = FALSE---------------------------------------------------------
ggplot2::theme_set(ggplot2::theme_bw())
## ---- out.width = "500px", echo = FALSE---------------------------------------
knitr::include_graphics("reversible-illness-death.png")
## -----------------------------------------------------------------------------
tmat <- rbind(c(NA, 1, 2),
c(3, NA, 4),
c(NA, NA, NA))
colnames(tmat) <- rownames(tmat) <- c("Healthy", "Sick", "Death")
print(tmat)
## ---- warning = FALSE, message = FALSE----------------------------------------
library("hesim")
library("data.table")
strategies <- data.table(strategy_id = c(1, 2),
strategy_name = c("SOC", "New"))
n_patients <- 1000
patients <- data.table(patient_id = 1:n_patients,
age = rnorm(n_patients, mean = 45, sd = 7),
female = rbinom(n_patients, size = 1, prob = .51))
states <- data.table(state_id = c(1, 2),
state_name = rownames(tmat)[1:2]) # Non-death health states
hesim_dat <- hesim_data(strategies = strategies,
patients = patients,
states = states)
## -----------------------------------------------------------------------------
labs <- get_labels(hesim_dat)
labs$transition_id <- c("Healthy-> Sick" = 1,
"Healthy -> Death" = 2,
"Sick -> Healthy" = 3,
"Sick -> Death" = 4)
print(labs)
## ---- warning = FALSE, message = FALSE----------------------------------------
library("flexsurv")
n_trans <- max(tmat, na.rm = TRUE) # Number of transitions
wei_fits_cr <- vector(length = n_trans, mode = "list")
for (i in 1:length(wei_fits_cr)){
wei_fits_cr[[i]] <- flexsurvreg(Surv(years, status) ~ factor(strategy_id),
data = mstate3_exdata$transitions,
subset = (trans == i) ,
dist = "weibull")
}
wei_fits_cr <- flexsurvreg_list(wei_fits_cr)
## -----------------------------------------------------------------------------
wei_fits_cf <- vector(length = n_trans, mode = "list")
for (i in 1:length(wei_fits_cf)){
wei_fits_cf[[i]] <- flexsurvreg(Surv(Tstart, Tstop, status) ~ factor(strategy_id),
data = mstate3_exdata$transitions,
subset = (trans == i) ,
dist = "weibull")
}
wei_fits_cf <- flexsurvreg_list(wei_fits_cf)
## -----------------------------------------------------------------------------
utility_tbl <- stateval_tbl(data.table(state_id = states$state_id,
mean = mstate3_exdata$utility$mean,
se = mstate3_exdata$utility$se),
dist = "beta")
head(utility_tbl)
## -----------------------------------------------------------------------------
drugcost_tbl <- stateval_tbl(data.table(strategy_id = strategies$strategy_id,
est = mstate3_exdata$costs$drugs$costs),
dist = "fixed")
medcost_tbl <- stateval_tbl(data.table(state_id = states$state_id,
mean = mstate3_exdata$costs$medical$mean,
se = mstate3_exdata$costs$medical$se),
dist = "gamma")
## -----------------------------------------------------------------------------
n_samples <- 1000
## -----------------------------------------------------------------------------
transmod_data <- expand(hesim_dat,
by = c("strategies", "patients"))
head(transmod_data)
## -----------------------------------------------------------------------------
transmod_cr <- create_IndivCtstmTrans(wei_fits_cr, transmod_data,
trans_mat = tmat, n = n_samples,
clock = "reset",
start_age = patients$age)
transmod_cf <- create_IndivCtstmTrans(wei_fits_cf, transmod_data,
trans_mat = tmat, n = n_samples,
clock = "forward",
start_age = patients$age)
## -----------------------------------------------------------------------------
# Predict hazard
transmod_data_pat1 <- transmod_data[patient_id == 1]
predict_haz <- function(fits, clock){
transmod_cr_pat1 <- create_IndivCtstmTrans(fits, transmod_data_pat1,
trans_mat = tmat,
clock = clock,
uncertainty = "none")
haz <- transmod_cr_pat1$hazard(t = seq(0, 20, 1))
title_clock <- paste(toupper(substr(clock, 1, 1)),
substr(clock, 2, nchar(clock)), sep="")
haz[, clock := title_clock]
return(haz[, ])
}
## ----predicted_hazards_plot, fig.width = 7, fig.height = 4--------------------
# Plot hazards
library("ggplot2")
haz <- rbind(predict_haz(wei_fits_cr, "reset"),
predict_haz(wei_fits_cf, "forward"))
set_labels(haz, labels = labs,
new_names = c("strategy_name", "trans_name"))
ggplot(haz[t > 0],
aes(x = t, y = hazard, col = clock, linetype = strategy_name)) +
geom_line() +
facet_wrap(~trans_name) +
xlab("Years") + ylab("Hazard") +
scale_linetype_discrete(name = "Strategy") +
scale_color_discrete(name = "Clock")
## -----------------------------------------------------------------------------
# Utility
utilitymod <- create_StateVals(utility_tbl, n = n_samples, hesim_data = hesim_dat)
# Costs
drugcostmod <- create_StateVals(drugcost_tbl, n = n_samples, hesim_data = hesim_dat)
medcostmod <- create_StateVals(medcost_tbl, n = n_samples, hesim_data = hesim_dat)
costmods <- list(Drug = drugcostmod,
Medical = medcostmod)
## -----------------------------------------------------------------------------
econmod_cr <- IndivCtstm$new(trans_model = transmod_cr,
utility_model = utilitymod,
cost_models = costmods)
econmod_cf <- IndivCtstm$new(trans_model = transmod_cf,
utility_model = utilitymod,
cost_models = costmods)
## ---- cache = FALSE-----------------------------------------------------------
# "Clock reset"
econmod_cr$sim_disease()
head(econmod_cr$disprog_)
# "Clock forward"
econmod_cf$sim_disease()
## -----------------------------------------------------------------------------
econmod_cr$sim_stateprobs(t = seq(0, 20 , 1/12))
## ----stprobs_by_strategy_plot, fig.width = 7, fig.height = 4------------------
# Short function to create state probability "dataset" for plotting
summarize_stprobs <- function(stateprobs){
x <- stateprobs[, .(prob_mean = mean(prob)),
by = c("strategy_id", "state_id", "t")]
set_labels(x, labels = labs, new_names = c("strategy_name", "state_name"))
}
# Plot of state probabilities
stprobs_cr <- summarize_stprobs(econmod_cr$stateprobs_)
ggplot(stprobs_cr, aes(x = t, y = prob_mean, col = strategy_name)) +
geom_line() + facet_wrap(~state_name) +
xlab("Years") + ylab("Probability in health state") +
scale_color_discrete(name = "Strategy") +
theme(legend.position = "bottom")
## ----stprobs_by_timescale_plot, fig.width = 7, fig.height = 4-----------------
econmod_cf$sim_stateprobs(t = seq(0, 20 , 1/12))
stprobs_cf <- summarize_stprobs(econmod_cf$stateprobs_)
# Compare "clock forward" and "clock reset" cases
stprobs <- rbind(data.table(stprobs_cf, clock = "Forward"),
data.table(stprobs_cr, clock = "Reset"))
ggplot(stprobs[strategy_id == 1],
aes(x = t, y = prob_mean, col = clock)) +
geom_line() + facet_wrap(~state_name) +
xlab("Years") + ylab("Probability in health state") +
scale_color_discrete(name = "Clock") +
theme(legend.position = "bottom")
## -----------------------------------------------------------------------------
econmod_cr$sim_qalys(dr = c(0,.03))
head(econmod_cr$qalys_)
## ----qalys_plot, fig.width = 6, fig.height = 4--------------------------------
qalys_summary <- econmod_cr$qalys_[, .(mean = mean(qalys)),
by = c("strategy_id", "state_id", "dr")]
set_labels(qalys_summary, labels = labs,
new_names = c("strategy_name", "state_name"))
ggplot(qalys_summary[dr == .03],
aes(x = strategy_name, y = mean, fill = state_name)) +
geom_bar(stat = "identity") +
scale_fill_discrete(name = "") +
xlab("Strategy") + ylab("Mean QALYs")
## -----------------------------------------------------------------------------
econmod_cr$sim_costs(dr = 0.03)
head(econmod_cr$costs_)
## -----------------------------------------------------------------------------
ce_sim <- econmod_cr$summarize()
format(summary(ce_sim, labels = labs))
cea_out <- cea(ce_sim, dr_qalys = .03, dr_costs = .03)
cea_pw_out <- cea_pw(ce_sim, comparator = 1, dr_qalys = .03, dr_costs = .03)
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.