Competing Risks
Rates, Probabilities, and Embedding
Learning Objectives
- Pitfalls of traditional transition probability conversion formulas.
- Discuss compound transitions and how to account for them.
- 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
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
Transition Probability Matrix
- Next, place these probabilities within the transition probability matrix array.
- Outer array dimensions goverened by the number of strategies under consideration.
Transition Probability Matrix
- 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.00000Construct the Markov Trace
sim_cohort(params)is a function that constructs a Markov trace from the supplied parameters.
Construct the Markov Trace
- Specify the initial state occupancy.
Construct the Markov Trace
- Iterate over the cycles in the model.
Construct the Markov Trace
- 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:
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
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.
3. Daily Cycle Length
$natural_history
to
from Healthy CVD
Healthy 0.99959 0.00041087
CVD 0.00000 1.000000003. Daily Cycle Length
| Healthy | CVD |
|---|---|
| 740.82 | 259.18 |
3. Daily Cycle Length
Summary
- After 2 years, 259 individuals are in the CVD health state.
- Because there is no death, after 100 cycles life expectancy is 100 years.
- 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.
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.0000000Construct the Markov Trace
Construct the Markov Trace
1. CVD State Occupancy at 2 Years
| Healthy | CVD | Dead |
|---|---|---|
| 723.8 | 244.5 | 31.7 |
2. Life Expectancy
3. Daily Cycle Length
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.0000000003. Daily Cycle Length
| Healthy | CVD | Dead |
|---|---|---|
| 726.1 | 231.5 | 42.4 |
3. Daily Cycle Length
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?
Parameterize the model
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.0000000Spot 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.0000000Recall 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.0000000Recall 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.000000Recall 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.000000Discretizing 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.000000Discretizing a Continuous Process
Discretizing a Continuous Process
Why Don’t the Formulas Work?
- When we use the standard rate-to-probability conversion formula (i.e., pH,CVD=1−exp(−rH,CVDt) ) 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, pH,CVD=1−exp(−rH,CVDt) captures the probability of following scenarios occurring in a cycle:
- Individual transitions from healthy to CVD.
- 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(θA)−f(θB), f(θA)f(θ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(θA)−f(θB), f(θA)f(θB)
- Net monetary benefit.
- Net health benefit.
- Incremental cost effectiveness ratio.
- With error it becomes (for NHB or NMB):
f(θA)+ϵA−(f(θB)+ϵB)
Specification Error
If ϵA−ϵB∼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.00What Is Commonly Taught
- Convert each rate to a probability for the time step using p=1−exp(−rt).
- Place the probabilities within a transition matrix P.
- Diagonal elements of P: 1 - (row sum of nondiagonal elements).
What We’ll Teach You
- Place each rate in a rate matrix R.
- Diagonal elements of R: -1 * (row sum of nondiagonal elements).
- Embed the transition probability matrix using matrix exponentiation P=exp(Rt).
What We’ll Teach You
- Place each rate in a rate matrix R.
- Diagonal elements of R: -1 * (row sum of nondiagonal elements).
- Embed the transition probability matrix using matrix exponentiation P=exp(Rt).
This process will ensure that the matrix captures “jumpover” states for compound transitions!
What We’ll Teach You
- Place each rate in a rate matrix R.
- Diagonal elements of R: -1 * (row sum of nondiagonal elements).
- Embed the transition probability matrix using matrix exponentiation P=exp(Rt).
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
- 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 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 R to the model cycle length (e.g., 1 year, 1 day, 1 month, etc.).
Transition Rate Matrix
Transition Probability Matrix
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.00Using 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.000Your Turn!
Case Study Overview
Learning Objectives
- Construct a natural history model for CVD.
- Background mortality based on cause-deleted life table data.
- CVD mortality based on cause-specific mortality.
- Markov model is properly embedded.
- Markov model results will closely match cause-deleted and cause-specific mortality in the US population.
Competing Risks Rates, Probabilities, and Embedding