M2_gmod_advanced_markov_tutorial.RmdThis vignette incorporates both cycle (age) dependency and state residence dependencies replicating our 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)NEW Syntax using layers = twigs
n_cycles <- 5
mygmod <- gmod() + # for illustration it is 75 in the tutorial
decisions("StandardOfCare", "StrategyA", "StrategyB", "StrategyAB") +
states(names=c("H", "S1", "S2", "D"),
init_probs=c(1,0,0,0),
max_cycle_in_states=c(1,n_cycles,1,1)) +
event(name = "die",
scenarios = c(T, F),
probs = c(pDie(state, cycle, cycle_in_state), 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, cycle_in_state),
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,0))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 cycle_in_state is a special keyward that consists of
the prefix cycle_in_ followed by one of the state names.
This keyword tells gmod which state to expand by adding additional
tunnel states. You can create other tunnels by adding
cycle_in_[State]. Note, adding tunnel states expands the
state space and can be computationally expensive.
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 at cycle 1 instead of 0 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
#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
#
##* Weibull parameters for state-residence-dependent transition probability of
##* becoming Sicker when Sick conditional on surviving
#r_S1S2_scale <- 0.08 # scale
#r_S1S2_shape <- 1.1 # shape
#
## 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 mortality
n_age_init <- 25 - 1 # age at baseline gmod starts at cycle 1 instead of 0 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
#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
#
##* Weibull parameters for state-residence-dependent transition probability of
##* becoming Sicker when Sick conditional on surviving
#r_S1S2_scale <- 0.08 # scale
#r_S1S2_shape <- 1.1 # shape
#
### Age-dependent mortality rates ----
#lt_usa_2015 <- read.csv("../data/LifeTable_USA_Mx_2015.csv")
#* Extract age-specific all-cause mortality for ages in model time horizon
v_r_mort_by_age0 <- lt_usa_2015 %>%
dplyr::filter(Age >= n_age_init & Age < n_age_max) %>%
dplyr::select(Total) %>%
as.matrix() # anyone above 100 have the same mortality
### 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 transitioning from Healthy to Sick
#ic_HS1 <- 1000 # increase in cost when transitioning from Healthy to Sick
#ic_D <- 2000 # increase in cost when dying
params <- list(
### 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
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
#* Weibull parameters for state-residence-dependent transition probability of
#* becoming Sicker when Sick conditional on surviving
r_S1S2_scale = 0.08, # scale
r_S1S2_shape = 1.1 , # shape
v_r_mort_by_age = v_r_mort_by_age0,
### 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 transitioning from Healthy to Sick
ic_HS1 = 1000, # increase in cost when transitioning from Healthy to Sick
ic_D = 2000 # increase in cost when dying
)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")
#
#pRecover("S2")probability 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_state){
if (state=="S1"){
hr_S1S2 <- if (decision %in% c("StrategyB", "StrategyAB")) hr_S1S2_trtB else 1
r_S1S2_tunnels <- (cycle_in_state*r_S1S2_scale)^r_S1S2_shape -
((cycle_in_state - 1)*r_S1S2_scale)^r_S1S2_shape
#p_S1S2_tunnels <- rate2prob(r_S1S2_tunnels)
#rProgress <- (state=="S1") * r_S1S2_tunnels * hr_S1S2
rProgress <- r_S1S2_tunnels * hr_S1S2
rate2prob(rProgress)
} else {
0
}
}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, cycle_in_state){
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)
}Then, we build the gmod Markov model by using the
gmod_build() function.
The output is the structure of the model including the formulae.
To evaluate the model, we run gmod_evaluate() which
returns the numerical outcomes.
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_struc
#> $is_cycle_dep
#> [1] TRUE
#>
#> $decisions
#> [1] "StandardOfCare" "StrategyA" "StrategyB" "StrategyAB"
#>
#> $n_decisions
#> [1] 4
#>
#> $tunnel_states
#> [1] "S1"
#>
#> $tunnel_lengths
#> S1
#> 5
#>
#> $states
#> [1] "H" "S1" "S2" "D"
#>
#> $n_states
#> [1] 4
#>
#> $states_expanded
#> [1] "H" "S1_tnl1" "S1_tnl2" "S1_tnl3" "S1_tnl4" "S1_tnl5" "S2"
#> [8] "D"
#>
#> $n_states_expanded
#> [1] 8
#>
#> $n_events
#> [1] 2
#>
#> $p0
#> $p0$StandardOfCare
#> H S1_tnl1 S1_tnl2 S1_tnl3 S1_tnl4 S1_tnl5 S2 D
#> 1 0 0 0 0 0 0 0
#>
#> $p0$StrategyA
#> H S1_tnl1 S1_tnl2 S1_tnl3 S1_tnl4 S1_tnl5 S2 D
#> 1 0 0 0 0 0 0 0
#>
#> $p0$StrategyB
#> H S1_tnl1 S1_tnl2 S1_tnl3 S1_tnl4 S1_tnl5 S2 D
#> 1 0 0 0 0 0 0 0
#>
#> $p0$StrategyAB
#> H S1_tnl1 S1_tnl2 S1_tnl3 S1_tnl4 S1_tnl5 S2 D
#> 1 0 0 0 0 0 0 0
#>
#>
#> $events
#> [1] "die" "get_event"
#>
#> $events_df
#> type event values
#> 1 event die TRUE
#> 2 event die FALSE
#> 3 event get_event recover
#> 4 event get_event getsick
#> 5 event get_event progress
#> 6 event get_event stay
#> probs
#> 1 pDie(state,cycle,cycle_in_state)
#> 2 1-(pDie(state,cycle,cycle_in_state))
#> 3 pRecover(state)
#> 4 pGetSick(state)
#> 5 pProgress(state,decision,cycle_in_state)
#> 6 1-(pRecover(state)+pGetSick(state)+pProgress(state,decision,cycle_in_state))
#> outcomes id
#> 1 D 1
#> 2 get_event 2
#> 3 H 3
#> 4 S1 4
#> 5 S2 5
#> 6 curr_state 6
#>
#> $payoffs
#> $payoffs$cost
#> cost(state, decision, get_event, die)
#>
#> $payoffs$utility
#> utility(state, decision, get_event)
#>
#>
#> $payoff_names
#> [1] "cost" "utility"
#>
#> $n_payoffs
#> [1] 2
#>
#> $discounts
#> c(cost = 0, utility = 0)
#>
#> $model_equations
#> # A tibble: 5 × 4
#> dest probs cost utility
#> <chr> <chr> <chr> <chr>
#> 1 D (pDie(state,cycle,cycle_in_state)) (pDi… (pDie(…
#> 2 H (1-(pDie(state,cycle,cycle_in_state)))*(pRecover(sta… (1-(… (1-(pD…
#> 3 S1 (1-(pDie(state,cycle,cycle_in_state)))*(pGetSick(sta… (1-(… (1-(pD…
#> 4 S2 (1-(pDie(state,cycle,cycle_in_state)))*(pProgress(st… (1-(… (1-(pD…
#> 5 curr_state (1-(pDie(state,cycle,cycle_in_state)))*(1-(pRecover(… (1-(… (1-(pD…
#>
#> attr(,"class")
#> [1] "gmod_markov"
#source("../R/myfun.R")
model_results <- my_markov_model(model_struc, params, return_trace = T)
dim(model_results$Trace)
#> [1] 5 8 4
model_results
#> $Trace
#> , , StandardOfCare
#>
#> H S1_tnl1 S1_tnl2 S1_tnl3 S1_tnl4 S1_tnl5 S2
#> 1 1.0000000 0.0000000 0.00000000 0.00000000 0.00000000 0 0.000000000
#> 2 0.8598357 0.1391509 0.00000000 0.00000000 0.00000000 0 0.000000000
#> 3 0.7939162 0.1196487 0.07578761 0.00000000 0.00000000 0 0.008359086
#> 4 0.7592510 0.1104679 0.06515199 0.04063795 0.00000000 0 0.020638371
#> 5 0.7375972 0.1056427 0.06014976 0.03493326 0.02164591 0 0.034426300
#> D
#> 1 0.000000000
#> 2 0.001013486
#> 3 0.002288444
#> 4 0.003852790
#> 5 0.005604840
#>
#> , , StrategyA
#>
#> H S1_tnl1 S1_tnl2 S1_tnl3 S1_tnl4 S1_tnl5 S2
#> 1 1.0000000 0.0000000 0.00000000 0.00000000 0.00000000 0 0.000000000
#> 2 0.8598357 0.1391509 0.00000000 0.00000000 0.00000000 0 0.000000000
#> 3 0.7939162 0.1196487 0.07578761 0.00000000 0.00000000 0 0.008359086
#> 4 0.7592510 0.1104679 0.06515199 0.04063795 0.00000000 0 0.020638371
#> 5 0.7375972 0.1056427 0.06014976 0.03493326 0.02164591 0 0.034426300
#> D
#> 1 0.000000000
#> 2 0.001013486
#> 3 0.002288444
#> 4 0.003852790
#> 5 0.005604840
#>
#> , , StrategyB
#>
#> H S1_tnl1 S1_tnl2 S1_tnl3 S1_tnl4 S1_tnl5 S2
#> 1 1.0000000 0.0000000 0.00000000 0.00000000 0.00000000 0 0.000000000
#> 2 0.8598357 0.1391509 0.00000000 0.00000000 0.00000000 0 0.000000000
#> 3 0.7939162 0.1196487 0.07906904 0.00000000 0.00000000 0 0.005077653
#> 4 0.7605380 0.1104679 0.06797292 0.04451414 0.00000000 0 0.012678655
#> 5 0.7413303 0.1058218 0.06275411 0.03826531 0.02496218 0 0.021348758
#> D
#> 1 0.000000000
#> 2 0.001013486
#> 3 0.002288444
#> 4 0.003828382
#> 5 0.005517580
#>
#> , , StrategyAB
#>
#> H S1_tnl1 S1_tnl2 S1_tnl3 S1_tnl4 S1_tnl5 S2
#> 1 1.0000000 0.0000000 0.00000000 0.00000000 0.00000000 0 0.000000000
#> 2 0.8598357 0.1391509 0.00000000 0.00000000 0.00000000 0 0.000000000
#> 3 0.7939162 0.1196487 0.07906904 0.00000000 0.00000000 0 0.005077653
#> 4 0.7605380 0.1104679 0.06797292 0.04451414 0.00000000 0 0.012678655
#> 5 0.7413303 0.1058218 0.06275411 0.03826531 0.02496218 0 0.021348758
#> D
#> 1 0.000000000
#> 2 0.001013486
#> 3 0.002288444
#> 4 0.003828382
#> 5 0.005517580
#>
#>
#> $summary_payoffs
#> cost utility
#> StandardOfCare 12963.21 4.756449
#> StrategyA 23002.89 4.911092
#> StrategyB 23498.72 4.763860
#> StrategyAB 33479.50 4.922386create datasets for the expanded variables
gmod_expand_functions(mygmod)
#> Note: The dataset df_cost created for function cost.
#> Note: The dataset df_pDie created for function pDie.
#> Note: The dataset df_pGetSick created for function pGetSick.
#> Note: The dataset df_pProgress created for function pProgress.
#> Note: The dataset df_pRecover created for function pRecover.
#> Note: The dataset df_utility created for function utility.removes probability paths that are 0, and simplifies those that are
1. Then, we build the gmod Markov model by using the
gmod_build() function.
The output is the structure of the model including the formulae.
To evaluate the model, we run gmod_evaluate() which
returns the numerical outcomes.
The 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:
Similarly, we can load the summary outcomes of the costs and utility of the strategies, such that: