Competing Risks

Rates, Probabilities, and Embedding

Learning Objectives

  1. Pitfalls of traditional transition probability conversion formulas.
  2. Discuss compound transitions and how to account for them.
  3. Demonstrate how to properly embed a transition probability matrix into a timestep.

Part 1: A World Without Death

  • We will first construct a very simple CVD model with only two states: Healthy and CVD.
  • Transitions are goverened by a single rate from Healthy CVD

  • We will first construct a very simple CVD model with only two states: Healthy and CVD.
  • Transitions are goverened by a single rate from Healthy CVD
G Healthy Healthy Healthy->Healthy CVD CVD Healthy->CVD r_H_CVD = 0.15 Dead Dead Healthy->Dead r_HD=0.006 CVD->CVD CVD->Dead hr_S * r_HD hr_S=10 Dead->Dead

Parameterize the model

params = 
  list(
    t_names = c("natural_history"),   # Strategy names. 
    n_treatments = 1,                 # Number of treatments
    s_names  = c("Healthy", "CVD"),   # State names
    n_states = 2,                     # 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
  )

Transition Probability Matrix

We’ll next construct a transition probability matrix using the standard rate-to-probability conversion formula:

\[ p = 1 - \exp(-rt) \] where \(r\) is the rate, \(t\) is the time-step (e.g., 1 for annual, 1/12 for monthly, etc.).

Transition Probability Matrix

params$mP <- 
  with(params,{
    tp_H_CVD = 1 - exp(-r_H_CVD * cycle) 
    array(data = c(1 - tp_H_CVD, 0, 
                 tp_H_CVD,1),
          dim = c(n_states, n_states, n_treatments),
                  dimnames = list(from = s_names,
                                  to = s_names,
                                  t_names))
  })

Transition Probability Matrix

params$mP <- 
  with(params,{
    tp_H_CVD = 1 - exp(-r_H_CVD * cycle) 
    array(data = c(1 - tp_H_CVD, 0, 
                 tp_H_CVD,1),
          dim = c(n_states, n_states, n_treatments),
                  dimnames = list(from = s_names,
                                  to = s_names,
                                  t_names))
  })
  • Next, place these probabilities within the transition probability matrix array.
  • Outer array dimensions goverened by the number of strategies under consideration.

Transition Probability Matrix

params$mP <- 
  with(params,{
    tp_H_CVD = 1 - exp(-r_H_CVD * cycle) 
    array(data = c(1 - tp_H_CVD, 0, 
                 tp_H_CVD,1),
          dim = c(n_states, n_states, n_treatments),
                  dimnames = list(from = s_names,
                                  to = s_names,
                                  t_names))
  })
  • Remaining code simply names the various dimensions of the array.

Transition Probability Matrix

$natural_history
         to
from      Healthy     CVD
  Healthy 0.86071 0.13929
  CVD     0.00000 1.00000

Construct the Markov Trace

sim_cohort <- function(params) {
  tr_ <- t(c("Healthy" = params$n_cohort, "CVD" = 0)) 
  trace <- 
    0:params$n_cycles %>% map_df(~({ 
        apply(params$mP,3, simplify=FALSE, function(tp) {  
          (tr_ %*% (tp %^% .x)) %>% data.frame()
        })
    })) %>% 
    as.matrix() 
}
  • sim_cohort(params) is a function that constructs a Markov trace from the supplied parameters.

Construct the Markov Trace

sim_cohort <- function(params) {
  tr_ <- t(c("Healthy" = params$n_cohort, "CVD" = 0)) 
  trace <- 
    0:params$n_cycles %>% map_df(~({  
        apply(params$mP,3, simplify=FALSE, function(tp) {  
          (tr_ %*% (tp %^% .x)) %>% data.frame()
        })
    })) %>% 
    as.matrix() 
}
  • Specify the initial state occupancy.

Construct the Markov Trace

sim_cohort <- function(params) {
  tr_ <- t(c("Healthy" = params$n_cohort, "CVD" = 0)) 
  trace <- 
    0:params$n_cycles %>% map_df(~({  
        apply(params$mP,3, simplify=FALSE, function(tp) {  
          (tr_ %*% (tp %^% .x)) %>% data.frame()
        })
    })) %>% 
    as.matrix() 
}
  • Iterate over the cycles in the model.

Construct the Markov Trace

sim_cohort <- function(params) {
  tr_ <- t(c("Healthy" = params$n_cohort, "CVD" = 0)) 
  trace <- 
    0:params$n_cycles %>% map_df(~({  
        apply(params$mP,3, simplify=FALSE, function(tp) {  
          (tr_ %*% (tp %^% .x)) %>% data.frame()
        })
    })) %>% 
    as.matrix() 
}
  • Within a cycle, iterate over the strategies under consideration.

Construct the Markov Trace

The markov trace tells us the fraction of the population in each state in any cycle:

trace <- sim_cohort(params)
trace[1:5,]
     natural_history.Healthy natural_history.CVD
[1,]                 1.00000             0.00000
[2,]                 0.86071             0.13929
[3,]                 0.74082             0.25918
[4,]                 0.63763             0.36237
[5,]                 0.54881             0.45119

Checkpoint 1: CVD State Occupancy at Two Years

1. CVD State Occupancy at Two Years

  • Suppose we start off with a cohort of 1000 healthy individuals.
  • How many are in the CVD health state at year 2?

1. CVD State Occupancy at Two Years

  • Suppose we start off with a cohort of 1000 healthy individuals.
  • How many are in the CVD health state at year 2?
  • Pin this number to your brain for a few minutes!
Healthy CVD
740.82 259.18

Checkpoint 2: Life Expectancy

2. Life Expectancy

  • We have a model with 100 annual cycles and no death.
  • Let’s verify that life expectancy is 100.

2. Life Expectancy

  • We have a model with 100 annual cycles and no death.
  • Let’s verify that life expectancy is 100
1payoff_life_exp = c("Healthy" = 1, "CVD" = 1)
2life_exp_cycle = trace[-1,] %*% payoff_life_exp
3total_life_exp = sum(life_exp_cycle)
1
Life expectancy payoff is 1 for every state where someone remains alive.
2
Total life-years alive in each state is the markov trace multiplied by the payoff vector.
3
Total life expectancy is the sum of cycle values.

2. Life Expectancy

  • We have a model with 100 annual cycles and no death.
  • Let’s verify that life expectancy is 100
payoff_life_exp = c("Healthy" = 1, "CVD" = 1)
life_exp_cycle = trace[-1,] %*% payoff_life_exp
total_life_exp = sum(life_exp_cycle)
total_life_exp * params$cycle
[1] 100

Checkpoint 3: Changing Cycle Lengths

3. Daily Cycle Length

  • We will next switch to a daily cycle length and verify that we can recapitulate these numbers.

3. Daily Cycle Length

  • We will next switch to a daily cycle length and verify that we can recapitulate these numbers.
params_daily <- modifyList(params,list(
  cycle = 1/365,
  n_cycles = params$n_cycles * 365))

3. Daily Cycle Length

params_daily$mP <- 
  with(params_daily,{
    tp_H_CVD = 1 - exp(-r_H_CVD * cycle) 
    array(data = c(1 - tp_H_CVD, 0, 
                 tp_H_CVD,1),
          dim = c(n_states, n_states, n_treatments),
                  dimnames = list(from = s_names,
                                  to = s_names,
                                  t_names))
  }) 
$natural_history
         to
from      Healthy        CVD
  Healthy 0.99959 0.00041087
  CVD     0.00000 1.00000000

3. Daily Cycle Length

# Markov Trace (Daily Cycle)
trace_daily <- 
  sim_cohort(params_daily)

# State Occupancy at 2 Years
1000 * trace_daily[365 * 2 + 1,]
Healthy CVD
740.82 259.18

3. Daily Cycle Length

# Life Expectancy
payoff_life_exp = c("Healthy" = 1, "CVD" = 1) 
life_exp_cycle_daily = trace_daily[-1,] %*% payoff_life_exp  
total_life_exp_daily = sum(life_exp_cycle_daily) 
total_life_exp_daily * params_daily$cycle
[1] 100

Summary

  1. After 2 years, 259 individuals are in the CVD health state.
  2. Because there is no death, after 100 cycles life expectancy is 100 years.
  3. We verified that the above results hold exactly when we convert to a daily cycle length.

Part 2: Include Background Mortality

  • We next include a background mortality rate that is increased by a hazard ratio (10.0) among people with CVD.

  • We next include a background mortality rate that is increased by a hazard ratio (10.0) among people with CVD.
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 the model

params_mort = 
  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 Probability Matrix

params_mort$mP <- 
  with(params_mort,{
    
    tp_H_CVD = 1 - exp(-r_H_CVD * cycle) 
    tp_H_D = 1 - exp(-r_H_D * cycle)
    tp_CVD_D <- 1 - exp(-(hr_CVD * r_H_D) * cycle)
    
    array(data = c(1 - tp_H_CVD - tp_H_D, 0, 0, 
                 tp_H_CVD,1-tp_CVD_D, 0,
                 tp_H_D, tp_CVD_D, 1),
          dim = c(n_states, n_states, n_treatments),
                  dimnames = list(from = s_names,
                                  to = s_names,
                                  t_names))
  })

Transition Probability Matrix

$natural_history
         to
from      Healthy     CVD      Dead
  Healthy 0.85076 0.13929 0.0099502
  CVD     0.00000 0.90484 0.0951626
  Dead    0.00000 0.00000 1.0000000

Construct the Markov Trace

sim_cohort <- function(params) {
  tr_ <- c(1,rep(0,params$n_states-1))
  trace <- 
    0:params$n_cycles %>% map_df(~({ 
        apply(params$mP,3, simplify=FALSE, function(tp) {  
          (tr_ %*% (tp %^% .x)) %>% data.frame()
        })
    })) %>% 
    as.matrix() 
}

trace_mort <- 
  sim_cohort(params_mort)

Construct the Markov Trace

trace_mort[1:5,]
     natural_history.Healthy natural_history.CVD natural_history.Dead
[1,]                 1.00000             0.00000            0.0000000
[2,]                 0.85076             0.13929            0.0099502
[3,]                 0.72379             0.24454            0.0316707
[4,]                 0.61577             0.32209            0.0621437
[5,]                 0.52387             0.37721            0.0989213

1. CVD State Occupancy at 2 Years

Healthy CVD Dead
723.8 244.5 31.7

2. Life Expectancy

payoff_life_exp_mort = c("Healthy" = 1, "CVD" = 1, "Dead" = 0) 
life_exp_cycle_mort = trace_mort[-1,] %*% payoff_life_exp_mort 
total_life_exp_mort = sum(life_exp_cycle_mort)
total_life_exp_mort
[1] 15.507

3. Daily Cycle Length

params_mort_daily <- modifyList(params_mort,list(
  cycle = 1/365,
  n_cycles = params$n_cycles * 365))

3. Daily Cycle Length

params_mort_daily$mP <- 
  with(params_mort_daily,{
    tp_H_CVD = 1 - exp(-r_H_CVD * cycle) 
    tp_H_D = 1 - exp(-r_H_D * cycle)
    tp_CVD_D <- 1 - exp(-(hr_CVD * r_H_D) * cycle)
    
    array(data = c(1 - tp_H_CVD - tp_H_D, 0, 0, 
                 tp_H_CVD,1-tp_CVD_D, 0,
                 tp_H_D, tp_CVD_D, 1),
          dim = c(n_states, n_states, n_treatments),
                  dimnames = list(from = s_names,
                                  to = s_names,
                                  t_names))
  })

3. Daily Cycle Length

$natural_history
         to
from      Healthy        CVD        Dead
  Healthy 0.99956 0.00041087 0.000027397
  CVD     0.00000 0.99972606 0.000273935
  Dead    0.00000 0.00000000 1.000000000

3. Daily Cycle Length

# Markov Trace (Daily Cycle)
trace_mort_daily <- sim_cohort(params_mort_daily)
Healthy CVD Dead
726.1 231.5 42.4

3. Daily Cycle Length

# Life Expectancy
payoff_mort_life_exp = c("Healthy" = 1, "CVD" = 1, "Death" = 0) 
life_exp_mort_daily = trace_mort_daily[-1,] %*% payoff_mort_life_exp  
total_life_exp_mort_daily = sum(life_exp_mort_daily) 
total_life_exp_mort_daily * params_mort_daily$cycle
[1] 15.624

Signs of Trouble

In the presence of competing events (mortality), we obtain different results for different cycle lengths!

Signs of Trouble

In the presence of competing events (mortality), we obtain different results for different cycle lengths!

Model With No Mortality Model With Mortality
Yearly Cycle Daily Cycle Yearly Cycle Daily Cycle
1. CVD State Occupancy at Two Years 259.18 259.18 244.540 231.488
2. Life Expectancy 100.00 100.00 15.507 15.624

Signs of Trouble

In the presence of competing events (mortality), we obtain different results for different cycle lengths! Notice how the state occupancy in CVD differs in the right panel of the figure below.

What Went Wrong?

What Went Wrong?

  • To gain intuition for what is going on, let’s split out death for those with vs. without CVD:

What Went Wrong?

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

Parameterize the model

params_mort2 = 
  list(
    t_names = c("natural_history"),                       
    n_treatments = 1,                                     
    s_names  = c("Healthy", "CVD", "CVDDeath", "Dead"),   
    n_states = 4,                                         
    n_cohort = 1,                                         
    n_cycles = 100,                                        
    cycle = 1,                                            
    initial_age = 55,                                     
    r_H_CVD = 0.15,                                       
    hr_CVD = 10,                                          
    r_H_D = 0.01                                          
  )

params_mort2$mP <- 
  with(params_mort2,{
    
    tp_H_CVD = 1 - exp(-r_H_CVD * cycle) 
    tp_H_D = 1 - exp(-r_H_D * cycle)
    tp_CVD_D <- 1 - exp(-(hr_CVD * r_H_D) * cycle)
    
    array(data = c(1 - tp_H_CVD - tp_H_D, 0, 0,0,  
                 tp_H_CVD,1-tp_CVD_D, 0,0,
                 0,tp_CVD_D,1,0,
                 tp_H_D, 0,0, 1),
          dim = c(n_states, n_states, n_treatments),
                  dimnames = list(from = s_names,
                                  to = s_names,
                                  t_names))
  })

Spot the Problem (1yr Cycle)

$natural_history
          to
from       Healthy     CVD CVDDeath      Dead
  Healthy  0.85076 0.13929 0.000000 0.0099502
  CVD      0.00000 0.90484 0.095163 0.0000000
  CVDDeath 0.00000 0.00000 1.000000 0.0000000
  Dead     0.00000 0.00000 0.000000 1.0000000

Spot the Problem (1yr Cycle)

$natural_history
          to
from       Healthy     CVD CVDDeath      Dead
  Healthy  0.85076 0.13929 0.000000 0.0099502
  CVD      0.00000 0.90484 0.095163 0.0000000
  CVDDeath 0.00000 0.00000 1.000000 0.0000000
  Dead     0.00000 0.00000 0.000000 1.0000000

Recall that in this model, CVD has a very high death rate—meaning that many people who develop CVD will die very soon after.

Spot the Problem (1yr Cycle)

$natural_history
          to
from       Healthy     CVD CVDDeath      Dead
  Healthy  0.85076 0.13929 0.000000 0.0099502
  CVD      0.00000 0.90484 0.095163 0.0000000
  CVDDeath 0.00000 0.00000 1.000000 0.0000000
  Dead     0.00000 0.00000 0.000000 1.0000000

Recall that in this model, CVD has a very high death rate—meaning that many people who develop CVD will die very soon after.

With this knowledge in mind, can you spot something wrong with the above (annual) transition matrix?

What if we change to a 10-year cycle?

  • See the problem now?
$natural_history
          to
from       Healthy     CVD CVDDeath     Dead
  Healthy  0.12797 0.77687  0.00000 0.095163
  CVD      0.00000 0.36788  0.63212 0.000000
  CVDDeath 0.00000 0.00000  1.00000 0.000000
  Dead     0.00000 0.00000  0.00000 1.000000

Recall that in this model, CVD has a very high death rate—meaning that many people who develop CVD will die very soon after.

With this knowledge in mind, can you spot something wrong with the above (annual) transition matrix?

Discretizing a Continuous Process

  • Our constructed model rules out compound transitions within a cycle.
  • Compound transition: Healthy CVD CVD Death in a single cycle.
$natural_history
          to
from       Healthy     CVD CVDDeath     Dead
  Healthy  0.12797 0.77687  0.00000 0.095163
  CVD      0.00000 0.36788  0.63212 0.000000
  CVDDeath 0.00000 0.00000  1.00000 0.000000
  Dead     0.00000 0.00000  0.00000 1.000000

Discretizing a Continuous Process

  • To accurately capture the underlying continuous time process, there should be a transition probability defined for a seemingly impossible transition: Healthy CVD Death.
$natural_history
          to
from       Healthy     CVD CVDDeath     Dead
  Healthy  0.12797 0.77687  0.00000 0.095163
  CVD      0.00000 0.36788  0.63212 0.000000
  CVDDeath 0.00000 0.00000  1.00000 0.000000
  Dead     0.00000 0.00000  0.00000 1.000000

Discretizing a Continuous Process

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

Discretizing a Continuous Process

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

Why Don’t the Formulas Work?

  • When we use the standard rate-to-probability conversion formula (i.e., \(p_{H,CVD}= 1 - \exp(-r_{H,CVD} t)\) ) we obtain the marginal transition probability.
  • This marginal probability reflects the union of all possible transitions that start Healthy and then transition to and from the CVD state in the cycle.

Why Don’t the Formulas Work?

That is, \(p_{H,CVD}= 1 - \exp(-r_{H,CVD} t)\) captures the probability of following scenarios occurring in a cycle:

  1. Individual transitions from healthy to CVD.
  2. Individual transitions from healthy to CVD to CVD Death.

10 Year Cycle

$natural_history
          to
from       Healthy     CVD CVDDeath     Dead
  Healthy  0.12797 0.77687  0.00000 0.095163
  CVD      0.00000 0.36788  0.63212 0.000000
  CVDDeath 0.00000 0.00000  1.00000 0.000000
  Dead     0.00000 0.00000  0.00000 1.000000
  • This matrix “captures” all transitions from Healthy CVD in a single transition probability (0.7769).

10 Year Cycle

$natural_history
          to
from       Healthy     CVD CVDDeath     Dead
  Healthy  0.12797 0.77687  0.00000 0.095163
  CVD      0.00000 0.36788  0.63212 0.000000
  CVDDeath 0.00000 0.00000  1.00000 0.000000
  Dead     0.00000 0.00000  0.00000 1.000000
  • But we actually need to distribute out the marginal probability so that we’re left only with the probability that someone transitions to CVD and remains there at the end of the cycle.
  • As constructed, we have ruled out the possibility of compound transitions (Healthy CVD CVDDeath) in a single cycle.

10 Year Cycle

$natural_history
         to
from      Healthy     CVD     Dead
  Healthy 0.12797 0.77687 0.095163
  CVD     0.00000 0.36788 0.632121
  Dead    0.00000 0.00000 1.000000
  • This problem isn’t fixed in the version of the model where death is a single health state.
  • The problem persists, but is “hidden” from view!

10 Year Cycle

$natural_history
          to
from       Healthy     CVD CVDDeath     Dead
  Healthy  0.12797 0.77687  0.00000 0.095163
  CVD      0.00000 0.36788  0.63212 0.000000
  CVDDeath 0.00000 0.00000  1.00000 0.000000
  Dead     0.00000 0.00000  0.00000 1.000000
  • This results in overcounting CVD state occupancy in some cases (people who develop CVD and die soon after) and undercounting CVD state occupancy in others (people who die of background causes before they would have developed CVD in the same cycle).

Discretizing a Continuous Process

Changing to shorter cycle does not get us out of this problem.

Quantified Absolute Error

  • Error after one cycle

Quantified Relative Error

  • For rates < 0.1 of the chosen time step the error is <5%.
  • This error will accumulate over each timestep.

Comparisons Across Strategies Bail You Out

  • Common CEA approaches compare strategies, i.e., \(f(\theta_A) - f(\theta_{B})\), \(\frac{f(\theta_A)}{ f(\theta_{B})}\)
    • Net monetary benefit.
    • Net health benefit.
    • Incremental cost effectiveness ratio.
  • With error it becomes (for NHB or NMB):

Comparisons Across Strategies Bail You Out

  • Common CEA approaches compare strategies, i.e., \(f(\theta_A) - f(\theta_{B})\), \(\frac{f(\theta_A)}{ f(\theta_{B})}\)
    • Net monetary benefit.
    • Net health benefit.
    • Incremental cost effectiveness ratio.
  • With error it becomes (for NHB or NMB):

\[f(\theta_A) + \epsilon_A - (f(\theta_{B}) + \epsilon_{B})\]

Specification Error

If \(\epsilon_A - \epsilon_{B} \sim 0\) the impact to the decision threshold is minimal.

  • CEA outcomes are thus semi-robust against this misspecification.
  • Restated: Single model runs are biased under misspecification, but comparison across strategies helps cancel this out.

Rate Misspecification Conclusions

In the presence of competing events (mortality), using standard conversion formulas will yield incorrect transition probability matrices.

Why Does This Matter?

  • Without properly embedding, you are no longer modeling a continuous time process.
  • So your model results will differ compared to a continuous time model (e.g., discrete event simulation)
  • Microsimulation suffers the same problem because the patient-level transition probability matrix is embedded in a time step and requires Markov Rate conversion.

Rate Misspecification Conclusions

  • Existing publications’ CEAs are generally okay if they misspecified rates.
  • Check the maximum rate of an event as a rough estimate of potential bias.
  • Rule of thumb, on the error plot use the estimated rate for full time scale of simulation, i.e. one step is entire simulation.

When Might it Matter?

  • Strategy compares strategies with very different rates of acute events or complications of treatment immediately after diagnosis (i.e., lots of compound transitions in a cycle!)
  • Particularly impactful in populations with high underlying death rates.

Rate Misspecification Error Decision Threshold

The amount of error for a decision threshold is the difference in contour. In this example: \(0.07 - 0.02 = 0.05\).

Additional Notes

  • Discrete Event Simulation automatically handles competing events.

  • Microsimulation misspecification error interacts with the discrete time step truncation and increases this error!

  • Published rates from clinical studies are usually adjusted for competing events (the purpose of survival studies).

So How Do We Obtain the Correct Transition Matrix?

  • Rather than use the conversion formulas, we need to properly “embed” the transition probability matrix.

  • The correct transition probability matrix should include a seemingly impossible transition: Healthy CVDDeath

  • Rather than use the conversion formulas, we need to properly “embed” the transition probability matrix.

  • The correct transition probability matrix should include a seemingly impossible transition: Healthy CVDDeath

         Healthy   CVD CVDDeath Dead
Healthy    0.202 0.415    0.333 0.05
CVD        0.000 0.368    0.632 0.00
CVDDeath   0.000 0.000    1.000 0.00
Dead       0.000 0.000    0.000 1.00

What Is Commonly Taught

  1. Convert each rate to a probability for the time step using \(p = 1-\exp(-rt)\).
  2. Place the probabilities within a transition matrix \(\mathbf{P}\).
  3. Diagonal elements of \(\mathbf{P}\): 1 - (row sum of nondiagonal elements).

What We’ll Teach You

  1. Place each rate in a rate matrix \(\mathbf{R}\).
  2. Diagonal elements of \(\mathbf{R}\): -1 * (row sum of nondiagonal elements).
  3. Embed the transition probability matrix using matrix exponentiation \(\mathbf{P} = \exp{(\mathbf{R}t)}\).

What We’ll Teach You

  1. Place each rate in a rate matrix \(\mathbf{R}\).
  2. Diagonal elements of \(\mathbf{R}\): -1 * (row sum of nondiagonal elements).
  3. Embed the transition probability matrix using matrix exponentiation \(\mathbf{P} = \exp{(\mathbf{R}t})\).

This process will ensure that the matrix captures “jumpover” states for compound transitions!

What We’ll Teach You

  1. Place each rate in a rate matrix \(\mathbf{R}\).
  2. Diagonal elements of \(\mathbf{R}\): -1 * (row sum of nondiagonal elements).
  3. Embed the transition probability matrix using matrix exponentiation \(\mathbf{P} = \exp{(\mathbf{R}t})\).

This process will ensure that the matrix captures “jumpover” states for compound transitions!

It also ensures that changing cycle lengths will not change your answers!

Steps to Embed a Transition Matrix

Parameterize the model

params_mort_embedded = 
  list(
    t_names = c("natural_history"),                       
    n_treatments = 1,                                     
    s_names  = c("Healthy", "CVD", "CVDDeath", "Dead"),   
    n_states = 4,                                         
    n_cohort = 1,                                         
    n_cycles = 100,                                        
    cycle = 1,                                            
    initial_age = 55,                                     
    r_H_CVD = 0.15,                                       
    hr_CVD = 10,                                          
    r_H_D = 0.01                                          
  )
  • Note: this is the same parameter list as earlier!

Transition Rate Matrix

params_mort_embedded$mR <- 
  with(params_mort2_emb,{

    r_CVD_D <- hr_CVD * r_H_D
    
    R_ <- 
      array(data = c(0, 0, 0,0,  
                 r_H_CVD,0, 0,0,
                 0,r_CVD_D,0,0,
                 r_H_D, 0,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_mort_embedded$mR <- 
  with(params_mort2_emb,{

    r_CVD_D <- hr_CVD * r_H_D
    
    R_ <- 
      array(data = c(0, 0, 0,0,  
                 r_H_CVD,0, 0,0,
                 0,r_CVD_D,0,0,
                 r_H_D, 0,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    
  })
  • Calculate the CVD Dead rate

Transition Rate Matrix

params_mort_embedded$mR <- 
  with(params_mort2_emb,{

    r_CVD_D <- hr_CVD * r_H_D
    
    R_ <- 
      array(data = c(0, 0, 0,0,  
                 r_H_CVD,0, 0,0,
                 0,r_CVD_D,0,0,
                 r_H_D, 0,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    
  })
  • Place the rates within the matrix. We’ll get the diagonal elements in a second to ensure things balance out.

Transition Rate Matrix

params_mort_embedded$mR <- 
  with(params_mort2_emb,{
    r_CVD_D <- hr_CVD * r_H_D
    R_ <- 
      array(data = c(0, 0, 0,0,  
                 r_H_CVD,0, 0,0,
                 0,r_CVD_D,0,0,
                 r_H_D, 0,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    
  })
  • Diagonal elements are the negative row sum (i.e., each row of \(\mathbf{R}\) sums to 0)

Transition Rate Matrix

params_mort_embedded$mR <- 
  with(params_mort2_emb,{
    r_CVD_D <- hr_CVD * r_H_D
    R_ <- 
      array(data = c(0, 0, 0,0,  
                 r_H_CVD,0, 0,0,
                 0,r_CVD_D,0,0,
                 r_H_D, 0,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    
  })
  • Convert \(\mathbf{R}\) to the model cycle length (e.g., 1 year, 1 day, 1 month, etc.).

Transition Rate Matrix

params_mort_embedded$mR
$natural_history
          to
from       Healthy   CVD CVDDeath Dead
  Healthy    -0.16  0.15      0.0 0.01
  CVD         0.00 -0.10      0.1 0.00
  CVDDeath    0.00  0.00      0.0 0.00
  Dead        0.00  0.00      0.0 0.00

Transition Probability Matrix

expm(params_mort_embedded$mR[["natural_history"]])
         Healthy     CVD  CVDDeath     Dead
Healthy  0.85214 0.13173 0.0068811 0.009241
CVD      0.00000 0.90484 0.0951626 0.000000
CVDDeath 0.00000 0.00000 1.0000000 0.000000
Dead     0.00000 0.00000 0.0000000 1.000000
  • Notice how the embedded matrix now has a transition probability from Healthy CVD Death!

Using Conversion Formulas:

          to
from       Healthy   CVD CVDDeath Dead
  Healthy    0.851 0.139    0.000 0.01
  CVD        0.000 0.905    0.095 0.00
  CVDDeath   0.000 0.000    1.000 0.00
  Dead       0.000 0.000    0.000 1.00

Using Rate Matrix Eponentiation:

         Healthy   CVD CVDDeath  Dead
Healthy    0.852 0.132    0.007 0.009
CVD        0.000 0.905    0.095 0.000
CVDDeath   0.000 0.000    1.000 0.000
Dead       0.000 0.000    0.000 1.000

Your Turn!

Case Study Overview

Learning Objectives

  1. Construct a natural history model for CVD.
  2. Background mortality based on cause-deleted life table data.
  3. CVD mortality based on cause-specific mortality.
  4. Markov model is properly embedded.
  5. Markov model results will closely match cause-deleted and cause-specific mortality in the US population.

Back to Website