M0 gmod: getting started with building Markov models • gmod

The purpose of this vignette is to demonstrate how gmod can be used to build common features of a Markov model. For decision modeling, please check the getting started with decision modeling vignette.

For a graphical interface to build the gmod syntax, please check https://www.dashlab.ca/projects/decision_twig/.

This vignette only explains the features one at a time to emphasize each individual feature. For examples on how to build a complete Markov model or a decision model, please refer to the detailed introductory and advanced Markov tutorial vignettes for Markov models, and the Herpes Simplex Virus (HSV) and the Doubilet vignettes for decision tree examples.

Overview of gmod:

gmod is built on the grammar of modeling principle which aims to describe a model in terms of its main components. The idea is to reduce the complexity of a decision/Markov model to it’s principal components (Twigs). gmod is inspired by the grammer of graphics (ggplot) and the grammar of data wranggling (dplyr) package for data wrangling. Similar to ggplot which builds a graph by combining layers, gmod then writes the model code for you!

gmod functions:

gmod has two main external functions. The first function is gmod() which instantiates the basic gmod object that can be supplemented with one or more twigs containing the components of the Markov model.

a gmod syntax for a Markov model currently consist of the following twigs:

decisions: a single twig that lists the decisions or strategies that we want to compare
states: a single twig that describes the Markov states, their initial probabitlities and the number of tunnel cycles (defaults to 1).
event: one or more twig describing, how events control transition among states. These transitions can be conditional on the decisions, Markov states, previous events, number of cycles, and number of cycles in a state (i.e., tunnels).
payoffs: a single payoff twig, such as costs, effectiveness, utility, … etc. Similar to the events, these can also depend on the decisions, Markov states, events, number of cycles, and number of cycles in a state (i.e., tunnels).

The second function is gmod_gen_model_function() which returns mainly creates a new function named my_markov_model that represents the model code in standard R that can be run independently of gmod. This function also returns other useful information, such as the model structure and the equations involved in the transition probabilities and payoffs …etc.

The my_markov_model() represents the entire model structure, and can be called directly. This function takes in a parameter list param which include the model parameter values, and returns the Markov results, such as the discounted lifetime costs and quality adjusted life-years (QALYs). It is important to pass the number of cycles n_cycles as a parameter as we discussed below. As described below, in addition to the main payoffs, this function can also return other intermediate outputs, such as the transitional probability matrices, the cohort’s Markov trace over time, and others.

We will start with building a minimal Markov model to illustrate the basic principles of gmod. Then, we will introduce how to add other functionality, such as cycle dependency, tunnels…etc. For each feature, we will only show how to add that feature to reduce complexity. For full models that contain many of these features in the same model, please refer to the introductory and advanced Markov tutorials vignettes.

A minimal example

We will adapt the sick-sicker model (ref) that we often use in our DARTH materials and teachings.

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)

Then, load the gmod library and the magrittr library needed for the pipe %>% operator.

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)

defining a basic Markov model using gmod() and Decision Twigs

We will start with a basic gmod() syntax looks like this:

mygmod <- gmod() + 
  decisions("StandardOfCare", "StrategyA") + 
  states(names = c("H", "S", "D"), init_probs = c(1,0,0)) + 
  event(name = "only_event",  
        scenarios = c("get_sick", "die", "stay"), 
        probs = c(pGetSick(state), pDie(decision, state), Inf), 
        outcomes = c("S", "D", "curr_state")) +  
  payoffs(cost = compute_cost(decision, state))

To have a graphical representation, we can copy this code into Decision Twigs to generate the twigs:

Here are the explanation of what each twig does:

decisions("StandardOfCare", "StrategyA"): define the first twig which is a decision with two strategies involved: StandardOfCare and StrategyA
states(...): defines the Markov twig. There must be only one of these per model. The states funciton takes in 3 arguments:
- names: H for healthy, S for sick and D for dead.
- init_probs for initial probabilities, and
- optional max_cycle_in_states which defines the maximum number of tunnel states per state. This option is omitted here, so the function uses the default of 1.
event(...): is the most powerful feature in gmod and it governs how the states are linked through events. In the Grammar of Modeling we distinguish between a state and an event. We use a noun for a state to indicate where the cohort stays a full cycle (e.g., one year), and use a verb to refer to an event, which in a Markov world is instantaneous and doesn’t take any time at all (e.g., get_sick for getting sick). There has to be at least one event twig, but there can be more events that can be chained together as we explain below. The event mapping takes in four arguments: name, scenarios, outcomes, and probs. We will explain each of these separately:
- event = "only_event": defines the event name. Each event must have a unique name that can be referenced by other events’s probability statements or payoff equations as described below.
- scenarios = c("get_sick", "die", "stay"): defines the event scenarios: here the event can only be one of these: getting sick, dying or staying in the current state.
- probs = c(pGetSick(state), pDie(decision, state), Inf)): these are the probabilities of the correspondign scenarios occuring. A powerful feature of gmod is that it allows these probabilities to be conditional on other model components. Here we define the probability of only_event==get_sick by calling pGetSick(state) which is dependent on the state. Similarly, the probability of only_event==die is defined by calling pDie(decision, state) which we tell gmod that is dependent on the state and the decision. pGetSick and pDie are user-defined functions that we describe below. The Inf is used to tell gmod to compute the complementary probability for the probability of staying in the same state: (i.e., 1-pGetSick(state)-pDie(state)) to reduce the amount of typing needed.
- outcomes = c("S", "D", "curr_state"): defines the outcomes of each one of the event values. getting sick ends up in the S state, dying ends up in the D state, and staying will return to whichever state the cohort was in using the special keyword "curr_state".
payoffs(cost = cost(decision, state)): defines a single payoff hers. Here we only have a cost payoff that is dependent on the state and decision. More than one payoff can be described here, but only a single payoff twig is allowed. And, similar to the probs argument these can be a function of other model components, such as decision, state, cycles, events, tunnels, …etc

User-defined functions for probabilities and payoffs

In this section, we define the three functions: pGetSick, pDie and cost that we use above in our event_mapping and the payoffs() layers above:

pGetSick <- function(state){
  if(state=="H") 0.1 else 0
}

pDie <- function(decision, state){
  switch(state, "H"=0.01, 
         "S"= if (decision=="StrategyA") 0.02 else 0.05, 
         "D"=0)
}

compute_cost <- function(decision, state){
  cost_state <- switch(state, "H"=1000, "S"=5000, "D"=0)
  cost_decision <- if(decision=="StrategyA" & state=="S") 2000 else 0
  cost_state+cost_decision
}

pGetSick is a simple function with a single if else statement. Probability of getting sick is 0.1 if one is Healthy, otherwise is 0.

pDie is a function of bothe the decision and state. Probability of death is 0.01 if Healthy, 0.02 if sick (+TreatmentA), 0.05 if Sick+StandardOfCare, and 0 if already dead.

compute_cost calculates the cost of the disease by combining the cost of the state $1000 if Healthy, $5000 if sick and $0 if Dead. In addition, it computes the cost of the decision as an additional $2000 if one receives TreatmentA compared to the standard of care.

As your functions get more complex, we recommended that you test them to make sure they return outcomes for all possible combinations of states and decisions and other input arguments. For example, the simple code below tests the functions for all possible combinations of decisions and states.

Evaluate the functions

gmod_expand_functions(mygmod)
#> Note: The dataset df_compute_cost created for function compute_cost.
#> Note: The dataset df_pDie created for function pDie.
#> Note: The dataset df_pGetSick created for function pGetSick.
head(df_pDie)
#>         decision state pDie
#> 1 StandardOfCare     H 0.01
#> 2      StrategyA     H 0.01
#> 3 StandardOfCare     S 0.05
#> 4      StrategyA     S 0.02
#> 5 StandardOfCare     D 0.00
#> 6      StrategyA     D 0.00

for (decision in c("StandardOfCare", "StrategyA")){
  for (state in c("H", "S", "D")){
    cat("decision:", decision, " state:", state, " pGetSick:", pGetSick(state), 
        " pDie:", pDie(decision, state), " cost:", compute_cost(decision, state), "\n")
  }
}
#> decision: StandardOfCare  state: H  pGetSick: 0.1  pDie: 0.01  cost: 1000 
#> decision: StandardOfCare  state: S  pGetSick: 0  pDie: 0.05  cost: 5000 
#> decision: StandardOfCare  state: D  pGetSick: 0  pDie: 0  cost: 0 
#> decision: StrategyA  state: H  pGetSick: 0.1  pDie: 0.01  cost: 1000 
#> decision: StrategyA  state: S  pGetSick: 0  pDie: 0.02  cost: 7000 
#> decision: StrategyA  state: D  pGetSick: 0  pDie: 0  cost: 0

The functions seem to be outputting the correct numbers. switch is similar to if else, but can be more compact than multiple if else statements.

Running the model

Now we have all the ingredients to run the model.

First, we need to generate the model function code

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)

This will create the model_struc which contains the model equations and formulae and can be examined, and as stated in the generated note, this function also produces a function representation of the Markov model. The default name of this function is my_markov_model.

Let’s examine the content of these two object. First, let’s check the model structure

model_struc # returns the model structure including the formulae before evaluation
#> $is_cycle_dep
#> [1] FALSE
#> 
#> $decisions
#> [1] "StandardOfCare" "StrategyA"     
#> 
#> $n_decisions
#> [1] 2
#> 
#> $states
#> [1] "H" "S" "D"
#> 
#> $n_states
#> [1] 3
#> 
#> $states_expanded
#> [1] "H" "S" "D"
#> 
#> $n_states_expanded
#> [1] 3
#> 
#> $n_events
#> [1] 1
#> 
#> $p0
#> $p0$StandardOfCare
#> H S D 
#> 1 0 0 
#> 
#> $p0$StrategyA
#> H S D 
#> 1 0 0 
#> 
#> 
#> $events
#> [1] "only_event"
#> 
#> $events_df
#>    type      event   values                                    probs   outcomes
#> 1 event only_event get_sick                          pGetSick(state)          S
#> 2 event only_event      die                     pDie(decision,state)          D
#> 3 event only_event     stay 1-(pGetSick(state)+pDie(decision,state)) curr_state
#>   id
#> 1  1
#> 2  2
#> 3  3
#> 
#> $payoffs
#> $payoffs$cost
#> compute_cost(decision, state)
#> 
#> 
#> $payoff_names
#> [1] "cost"
#> 
#> $n_payoffs
#> [1] 1
#> 
#> $discounts
#> cost 
#>    0 
#> 
#> $model_equations
#> # A tibble: 3 × 3
#>   dest       probs                                      cost                    
#>   <chr>      <chr>                                      <chr>                   
#> 1 D          (pDie(decision,state))                     (pDie(decision,state))*…
#> 2 S          (pGetSick(state))                          (pGetSick(state))*compu…
#> 3 curr_state (1-(pGetSick(state)+pDie(decision,state))) (1-(pGetSick(state)+pDi…
#> 
#> attr(,"class")
#> [1] "gmod_markov"

examining the content of this structure can be helpful to debug a model by examining the individual layers and formulae. For example $model_equations returns the probability and cost formulae as plain text. These can be evaluated by calling them directly if needed. For example, if we want to know what is the probability of staying at the Healthy state we can just copy the content of the probs column corresponding to curr_state and evaluate it with the corresponding parameters:

(1-(pGetSick(state="H")+pDie(decision="StandardOfCare",state="H"))) 
#> [1] 0.89

Next, we can also examine the contents of the my_markov_model function which generates the text.

my_markov_model #without parenthesis returns the function text
#> function (model_struc, params = NULL, return_transition_prob = FALSE, 
#>     return_state_payoffs = FALSE, return_trace = FALSE, return_cycle_payoffs = FALSE, 
#>     return_payoff_summary = TRUE) 
#> {
#>     if (!is.null(params)) 
#>         list2env(params, envir = .GlobalEnv)
#>     decisions <- model_struc$decisions
#>     tunnel_lengths <- model_struc$tunnel_lengths
#>     tunnel_states <- model_struc$tunnel_states
#>     tunnel_lengths[is.infinite(tunnel_lengths)] <- n_cycles
#>     n_decisions <- model_struc$n_decisions
#>     cycles <- 1:n_cycles
#>     states_expanded <- model_struc$states_expanded
#>     n_states_expanded <- model_struc$n_states_expanded
#>     payoff_names <- model_struc$payoff_names
#>     n_payoffs <- model_struc$n_payoffs
#>     discounts <- eval(model_struc$discounts)
#>     P <- array(0, dim = c(n_states_expanded, n_states_expanded, 
#>         n_decisions), dimnames = list(states_expanded, states_expanded, 
#>         decisions))
#>     state_payoffs <- array(0, dim = c(n_states_expanded, n_decisions, 
#>         n_payoffs), dimnames = list(states_expanded, decisions, 
#>         payoff_names))
#>     dest <- "D"
#>     for (state_expanded in states_expanded) {
#>         state <- tunnel2state(state_expanded)$state
#>         for (decision in decisions) {
#>             P[state_expanded, dest, decision] <- (pDie(decision, 
#>                 state))
#>         }
#>     }
#>     dest <- "S"
#>     for (state_expanded in states_expanded) {
#>         state <- tunnel2state(state_expanded)$state
#>         P[state_expanded, dest, ] <- (pGetSick(state))
#>     }
#>     dest <- "curr_state"
#>     for (state_expanded in states_expanded) {
#>         state <- tunnel2state(state_expanded)$state
#>         for (decision in decisions) {
#>             dest <- state_expanded
#>             P[state_expanded, dest, decision] <- P[state_expanded, 
#>                 dest, decision] + (1 - (pGetSick(state) + pDie(decision, 
#>                 state)))
#>         }
#>     }
#>     for (state_expanded in states_expanded) {
#>         state <- state_expanded
#>         for (decision in decisions) {
#>             state_payoffs[state_expanded, decision, "cost"] <- (pDie(decision, 
#>                 state)) * compute_cost(decision, state) + (pGetSick(state)) * 
#>                 compute_cost(decision, state) + (1 - (pGetSick(state) + 
#>                 pDie(decision, state))) * compute_cost(decision, 
#>                 state)
#>         }
#>     }
#>     Trace <- array(NA, dim = c(n_cycles, n_states_expanded, n_decisions), 
#>         dimnames = list(cycles, states_expanded, decisions))
#>     cycle_payoffs <- array(0, dim = c(n_cycles, n_states_expanded, 
#>         n_decisions, n_payoffs), dimnames = list(cycles, states_expanded, 
#>         decisions, payoff_names))
#>     for (decision in decisions) {
#>         Trace[1, , decision] <- model_struc$p0[[decision]]
#>         for (i in 2:n_cycles) {
#>             P_temp <- P[, , decision]
#>             Trace[i, , decision] <- Trace[i - 1, , decision] %*% 
#>                 P_temp
#>         }
#>     }
#>     cycle_ones <- matrix(1, nrow = n_cycles, ncol = 1)
#>     state_ones <- matrix(1, nrow = 1, ncol = n_states_expanded)
#>     summary_payoffs <- matrix(0, nrow = n_decisions, ncol = n_payoffs, 
#>         dimnames = list(decisions, payoff_names))
#>     for (decision in decisions) {
#>         for (payoff_name in payoff_names) {
#>             mat_discounts <- matrix(1/(1 + discounts[payoff_name])^(cycles - 
#>                 1), ncol = 1) %*% state_ones
#>             Payoff_temp <- cycle_ones %*% state_payoffs[, decision, 
#>                 payoff_name]
#>             cycle_payoffs[, , decision, payoff_name] <- Trace[, 
#>                 , decision] * Payoff_temp * mat_discounts
#>             summary_payoffs[decision, payoff_name] <- sum(cycle_payoffs[, 
#>                 , decision, payoff_name])
#>         }
#>     }
#>     model_results <- list()
#>     if (return_transition_prob) 
#>         model_results$P <- P
#>     if (return_state_payoffs) 
#>         model_results$state_payoffs <- state_payoffs
#>     if (return_trace) 
#>         model_results$Trace <- Trace
#>     if (return_cycle_payoffs) 
#>         model_results$cycle_payoffs <- cycle_payoffs
#>     if (return_payoff_summary) 
#>         model_results$summary_payoffs <- summary_payoffs
#>     return(model_results)
#> }
#> <environment: 0x126f50e18>

The generated code is build on our introductory and advanced Markov tutorials in Medical Decision Making. For further information, please refer to these tutorials.

Now, let’s run the model function and generate our summary results of running the model for 5 cycles and computing the cumulative cost for our cohort.

n_cycles <- 10
my_markov_model(model_struc)
#> $summary_payoffs
#>                    cost
#> StandardOfCare 20998.32
#> StrategyA      28735.17

You can also retrieve the other components of the Markov model run, such as the Markov trace which shows the proportion of the cohort in each of the states for each of the 5 cycles.

model_results <- my_markov_model(model_struc, return_trace = T)
model_results$Trace
#> , , StandardOfCare
#> 
#>            H         S          D
#> 1  1.0000000 0.0000000 0.00000000
#> 2  0.8900000 0.1000000 0.01000000
#> 3  0.7921000 0.1840000 0.02390000
#> 4  0.7049690 0.2540100 0.04102100
#> 5  0.6274224 0.3118064 0.06077119
#> 6  0.5584059 0.3589583 0.08263573
#> 7  0.4969813 0.3968510 0.10616771
#> 8  0.4423133 0.4267066 0.13098007
#> 9  0.3936589 0.4496026 0.15673853
#> 10 0.3503564 0.4664883 0.18315525
#> 
#> , , StrategyA
#> 
#>            H         S          D
#> 1  1.0000000 0.0000000 0.00000000
#> 2  0.8900000 0.1000000 0.01000000
#> 3  0.7921000 0.1870000 0.02090000
#> 4  0.7049690 0.2624700 0.03256100
#> 5  0.6274224 0.3277175 0.04486009
#> 6  0.5584059 0.3839054 0.05768866
#> 7  0.4969813 0.4320679 0.07095083
#> 8  0.4423133 0.4731246 0.08456200
#> 9  0.3936589 0.5078935 0.09844763
#> 10 0.3503564 0.5371015 0.11254209

changing basic features of our Markov model:

We start by modifying some of the basic features of our markov model. We will explore how to change the number of Markov cycles, the initial probabilities, and the discount rates.

changing the number of cycles:

Changing the number of Markov cycles is straightforward by changing the n_cycles parameter in the gmod function. For example, below we change the number of cycles from 5 to 2:

mygmod <- gmod()

changing initial probabilities:

We can set the initial probabilities in various ways. Most of the time we start the cohort in the same health state. Above we showed how to set it up when everyone starts at the Healthy state.

states(names = c("H","S","D"), init_probs = c(1,0,0))
#> [[1]]
#> [[1]]$type
#> [1] "states"
#> 
#> [[1]]$states
#> [1] "H" "S" "D"
#> 
#> 
#> [[2]]
#> [[2]]$type
#> [1] "initial_prob"
#> 
#> [[2]]$states
#> [1] "H" "S" "D"
#> 
#> [[2]]$probs
#> [1] 1 0 0

However, we can do any distribution we want. For example, we can define the starting probabilities for all the states.

states(names = c("H","S","D"), init_probs = c(0.7,0.2,0.1))
#> [[1]]
#> [[1]]$type
#> [1] "states"
#> 
#> [[1]]$states
#> [1] "H" "S" "D"
#> 
#> 
#> [[2]]
#> [[2]]$type
#> [1] "initial_prob"
#> 
#> [[2]]$states
#> [1] "H" "S" "D"
#> 
#> [[2]]$probs
#> [1] 0.7 0.2 0.1

Or, we can even use a combination of variables and Inf as a placeholder for the complementary probabilities.

pH0 <- 0.7
pS0 <- 0.2
states(names = c("H","S","D"), init_probs = c(pH0, pS0, Inf))
#> [[1]]
#> [[1]]$type
#> [1] "states"
#> 
#> [[1]]$states
#> [1] "H" "S" "D"
#> 
#> 
#> [[2]]
#> [[2]]$type
#> [1] "initial_prob"
#> 
#> [[2]]$states
#> [1] "H" "S" "D"
#> 
#> [[2]]$probs
#> [1] 0.7 0.2 Inf

To run the model with this updated probabilities, we need to add it as a new init_prob layer to the gmod code.

mygmod <- gmod() + 
  decisions("StandardOfCare", "StrategyA") + 
  states(names = c("H","S","D"), init_probs = c(pH0, pS0, Inf)) + 
  event(name = "only_event",  
              scenarios = c("get_sick", "die", "stay"), 
              probs = c(pGetSick(state), pDie(decision, state), Inf), 
              outcomes = c("S", "D", "curr_state")) +  
  payoffs(cost = compute_cost(decision, state))
  
# creating the updated function and evaluating it
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_results <- my_markov_model(model_struc, return_trace = T)
model_results$Trace
#> , , StandardOfCare
#> 
#>            H         S         D
#> 1  0.7000000 0.2000000 0.1000000
#> 2  0.6230000 0.2600000 0.1170000
#> 3  0.5544700 0.3093000 0.1362300
#> 4  0.4934783 0.3492820 0.1572397
#> 5  0.4391957 0.3811657 0.1796386
#> 6  0.3908842 0.4060270 0.2030888
#> 7  0.3478869 0.4248141 0.2272990
#> 8  0.3096193 0.4383621 0.2520186
#> 9  0.2755612 0.4474059 0.2770329
#> 10 0.2452495 0.4525917 0.3021588
#> 
#> , , StrategyA
#> 
#>            H         S         D
#> 1  0.7000000 0.2000000 0.1000000
#> 2  0.6230000 0.2660000 0.1110000
#> 3  0.5544700 0.3229800 0.1225500
#> 4  0.4934783 0.3719674 0.1345543
#> 5  0.4391957 0.4138759 0.1469284
#> 6  0.3908842 0.4495179 0.1595979
#> 7  0.3478869 0.4796160 0.1724971
#> 8  0.3096193 0.5048124 0.1855683
#> 9  0.2755612 0.5256780 0.1987607
#> 10 0.2452495 0.5427206 0.2120299

We can see from the first line of the Trace that H=0.7, S=0.2 and D=0.1.

gmod tries to check the composition of the probability vector and returns an error if there are any inconsistencies.

changing the discount rates

Changing the discount rate is similar to setting the initial probabilities. These can be set for all the payoffs. To illustrate, we will add an additional payoff (effectiveness), and use a discount rate of 3% per cycle for both costs and utilities.

mygmod <- gmod() + 
  decisions("StandardOfCare", "StrategyA") + 
  states(names = c("H", "S", "D"), init_probs = c(1,0,0)) +
  event(name = "only_event",  
        scenarios = c("get_sick", "die", "stay"), 
        probs = c(pGetSick(state), pDie(decision, state), Inf), 
        outcomes = c("S", "D", "curr_state")) +  
  payoffs(cost = compute_cost(decision, state), 
          effectiveness = compute_effectiveness(state), 
          discount_rates = c(0.03, 0.03))

# define the effectiveness function, the rest of the function will be just like before
compute_effectiveness <- function(state){
  switch(state, "H"=1, "S"=0.7, "D"=0)
}

# creating the updated function and evaluating it
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)
# we look at the cycle payoffs which apply the discounts to the cycle payoffs.
my_markov_model(model_struc, return_cycle_payoffs=TRUE)
#> $cycle_payoffs
#> , , StandardOfCare, cost
#> 
#>            H         S D
#> 1  1000.0000    0.0000 0
#> 2   864.0777  485.4369 0
#> 3   746.6302  867.1882 0
#> 4   645.1465 1162.2757 0
#> 5   557.4567 1385.1797 0
#> 6   481.6859 1548.2030 0
#> 7   416.2140 1661.7823 0
#> 8   359.6412 1734.7575 0
#> 9   310.7580 1774.6022 0
#> 10  268.5190 1787.6224 0
#> 
#> , , StrategyA, cost
#> 
#>            H         S D
#> 1  1000.0000    0.0000 0
#> 2   864.0777  679.6117 0
#> 3   746.6302 1233.8580 0
#> 4   645.1465 1681.3806 0
#> 5   557.4567 2038.2093 0
#> 6   481.6859 2318.1211 0
#> 7   416.2140 2532.9503 0
#> 8   359.6412 2692.8555 0
#> 9   310.7580 2806.5507 0
#> 10  268.5190 2881.5051 0
#> 
#> , , StandardOfCare, effectiveness
#> 
#>            H          S D
#> 1  1.0000000 0.00000000 0
#> 2  0.8640777 0.06796117 0
#> 3  0.7466302 0.12140635 0
#> 4  0.6451465 0.16271859 0
#> 5  0.5574567 0.19392516 0
#> 6  0.4816859 0.21674842 0
#> 7  0.4162140 0.23264952 0
#> 8  0.3596412 0.24286605 0
#> 9  0.3107580 0.24844430 0
#> 10 0.2685190 0.25026713 0
#> 
#> , , StrategyA, effectiveness
#> 
#>            H          S D
#> 1  1.0000000 0.00000000 0
#> 2  0.8640777 0.06796117 0
#> 3  0.7466302 0.12338580 0
#> 4  0.6451465 0.16813806 0
#> 5  0.5574567 0.20382093 0
#> 6  0.4816859 0.23181211 0
#> 7  0.4162140 0.25329503 0
#> 8  0.3596412 0.26928555 0
#> 9  0.3107580 0.28065507 0
#> 10 0.2685190 0.28815051 0
#> 
#> 
#> $summary_payoffs
#>                    cost effectiveness
#> StandardOfCare 18057.18      7.387116
#> StrategyA      24515.17      7.536633

Transition probability and payoffs conditional on Markov cycle:

Sometimes we want a probability of transition from one state to another, or a payoff’s value to depend on the cohort’s age. We can achieve this by using the special variable cycle for both transition probabilities and payoffs.

The code below shows the previous code in a single block where we modify the probability of death and cost payoff to depend on cycle. See if you can spot the difference:

mygmod <- gmod() + 
  decisions("StandardOfCare", "StrategyA") + 
  states(names = c("H", "S", "D"), init_probs = c(1,0,0)) +
  event(name = "only_event",  
              scenarios = c("get_sick", "die", "stay"), 
              probs = c(pGetSick(state), pDie(decision, state, cycle), Inf), 
              outcomes = c("S", "D", "curr_state")) +  
  payoffs(cost = compute_cost(decision, state, cycle))

pGetSick <- function(state){
  if(state=="H") 0.1 else 0
}

pDie <- function(decision, state, cycle){
  switch(state, "H"=0.01, 
       "S"= if (decision=="StrategyA") 0.02 else 0.05, 
       "D"=0) * cycle
}

compute_cost <- function(decision, state, cycle){
  cost_state <- switch(state, "H"=1000, "S"=5000, "D"=0) * cycle
  cost_decision <- if(decision=="StrategyA" & state=="S") 2000 else 0
  cost_state+cost_decision
}

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)

Of course, we only need to add the variable cycle to the function definitions pDie and compute_cost and also remember to add them to the function calls in the gmod function. For illustration we assume that pDie and the cost of the states both increased as a linear function of the Markov cycles, 1X, 2X, …, 5X for cycles 1 through 5, respectively. Let’s check first to see that our functions are correct, and then we will run the model and inspect the outputs.

for (decision in c("StandardOfCare", "StrategyA")){
  for (state in c("H", "S", "D")){
    for (cycle in 1:2){ # just print the first two cycles
          cat("decision:", decision, " state:", state, " cycle:", cycle, " pGetSick:", pGetSick(state), 
        " pDie:", pDie(decision, state, cycle), " cost:", compute_cost(decision, state, cycle), "\n")
    }
  }
}
#> decision: StandardOfCare  state: H  cycle: 1  pGetSick: 0.1  pDie: 0.01  cost: 1000 
#> decision: StandardOfCare  state: H  cycle: 2  pGetSick: 0.1  pDie: 0.02  cost: 2000 
#> decision: StandardOfCare  state: S  cycle: 1  pGetSick: 0  pDie: 0.05  cost: 5000 
#> decision: StandardOfCare  state: S  cycle: 2  pGetSick: 0  pDie: 0.1  cost: 10000 
#> decision: StandardOfCare  state: D  cycle: 1  pGetSick: 0  pDie: 0  cost: 0 
#> decision: StandardOfCare  state: D  cycle: 2  pGetSick: 0  pDie: 0  cost: 0 
#> decision: StrategyA  state: H  cycle: 1  pGetSick: 0.1  pDie: 0.01  cost: 1000 
#> decision: StrategyA  state: H  cycle: 2  pGetSick: 0.1  pDie: 0.02  cost: 2000 
#> decision: StrategyA  state: S  cycle: 1  pGetSick: 0  pDie: 0.02  cost: 7000 
#> decision: StrategyA  state: S  cycle: 2  pGetSick: 0  pDie: 0.04  cost: 12000 
#> decision: StrategyA  state: D  cycle: 1  pGetSick: 0  pDie: 0  cost: 0 
#> decision: StrategyA  state: D  cycle: 2  pGetSick: 0  pDie: 0  cost: 0

Now we can examine the impact of this on the model results

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)
my_markov_model(model_struc)
#> $summary_payoffs
#>                     cost
#> StandardOfCare  64950.86
#> StrategyA      100340.56

Transition probabilty and payoff dependency on the number of cycles in a state (tunnel states or state residency):

Sometimes we have a different situation where we want to keep track of how long the cohort has been in a particular state to model history of a disease. We usually do this by creating what is called tunnel states and expanding the number of states for example S1, S2, …, etc. where a transition to S2 is only possible from S1 and to S3 from S2, … etc.

Here we show how to do this again for pDie and compute_cost. We are only interested of a history of 3 years. Again - can you spot the difference compared to the minimal model?

mygmod <- gmod() + 
  decisions("StandardOfCare", "StrategyA") +
  states(names = c("H", "S", "D"), init_probs = c(1,0,0), max_cycle_in_states = c(1,3,1)) +  #
  event(name = "only_event",  
              scenarios = c("get_sick", "die", "stay"), 
              probs = c(pGetSick(state), pDie(decision, state, cycle_in_state), Inf), 
              outcomes = c("S", "D", "curr_state")) + #
  payoffs(cost = compute_cost(decision, state, cycle_in_state)) #

pGetSick <- function(state){
  if(state=="H") 0.1 else 0
}

pDie <- function(decision, state, cycle_in_state){ #
  switch(state, "H"=0.01, 
       "S"= cycle_in_state * if (decision=="StrategyA") 0.02 else 0.05, 
       "D"=0) 
}

compute_cost <- function(decision, state, cycle_in_state){ #
  cost_state <- switch(state, "H"=1000, "S"=5000 * cycle_in_state, "D"=0)  #
  cost_decision <- if(decision=="StrategyA" & state=="S") 2000 else 0
  cost_state+cost_decision
}

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)

It is quite simple to add tunnel dependency by adding the cycle_in_state variable.

IMPORTANT: when a model has at least on state that has time dependency, the cycle_in_state variable becomes avaialble for all the states. States that are not defined in the tunnels() layer will always return cycle_in_state=1

Similar to above, we assume that the probability of death is a linear function of the number of cycles in S. For the cost, we also assume that the longer you stay in S the higher the cost will be - also linearly related.

The other important change is that we have to add a tunnels() layer specifying the states that have residency time history dependency. Here we assume only S has time dependency for 3 cycles only. This can be expanded to multiple states with various tunnel lengths.

Let’s check first to see that our functions are correct, and then we will run the model and inspect the outputs.

for (decision in c("StandardOfCare", "StrategyA")){
  for (state in c("H", "S", "D")){
    for (cycle_in_state in 1:2){ # just print the first two cycles
          cat("decision:", decision, " state:", state, " cycle_in_state:", cycle_in_state, " pGetSick:", pGetSick(state), 
        " pDie:", pDie(decision, state, cycle_in_state), " cost:", compute_cost(decision, state, cycle_in_state), "\n")
    }
  }
}
#> decision: StandardOfCare  state: H  cycle_in_state: 1  pGetSick: 0.1  pDie: 0.01  cost: 1000 
#> decision: StandardOfCare  state: H  cycle_in_state: 2  pGetSick: 0.1  pDie: 0.01  cost: 1000 
#> decision: StandardOfCare  state: S  cycle_in_state: 1  pGetSick: 0  pDie: 0.05  cost: 5000 
#> decision: StandardOfCare  state: S  cycle_in_state: 2  pGetSick: 0  pDie: 0.1  cost: 10000 
#> decision: StandardOfCare  state: D  cycle_in_state: 1  pGetSick: 0  pDie: 0  cost: 0 
#> decision: StandardOfCare  state: D  cycle_in_state: 2  pGetSick: 0  pDie: 0  cost: 0 
#> decision: StrategyA  state: H  cycle_in_state: 1  pGetSick: 0.1  pDie: 0.01  cost: 1000 
#> decision: StrategyA  state: H  cycle_in_state: 2  pGetSick: 0.1  pDie: 0.01  cost: 1000 
#> decision: StrategyA  state: S  cycle_in_state: 1  pGetSick: 0  pDie: 0.02  cost: 7000 
#> decision: StrategyA  state: S  cycle_in_state: 2  pGetSick: 0  pDie: 0.04  cost: 12000 
#> decision: StrategyA  state: D  cycle_in_state: 1  pGetSick: 0  pDie: 0  cost: 0 
#> decision: StrategyA  state: D  cycle_in_state: 2  pGetSick: 0  pDie: 0  cost: 0

Now we can examine the impact of this on the model results

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_results <- my_markov_model(model_struc, return_transition_prob = T, return_state_payoffs = T)

Let’s examine the transition probability matrix structure with the tunnel state

model_results$P # shows transition probabilities to Death
#> , , StandardOfCare
#> 
#>           H S_tnl1 S_tnl2 S_tnl3    D
#> H      0.89    0.1   0.00   0.00 0.01
#> S_tnl1 0.00    0.0   0.95   0.00 0.05
#> S_tnl2 0.00    0.0   0.00   0.90 0.10
#> S_tnl3 0.00    0.0   0.00   0.85 0.15
#> D      0.00    0.0   0.00   0.00 1.00
#> 
#> , , StrategyA
#> 
#>           H S_tnl1 S_tnl2 S_tnl3    D
#> H      0.89    0.1   0.00   0.00 0.01
#> S_tnl1 0.00    0.0   0.98   0.00 0.02
#> S_tnl2 0.00    0.0   0.00   0.96 0.04
#> S_tnl3 0.00    0.0   0.00   0.94 0.06
#> D      0.00    0.0   0.00   0.00 1.00

This returns a list of two 5x5 matrices for the transition probabilities of the two strategies. These matrices describe the probability of transition from an origin state (row) to a destination state (column). Now notice that S is being replaced by S_tnl1, S_tnl2 and S_tnl3 indicating the three cycles in S as separate states. Also notice the special way transitions occur among these states. i.e., each tunnel state goes to the next one, but never back or stay in the same state except for the last state in the tunnel, here S_tnl3. Since we made the probability of death a function of the tunnel state, we can see that the probabiltiy of death for StandardOfCare (last column) for S_tnl1 = 0.05, S_tnl2=0.10 and S_tnl3=0.15 as we specified in our function. Similarly, the probability of death under StrategyA also incrases linearly from 0.02 for S_tnl1, to 0.04 for S_tnl2 and 0.06 for S_tnl3.

And, let’s examine the cycle payoffs.

model_results$state_payoffs # shows cycle payoffs of the Sick state
#> , , cost
#> 
#>        StandardOfCare StrategyA
#> H                1000      1000
#> S_tnl1           5000      7000
#> S_tnl2          10000     12000
#> S_tnl3          15000     17000
#> D                   0         0

Here you can see that for each year in S, the cost increases by $5000.

Multiple Events.

It is relatively easy to add additional events to a Markov model. Please keep in mind that Events are quite flexible and are central to the Grammar of Modeling framework. While users are free to make branching events, we recommend to use a single linear sequence of events, and call previous event values to make the same sequence apply to all decisions and states. For example, if an event doesn’t happen for a particular decision (e.g., getting a side effect), then one can make the get_side_effect event a function of the decision, and only have a non-zero probability for that particular treatment. The Grammer of Modeling is designed aroudn the concept of a single sequence of events as it makes the code more concise and easier to manage.

Below we slightly modify the minimal example by converting the Sick state to an event. Assume that Sick is an acute disease that only takes a few days rather than a full cycle. We will have a sequence of events get_sick followed by die. Because Sick is no longer a state, we will use these two health states Alive and Dead and show how the probability of death can be a function of the event get_sick occuring before the event die. In addition, we will show how the cost can also be a function of the value of the sick_event. We will also look at one strategy “StandardOfCare to simplify the illustration.

mygmod <- gmod() + 
  decisions("StandardOfCare") + 
  states(names = c("Alive", "Dead"), init_probs = c(1,0)) + 
  event(name = "get_sick",  
              scenarios = c(TRUE, FALSE), 
              probs = c(0.05, Inf), 
              outcomes = c("die", "die")) +  
  event(name = "die",  
              scenarios = c(TRUE, FALSE), 
              probs = c(pDie(state), Inf), 
              outcomes = c("Dead", "curr_state")) +  
  payoffs(cost = compute_cost(state))

pDie <- function(state){
  # if the cohort gets sick, the probability of death = 0.01 otherwise it is 0.1
  switch(state, "Alive"=0.01,  
                "Dead"=0)
}

compute_cost <- function(state){
  cost_state <- switch(state, "Alive"=1000, "Dead"=0)
  # add 2000 if someone in the Alive state gets sick
  cost_state
}

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)
my_markov_model(model_struc,return_transition_prob=TRUE,return_state_payoffs=TRUE,return_trace=TRUE,return_cycle_payoffs=TRUE,return_payoff_summary=TRUE)
#> $P
#> , , StandardOfCare
#> 
#>       Alive Dead
#> Alive  0.99 0.01
#> Dead   0.00 1.00
#> 
#> 
#> $state_payoffs
#> , , cost
#> 
#>       StandardOfCare
#> Alive           1000
#> Dead               0
#> 
#> 
#> $Trace
#> , , StandardOfCare
#> 
#>        Alive       Dead
#> 1  1.0000000 0.00000000
#> 2  0.9900000 0.01000000
#> 3  0.9801000 0.01990000
#> 4  0.9702990 0.02970100
#> 5  0.9605960 0.03940399
#> 6  0.9509900 0.04900995
#> 7  0.9414801 0.05851985
#> 8  0.9320653 0.06793465
#> 9  0.9227447 0.07725531
#> 10 0.9135172 0.08648275
#> 
#> 
#> $cycle_payoffs
#> , , StandardOfCare, cost
#> 
#>        Alive Dead
#> 1  1000.0000    0
#> 2   990.0000    0
#> 3   980.1000    0
#> 4   970.2990    0
#> 5   960.5960    0
#> 6   950.9900    0
#> 7   941.4801    0
#> 8   932.0653    0
#> 9   922.7447    0
#> 10  913.5172    0
#> 
#> 
#> $summary_payoffs
#>                    cost
#> StandardOfCare 9561.792

Notice that now we have two sequential events. The first event is get_sick, the second one is die. Here get_sick doesn’t do much because if the cohort gets sick or not, they face the next event die. Next, we will show how to make the probabilit of death a function of the event get_sick and also make the cost also a function of this event.

Transition and Payoff dependent on prior events

In the previous section, we illustrated how to add multiple sequential events. gmod can also use these prior events to influence the transition probabilities downstream and the payoffs. This can be particularly important to put all the events on a single sequence across decisions and states. Traditional models often duplicate (copy or clone) these events which creates a nested branching problem that can be difficult to manage and troubleshoot.

The other advantage of being able to refer to any event value is the ability to associate payoffs to events. This is often referred to as transitional rewards or payoffs because these payoffs occur when an event occurs. These payoffs are quite different from state payoffs as state payoffs are associated with spending time in states compared to event (transitional) payoffs which are associated with the events that are considered instantaneous and only happens for a portion of the cohort residing in a state.

Below we modify the previous example by making the probability of death a function of the event get_sick occurring upstream. In addition, we will show how the cost can also be a function of the value of the sick_event. We will also look at one strategy “StandardOfCare to simplify the illustration.

mygmod <- gmod() + 
  decisions("StandardOfCare") + 
  states(names = c("Alive", "Dead"), init_probs = c(1,0)) + 
  event(name = "get_sick",  
              scenarios = c(TRUE, FALSE), 
              probs = c(0.05, Inf), 
              outcomes = c("die", "die")) +  
    event(name = "die",  
              scenarios = c(TRUE, FALSE), 
              probs = c(pDie(state, get_sick), Inf), 
              outcomes = c("Dead", "curr_state")) +  
  payoffs(cost = compute_cost(state, get_sick))

pDie <- function(state, get_sick){
  # if the cohort gets sick, the probability of death = 0.01 otherwise it is 0.1
  switch(state, "Alive"= if(get_sick==TRUE) 0.1 else 0.01,  
                "Dead"=0)
}

compute_cost <- function(state, get_sick){
  cost_state <- switch(state, "Alive"=1000, "Dead"=0)
  # add 2000 if someone in the Alive state gets sick
  cost_get_sick <- if(get_sick==TRUE & state=="Alive") 2000 else 0
  cost_state+cost_get_sick
}

Notice that now we have two sequential events. The first event is die, the second one is get_sick. Now, if die=TRUE, then the cohort moves to the state D, but if D=FALSE, then the cohort moves to the next event get_sick. Then, if get_sick=TRUE, then the cohort moves to the state S, but if not, they will stay in their current state.

The rest of the user-functions can stay the same. Below we show how to use prior events to influence the probability of transition and the payoffs.

 for (state in c("Alive", "Dead")){
   for (get_sick in c(TRUE, FALSE)){ # just print the first two cycles
         cat(" state:", state, " get_sick:", get_sick, #
       " pDie:", pDie(state, get_sick), " cost:", compute_cost(state, get_sick), "\n")
   }
 }
#>  state: Alive  get_sick: TRUE  pDie: 0.1  cost: 3000 
#>  state: Alive  get_sick: FALSE  pDie: 0.01  cost: 1000 
#>  state: Dead  get_sick: TRUE  pDie: 0  cost: 0 
#>  state: Dead  get_sick: FALSE  pDie: 0  cost: 0

Now we can examine the impact of this on the model results

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)
my_markov_model(model_struc,return_transition_prob=TRUE,return_state_payoffs=TRUE,return_trace=TRUE,return_cycle_payoffs=TRUE,return_payoff_summary=TRUE)
#> $P
#> , , StandardOfCare
#> 
#>        Alive   Dead
#> Alive 0.9855 0.0145
#> Dead  0.0000 1.0000
#> 
#> 
#> $state_payoffs
#> , , cost
#> 
#>       StandardOfCare
#> Alive           1100
#> Dead               0
#> 
#> 
#> $Trace
#> , , StandardOfCare
#> 
#>        Alive       Dead
#> 1  1.0000000 0.00000000
#> 2  0.9855000 0.01450000
#> 3  0.9712103 0.02878975
#> 4  0.9571277 0.04287230
#> 5  0.9432493 0.05675065
#> 6  0.9295722 0.07042777
#> 7  0.9160934 0.08390656
#> 8  0.9028101 0.09718992
#> 9  0.8897193 0.11028066
#> 10 0.8768184 0.12318159
#> 
#> 
#> $cycle_payoffs
#> , , StandardOfCare, cost
#> 
#>        Alive Dead
#> 1  1100.0000    0
#> 2  1084.0500    0
#> 3  1068.3313    0
#> 4  1052.8405    0
#> 5  1037.5743    0
#> 6  1022.5295    0
#> 7  1007.7028    0
#> 8   993.0911    0
#> 9   978.6913    0
#> 10  964.5002    0
#> 
#> 
#> $summary_payoffs
#>                    cost
#> StandardOfCare 10309.31

Running our model with other parameter values

gmod allows using variables inside the custom functions. These values must be defined in the global environment so that they are accessible within all the functions. To illustrate the use of parameters, we can begin by creating a new list (params) that contains 4 parameters as shown below. In the functions below we replace the values with the variable names we use in the params list. For example, we will use pDieSick in the params = 0.1, and replace the 0.1 we used in the pDie function with the variable name pDieSick. We do the same for the other variables.

params <- list(pDieSick = 0.1,
               pDieNotSick = 0.01,
               costAlive = 1000,
               costGetSick = 2000)

mygmod <- gmod() + 
  decisions("StandardOfCare") + 
  states(names = c("Alive", "Dead"), init_probs = c(1,0)) + 
  event(name = "get_sick",  
              scenarios = c(TRUE, FALSE), 
              probs = c(0.05, Inf), 
              outcomes = c("die", "die")) +  
    event(name = "die",  
              scenarios = c(TRUE, FALSE), 
              probs = c(pDie(state, get_sick), Inf), 
              outcomes = c("Dead", "curr_state")) +  
  payoffs(cost = compute_cost(state, get_sick))

pDie <- function(state, get_sick){
  # if the cohort gets sick, the probability of death = 0.01 otherwise it is 0.1
  switch(state, "Alive"= if(get_sick==TRUE) pDieSick else pDieNotSick,  
                "Dead"=0)
}

compute_cost <- function(state, get_sick){
  cost_state <- switch(state, "Alive"=costAlive, "Dead"=0)
  # add 2000 if someone in the Alive state gets sick
  cost_get_sick <- if(get_sick==TRUE & state=="Alive") costGetSick else 0
  cost_state+cost_get_sick
}

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)
my_markov_model(model_struc,params=params)
#> $summary_payoffs
#>                    cost
#> StandardOfCare 10309.31

gmod allows passing parameter values as a list. The values passed in the list replace the parameter values in the global environment.

Rerunnign the model with a different set of parameter values

Now, we can simply change the values of the parameter list and rerun the model by simply calling it with a different parameter values. For example, we change the value of pDieSick to 0.3 from 0.1, and change the costGetSick to 10000 from 2000.

params <- list(pDieSick = 0.3, # was 0.1
               pDieNotSick = 0.01,
               costAlive = 1000, 
               costGetSick = 10000) # was 2000
my_markov_model(model_struc,params=params)
#> $summary_payoffs
#>                   cost
#> StandardOfCare 13449.8

Being able to pass parameters as a list can be helpful when we conduct analyses using the gmod structure. For example, if we want to run a probabilistic sensitivity analysis, deterministic analyses, do value of information, or calibrate our model.

Summary:

In this Vignette we went over the basic features of the gmod library for building and running Markov models. We illustrated these features one at a time, and show how to define decisions, states, events, payoffs. In addition, we show how to change the number of cycles, the discounting rates. Furthermore, we show how to make transition probabilities and payoffs a function of decisions, states cycles, number of cycles in a state (tunnel) and previous events.

For further details on gmod, please refer to the other vignettes, mainly the introductory and advanced Markov model tutorial vignettes that replicate our anlayses published in Medical Decision Making.