/
a_introduction.Rmd
194 lines (134 loc) · 6.87 KB
/
a_introduction.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
---
title: "An Introduction to `heemod`"
date: "`r Sys.Date()`"
output:
rmarkdown::html_vignette:
toc: true
vignette: >
%\VignetteIndexEntry{An Introduction to `heemod`}
%\VignetteEngine{knitr::rmarkdown}
\usepackage[utf8]{inputenc}
---
```{r, echo=FALSE, include=FALSE}
library(heemod)
```
This document is a presentation of the basic steps to define and run a model in `heemod`. Note that decision trees are actually a subset even of *Markov model*, and thus can be specified easily with this package.
When building a Markov model the following steps must be performed:
* Specify transition probabilities between states.
* Specify values attached to states (costs, utilities, etc.).
* Combine this information and run the model.
## Transition probabilities
The probability to transition from one state to another during a time period is called a *transition probability*. Said time period is called a *cycle*.
Transition probabilities between states can be described through a 2-way table where the lines correspond to the states at the beginning of a cycle and the columns to the states at the end of a cycle. Consider a model with 2 states `A` and `B`:
| | A | B |
|-------|---|---|
| **A** | 1 | 2 |
| **B** | 3 | 4 |
When starting a cycle in state **A** (row **A**), the probability to still be in state **A** at the end of the cycle is found in colunm **A** (cell **1**) and the probability to change to state **B** is found in column **B** (cell **2**).
Similarly, when starting a cycle from state **B** (row **B**), the probability to be in state **A** or **B** at the end of the cycle are found in cells **3** or **4** respectively.
In the context of Markov models, this 2-way table is called a *transition matrix*. A transition matrix can be defined easily in `heemod` with the `define_transition()` function. If we consider the previous example, where cell values have been replaced by actual probabilities:
| | A | B |
|-------|-----|-----|
| **A** | 0.9 | 0.1 |
| **B** | 0.2 | 0.8 |
That transition matrix can be defined with the following command:
```{r}
mat_trans <- define_transition(
.9, .1,
.2, .8
)
mat_trans
```
## Attach values to states
Values are attached to states. Cost and utility are classical examples of such values. To continue with the previous example, the following values can be attached to state **A** and **B**:
* State **A** has a cost of *1234* per cycle and an utility of *0.85*.
* State **B** has a cost of *4321* per cycle and an utility of *0.50*.
A state and its values can be defined with `define_state()`:
```{r}
state_A <- define_state(
cost = 1234,
utility = 0.85
)
state_A
state_B <- define_state(
cost = 4321,
utility = 0.50
)
state_B
```
## Combine information in a model
Now that the transition matrix and the state values are defined for a given strategy, we can combine them with `define_strategy()`:
```{r}
strat <- define_strategy(
transition = mat_trans,
state_A,
state_B
)
strat
```
## Run the model
A model is a made of one or more strategy, run with `run_model()` for a given number of cycles. Here we run only one strategy, so no comparison of the cost-effectiveness of the different strategies will be made. The variables corresponding to valuation of cost and effect must be given at that point.
```{r}
res_mod <- run_model(
strat,
cycles = 10,
cost = cost,
effect = utility
)
res_mod
```
By default the model is run for *1000* persons starting in the first state (here state **A**).
## Analyze results
We can plot the state membership counts over time. Other plot types are available.
```{r, fig.width = 6, fig.height=4, fig.align='center'}
plot(res_mod)
```
Plots can be modified using `ggplot2` syntax.
```{r, fig.align='center', fig.height=4, fig.width=6, message=FALSE}
library(ggplot2)
plot(res_mod) +
xlab("Time") +
ylab("N") +
theme_minimal() +
scale_color_brewer(
name = "State",
palette = "Set1"
)
```
And black & white plots for publication are available with the `bw` plot option
```{r, fig.align='center', fig.height=4, fig.width=6, message=FALSE}
library(ggplot2)
plot(res_mod, bw = TRUE)
```
The state membership counts and the values can be accessed with `get_counts()` and `get_values()` respectively.
```{r}
head(get_counts(res_mod))
head(get_values(res_mod))
```
## Convenience functions
Convenience functions are available to easily compute transition probabilities from incidence rates, OR, RR, or probabilities estimated on a different timeframe.
Example : convert an incidence rate of *162* cases per *1,000* person-years to a 5-year probability.
```{r}
rate_to_prob(r = 162, per = 1000, to = 5)
```
See `?probability` to see a list of the convenience functions available.
## External data
Mortality rates by age and sex (often used as transition probabilities) can be downloaded from the WHO online database with the function `get_who_mr()`.
External data contained in user-defined data frames can be referenced in a model with the function `look_up()`.
## Compare strategies
In order to compare different strategies it is possible to run a model with multiple strategies in parallel, examples are provided in `vignette("c-homogeneous", "heemod")` or `vignette("d-non-homogeneous", "heemod")`.
## Time-variability
Time varying transition probabilities and state values are available through the `model_time` and `state_time` variables. Time-varying probabilities and values are explained in `vignette("b-time-dependency", "heemod")`.
## Transition costs
Transition costs between states can be defined by adding an additional state called a transition state. This is a special case of a tunnel state. A situation where the transition from `A` to `B` (noted `A->B`) costs \$100 and reduces utility by 0.1 points compared to the usual values of `B` can be modelled by adding a state `B_trans`.
The cost and utility of `B_trans` are the same as those of `B` +\$100 and -0.1 QALY respectively.
`A->B_trans` is equal to the former value of `A->B`. `A->B` is set to 0 (all transitions from `A` to `B` must pass through `B_trans` from now on).
The probability to stay in `B_trans` is 0, `B_trans->B` is equal to `B->B` and similarly all `B_trans->*` are equal to `B->*`.
## Parallel computing
Some operation such as PSA can be significantly sped up using parallel computing. This can be done in the following way:
* Define a cluster with the `use_cluster()` functions (i.e. `use_cluster(4)` to use 4 cores).
* Run the analysis as usual.
* To stop using parallel computing use the `close_cluster()` function.
Results may vary depending on the machine, but we found speed gains to be quite limited beyond 4 cores.
## Going further
Probabilistic uncertainty analysis `vignette("e-probabilistic", "heemod")` and deterministic sensitivity analysis `vignette("f-sensitivity", "heemod")` can be performed. Population-level estimates and heterogeneity can be computed : `vignette("g-heterogeneity", "heemod")`.