-
Notifications
You must be signed in to change notification settings - Fork 25
/
c_homogeneous.Rmd
226 lines (177 loc) · 6.45 KB
/
c_homogeneous.Rmd
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
---
title: "Simple Markov Models (Homogeneous)"
date: "`r Sys.Date()`"
output:
rmarkdown::html_vignette:
toc: true
vignette: >
%\VignetteIndexEntry{Simple Markov Models (Homogeneous)}
%\VignetteEngine{knitr::rmarkdown}
\usepackage[utf8]{inputenc}
---
```{r, echo=FALSE, include=FALSE}
library(heemod)
library(ggplot2)
```
The most simple Markov models in health economic evaluation are models were transition probabilities between states do not change with time. Those are called *homogeneous* or *time-homogeneous* Markov models.
## Model description
In this example we will model the cost effectiveness of lamivudine/zidovudine combination therapy in HIV infection ([Chancellor, 1997](https://pubmed.ncbi.nlm.nih.gov/10169387/) further described in [Decision Modelling for Health Economic Evaluation](https://global.oup.com/academic/product/decision-modelling-for-health-economic-evaluation-9780198526629/), page 32. For the sake of simplicity we will not reproduce exactly the analysis from the book. See vignette `vignette("i-reproduction", "heemod")` for an exact reproduction of the analysis.
This model aims to compare costs and utilities of two treatment strategies, *monotherapy* and *combined therapy*.
Four states are described, from best to worst health-wise:
* __A__: CD4 cells > 200 and < 500 cells/mm3;
* __B__: CD4 < 200 cells/mm3, non-AIDS;
* __C__: AIDS;
* __D__: Death.
## Transition probabilities
Transition probabilities for the monotherapy study group are rather simple to implement with `define_transition()`:
```{r}
mat_mono <- define_transition(
.721, .202, .067, .010,
0, .581, .407, .012,
0, 0, .750, .250,
0, 0, 0, 1
)
mat_mono
```
The combined therapy group has its transition probabilities multiplied by $rr = 0.509$, the relative risk of event for the population treated by combined therapy. Since $rr < 1$, the combined therapy group has less chance to transition to worst health states.
The probabilities to stay in the same state are equal to $1 - \sum P_{trans}$ where $P_{trans}$ are the probabilities to change to another state (because all transition probabilities from a given state must sum to 1).
We use the alias `C` as a convenient way to specify the probability complement, equal to $1 - \sum P_{trans}$.
```{r}
rr <- .509
mat_comb <- define_transition(
C, .202*rr, .067*rr, .010*rr,
0, C, .407*rr, .012*rr,
0, 0, C, .250*rr,
0, 0, 0, 1
)
mat_comb
```
We can plot the transition matrix for the monotherapy group:
```{r, fig.width = 6, fig.height=6, fig.align='center'}
plot(mat_mono)
```
And the combined therapy group:
```{r, fig.width = 6, fig.height=6, fig.align='center'}
plot(mat_comb)
```
## State values
The costs of lamivudine and zidovudine are defined:
```{r}
cost_zido <- 2278
cost_lami <- 2086
```
In addition to drugs costs (called `cost_drugs` in the model), each state is associated to healthcare costs (called `cost_health`). Cost are discounted at a 6% rate with the `discount` function.
Efficacy in this study is measured in terms of life expectancy (called `life_year` in the model). Each state thus has a value of 1 life year per year, except death who has a value of 0. Life-years are not discounted in this example.
Only `cost_drug` differs between the monotherapy and the combined therapy treatment groups, the function `dispatch_strategy()` can be used to account for that. For example state A can be defined with `define_state()`:
```{r}
state_A <- define_state(
cost_health = discount(2756, .06),
cost_drugs = discount(dispatch_strategy(
mono = cost_zido,
comb = cost_zido + cost_lami
), .06),
cost_total = cost_health + cost_drugs,
life_year = 1
)
state_A
```
The other states for the monotherapy treatment group can be specified in the same way:
```{r}
state_B <- define_state(
cost_health = discount(3052, .06),
cost_drugs = discount(dispatch_strategy(
mono = cost_zido,
comb = cost_zido + cost_lami
), .06),
cost_total = cost_health + cost_drugs,
life_year = 1
)
state_C <- define_state(
cost_health = discount(9007, .06),
cost_drugs = discount(dispatch_strategy(
mono = cost_zido,
comb = cost_zido + cost_lami
), .06),
cost_total = cost_health + cost_drugs,
life_year = 1
)
state_D <- define_state(
cost_health = 0,
cost_drugs = 0,
cost_total = 0,
life_year = 0
)
```
## Strategy definitions
Strategies can now be defined by combining a transition matrix and a state list with `define_strategy()`:
```{r}
strat_mono <- define_strategy(
transition = mat_mono,
state_A,
state_B,
state_C,
state_D
)
strat_mono
```
For the combined therapy model:
```{r}
strat_comb <- define_strategy(
transition = mat_comb,
state_A,
state_B,
state_C,
state_D
)
```
## Running the model
Both strategies can then be combined in a model and run for 50 years with `run_model()`. Strategies are given names (`mono` and `comb`) in order to facilitate result interpretation.
```{r}
res_mod <- run_model(
mono = strat_mono,
comb = strat_comb,
cycles = 50,
cost = cost_total,
effect = life_year
)
```
By default models are run for 1000 persons starting in the first state (here state **A**).
## Result interpretation
Strategy values can then be compared with `summary()` (optionally net monetary benefits can be calculated with the `threshold` option):
```{r}
summary(res_mod,
threshold = c(1000, 5000, 6000, 1e4))
```
The incremental cost-effectiveness ratio of the combined therapy strategy is thus £`r round(summary(res_mod)$res_comp$.icer[2])` per life-year gained.
The counts per state can be plotted by model:
```{r, fig.align='center', fig.width=6, fig.height=6, message=FALSE}
plot(res_mod, type = "counts", panel = "by_strategy") +
xlab("Time") +
theme_bw() +
scale_color_brewer(
name = "State",
palette = "Set1"
)
```
Or by state:
```{r, fig.align='center', fig.width=6, fig.height=8, message=FALSE}
plot(res_mod, type = "counts", panel = "by_state") +
xlab("Time") +
theme_bw() +
scale_color_brewer(
name = "Strategy",
palette = "Set1"
)
```
The values can also be represented:
```{r, fig.align='center', fig.width=6, fig.height=8, message=FALSE}
plot(res_mod, type = "values", panel = "by_value",
free_y = TRUE) +
xlab("Time") +
theme_bw() +
scale_color_brewer(
name = "Strategy",
palette = "Set1"
)
```
Note that classic `ggplot2` syntax can be used to modifiy plot appearance.