M3_Markov_age_dep.RmdThis vignette replicates the age dependent part of the Advanced tutorial for Markov models 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.
This vignette builds on the previous introductory tutorial. If you haven’t already, please review the introductory Markov tutorial first.
The outcomes may be slightly different from the tutorial. The reason is that the tutorial applies a within cycle correction and also assumes that transition rewards happen after the transition. gmod links the transition rewards to the events so they occur during each cycle.
First, install gmod using devtools. Also, make sure to enable use_vignettes so that vignettes render correctly.
#library(devtools)
#install_github("hjalal/gmod", use_vignettes = TRUE, force = TRUE)
rm(list = ls())
library(gmod)
#> Please note that active development of gmod is now moved to twig. Please install twig https://hjalal.github.io/twig/ and review the vignettes for the latest features.
library(magrittr)Define the model structure. Note that the model is structured in a way that the event “die” happens first and if one survives they go to the rest of the evnets, whether recover, getsick progress or stay in their health state depending on where they were and the decision taken
mygmod <- gmod() +
decisions("StandardOfCare", "StrategyA", "StrategyB", "StrategyAB") +
states(names = c("H", "S1", "S2", "D"),
init_probs = c(1,0,0,0)) +
event(name = "die",
scenarios = c(T, F),
probs = c(pDie(state, cycle), Inf),
outcomes = c("D", "get_event")) +
event(name = "get_event",
scenarios = c("recover", "getsick", "progress", "stay"),
probs = c(pRecover(state),
pGetSick(state),
pProgress(state, decision),
Inf),
outcomes = c("H", "S1", "S2", "curr_state")) +
payoffs(cost = cost(state, decision, get_event, die),
utility = utility(state, decision, get_event),
discount_rates = c(0.03, 0.03))The gmod() function defines the model type as a Markov
model and number of cycles as 75. The decision() function
defines the treatment strategies. The initial_probs()
defines the initial probabilities. Here we only define the “Moderate”
states as the starting state with probability equals 1. The other states
will be set to 0.
Since the tutorial assumes that all the events happen simultaneously
(i.e., no sequence), we add a single event_mapping() layer
we model the individual final_outcomes separately. We can interpret the
event_mapping as the following: The event get_event can
have 5 values: recover, getsick, progress, die, or stay in the same
state. These values will outcome in the following states: H, S1, S2, D
and staying in the current state, respectively. The prob argument
defines the probability of these final_outcomes. For example, pRecover()
defins the probability of recovering and it is a function of state
because this probability is only applied if the state is “S1”, and is 0
otherwise. Note here we use Inf here is as a placeholder to
reduce the amount of typing required. Then, we define the payoffs of
cost and utility and define each as a function of the state and the
decision taken.
The "curr_state" is a build-in function that instructs
gmod to compute the probability of a state returning to itself.
Before we can evaluate the model, first we need to pass the parameter
values. Here we use the global environment to define those. Later, we
show how to pass these as a parameter list param which will
be helpful for analyses such as probabilistic sensitivity analyses.
n_age_init <- 25 - 1 # age at baseline - gmod starts from cycle 1 instead of 0 as in the tutorial.
n_age_max <- 100 # maximum age of follow up
### Transition rates (annual), and hazard ratios (HRs) ----
r_HS1 <- 0.15 # constant annual rate of becoming Sick when Healthy
r_S1H <- 0.5 # constant annual rate of becoming Healthy when Sick
r_S1S2 <- 0.105 # constant annual rate of becoming Sicker when Sick
hr_S1 <- 3 # hazard ratio of death in Sick vs Healthy
hr_S2 <- 10 # hazard ratio of death in Sicker vs Healthy
### Effectiveness of treatment B ----
hr_S1S2_trtB <- 0.6 # hazard ratio of becoming Sicker when Sick under treatment B
## Age-dependent mortality rates ----
lt_usa_2015 <- read.csv("../inst/extdata/LifeTable_USA_Mx_2015.csv")
#* Extract age-specific all-cause mortality for ages in model time horizon
v_r_mort_by_age <- lt_usa_2015 %>%
dplyr::filter(Age >= n_age_init & Age < n_age_max) %>%
dplyr::select(Total) %>%
as.matrix() # anyone above 100 have the same mortalityDefine the cost and utilities
### State rewards ----
#### Costs ----
c_H <- 2000 # annual cost of being Healthy
c_S1 <- 4000 # annual cost of being Sick
c_S2 <- 15000 # annual cost of being Sicker
c_D <- 0 # annual cost of being dead
c_trtA <- 12000 # annual cost of receiving treatment A
c_trtB <- 13000 # annual cost of receiving treatment B
#### Utilities ----
u_H <- 1 # annual utility of being Healthy
u_S1 <- 0.75 # annual utility of being Sick
u_S2 <- 0.5 # annual utility of being Sicker
u_D <- 0 # annual utility of being dead
u_trtA <- 0.95 # annual utility when receiving treatment A
### Transition rewards ----
du_HS1 <- 0.01 # disutility when the event get sick happens which is conditional on not dying.
ic_HS1 <- 1000 # increase in cost when the event get sick occurs conditional on not dying.
ic_D <- 2000 # increase in cost when the die event occurs.The model has six user-defined functions: four for computing
probabilities pRecover(), pGetSick(),
pProgress() and pDie(), and 2 for payoffs
cost() and utility().
probability of recover is only dependent on state being S1 otherwise is 0. Following the tutorial, we convert the rate of recover r_S1H to a probability by using the function rate2prob().
pRecover <- function(state){
rRecover <- if (state=="S1") r_S1H else 0
rate2prob(rRecover)
}This is a standard function and can be evaluated like any other function. For example, we can type
pRecover("S1")
#> [1] 0.3934693
pRecover("S2")
#> [1] 0probability of getting sick is also a function of a single state “H” otherwise is 0.
pGetSick <- function(state){
rGetSick <- if (state=="H") r_HS1 else 0
rate2prob(rGetSick)
}probability of progress is slightly more complex as it is a function of both the state and the decision.
pProgress <- function(state, decision){ #, cycle_in_S1){
hr_S1S2 <- if (decision %in% c("StrategyB","StrategyAB")) hr_S1S2_trtB else 1
rProgress <- (state=="S1") * r_S1S2 * hr_S1S2
rate2prob(rProgress)
}hr_S1S2 is a function of the decision if it contains “B” then it is hr_S1S2_trtB otherwise is 1. then the rate of progression is only conditional on state being S1 which gets multiplied by the r_S1S2 and the hr_S1S2.
Using expressions such as (state==“S1”) are very helpful as they can reduce the typing required and are easy to read once one gets using to them. This expression is shorthand of a full if…then statement:
state <- "S1"
# these two statements are equivalent and both return TRUE:
(state=="S1")
#> [1] TRUE
if (state=="S1") TRUE else FALSE
#> [1] TRUEProbability of death is a function of the state. For states S1 and S2 the rate of mortality is relative to the healthy rate multiplied by a hazard ratio for each state.
pDie <- function(state, cycle){
r_HD <- v_r_mort_by_age[cycle]
rDie <- switch(state,
"H" = r_HD,
"S1" = r_HD*hr_S1,
"S2" = r_HD*hr_S2,
"D" = Inf) #translates to probability of 1
rate2prob(rDie)
}The payoffs are a cost and utility function. Both are dependent on both the state and the decision. For example only those in S1 and S2 receive treatment.
cost <- function(state, decision, get_event, die){
# cost of decision is only applied if the state is either S1 or S2
trans_cost_getting_sick <- (get_event=="getsick")*ic_HS1 # increase in cost when transitioning from Healthy to Sick
trans_cost_dying <- (die==TRUE)*ic_D # increase in cost when dying
if (state %in% c("S1","S2")){
c_decision <- switch(decision,
"StandardOfCare" = 0,
"StrategyA" = c_trtA ,
"StrategyB" = c_trtB,
"StrategyAB" = c_trtA + c_trtB)}
else {
c_decision <- 0
}
# cost of the state is a function of the state
c_state <- switch(state,
"H" = c_H,
"S1" = c_S1,
"S2" = c_S2,
"D" = c_D)
# combine both
return(c_decision + c_state + trans_cost_getting_sick + trans_cost_dying)
}Similarly the utility function is a function of the state and decision. For H, S2 and D, the utility is independent of the decision. But, Strategies involving “A” improve the utility of S1 from u_S1 to u_trtA. One way to code this is:
utility <- function(state, decision, get_event){
trans_util_getting_sick <- -du_HS1*(get_event=="getsick")
u_state <- switch(state,
"H" = u_H,
"S1" = if (decision %in% c("StrategyA", "StrategyAB")) u_trtA else u_S1,
"S2" = u_S2,
"D" = u_D)
return(u_state + trans_util_getting_sick)
}The output is the structure of the model including the formulae using
gmod_gen_model_function() and then calling the generated
function my_markov_model().
To evaluate the model, we run gmod_evaluate() which returns
the numerical outcomes.
n_cycles <- 10
model_struc <- gmod_gen_model_function(mygmod)
#> Warning in gmod_gen_model_function.gmod_markov(mygmod): Model function my_markov_model is generated. It can be run by calling it directly, for example this function returns the summary results:
#> my_markov_model(model_struc,params,return_transition_prob=FALSE,return_state_payoffs=FALSE,return_trace=FALSE,return_cycle_payoffs=FALSE,return_payoff_summary=TRUE)
model_res <- my_markov_model(model_struc, return_trace = T)
#model_resThe warning just indicates that no parameters were passed, so gmod
uses the values in the global environment. For analyses, such as
probabilistic sensitivity analysis (PSA), we can pass the parameters as
a list param.
We can load individual outcomes here that we are interested in. For example, if we are interested in the Markov Trace of the Standard of Care treatment, we can type:
model_res$Trace
#> , , StandardOfCare
#>
#> H S1 S2 D
#> 1 1.0000000 0.0000000 0.00000000 0.000000000
#> 2 0.8598357 0.1391509 0.00000000 0.001013486
#> 3 0.7939162 0.1899669 0.01382842 0.002288444
#> 4 0.7571058 0.2064451 0.03255561 0.003893471
#> 5 0.7319044 0.2096413 0.05271416 0.005740082
#> 6 0.7114248 0.2077181 0.07292855 0.007928592
#> 7 0.6930599 0.2038953 0.09271446 0.010330419
#> 8 0.6757341 0.1993907 0.11184616 0.013029051
#> 9 0.6590119 0.1946753 0.13021137 0.016101511
#> 10 0.6427558 0.1899520 0.14781676 0.019475491
#>
#> , , StrategyA
#>
#> H S1 S2 D
#> 1 1.0000000 0.0000000 0.00000000 0.000000000
#> 2 0.8598357 0.1391509 0.00000000 0.001013486
#> 3 0.7939162 0.1899669 0.01382842 0.002288444
#> 4 0.7571058 0.2064451 0.03255561 0.003893471
#> 5 0.7319044 0.2096413 0.05271416 0.005740082
#> 6 0.7114248 0.2077181 0.07292855 0.007928592
#> 7 0.6930599 0.2038953 0.09271446 0.010330419
#> 8 0.6757341 0.1993907 0.11184616 0.013029051
#> 9 0.6590119 0.1946753 0.13021137 0.016101511
#> 10 0.6427558 0.1899520 0.14781676 0.019475491
#>
#> , , StrategyB
#>
#> H S1 S2 D
#> 1 1.0000000 0.0000000 0.000000000 0.000000000
#> 2 0.8598357 0.1391509 0.000000000 0.001013486
#> 3 0.7939162 0.1953247 0.008470643 0.002288444
#> 4 0.7592072 0.2166710 0.020268118 0.003853620
#> 5 0.7377216 0.2234403 0.033235108 0.005602964
#> 6 0.7218365 0.2240962 0.046446164 0.007621213
#> 7 0.7084326 0.2222402 0.059542002 0.009785170
#> 8 0.6961415 0.2193453 0.072344049 0.012169122
#> 9 0.6843753 0.2160275 0.084758295 0.014838971
#> 10 0.6729252 0.2125690 0.096774962 0.017730805
#>
#> , , StrategyAB
#>
#> H S1 S2 D
#> 1 1.0000000 0.0000000 0.000000000 0.000000000
#> 2 0.8598357 0.1391509 0.000000000 0.001013486
#> 3 0.7939162 0.1953247 0.008470643 0.002288444
#> 4 0.7592072 0.2166710 0.020268118 0.003853620
#> 5 0.7377216 0.2234403 0.033235108 0.005602964
#> 6 0.7218365 0.2240962 0.046446164 0.007621213
#> 7 0.7084326 0.2222402 0.059542002 0.009785170
#> 8 0.6961415 0.2193453 0.072344049 0.012169122
#> 9 0.6843753 0.2160275 0.084758295 0.014838971
#> 10 0.6729252 0.2125690 0.096774962 0.017730805Similarly, we can load the summary outcomes of the costs and utility of the strategies, such that:
model_res$summary_payoffs
#> cost utility
#> StandardOfCare 28525.55 8.067497
#> StrategyA 52980.50 8.367632
#> StrategyB 51650.73 8.140316
#> StrategyAB 75088.41 8.461766