Health states and transition probabilities

The model [usually] starts with everyone in the healthiest state

At each cycle:

• the populations in each state are updated
• the total costs and effects over that cycle are calculated

Transition matrix

We can represent the model graphically using the information in the transition matrix

Our simple model.

##         healthy  ill dead
## healthy    0.97 0.02 0.01
## ill        0.00 0.96 0.04
## dead       0.00 0.00 1.00

The transition probabilities do not reflect how long someone has been in their state, but only which state they are in.

This memoryless behaviour is known as the Markov property.

A cohort simulation

We can carry this on for many cycles, each time multiplying the state population vector from the previous time by the transition matrix.

In Markov models, a cycle is usually a fairly long period of time, commonly a year.

##      healthy       ill      dead
## 1  1000.0000   0.00000   0.00000
## 2   970.0000  20.00000  10.00000
## 3   940.9000  38.60000  20.50000
## 4   912.6730  55.87400  31.45300
## 5   885.2928  71.89250  42.81469
## 6   858.7340  86.72266  54.54332
## 7   832.9720 100.42843  66.59956
## 8   807.9828 113.07073  78.94642
## 9   783.7434 124.70756  91.54908
## 10  760.2311 135.39413 104.37482

How many cycles?

Because there is a non-zero probability of remaining in each state, the system will never reach \(\left(0,0,N\right)\)

We therefore need to employ a stopping criterion to choose the number of cycles

• We will stop when the number of people in the ‘dead’ state rounds up to \(N\)
• Other stopping criteria are available - see Sonnenberg and Beck (1993).

Costs and QALYs

For our example, let’s say our costs and QALYs are

##     State Cost QALY
## 1 healthy   40 0.95
## 2     ill 1000 0.60
## 3    dead    0 0.00

For the first cycle we have 1000 people all in the ‘healthy’ state.

Therefore over that cycle the cohort generate

• a cost of \(1000 \times{40} = 40000\) (in some monetary unit)
• an effect of \(1000\times{0.95}=950\) QALYs.

##   costs QALYs
## 1 40000   950

Cost-effectiveness measures

By

1. summing each column and
2. diving by the number of people in the cohort (1000 in this example)

we can find the

• average cost per person
• average QALYs accrued per person
• average cost per QALY

heemod: running the model

We now collect the information into a strategy and run the model.

The default is to have 1000 people starting in the first state (‘healthy’ in our case):

## the strategy combines the transition matrix with the information about each state

strat_simp = define_strategy(
  transition = mm_simp_heemod,
  healthy = healthy_state, #[<name in transition matrix> = <defined state>]
  ill=ill_state, 
  dead=dead_state
)

We can now run this model with the information in the strategy for some number of cycles

We link the effect and cost arguments to the relevant parts of our state definitions.

res_simp = run_model(
  sim = strat_simp,
  cycles = 87,
  cost = cost,
  effect= QALYs
)

heemod: output

The output is formatted long, so it isn’t so easy to read (but it is good for ggplot2)

## See the first few rows of the counts in each state
head(get_counts(res_simp))

##   .strategy_names markov_cycle state_names    count
## 1             sim            1     healthy 950.0000
## 2             sim            2     healthy 855.0000
## 3             sim            3     healthy 769.5000
## 4             sim            4     healthy 692.5500
## 5             sim            5     healthy 623.2950
## 6             sim            6     healthy 560.9655

## see the first few rows of the costs and QALYs for each state
head(get_values(res_simp))

##   markov_cycle .strategy_names value_names    value
## 1            1             sim        cost  78000.0
## 2            2             sim        cost 145400.0
## 3            3             sim        cost 197036.0
## 4            4             sim        cost 235567.3
## 5            5             sim        cost 263257.2
## 6            6             sim        cost 282028.6

heemod: plotting output

We can plot the state populations and the costs and effects over time

# Plot the populations of each state over time
plot(res_simp, type="counts")

# Plot the populations of each state over time
plot(res_simp, type="value", panel = "by_value", free_y = TRUE)

Example - two drugs

Suppose we are performing a cost-effectiveness analysis for two treatments of a chronic disease

• ‘drug A’: the current standard

• ‘drug B’: a newer, more effective but more costly treatment
  
  Our model has three states:

• ‘mild’, ‘severe’ and ‘dead’

• there is no returning from ‘severe’ to ‘mild’.

Setting up the model

The transition matrix for drug A is

\[\begin{pmatrix} 0.75& 0.23 & 0.02\\ 0 & 0.9 & 0.1\\ 0 & 0 & 1 \end{pmatrix}\]

We can specify this in R as follows:

tm_drugA = define_transition(
  state_names = c("mild", "severe", "dead"),
  0.75, 0.23, 0.02,
  0, 0.9, 0.1, 
  0, 0, 1
)
tm_drugA

## A transition matrix, 3 states.
## 
##        mild severe dead
## mild   0.75 0.23   0.02
## severe      0.9    0.1 
## dead               1

Setting up the model

It is thought that while on drug B, the probability of moving to a worse state is generally less.

We represent this by a parameter

Relative risk

The probability of transitioning from one state to a worse state with drug B is given by the relative risk \((rr)\) multiplied by the probability of the same transition with drug A.

So, in our example:

Where \(C\) is the complement - the value necessary for the probabilities to sum to one.

Relative risk in heemod

In the code below

• rr is our relative risk
• C stands for complement, and is one minus the sum of the other probabilities in that row (ensuring the row sums to one).

rr = 0.7  ## Our relative risk
tm_drugB = define_transition(
  state_names = c("mild", "severe", "dead"),
  C, 0.23*rr, 0.02*rr,
  0, C, 0.1*rr,
  0, 0, 1
)

tm_drugB

## A transition matrix, 3 states.
## 
##        mild severe    dead     
## mild   C    0.23 * rr 0.02 * rr
## severe      C         0.1 * rr 
## dead                  1

Splitting costs

In this example we will split the costs of drugs and healthcare.
We will have to specify costs for each drug (strategy)

The function dispatch_strategy specifies a different return value for each strategy (each drug, in this example).

It can be used in any place where the value returned depends on the strategy.
The name used needs to be the same as in run_model.

Dispatch strategy

# setting up some initial cost variables
cost_drugA = 2000
cost_drugB = 2700

state_mild= define_state(

# cost of healthcare doesn't depend on strategy (drug)
  cost_health = 1500,     
# cost of drug depends on strategy
  cost_drugs = dispatch_strategy(
    drugA = cost_drugA, 
    drugB = cost_drugB
  ), 

  cost_total = cost_health + cost_drugs,
  QALYs = 0.9  # QALYs is our choice of effect here. 
  ## We could use a different measure, so long as we linked this to the 
  ## 'effect' argument in run_model.
)

Note that we now specify cost_health, cost_drugs and cost_total.

Comparing treatments: strategies

In this model, the different treatment strategies are distinguished by:

• Difference in probabilities of transitioning between states (seen in the transition matrix)
• Cost of drug (seen in the state definitions)

We now combine this information together to generate our two Strategies.

strat_drugA = define_strategy(
  transition = tm_drugA,
  mild = state_mild,        # the state definitions distinguish between strategies
  severe = state_severe,    # via `dispatch_strategy`
  dead = state_dead
)

strat_drugB = define_strategy(
  transition = tm_drugB,
  mild = state_mild,
  severe = state_severe,
  dead = state_dead
)

Comparing treatments: running the model

We include both strategies in the run_model function.

res_drugsAB <- run_model(
  drugA = strat_drugA,    # the name here is what's used in `dispatch_strategy`
  drugB = strat_drugB,
  cycles = 143,            # How many cycles do we want to run for
  cost = cost_total,
  effect = QALYs
)

Comparing treatments: results

summary(res_drugsAB)

## 2 strategies run for 143 cycles.
## 
## Initial state counts:
## 
## mild = 1000L
## severe = 0L
## dead = 0L
## 
## Counting method: 'life-table'.
## 
##  
## 
## Counting method: 'beginning'.
## 
##  
## 
## Counting method: 'end'.
## 
## Values:
## 
##       cost_health cost_drugs cost_total     QALYs
## drugA    24645721   19649231   44294952  8669.997
## drugB    31508935   34547463   66056398 12578.177
## 
## Efficiency frontier:
## 
## drugA -> drugB
## 
## Differences:
## 
##       Cost Diff. Effect Diff.    ICER  Ref.
## drugB   21761.45     3.908179 5568.18 drugA

Comparing treatments: Plotting

Costs and effects over time:

plot(res_drugsAB, type = "values", panel = "by_value",
     free_y = TRUE)

These are all built in ggplot2, so if you would like to change the appearance in any way (eg. axis labels, colours) you can use the usual ggplot2 syntax.