Transition Rate Matrix

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

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

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

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

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.

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

• For now, just include zeros.

Non-Markovian Accumulators

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

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

• The transition matrix now includes a Markovian submatrix (green) and a non-markovian accumulator component (purple).

expm(R)  

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!

How Many with CVD After 1 Years?

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

Transition States

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

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

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

• 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?

• Note how there are some deaths in the first cycle, despite everyone starting from healthy!

Incorporating a one-time 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

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

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

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

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")]

Tunnel State

# It is possible to exit the tunnel state to death. 
P["tunCVDy1","Dead"] = P["CVD","Dead"]

Tunnel State

# It is possible to exit the tunnel state to death. 
P["tunCVDy1","Dead"] = P["CVD","Dead"]

Tunnel State

P["tunCVDy1","tunCVDy1"] = 0
P["tunCVDy1","tunCVDy2"] = 1 - P["CVD","Dead"]

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"]