Advanced Bookkeeping

Non-Markovian Accumulators, Transition States, and Backwards Conversion

Learning Objectives

  • Establish the transition rate matrix as the central “hub” of the modeling process.
  • Include transition states and accumulators to accurately track and count events, costs, and QALYs.
  • Demonstrate how to include tunnel states.

Review of Key Concepts

Transition Rate Matrix

  • The central “hub” of a Markov model.
  • Straightforward to convert rate matrix into a transition probability matrix.
  • Facilitates modeling using alternative techniques:
    • Continuous time Markov
    • Discrete event simulation

Transition Rate Matrix

  • The “workshop” where you can include additional health states, new evidence, change cycle lengths, etc.
  • Ensures that results obtained via discrete time Markov model will match those under an identical model but different method (e.g., DES).

What they didn’t teach you

  • The transition rate matrix is also essential for accurate and transparent accounting of costs and QALYs.

What we’ll teach you

  • We’ll demonstrate how to augment the matrix with non-Markovian elements using the basic CVD model from the last session.
  • We’ll show you how to backwards convert an existing transition probability matrix into a continuous generator (rate) matrix.
  • We’ll then replicate and extend a recent didactic model (Green et al. 2023) in an applied case study using these methods.

CVD Model

CVD Model

  • Let’s quickly build up the basic CVD model again.

CVD Model

G Healthy Healthy Healthy->Healthy CVD CVD Healthy->CVD r_H_CVD = 0.15 Dead Dead Healthy->Dead r_HD=0.01 CVD->CVD CVD->Dead hr_CVD * r_HD hr_CVD=10 Dead->Dead

Parameterize

params = 
  list(
    t_names = c("natural_history"),           # Strategy names. 
    n_treatments = 1,                         # Number of treatments
    s_names  = c("Healthy", "CVD", "Dead"),   # State names
    n_states = 3,                             # Number of states
    n_cohort = 1,                             # Cohort size
    n_cycles = 100,                           # Number of cycles in model.  
    cycle = 1,                                # Cycle length
    initial_age = 55,                         # Cohort starting age
    r_H_CVD = 0.15,                           # Rate of healthy -> CVD
    hr_CVD = 10,                              # Hazard Ratio: CVD Death
    r_H_D = 0.01                              # Rate of healthy -> dead
  )

Transition Rate Matrix

params$mR <- 
  with(params,{
    r_CVD_D <- hr_CVD * r_H_D
    R_ <- 
      array(data = c(0, 0, 0,  
                 r_H_CVD,0, 0,
                 r_H_D,r_H_D + r_CVD_D,0,
                 r_H_D, 0,0),
          dim = c(n_states, n_states, n_treatments),
                  dimnames = list(from = s_names,
                                  to = s_names,
                                  t_names))
    R <- apply(R_,3, simplify=FALSE,function(x) {
      diag(x) = -rowSums(x)
      x * cycle
    })
  R    
  })

Transition Rate Matrix

params$mR
$natural_history
         to
from      Healthy   CVD Dead
  Healthy   -0.16  0.15 0.01
  CVD        0.00 -0.11 0.11
  Dead       0.00  0.00 0.00

Transition Probability Matrix

expm(params$mR[["natural_history"]])
        Healthy     CVD     Dead
Healthy 0.85214 0.13107 0.016785
CVD     0.00000 0.89583 0.104166
Dead    0.00000 0.00000 1.000000

Advanced Bookkeeping

  1. Given compound transitions, how can we track the total number of people who develop CVD?
  2. What if there is a cost associated with a single transition type (e.g., CVD CVD death)?
  3. What if we want to include a tunnel state to capture transient utility changes after contracting CVD or experiencing an acute event?

Advanced Bookkeeping

Augmenting the transition matrix with non-markovian elements can address each question:

  1. Accumulator states to count total number of people who enter a health state.
  2. Transition states to count the number of people who transition into a state in a given cycle.
  3. Tunnel states to capture transient events (e.g., 2 cycle utility decrement, etc.)

Markovian Transition Matrices

  • Traditional Markov modeling approaches restrict to a transition probability matrix.
  • Row sums all equal 1.0.

Markovian Transition Matrices

  • Traditional Markov modeling approaches restrict to a transition probability matrix.
  • Row sums all equal 1.0.
Healthy CVD Dead
Healthy 0.852 0.131 0.017
CVD 0.000 0.896 0.104
Dead 0.000 0.000 1.000

Non-Markovian Elements

  • Augmenting the transition matrix with non-markovian rows and columns can facilitate accurate counting and bookkeeping.
  • Depending on the objective, can include non-markovian accumulators and/or transition states.
  • Tunnel states are also non-Markovian, but occupy a middle ground between transition states and accumulators.

Non-Markovian Elements

  • Accumulator: tracks the total number of individuals who have entered a given state up until a given cycle (even if they moved out of the state later).

Non-Markovian Elements

  • Transition state: tracks the total number of individuals who enter a given state in a given cycle.
    • Can construct as a tunnel state if you move them through to the next state (or exit for some other reason like background mortality).
    • But you don’t have to do this. Depending on objective, you could just count the total number who enter the state to count up transitory costs, etc.

Where we’re headed

Healthy CVD Dead trCVD accCVD
Healthy 0.852 0.131 0.017 . .
CVD 0 0.896 0.104 . .
Dead 0 0 1 . .
trCVD . . . . .
accCVD . . . . .

Augmenting the Transition Matrix

Augmenting the Transition Matrix

  • As noted earlier, all roads begin from the transition rate matrix.
  • The way in which we construct the Markovian components is the same.
    • Place rates as appropriate.
    • Diagonal of Markovian elements is the negative sum of the other row elements.

Augmenting the Transition Matrix

  • Rate matrix for basic CVD model
Healthy CVD Dead
Healthy -0.16 0.15 0.01
CVD 0.00 -0.11 0.11
Dead 0.00 0.00 0.00

Total CVD Cases

  • Suppose we wanted to track the total number of people who develop CVD over the simulated time horizon.
  • We can’t do this accurately using CVD state occupancy in the Markov trace due to compound transitions.

Total CVD Cases

  • We need some way to track everyone who moves into the CVD state.
  • We can do this by adding a non-markovian accumulator to the rate matrix.

1. Accumulators

Non-Markovian Accumulators

We’ll start with the original transition rate matrix:

R_ <- params$mR[["natural_history"]]
R_
         to
from      Healthy   CVD Dead
  Healthy   -0.16  0.15 0.01
  CVD        0.00 -0.11 0.11
  Dead       0.00  0.00 0.00

Non-Markovian Accumulators

  • Next, we will add both a column and a row for the accumulator that tracks movement into the CVD health state.
  • For now, just include zeros.

Non-Markovian Accumulators

  • Next, we will add both a column and a row for the accumulator that tracks movement into the CVD health state.

  • For now, just include zeros.

R <- 
  cbind(R_,"accCVD" = c(0,0,0)) %>% 
  rbind(.,"accCVD" = c(0,0,0,0))  
R
        Healthy   CVD Dead accCVD
Healthy   -0.16  0.15 0.01      0
CVD        0.00 -0.11 0.11      0
Dead       0.00  0.00 0.00      0
accCVD     0.00  0.00 0.00      0

Non-Markovian Accumulators

  • Next fill in the appropriate transition rate.
R["Healthy","accCVD"] = params$r_H_CVD
R
        Healthy   CVD Dead accCVD
Healthy   -0.16  0.15 0.01   0.15
CVD        0.00 -0.11 0.11   0.00
Dead       0.00  0.00 0.00   0.00
accCVD     0.00  0.00 0.00   0.00

Non-Markovian Accumulators

  • Embed using matrix exponentiation to obtain the transition probability matrix.
  • The transition matrix now includes a Markovian submatrix (green) and a non-markovian accumulator component (purple).
expm(R)
Healthy CVD Dead accCVD
Healthy 0.852 0.131 0.017 0.139
CVD 0 0.896 0.104 0
Dead 0 0 1 0
accCVD 0 0 0 1

Non-Markovian Accumulators

  • Note that the Markovian submatrix still has rows that sum to one.
  • Also note that the accumulator captures the same transition probability that was (incorrectly) placed in the Healthy CVD cell in lecture 2!
    • Recall that this probability is the marginal probability of transition (i.e., all transitions into CVD state). So it’s exactly what we want to accumulate!

How Many with CVD After 1 Years?

  • If we simply went based on state occupancy after one cycle, we’d understate the number of people who develop CVD.
    • This is the number who remain in the CVD health state.
    • Some had a compound transition and died within the cycle after they developed CVD!
cycle Healthy CVD Dead accCVD
0 100000 0 0.0 0
1 85214 13107 1678.5 13862

How Many with CVD After 1 Years?

  • The accumulator column tells us how many total people developed CVD.
cycle Healthy CVD Dead accCVD
0 100000 0 0.0 0
1 85214 13107 1678.5 13862

2. Transition States

Transition States

  • What if we need to count how many people transition to a given state in a cycle?
  • For example, we may want to apply a one-time cost to a transition such as CVD CVD death.
  • This type of transition is not immediately captured in the Markov trace.
    • Observed cycle-to-cycle change in death counts capture CVD CVD Death as well as death from background causes.

Transition States

  • Transition states are straightforward to capture as another element in the transition matrix.
  • Similar to an accumulator state, but we require state exit after one cycle.
  • Amounts so simply zeroing out the accCVD accCVD transition probability.

Transition States

We’ll start with the original transition rate matrix:

R_ <- params$mR[["natural_history"]]
R_
         to
from      Healthy   CVD Dead
  Healthy   -0.16  0.15 0.01
  CVD        0.00 -0.11 0.11
  Dead       0.00  0.00 0.00

Transition States

  • Next, we will add both a column and a row for the transition state that tracks movement into the CVD Death state.

Transition States

  • Next, we will add both a column and a row for the transition state that tracks movement into the CVD Death state.
  • For now, just include zeros.
R <- 
  cbind(R_,"trCVDDeath" = c(0,0,0)) %>% 
  rbind(.,"trCVDDeath" = c(0,0,0,0))  
R
           Healthy   CVD Dead trCVDDeath
Healthy      -0.16  0.15 0.01          0
CVD           0.00 -0.11 0.11          0
Dead          0.00  0.00 0.00          0
trCVDDeath    0.00  0.00 0.00          0

Transition States

  • We then fill in the appropriate transition rate.
R["CVD","trCVDDeath"] = params$hr_CVD * params$r_H_D
R
           Healthy   CVD Dead trCVDDeath
Healthy      -0.16  0.15 0.01        0.0
CVD           0.00 -0.11 0.11        0.1
Dead          0.00  0.00 0.00        0.0
trCVDDeath    0.00  0.00 0.00        0.0

Transition States

  • Embed using matrix exponentiation to obtain the transition probability matrix.
  • Note that we need to manually zero out the transition probability that keeps people in the transition state.
P = expm(R)  
P["trCVDDeath","trCVDDeath"] = 0
P
           Healthy     CVD     Dead trCVDDeath
Healthy    0.85214 0.13107 0.016785  0.0068583
CVD        0.00000 0.89583 0.104166  0.0946962
Dead       0.00000 0.00000 1.000000  0.0000000
trCVDDeath 0.00000 0.00000 0.000000  0.0000000

Transition States

Healthy CVD Dead trCVDDeath
Healthy 0.852 0.131 0.017 0.007
CVD 0.000 0.896 0.104 0.095
Dead 0.000 0.000 1.000 0.000
trCVDDeath 0.000 0.000 0.000 0.000
  • Note that the Markovian submatrix still has rows that sum to one.
  • Also note that the transition state column captures transitions to death that occur from both the Healthy state (compound transitions) and the CVD state!

How Many CVD Death Transitions Are There?

  • Hypothetical cohort of 100,000.
  • Note how there are some deaths in the first cycle, despite everyone starting from healthy!
cycle Healthy CVD Dead trCVDDeath
0 100000 0 0.0 0.00
1 85214 13107 1678.5 685.83

Incorporating a one-time cost

  • Can define a cost “payoff” for trCVDDeath and apply as usual to calculate accumulated costs from a one-time CVD death cost.
# Define the cost vector
cost_payoff = c("Healthy" = 0, "CVD" = 500, "Dead" = 0, "trCVDDeath" = 2000) 

# Sum up total costs in each cycle
costs_cycle = trace[["natural_history"]][-1,] %*% cost_payoff

# Function for cycle correction (alternative Simpson's method)
alt_simp_coef <- function(i) c(17, 59, 43, 49, rep(48, i-8), 49, 43, 59, 17) / 48
cycle_adj      <- function(x) sum(alt_simp_coef(length(x)) * x)

total_costs = cycle_adj(costs_cycle) 
total_costs
[1] 5926.6

3. Tunnel States

Tunnel States

  • Tunnel states are structured similarly to accumulators and transition states.
  • One additional complication: need to allow for tunnel state exit to background mortality (or other exit states).

Accumulators and Transition States

  • Occupy both ends of a spectrum.
  • Accumulators: probability of remaining in accumulator state = 1.0.
  • Transition States: probability of remaining in transition state = 0.0.

Tunnel States

  • Probability of remaining in tunnel: 0 < p < 1.
  • Reflects the idea that most will remain in the tunnel for the cycle, but some might leave due to secular death, transition back to a less severe disease state (e.g., remission), etc.

Tunnel States

  • Because we hold individuals “captive” in the tunnel state (unless they exit to background death), tunnel states are non-Markovian.
  • If we built them into the rate matrix, then embedded that into a transition probability matrix, we’d have compound transitions through the tunnel!
  • We therefore must define tunnel states in the probability space, not the rate space.

Tunnel States

  • Suppose we want to record a temporary 2 year (cycle) initial treatment cost once someone develops CVD.
  • We can set up a tunnel state to capture this.

Tunnel State

We’ll start with the original transition rate matrix:

R_ <- params$mR[["natural_history"]]
R_
         to
from      Healthy   CVD Dead
  Healthy   -0.16  0.15 0.01
  CVD        0.00 -0.11 0.11
  Dead       0.00  0.00 0.00

Tunnel State

  • Next, we will add columns and rows for the tunnel states.
  • Given that our model has an annual cycle, and the tunnel lasts two years, we need to add two columns and two rows.
  • If we had a monthly time cycle, we’d need to add 24 …

Tunnel States

  • For now, just include zeros.
R <- 
    cbind(R_,"tunCVDy1" = c(0,0,0), "tunCVDy2" = c(0,0,0)) %>% 
    rbind(.,"tunCVDy1" = c(0,0,0,0,0,0), "tunCVDy2" = c(0,0,0,0,0,0))  
R
         Healthy   CVD Dead tunCVDy1 tunCVDy2
Healthy    -0.16  0.15 0.01        0        0
CVD         0.00 -0.11 0.11        0        0
Dead        0.00  0.00 0.00        0        0
tunCVDy1    0.00  0.00 0.00        0        0
tunCVDy2    0.00  0.00 0.00        0        0

Tunnel State

  • We then fill in the appropriate transition rate for entry into the tunnel.
  • We zero out the transition from healthy to CVD since this can only occur after 2 years.
R["Healthy","tunCVDy1"] = R["Healthy","CVD"]
R["Healthy","CVD"] = 0
R
         Healthy   CVD Dead tunCVDy1 tunCVDy2
Healthy    -0.16  0.00 0.01     0.15        0
CVD         0.00 -0.11 0.11     0.00        0
Dead        0.00  0.00 0.00     0.00        0
tunCVDy1    0.00  0.00 0.00     0.00        0
tunCVDy2    0.00  0.00 0.00     0.00        0

Tunnel State

  • We now need to embed the transition probability matrix because tunnel states can only be defined on the probability scale.
  • Intuition: if we defined the tunnel state in the rate matrix, when we embed we’d end up with compound transitions that the tunnel state rules out (i.e., Healthy to tunCVDy1)

Tunnel State

  • We embed and rearrange the initial matrix to reflect the time sequencing through the tunnel state.
P_ = expm(R)  
P = P_[c("Healthy","tunCVDy1","tunCVDy2","CVD","Dead"),c("Healthy","tunCVDy1","tunCVDy2","CVD","Dead")]
Healthy tunCVDy1 tunCVDy2 CVD Dead
Healthy 0.852 0.139 0 0 0.009
tunCVDy1 0 1 0 0 0
tunCVDy2 0 0 1 0 0
CVD 0 0 0 0.896 0.104
Dead 0 0 0 0 1

Tunnel State

  • Unlike a transition state or accumulator, we need to allow for exit from the tunnel to death in the cycle.
# It is possible to exit the tunnel state to death. 
P["tunCVDy1","Dead"] = P["CVD","Dead"]

Tunnel State

  • Unlike a transition state or accumulator, we need to allow for exit from the tunnel to death in the cycle.
# It is possible to exit the tunnel state to death. 
P["tunCVDy1","Dead"] = P["CVD","Dead"]
Healthy tunCVDy1 tunCVDy2 CVD Dead
Healthy 0.852 0.139 0 0 0.009
tunCVDy1 0 1 0 0 0.104
tunCVDy2 0 0 1 0 0
CVD 0 0 0 0.896 0.104
Dead 0 0 0 0 1

Tunnel State

  • We also must force exit into the next tunnel state.
P["tunCVDy1","tunCVDy1"] = 0
P["tunCVDy1","tunCVDy2"] = 1 - P["CVD","Dead"]
Healthy tunCVDy1 tunCVDy2 CVD Dead
Healthy 0.852 0.139 0 0 0.009
tunCVDy1 0 0 0.896 0 0.104
tunCVDy2 0 0 1 0 0
CVD 0 0 0 0.896 0.104
Dead 0 0 0 0 1

Tunnel State

  • We repeat a similar exercise for the 2nd year in the tunnel.
    • Allow for exit to death.
    • Force entry into the CVD state after the cycle.

Tunnel State

P["tunCVDy2","Dead"] = P["CVD","Dead"]
P["tunCVDy2","tunCVDy2"] = 0
P["tunCVDy2","CVD"] = 1 - P["CVD","Dead"]
Healthy tunCVDy1 tunCVDy2 CVD Dead
Healthy 0.852 0.139 0 0 0.009
tunCVDy1 0 0 0.896 0 0.104
tunCVDy2 0 0 0 0.896 0.104
CVD 0 0 0 0.896 0.104
Dead 0 0 0 0 1

Your Turn!

Case Study Overview

G Asymptomatic Asymptomatic Asymptomatic->Asymptomatic Progressive Progressive Disease Asymptomatic->Progressive tpProg * (1-effect) Dead Dead Asymptomatic->Dead tpDn Progressive->Progressive Progressive->Dead   tpDcm+tpDn Dead->Dead

Source: Green et al. (2023)

Back to Website

Green, Nathan, Felicity Lamrock, Nichola Naylor, Jack Williams, and Andrew Briggs. 2023. “Health Economic Evaluation Using Markov Models in r for Microsoft Excel Users: A Tutorial.” PharmacoEconomics 41 (1): 5–19.