--- title: "A Tutorial for Time Dependent Markov Model" output: rmarkdown::html_vignette vignette: > %\VignetteIndexEntry{A Tutorial for Time Dependent Markov Model} %\VignetteEngine{knitr::rmarkdown} %\VignetteEncoding{UTF-8} --- ```{r, include = FALSE} knitr::opts_chunk$set( collapse = TRUE, comment = "#>" ) ``` This vignette is based on our tutorial for time-dependent Markov models in R published in *Medical Decision Making*: *Alarid-Escudero, F., Krijkamp, E., Enns, E. A., Yang, A., Hunink, M. M., Pechlivanoglou, P., & Jalal, H. (2023). A tutorial on time-dependent cohort state-transition models in r using a cost-effectiveness analysis example. Medical Decision Making, 43(1), 21-41.* First, we load `twig` library ```{r setup} library(twig) ``` # Twig ## model definition Then, we define the model's structure `twig` ```{r} n_cycles <- 25 # number of cycles mytwig <- twig() + decisions(names = c(StandardOfCare, StrategyA, StrategyB, StrategyAB)) + # define decisions states(names = c(H, S1, S2, D), # Markov state names init_probs = c(1,0,0,0), # everyone starts at H max_cycles = c(1,n_cycles, 1, 1)) + # the cohort can stay in S1 for n_cycles event(name = death_event, # first event is death options = c(yes,none), # which 2 options probs = c(pDie, leftover), # probability function name and its complement transitions = c(D, second_event)) + # if death occurs go to D, otherwise, go to the next event (second_event) event(name = second_event, # the second event options = c(recover, getsick, progress, none), # has 4 options probs = c(pRecover, pGetSick, pProgress, leftover), # and 3 named probabilities and a complement transitions = c(H, S1, S2, stay)) + # resulting in transitions to H, S1, S2 or else staying in the original state payoffs(names = c(cost, utility), # payoff names discount_rates = c(0.03, 0.03)) # payoff discount rates ``` ## DecisionTwig In [DecisionTwig](https://www.dashlab.ca/projects/decision_twig/), this Markov model looks like this: ![Markov](../man/figures/markov_twig.png){width=700px} Next, we create a data frame of random samples from the model input parameters' distributions ``` {r} n_sims <- 1000 # Create the data.table with n_sim rows of random samples params <- data.frame( r_HS1 = rbeta(n_sims, 2, 10), # Transition rate with beta distribution r_S1H = rbeta(n_sims, 5, 5), # Another transition rate with a different shape hr_S1 = rlnorm(n_sims, log(3), 0.2), # Hazard ratio, log-normal to allow skewness hr_S2 = rlnorm(n_sims, log(10), 0.2), # Higher hazard ratio, same distribution hr_S1S2_trtB = rbeta(n_sims, 6, 4), # Hazard ratio under treatment with beta distribution r_S1S2_scale = rgamma(n_sims, shape = 2, rate = 25), # Scale parameter, gamma distribution r_S1S2_shape = rgamma(n_sims, shape = 3, rate = 3), # Shape parameter, gamma distribution c_H = rnorm(n_sims, mean = 2000, sd = 50), # Annual cost, slight variation for simulation c_S1 = rnorm(n_sims, mean = 4000, sd = 100), # Higher annual cost, slightly varied c_S2 = rnorm(n_sims, mean = 15000, sd = 500), # Large cost with moderate variation c_D = 0, # Constant, no variation c_trtA = rnorm(n_sims, mean = 12000, sd = 200), # Cost of treatment A with small variation c_trtB = rnorm(n_sims, mean = 13000, sd = 200), # Cost of treatment B u_H = rbeta(n_sims, 10, 1), # Utility close to 1 for Healthy u_S1 = rbeta(n_sims, 7.5, 2.5), # Utility less than Healthy, beta distribution u_S2 = rbeta(n_sims, 5, 5), # Utility for Sicker u_D = 0, # Utility for Dead is constant u_trtA = rbeta(n_sims, 9.5, 1), # Utility with treatment A, close to Healthy du_HS1 = rnorm(n_sims, mean = 0.01, sd = 0.005), # Disutility with slight variation ic_HS1 = rnorm(n_sims, mean = 1000, sd = 100), # Cost increase with transition ic_D = rnorm(n_sims, mean = 2000, sd = 100), # Cost increase when dying p0_H = rbeta(n_sims, 1, 9) # Initial probability of being Healthy ) ``` # Probability and Payoff Functions Then, we define the probability and payoff functions used in the `twig` above: ## Probability of recovery (pRecover) Below, we define the probability of recovery `pRecover` which returns the probability of recovery based on the state and the recovery rate `r_S1H`. However, it is important to first explain logical conditions in R which we will use as shorthand conditions replacing `ifelse` statements in the vectorized function definitions that we use in `twig`. ### Using Logical Conditions in R In R, a common way to assign values based on conditions is to use logical conditions. For example, the following function `pRecover` assigns a recovery rate `r_S1H` if the state is "S1" and 0 otherwise using a logical condition: ``` {r} pRecover <- function(state, r_S1H){ rRecover <- r_S1H * (state=="S1") # Assigns r_S1H if in S1, otherwise 0 rate2prob(rRecover) } ``` The shorthand notation `(state == "S1")` is a concise and efficient way to create a logical vector that evaluates to 1 (`TRUE`) when state is "S1" and 0 (`FALSE`) otherwise. This is particularly useful and easier to read than using ifelse statements which can become difficult to read especially when mutliple conditions are nested. In this example an equivalent function using `ifelse` statement will look like this ``` {r} pRecover_ifelse <- function(state, r_S1H){ rRecover <- ifelse(state == "S1", r_S1H, 0) # Assigns r_S1H if in S1, otherwise 0 rate2prob(rRecover) } ``` Although, both of these functions are equivalent in `twig`, the shorthand notation using conditional logic is more concise and easier to read and we will use it throughout this example and the rest of the vignettes. ## Probability of getting sick (pGetSick) Only those who are H can get sick ``` {r} pGetSick <- function(state, r_HS1){ rGetSick <- r_HS1 * (state=="H") rate2prob(rGetSick) } ``` ## Probability of progressing (pProgress) This depends on the state (only those who are in S1 can progress), the decision (different rates for different decisions), and the cycle_in_state (number of cycles spent in a state - tunnel state) ``` {r} pProgress <- function(state, decision, cycle_in_state, hr_S1S2_trtB, r_S1S2_scale, r_S1S2_shape){ hr_S1S2 <- hr_S1S2_trtB ^ (decision %in% c("StrategyB", "StrategyAB")) # hazard rate of progression for B or 1 otherwise r_S1S2_tunnels <- ((cycle_in_state*r_S1S2_scale)^r_S1S2_shape - ((cycle_in_state - 1)*r_S1S2_scale)^r_S1S2_shape) # hazard rate based on cycle_in_state (tunnel) which follows a weibull distribution # only those who are at S1 can progress rProgress <- r_S1S2_tunnels * (state=="S1") * hr_S1S2 rate2prob(rProgress) } ``` ## Probability of dying (pDie) Probabilty of dying depends on age. So, we define age-specific mortality data from the 2015 US Life Tables for ages 24 through 1000. Then, we define the `pDie` as a function of the state (different states have different rates of death), and the cycles since the simulation start which reflects the cohort's age. ``` {r} v_r_mort_by_age <- vector <- c( 0.000979, 0.001014, 0.000999, 0.001070, 0.001087, 0.001162, 0.001167, 0.001213, 0.001289, 0.001331, 0.001375, 0.001420, 0.001490, 0.001550, 0.001616, 0.001657, 0.001747, 0.001902, 0.002052, 0.002173, 0.002395, 0.002559, 0.002807, 0.003023, 0.003349, 0.003712, 0.004085, 0.004490, 0.004905, 0.005364, 0.005806, 0.006253, 0.006775, 0.007395, 0.007895, 0.008418, 0.008974, 0.009666, 0.010456, 0.011384, 0.011838, 0.012667, 0.013593, 0.014700, 0.015732, 0.017340, 0.018758, 0.020967, 0.022917, 0.024913, 0.026767, 0.029707, 0.032412, 0.035982, 0.039238, 0.043595, 0.048727, 0.053735, 0.059911, 0.066618, 0.074051, 0.082190, 0.090754, 0.103968, 0.115093, 0.124341, 0.137872, 0.154177, 0.172393, 0.194100, 0.212654, 0.243752, 0.259087, 0.287781, 0.316429, 0.339149 ) # death depends on the state and age. pDie <- function(state, cycle, hr_S1, hr_S2){ r_HD <- v_r_mort_by_age[cycle] # get age-specific mortality rDie <- r_HD * (state=="H") + # baseline mortality if healthy r_HD*hr_S1 * (state=="S1") + # multiplied by a hazard rate if S1 or r_HD*hr_S2 * (state=="S2") # S2 # else 0 rate2prob(rDie) } ``` ## Cost Cost is a function of the state, decision and whether either event have occured. The events capture transition costs. ``` {r} cost <- function(state, decision, second_event, death_event, ic_HS1, ic_D, c_trtA, c_trtB, c_H, c_S1, c_S2, c_D){ # cost of decision is only applied if the state is either S1 or S2 trans_cost_getting_sick <- ic_HS1 * (second_event=="getsick") # increase in cost when transitioning from Healthy to Sick trans_cost_dying <- ic_D * (death_event=="yes") # increase in cost when dying c_decision <- (state %in% c("S1","S2")) * ( c_trtA * (decision=="StrategyA") + c_trtB * (decision=="StrategyB") + (c_trtA + c_trtB) * (decision=="StrategyAB") ) # cost of the state is a function of the state c_state <- c_H * (state=="H") + c_S1 * (state=="S1") + c_S2 * (state=="S2") + c_D * (state=="D") # combine both return(c_decision + c_state + trans_cost_getting_sick + trans_cost_dying) } ``` ## Utility Similarly, utility depends on the state, decision and only if the cohort gotsick to apply a transition utility discount for those who make that transition. ``` {r} utility <- function(state, decision, second_event, du_HS1, u_H, u_trtA, u_S1, u_S2, u_D){ trans_util_getting_sick <- -du_HS1 * (second_event=="getsick") # apply a utility discount for those who get sick. # calcualte state utilities. note that S1 will have utility u_trtA if the decision involves A, and another utility if the decision does not involve A. u_state <- u_H * (state=="H") + u_trtA * (state=="S1" & decision %in% c("StrategyA", "StrategyAB")) + u_S1 * (state=="S1" & decision %out% c("StrategyA", "StrategyAB")) + u_S2 * (state=="S2") + u_D * (state=="D") # combine the two utilities. return(u_state + trans_util_getting_sick) } ``` # Running the `twig` We can run the model and check the summary results under `results$mean_ev` of average costs and utilities across all probabilistic analyses. The results also returns the individual simulation results under `results$sim_ev` ``` {r} results <- run_twig(twig_obj = mytwig, params = params, n_cycles = n_cycles, progress_bar = FALSE) results$mean_ev ``` ## Parallelization You can run the same command with parallelization, by setting `parallel = TRUE`. This will speed up the simulation by running each simulation in parallel. [It is commented out here because testing of the `twig` package only runs on a single core.] ``` {r} # results <- run_twig(twig_obj = mytwig, params = params, n_cycles = n_cycles, parallel = TRUE) ``` ## Incremental Cost-Effectiveness Ratio (ICER) using the `calculate_icers` function adapted from the dampack package, we can retrieve ICER table by passing the payoffs summary table. ``` {R} calculate_icers(results$mean_ev) ``` ## Cost-Effectiveness Acceptibility Curve (CEAC) Similarly, we can produce the CEAC curve from the simulation results and specifying a series of willingness to pay (WTP) values. ```{r fig.width=7, fig.height=5} plot_ceac(results$sim_ev, wtp_range = seq(0, 100000, by = 1000)) ``` ## Summary This example illustrated the following features of `twig` * Markov model * probabilistic anslysis * Multiple decisions * Multiple states * Multiple sequential events * tunnel states * cycle-dependency * dependency on prior events in a cycle * transition payoffs * using parallel computations * expected payoffs across probabilistic analyses * ICER * CEAC In addition, various intermediate objects and computations can be returned by enabling `verbose`. To reduce the size of the returned results, this will only use the first PSA sample. For details on these objects check the documentations.