Markov Cohort Model
DW Richardson
2022-08-20
Introduction
Markov Cohort Models area common practice in health economics. The general notion is to follow a cohort across time periods as they transition to various stages of health which could range from having a headache or death. The states of health are arbitrary. The model relies on transition probabilities which are used in the simulation to determine the flow of cohort members among states between time periods or cycles. Markov Cohort Models also tend to utilize a costs and benefits associated with different states. This typically comes in the form of a payoff matrix in a decision model.
Defining Base Parameters
The model here is quite generic with only three states Healthy, Sick, and Dead denoted as n_s and v_state_names. The number of cycles or time periods is n_t. While the number of individuals in this cohort is 1000.
#Create number of cycles
n_t <- 40
#States
n_s <- 3
#Cohort
n_c <- 1000
#State names
v_state_names <- c("Healthy", "Sick", "Dead")Transistion Probability Matrix
From an R perspective an actual matrix is not being setup, rather an array is, but for functional purposes we are creating a matrix of probabilities associated with either remaining in the current state or transitioning to another state with each cycle. This is setup as an array because it has three dimensions instead of two.
After creating the array full of NA’s, the array is filled in by running each of the probability setting created. The first three probabilities are more like rules. In this instance we have set the rules so no one in cohort can revert to a better state. The remaining probabilities are defining the odds of transitioning for the state first cited to the second state cited in the call. Notice several of the probabilities ignore a third argument as they are applied equally across all cycles. However, transitions from Healthy to Dead include a third argument allowing for the probability of death to increase with age. The final probability being defined is leveraging R to give us the probability of remaining healthy by simply subtracting the other probabilities of transition from 100%.
Note that row sums should equate to 1.0.
# Matrix Probabilities with transition periods
trans_mat <- array(NA_real_,
dim = c(n_s, n_s, n_t),
dimnames = list(from = v_state_names,
to = v_state_names,
cycle = 1:n_t))
#Rules (like dead can't go healthy)
trans_mat[2,1, ] <- 0
trans_mat[3,1, ] <- 0
trans_mat[3,2, ] <- 0
#Base transition probabilities
trans_mat[1,2, ] <- 0.03
trans_mat[3,3, ] <- 1
#Vary transition probabilities
trans_mat[1,3, 1:10] <- 0.01
trans_mat[1,3, 11:20 ] <- 0.02
trans_mat[1,3, 21:30] <- 0.04
trans_mat[1,3, 31:40] <- 0.08
trans_mat[2,3, ] <- trans_mat[1,3, ] + 0.04
trans_mat[2, 2, ] <- 1 - trans_mat[2, 3, ]
trans_mat[1,1,] <- 1 - apply(trans_mat[1, , ], 2, sum, na.rm = TRUE)
head(trans_mat, 5)## , , cycle = 1
##
## to
## from Healthy Sick Dead
## Healthy 0.96 0.03 0.01
## Sick 0.00 0.95 0.05
## Dead 0.00 0.00 1.00
##
## , , cycle = 2
##
## to
## from Healthy Sick Dead
## Healthy 0.96 0.03 0.01
## Sick 0.00 0.95 0.05
## Dead 0.00 0.00 1.00
##
## , , cycle = 3
##
## to
## from Healthy Sick Dead
## Healthy 0.96 0.03 0.01
## Sick 0.00 0.95 0.05
## Dead 0.00 0.00 1.00
##
## , , cycle = 4
##
## to
## from Healthy Sick Dead
## Healthy 0.96 0.03 0.01
## Sick 0.00 0.95 0.05
## Dead 0.00 0.00 1.00
##
## , , cycle = 5
##
## to
## from Healthy Sick Dead
## Healthy 0.96 0.03 0.01
## Sick 0.00 0.95 0.05
## Dead 0.00 0.00 1.00
##
## , , cycle = 6
##
## to
## from Healthy Sick Dead
## Healthy 0.96 0.03 0.01
## Sick 0.00 0.95 0.05
## Dead 0.00 0.00 1.00
##
## , , cycle = 7
##
## to
## from Healthy Sick Dead
## Healthy 0.96 0.03 0.01
## Sick 0.00 0.95 0.05
## Dead 0.00 0.00 1.00
##
## , , cycle = 8
##
## to
## from Healthy Sick Dead
## Healthy 0.96 0.03 0.01
## Sick 0.00 0.95 0.05
## Dead 0.00 0.00 1.00
##
## , , cycle = 9
##
## to
## from Healthy Sick Dead
## Healthy 0.96 0.03 0.01
## Sick 0.00 0.95 0.05
## Dead 0.00 0.00 1.00
##
## , , cycle = 10
##
## to
## from Healthy Sick Dead
## Healthy 0.96 0.03 0.01
## Sick 0.00 0.95 0.05
## Dead 0.00 0.00 1.00
##
## , , cycle = 11
##
## to
## from Healthy Sick Dead
## Healthy 0.95 0.03 0.02
## Sick 0.00 0.94 0.06
## Dead 0.00 0.00 1.00
##
## , , cycle = 12
##
## to
## from Healthy Sick Dead
## Healthy 0.95 0.03 0.02
## Sick 0.00 0.94 0.06
## Dead 0.00 0.00 1.00
##
## , , cycle = 13
##
## to
## from Healthy Sick Dead
## Healthy 0.95 0.03 0.02
## Sick 0.00 0.94 0.06
## Dead 0.00 0.00 1.00
##
## , , cycle = 14
##
## to
## from Healthy Sick Dead
## Healthy 0.95 0.03 0.02
## Sick 0.00 0.94 0.06
## Dead 0.00 0.00 1.00
##
## , , cycle = 15
##
## to
## from Healthy Sick Dead
## Healthy 0.95 0.03 0.02
## Sick 0.00 0.94 0.06
## Dead 0.00 0.00 1.00
##
## , , cycle = 16
##
## to
## from Healthy Sick Dead
## Healthy 0.95 0.03 0.02
## Sick 0.00 0.94 0.06
## Dead 0.00 0.00 1.00
##
## , , cycle = 17
##
## to
## from Healthy Sick Dead
## Healthy 0.95 0.03 0.02
## Sick 0.00 0.94 0.06
## Dead 0.00 0.00 1.00
##
## , , cycle = 18
##
## to
## from Healthy Sick Dead
## Healthy 0.95 0.03 0.02
## Sick 0.00 0.94 0.06
## Dead 0.00 0.00 1.00
##
## , , cycle = 19
##
## to
## from Healthy Sick Dead
## Healthy 0.95 0.03 0.02
## Sick 0.00 0.94 0.06
## Dead 0.00 0.00 1.00
##
## , , cycle = 20
##
## to
## from Healthy Sick Dead
## Healthy 0.95 0.03 0.02
## Sick 0.00 0.94 0.06
## Dead 0.00 0.00 1.00
##
## , , cycle = 21
##
## to
## from Healthy Sick Dead
## Healthy 0.93 0.03 0.04
## Sick 0.00 0.92 0.08
## Dead 0.00 0.00 1.00
##
## , , cycle = 22
##
## to
## from Healthy Sick Dead
## Healthy 0.93 0.03 0.04
## Sick 0.00 0.92 0.08
## Dead 0.00 0.00 1.00
##
## , , cycle = 23
##
## to
## from Healthy Sick Dead
## Healthy 0.93 0.03 0.04
## Sick 0.00 0.92 0.08
## Dead 0.00 0.00 1.00
##
## , , cycle = 24
##
## to
## from Healthy Sick Dead
## Healthy 0.93 0.03 0.04
## Sick 0.00 0.92 0.08
## Dead 0.00 0.00 1.00
##
## , , cycle = 25
##
## to
## from Healthy Sick Dead
## Healthy 0.93 0.03 0.04
## Sick 0.00 0.92 0.08
## Dead 0.00 0.00 1.00
##
## , , cycle = 26
##
## to
## from Healthy Sick Dead
## Healthy 0.93 0.03 0.04
## Sick 0.00 0.92 0.08
## Dead 0.00 0.00 1.00
##
## , , cycle = 27
##
## to
## from Healthy Sick Dead
## Healthy 0.93 0.03 0.04
## Sick 0.00 0.92 0.08
## Dead 0.00 0.00 1.00
##
## , , cycle = 28
##
## to
## from Healthy Sick Dead
## Healthy 0.93 0.03 0.04
## Sick 0.00 0.92 0.08
## Dead 0.00 0.00 1.00
##
## , , cycle = 29
##
## to
## from Healthy Sick Dead
## Healthy 0.93 0.03 0.04
## Sick 0.00 0.92 0.08
## Dead 0.00 0.00 1.00
##
## , , cycle = 30
##
## to
## from Healthy Sick Dead
## Healthy 0.93 0.03 0.04
## Sick 0.00 0.92 0.08
## Dead 0.00 0.00 1.00
##
## , , cycle = 31
##
## to
## from Healthy Sick Dead
## Healthy 0.89 0.03 0.08
## Sick 0.00 0.88 0.12
## Dead 0.00 0.00 1.00
##
## , , cycle = 32
##
## to
## from Healthy Sick Dead
## Healthy 0.89 0.03 0.08
## Sick 0.00 0.88 0.12
## Dead 0.00 0.00 1.00
##
## , , cycle = 33
##
## to
## from Healthy Sick Dead
## Healthy 0.89 0.03 0.08
## Sick 0.00 0.88 0.12
## Dead 0.00 0.00 1.00
##
## , , cycle = 34
##
## to
## from Healthy Sick Dead
## Healthy 0.89 0.03 0.08
## Sick 0.00 0.88 0.12
## Dead 0.00 0.00 1.00
##
## , , cycle = 35
##
## to
## from Healthy Sick Dead
## Healthy 0.89 0.03 0.08
## Sick 0.00 0.88 0.12
## Dead 0.00 0.00 1.00
##
## , , cycle = 36
##
## to
## from Healthy Sick Dead
## Healthy 0.89 0.03 0.08
## Sick 0.00 0.88 0.12
## Dead 0.00 0.00 1.00
##
## , , cycle = 37
##
## to
## from Healthy Sick Dead
## Healthy 0.89 0.03 0.08
## Sick 0.00 0.88 0.12
## Dead 0.00 0.00 1.00
##
## , , cycle = 38
##
## to
## from Healthy Sick Dead
## Healthy 0.89 0.03 0.08
## Sick 0.00 0.88 0.12
## Dead 0.00 0.00 1.00
##
## , , cycle = 39
##
## to
## from Healthy Sick Dead
## Healthy 0.89 0.03 0.08
## Sick 0.00 0.88 0.12
## Dead 0.00 0.00 1.00
##
## , , cycle = 40
##
## to
## from Healthy Sick Dead
## Healthy 0.89 0.03 0.08
## Sick 0.00 0.88 0.12
## Dead 0.00 0.00 1.00Determining State Membership
Similar to the transition probability matrix this is technically setup as an array. Initially the array is filled with NA’s and the initial state is defined. In this example the entire cohort is initially placed in the Healthy state. A loop is then run to replace the NA’s. Each iteration of the loop is matrix multiplying our transition probability matrix by the previous cycle’s state membership until completion.
#Create parking lot for state membership
state_membership <- array(NA_real_, dim = c(n_t, n_s),
dimnames = list(cycle = 1:n_t, state = v_state_names))
#Create initial state
state_membership[1, ] <- c(n_c, 0, 0)
#State membership iteration loop
for (i in 2:n_t) { state_membership[i, ] <- state_membership[i-1, ] %*%
trans_mat[ , , i - 1]}
head(state_membership, 5)## state
## cycle Healthy Sick Dead
## 1 1000.0000 0.0000 0.00000
## 2 960.0000 30.0000 10.00000
## 3 921.6000 57.3000 21.10000
## 4 884.7360 82.0830 33.18100
## 5 849.3466 104.5209 46.13251Defining Payoff Matrix
For our model each state in each cycle has a cost/benefit component. The cost is clear denoted as cost this can be framed as the literal health care cost of dealing with cohort members. The benefit is less clear but QALY can be thought of as a rather all-encompassing quantification of quality of life.
#Create payoff matrix
payoffs <- array(NA_real_,
dim = c(n_s, 2, n_t),
dimnames = list(state = v_state_names,
payoff = c("Cost", "QALY"),
cycle = 1:n_t))
payoffs[, , 1:10] <- c(10, 800, 0, 0.95, 0.65, 0)
payoffs[, , 11:20] <- c(25, 1000, 0, 0.92, 0.60, 0)
payoffs[, , 21:30] <- c(40, 1200, 0, 0.88, 0.55, 0)
payoffs[, , 31:40] <- c(80, 1000, 0, 0.85, 0.50, 0)
head(payoffs, 5)## , , cycle = 1
##
## payoff
## state Cost QALY
## Healthy 10 0.95
## Sick 800 0.65
## Dead 0 0.00
##
## , , cycle = 2
##
## payoff
## state Cost QALY
## Healthy 10 0.95
## Sick 800 0.65
## Dead 0 0.00
##
## , , cycle = 3
##
## payoff
## state Cost QALY
## Healthy 10 0.95
## Sick 800 0.65
## Dead 0 0.00
##
## , , cycle = 4
##
## payoff
## state Cost QALY
## Healthy 10 0.95
## Sick 800 0.65
## Dead 0 0.00
##
## , , cycle = 5
##
## payoff
## state Cost QALY
## Healthy 10 0.95
## Sick 800 0.65
## Dead 0 0.00
##
## , , cycle = 6
##
## payoff
## state Cost QALY
## Healthy 10 0.95
## Sick 800 0.65
## Dead 0 0.00
##
## , , cycle = 7
##
## payoff
## state Cost QALY
## Healthy 10 0.95
## Sick 800 0.65
## Dead 0 0.00
##
## , , cycle = 8
##
## payoff
## state Cost QALY
## Healthy 10 0.95
## Sick 800 0.65
## Dead 0 0.00
##
## , , cycle = 9
##
## payoff
## state Cost QALY
## Healthy 10 0.95
## Sick 800 0.65
## Dead 0 0.00
##
## , , cycle = 10
##
## payoff
## state Cost QALY
## Healthy 10 0.95
## Sick 800 0.65
## Dead 0 0.00
##
## , , cycle = 11
##
## payoff
## state Cost QALY
## Healthy 25 0.92
## Sick 1000 0.60
## Dead 0 0.00
##
## , , cycle = 12
##
## payoff
## state Cost QALY
## Healthy 25 0.92
## Sick 1000 0.60
## Dead 0 0.00
##
## , , cycle = 13
##
## payoff
## state Cost QALY
## Healthy 25 0.92
## Sick 1000 0.60
## Dead 0 0.00
##
## , , cycle = 14
##
## payoff
## state Cost QALY
## Healthy 25 0.92
## Sick 1000 0.60
## Dead 0 0.00
##
## , , cycle = 15
##
## payoff
## state Cost QALY
## Healthy 25 0.92
## Sick 1000 0.60
## Dead 0 0.00
##
## , , cycle = 16
##
## payoff
## state Cost QALY
## Healthy 25 0.92
## Sick 1000 0.60
## Dead 0 0.00
##
## , , cycle = 17
##
## payoff
## state Cost QALY
## Healthy 25 0.92
## Sick 1000 0.60
## Dead 0 0.00
##
## , , cycle = 18
##
## payoff
## state Cost QALY
## Healthy 25 0.92
## Sick 1000 0.60
## Dead 0 0.00
##
## , , cycle = 19
##
## payoff
## state Cost QALY
## Healthy 25 0.92
## Sick 1000 0.60
## Dead 0 0.00
##
## , , cycle = 20
##
## payoff
## state Cost QALY
## Healthy 25 0.92
## Sick 1000 0.60
## Dead 0 0.00
##
## , , cycle = 21
##
## payoff
## state Cost QALY
## Healthy 40 0.88
## Sick 1200 0.55
## Dead 0 0.00
##
## , , cycle = 22
##
## payoff
## state Cost QALY
## Healthy 40 0.88
## Sick 1200 0.55
## Dead 0 0.00
##
## , , cycle = 23
##
## payoff
## state Cost QALY
## Healthy 40 0.88
## Sick 1200 0.55
## Dead 0 0.00
##
## , , cycle = 24
##
## payoff
## state Cost QALY
## Healthy 40 0.88
## Sick 1200 0.55
## Dead 0 0.00
##
## , , cycle = 25
##
## payoff
## state Cost QALY
## Healthy 40 0.88
## Sick 1200 0.55
## Dead 0 0.00
##
## , , cycle = 26
##
## payoff
## state Cost QALY
## Healthy 40 0.88
## Sick 1200 0.55
## Dead 0 0.00
##
## , , cycle = 27
##
## payoff
## state Cost QALY
## Healthy 40 0.88
## Sick 1200 0.55
## Dead 0 0.00
##
## , , cycle = 28
##
## payoff
## state Cost QALY
## Healthy 40 0.88
## Sick 1200 0.55
## Dead 0 0.00
##
## , , cycle = 29
##
## payoff
## state Cost QALY
## Healthy 40 0.88
## Sick 1200 0.55
## Dead 0 0.00
##
## , , cycle = 30
##
## payoff
## state Cost QALY
## Healthy 40 0.88
## Sick 1200 0.55
## Dead 0 0.00
##
## , , cycle = 31
##
## payoff
## state Cost QALY
## Healthy 80 0.85
## Sick 1000 0.50
## Dead 0 0.00
##
## , , cycle = 32
##
## payoff
## state Cost QALY
## Healthy 80 0.85
## Sick 1000 0.50
## Dead 0 0.00
##
## , , cycle = 33
##
## payoff
## state Cost QALY
## Healthy 80 0.85
## Sick 1000 0.50
## Dead 0 0.00
##
## , , cycle = 34
##
## payoff
## state Cost QALY
## Healthy 80 0.85
## Sick 1000 0.50
## Dead 0 0.00
##
## , , cycle = 35
##
## payoff
## state Cost QALY
## Healthy 80 0.85
## Sick 1000 0.50
## Dead 0 0.00
##
## , , cycle = 36
##
## payoff
## state Cost QALY
## Healthy 80 0.85
## Sick 1000 0.50
## Dead 0 0.00
##
## , , cycle = 37
##
## payoff
## state Cost QALY
## Healthy 80 0.85
## Sick 1000 0.50
## Dead 0 0.00
##
## , , cycle = 38
##
## payoff
## state Cost QALY
## Healthy 80 0.85
## Sick 1000 0.50
## Dead 0 0.00
##
## , , cycle = 39
##
## payoff
## state Cost QALY
## Healthy 80 0.85
## Sick 1000 0.50
## Dead 0 0.00
##
## , , cycle = 40
##
## payoff
## state Cost QALY
## Healthy 80 0.85
## Sick 1000 0.50
## Dead 0 0.00Determining Payoffs
Again an empty array of NA’s is initially setup, but then we use matrix multiplication of the payoff matrix with state membership matrix via a loop. The result is a calculation of cost and QALY for each cycle.
payoff_trace <- array(NA_real_,
dim = c(n_t, 2),
dimnames = list(cycle = 1:n_t, payoff = c("Cost", "QALY")))
for(i in 1:n_t) {
payoff_trace[i, ] <- state_membership[i, ] %*% payoffs[, , i]
}
head(payoff_trace, 5)## payoff
## cycle Cost QALY
## 1 10000.00 950.0000
## 2 33600.00 931.5000
## 3 55056.00 912.7650
## 4 74513.76 893.8531
## 5 92110.21 874.8178Findings
The average is calculated for both cost and QALY. These can serve as a great point of reference to see how model performs when other parameters are changed. The first plot is showing trends for state as expected the Healthy cohort decline while the Sick rises. The Death state accelerates upward before dipping. The second plot shows the payoff trace.
#Get Avg
colSums(payoff_trace)/n_c## Cost QALY
## 6918.99231 20.10628#Plotting for state trends
matplot(1:n_t, state_membership, type = 'l')
matplot(1:n_t, payoff_trace, type = 'l')