Goal

A Novel Decision Model for Cardiovascular Disease (CVD) Prevention and Treatment in South Africa

Goal

An Adapted Decision Model for Cardiovascular Disease (CVD) Prevention and Treatment in South Africa

What Never Happens

Solution

Solution: New Model Creation

Solution: New Model Creation

Solution: New Model Creation

Solution: Model Adaptation

Transition Rate Matrix

Transition Rate Matrix

Learning Objectives

Rates vs. Probabilities

Rates vs. Probabilities

Rates

Rates

Probabilities

Rates vs. Probabilities: Example¹

Rates vs. Probabilities: Example¹

Rates vs. Probabilities: Example¹

Rates vs. Probabilities: Example¹

Summary: Rates vs. Probabilities

What You Were Taught

• You were likely taught to convert each rate to a probability, then insert that probability into a transition matrix.
• This approach will yield the correct transition probabilities only if there is a single, non competing state you can transition to.

What We’ll Teach You

• We’ll teach you how to “embed” a transition probability matrix from a transition rate matrix.
• Essentially, fill out the transition matrix in terms of rates, then do a simple math operation to the entire matrix.
• This process is the matrix analogue to the usual rate-to-probability conversion formula.

From the rate matrix we can…

1. Swap in country-specific disease incidence and death rates, non-cause (“background”) mortality rates, etc.
2. Add non-markovian transition states to calculate years of life lost to disease (YLL), a key input into DALYs.
3. Add additional health states or strategies based on new evidence or country-specific context.
4. Change the model time step (e.g., convert from monthly to yearly cycles).

Embedding 101

• “Embedding” refers to the process of converting a transition rate matrix to a transition probability matrix.
• We embed a transition probability matrix to using a four-step process.

Embedding: A Four-Step Process

1. Place rates in a rate matrix.

1. Place rates in a rate matrix.

1. Place rates in a rate matrix.

1. Place rates in a rate matrix.

1. Place rates in a rate matrix.

1. Place rates in a rate matrix.

1. Place rates in a rate matrix.

1. Place rates in a rate matrix.

For example, we would use the rate matrix to “swap in” a different background mortality rate, or a different disease incidence rate from another country, region, or population of interest.

1. Place rates in a rate matrix

2. Make adjustments to the rate matrix as needed.

Suppose we want to:

2. Make adjustments to the rate matrix as needed.

• Think of the rate matrix as an automobile body shop: it’s where you go to make structural changes to your model.
• For now, we’ll assume you’re building a model from scratch.
• We’ll cover how to backwards-convert an existing model (based on a transition probability matrix) to a rate matrix later in this lecture.

2. Make adjustments to the rate matrix as needed.

• Suppose we want to model different strategies under consideration.
• We’ll do this by adapting a new rate matrix for each strategy.
• Note that not all adaptions can occur on the rate matrix; some may require adaptation of the transition probability matrix instead.

2. Make adjustments to the rate matrix as needed.

3. “Embed” the transition probability matrix using the rate matrix

3. “Embed” the transition probability matrix using the rate matrix

3. “Embed” the transition probability matrix using the rate matrix

3. “Embed” the transition probability matrix using the rate matrix

\[ p = 1 - \exp(-rt) \] where \(r\) is the rate and \(t\) is the time-step.

3. “Embed” the transition probability matrix using the rate matrix

3. “Embed” the transition probability matrix using the rate matrix

3. “Embed” the transition probability matrix using the rate matrix

\[ \mathbf{P} = e^{{\mathbf{R}t}} \]

3. “Embed” the transition probability matrix using the rate matrix

\[ \mathbf{P} = e^{{\mathbf{R}t}} \]

3. Calculate the Exponential of the Rate Matrix

\[ \mathbf{P} = e^{\mathbf{R}t} \]

3. Calculate the Exponential of the Rate Matrix

\[ \mathbf{P} = e^{\mathbf{R}t} \]

Pros:

3. Calculate the Exponential of the Rate Matrix

Cons:

3. Calculate the Exponential of the Rate Matrix

Cons:

library(expm)
P = expm(R)

3. Calculate the Exponential of the Rate Matrix

For our CVD example, the matrix exponential of the transition rate matrix yields the correct transition probability matrix:

3. “Embed” the transition probability matrix using the rate matrix

Why do we obtain this probability?

So Where Are We Now?

A Note of Caution

4. Make adjustments to \(\mathbf{P}\) as needed.

4. Make adjustments to \(\mathbf{P}\) as needed.

4. Make adjustments to \(\mathbf{P}\) as needed.

Full Process

Goal

An Adapted Decision Model for HIV Prevention and Treatment in South Africa

III. Backwards Conversion

• What if we have a model that is already defined in terms of a transition probability matrix?
• How can we convert to a new country or setting?
• As shown earlier, we need the transition rate matrix to swap in new background mortality, new strategies, country-specific incidence, etc.

III. Backwards Conversion

III. Backwards Conversion

• We will next show you several ways to work backwards.
• Boils down to solving for the continuous “generator matrix” for the observed transition probability matrix.
• We’ll show you one approach to do this, which can be done in Excel but is technically incorrect (it will hopefully get you close, however!)

Working Example: HIV Model

Transition probability matrix:

Three Methods to Solve for the Generator

A. Cellwise Conversion Using Rate-to-Probability Formula

A. Cellwise Conversion Using Rate-to-Probability Formula

\[ r = -ln(1-p) \]

A. Cellwise Conversion Using Rate-to-Probability Formula

\[ r = -ln(1-p) \]

A. Cellwise Conversion Using Rate-to-Probability Formula

Converted rate matrix:

Alternative (Advanced Approaches)

• The two additional approaches (B. matrix logarithm; and C. maximum likelihood approach) will get you closer to a more accurate rate matrix.
• These approaches are beyond the capabilities of Excel.
• You will need to import the transition probability matrix into a statistical software package (e.g., R) to do them.
• We provide the process and code here, but will not cover it in person.

A Word of Warning

• As you’ll see, the next two approaches are not always going to work perfectly. Or they may require you to impose some structural assumptions on your model to make them work.
• If the original transition probability matrix was defined incorrectly (e.g., no jumpover probabilities), the generator matrix may not exist.
• Identifiability is a more general issue, however.

B. Matrix Logarithm of P

• In mathematical terms, the generator matrix is the matrix logarithm of the transition probability matrix.
• A matrix has a logarithm if and only if if it is invertible.

B. Matrix Logarithm of P

\[ R = \log P \]

B. Matrix Logarithm of P

• The \(\log\) can be found using spectral or eigenvalue decomposition.
• This is beyond the capabilities of Excel; what follows is code for the open-source R language.
• If \(V\) is a matrix where each column is an eigenvector of \(P\), then

B. Matrix Logarithm of P

\[A' = V^{-1} A V\]

B. Matrix Logarithm of P

\[A' = V^{-1} A V\]

\[\log P = V (\log A') V^{-1}\]

B. Matrix Logarithm of P

# read the transition probability matrix as stored in an
# excel file 
mP <- readxl::read_xlsx("path-excel-file.xls")

# eigenvalue decomposition 
V  <- eigen(mP)$vectors
iV <- solve(V)
Ap <- iV %*% mP %*% V
lAp <- diag(log(diag(Ap)), nrow(Ap), ncol(Ap))
R  <- V %*% lAp %*% iV
R[abs(R) < 1e-6 ] <- 0
dimnames(R) = dimnames(mP)
round(R,3)

B. Matrix Logarithm of P

V  <- eigen(mP)$vectors
iV <- solve(V)
Ap <- iV %*% mP %*% V
lAp <- diag(log(diag(Ap)), nrow(Ap), ncol(Ap))
R  <- V %*% lAp %*% iV
R[abs(R) < 1e-6 ] <- 0
dimnames(R) = dimnames(mP)
round(R,3)

        Healthy LowCD4   AIDS  Death
Healthy  -0.327  0.311  0.002  0.013
LowCD4    0.000 -0.543  0.615 -0.072
AIDS      0.000  0.000 -0.288  0.288
Death     0.000  0.000  0.000  0.000

B. Matrix Logarithm of P

• Note there is a negative rate from LowCD4 Death!
• This is a transition from Death LowCD4 and should be moved to the other side of the matrix.

B. Matrix Logarithm of P

• The negative rate implies an identifiability issue.
• Essentially, the math only works out if we allow for negative rates—in this case, bringing people back alive from the death state!

B. Matrix Logarithm of P

• expm::logm Higham (2008) method returns same as eigenvalue method.
• Tweaking rates to get a proper model is ad-hoc and difficult.

B. Matrix Logarithm of P

R = expm::logm(mP)

        Healthy LowCD4   AIDS  Death
Healthy  -0.327  0.311  0.002  0.013
LowCD4    0.000 -0.543  0.615 -0.072
AIDS      0.000  0.000 -0.288  0.288
Death     0.000  0.000  0.000  0.000

C. Multi-State Model Based Approach

• Imposition of structural assumptions and fitting to pseudo-data derived from original Markov model can result in reasonable rate estimates.
• We could assume that a patient with HIV at this point in time only gets sicker and external causes of death are negligible. Healthy LowCD4 AIDS Death constrains model.
• We use the reported data from the original model at year 1.

C. Multi-State Model Based Approach

Steps:

1. Run the original model for 1+ cycles to obtain the Markov trace.
2. Construct psuedo-data from the resulting trace based on a cohort with reasonable size (e.g., 1,000 patients)
3. Estimate transition hazards in the pseudo-data based on a multistate model.¹
4. Use the estimated transition hazards to construct the rate matrix.

1. Run the original model

tr0 <- # Initial state occupancy.
  c("Healthy" = 1000,"LowCD4" = 0, "AIDS" = 0, "Death" = 0)
tr1 <- # State occupancy after one cycle.
  tr0 %*% mP 

tr <- # Bind together into a data frame. 
    rbind.data.frame(tr0, tr1) %>%
    mutate(cycle = c(0, 1)) %>% 
    select(cycle,everything())
tr

  cycle Healthy LowCD4 AIDS Death
1     0    1000      0    0     0
2     1     721    202   67    10

2. Construct psuedo-data

• Idea: create data for 1,000 simulated “patients” followed for two cycles.
• Counts in the Markov trace govern state occupancy in each cycle.
  • In first cycle, all 1,000 are in Healthy state.
  • In second cycle, 721 remain in Healthy while 202 are in LowCD4, etc.

2. Construct psuedo-data

tr <- 
    rbind.data.frame(tr0, tr1) %>%
    mutate(cycle = c(0, 1)) %>% 
    gather(state,count,-cycle) %>% 
    mutate(count = round(count,0)) %>% 
    arrange(cycle) %>% 
    lapply(.,rep,.$count) %>% 
    cbind.data.frame() %>% 
    as_tibble() %>% 
    mutate(state = factor(state,
      levels = c("Healthy","LowCD4","AIDS","Death"))) %>% 
    arrange(cycle,state) %>% 
    mutate(ptnum = rep(1:1000,2)) %>% 
    select(-count) %>% 
    arrange(ptnum,cycle) %>% 
    mutate(state = as.numeric(state)) %>% 
    select(ptnum,cycle,state)

# A tibble: 2,000 × 3
   ptnum cycle state
   <int> <dbl> <dbl>
 1     1     0     1
 2     1     1     1
 3     2     0     1
 4     2     1     1
 5     3     0     1
 6     3     1     1
 7     4     0     1
 8     4     1     1
 9     5     0     1
10     5     1     1
# ℹ 1,990 more rows

3. Fit Multistate Model

• Next, fit a multistate model to the pseudo-data.
• The multistate model will use likelihood-based methods to estimate transition intensities among health states in the data.
• These transition intensities are the rates that you can use in the rate matrix!

3. Fit Multistate Model

library(msm)
statetable.msm(state, ptnum, data=tr)

    to
from   1   2   3   4
   1 721 202  67  10

# Markov trace
rbind.data.frame(tr0, tr1) %>%
    mutate(cycle = c(0, 1)) %>% 
    select(cycle,everything())

  cycle Healthy LowCD4 AIDS Death
1     0    1000      0    0     0
2     1     721    202   67    10

3. Fit Multistate Model

Q.init <- rbind(rbind(c(0, 0.3, 0,   .0001),    
                      c(0, 0,   0.6, 0.01),    
                      c(0, 0,   0,   0.3),  
                      c(0, 0,   0,   0)))
dimnames(Q.init) = dimnames(mP)

hiv.msm <- 
  msm(state ~ cycle, subject = ptnum, data = tr, qmatrix = Q.init)

3. Fit Multistate Model

hiv.msm

Call:
msm(formula = state ~ cycle, subject = ptnum, data = tr, qmatrix = Q.init)

Maximum likelihood estimates

Transition intensities
                  Baseline                          
Healthy - Healthy -0.32712 ( -3.680e-01, -2.907e-01)
Healthy - LowCD4   0.32702 (  1.354e-01,  7.898e-01)
Healthy - Death    0.00010 (  0.000e+00,        Inf)
LowCD4 - LowCD4   -0.64477 ( -1.139e+01, -3.650e-02)
LowCD4 - AIDS      0.63465 (  9.115e-06,  4.419e+04)
LowCD4 - Death     0.01012 ( 8.238e-299, 1.242e+294)
AIDS - AIDS       -0.34837 ( -3.805e+40, -3.189e-42)
AIDS - Death       0.34837 (  3.189e-42,  3.805e+40)

-2 * log-likelihood:  1572.2 

4. Construct the Rate Matrix

Rmsm <- msm::qmatrix.msm(hiv.msm, ci = "none")
round(Rmsm,3)

        Healthy LowCD4   AIDS Death
Healthy  -0.327  0.327  0.000 0.000
LowCD4    0.000 -0.645  0.635 0.010
AIDS      0.000  0.000 -0.348 0.348
Death     0.000  0.000  0.000 0.000

5. Embed the Transition Probability Matrix

round(expm(Rmsm),3)

        Healthy LowCD4  AIDS Death
Healthy   0.721  0.202 0.067 0.010
LowCD4    0.000  0.525 0.388 0.088
AIDS      0.000  0.000 0.706 0.294
Death     0.000  0.000 0.000 1.000

round(mP,3)

        Healthy LowCD4  AIDS Death
Healthy   0.721  0.202 0.067 0.010
LowCD4    0.000  0.581 0.407 0.012
AIDS      0.000  0.000 0.750 0.250
Death     0.000  0.000 0.000 1.000

Markov Trace Comparison

# Original Model
c(1000, 0, 0, 0) %*% mP

     Healthy LowCD4 AIDS Death
[1,]     721    202   67    10

# MSM-Based Model
c(1000, 0, 0, 0) %*% pmatrix.msm(hiv.msm, t=1)

     Healthy LowCD4 AIDS Death
[1,]     721    202   67    10

• Resulting expectation after 1 cycle is nearly identical.
• The msm method has bits down at 4th decimal. Round off made them identical.
• Only a tiny bit of error (721 vs. 721.0005) makes log numerically unstable converting from probability to rate.

Key Takeaways

• By imposing some constraints on the underlying transitions, we were able to yield a generator matrix that makes sense!
• Critically, our multistate model-based generator matrix closely approximates the original transition probability matrix, but does not imply negative rates that bring people back from death.

Key Takeaways

round(logm(mP),3)

       [,1]   [,2]   [,3]   [,4]
[1,] -0.327  0.311  0.002  0.013
[2,]  0.000 -0.543  0.615 -0.072
[3,]  0.000  0.000 -0.288  0.288
[4,]  0.000  0.000  0.000  0.000

round(Rmsm,3)

        Healthy LowCD4   AIDS Death
Healthy  -0.327  0.327  0.000 0.000
LowCD4    0.000 -0.645  0.635 0.010
AIDS      0.000  0.000 -0.348 0.348
Death     0.000  0.000  0.000 0.000

Key Takeaways

• With a reasonable generator matrix defined, we can now augment the model as we see fit:

  • Add additional health states.
  • Add new strategies with evidence on efficiacy drawn from the literature (e.g., hazard rates).
  • Add accumulators and transition states.
  • Change the cycle length.

• Once done with the above, just embed the new matrix.