-
Notifications
You must be signed in to change notification settings - Fork 4
/
foi_models.Rmd
251 lines (215 loc) · 13.2 KB
/
foi_models.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
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
---
title: "An introduction to Force-of-Infection (FoI) Models"
output: rmarkdown::html_vignette
bibliography: references.bib
link-citations: true
vignette: >
%\VignetteIndexEntry{An introduction to Force-of-Infection (FoI) Models}
%\VignetteEngine{knitr::rmarkdown}
%\VignetteEncoding{UTF-8}
---
```{r, include = FALSE}
knitr::opts_chunk$set(
collapse = TRUE,
comment = "#>"
)
```
```{r cleaning, include = FALSE, echo = TRUE}
library(serofoi)
```
The current version of ***serofoi*** supports three different models for estimating the *Force-of-Infection (FoI)*, including constant and time-varying trajectories. For fitting the model to the seroprevalence data we use a suit of bayesian models that include prior and upper prior distributions
## What is the Force-of-Infection
The force of infection, also known as the hazard rate or the infection pressure, is a key concept in mathematical modelling of infectious diseases. It represents the rate at which susceptible individuals become infected, given their exposure to a pathogen. In simple terms, the force of infection quantifies the risk of a susceptible individual becoming infected over a period of time. It is usually expressed as a rate per unit of time (e.g., per day or per year).
## Constant vs Time-varying FoI
The *FoI* is one of the most important parameters in epidemiology, but it is often incorrectly assumed to be constant over time. Identifying whether the *FoI* follows a constant or a time-varying trend can be important in the identification and characterization of the spread of disease. In Table 1 there is a summary of the models currently supported by serofoi.
::: l-body-outset
| Model Option | Description and usage |
|-----------------|:-----------------------------------------------------------:|
| `constant` | Constant FoI |
| `tv_normal` | Time-varying normal FoI: slow change in FoI |
| `tv_normal_log` | Time-varying normal-log FoI: fast epidemic change in FoI |
:::
Table 1. Model options and descriptions.
## Model 1. Constant Force-of-Infection (endemic model)
The *endemic constant model* is a simple mathematical model used in epidemiology to describe the seroprevalence of an infectious disease within a population, as a product of a long-term transmission.
For a constant FoI endemic model, the rate of infection acquisition $\lambda$ is constant over time for each trajectory, and the seroprevalence $P$ behaves as a cumulative process increasing monotonically with age. For the seroprevalence at age $a$ and time $t$, we have:
$$
P(a,t) = 1-\exp\left(-\lambda a\right)
$$
The number of positive cases follows a binomial distribution, where $n$ is the number of trials (size of the age group) and $P$ is the probability of successes (seroprevalence) for a certain age group:
$$
p(a,t) \sim binom(n(a,t), P(a,t))
$$
In ***serofoi***, for the constant model, the *FoI* ($\lambda$) is modelled within a Bayesian framework using a uniform prior distribution $\sim U(0,2)$. Future versions of the package may allow to choose different default distributions. This model can be implemented for the previously prepared dataset `data_test` by means of the `fit_seromodel` function specifying `fit_seromodel="constant"`.
The object `simdata_constant` contains a minimal simulated dataset that emulates an hypothetical endemic situation where the *FoI* is constant with value 0.2 and includes data for 250 samples of individuals between 2 and 47 years old with a number of trials $n=5$. The following code shows how to implement the constant model to this simulated serosurvey:
```{r model_1, include = TRUE, echo = TRUE, results="hide", errors = FALSE, warning = FALSE, message = FALSE, fig.width=4, fig.asp=1.5, fig.align="center", out.width ="50%", fig.keep="none"}
data("simdata_constant")
serodata_constant <- prepare_serodata(simdata_constant)
model_1 <- fit_seromodel(
serodata = serodata_constant,
foi_model = "constant",
iter = 800
)
plot_seromodel(
model_1,
serodata = serodata_constant,
size_text = 6
)
```
```{r model_1_plot, include = TRUE, echo = FALSE, results="hide", errors = FALSE, warning = FALSE, message = FALSE, fig.width=4, fig.asp=1.5, fig.align="center", out.width ="50%", fig.keep="all"}
foi_sim_constant <- rep(0.02, 50)
model_1_plot <- plot_seromodel(
model_1,
serodata = serodata_constant,
foi_sim = foi_sim_constant,
size_text = 6
)
plot(model_1_plot)
```
Figure 1. Constant serofoi model plot. Simulated (red) vs modelled (blue) *FoI*.
In this case, 800 iterations are enough to ensure convergence. The `plot_seromodel` method provides a visualisation of the results, including a summary where the expected log pointwise predictive density (`elpd`) and its standard error (`se`) are shown. We say that a model converges if all the R-hat estimates are below 1.1.
## Time-varying FoI models
For the time-varying *FoI* models, the probability for a case to be positive at age a at time $t$ also follows a binomial distribution, as described above. However, the seroprevalence is obtained from a cumulative of the yearly-varying values of the *FoI* over time:
$$
P(a,t) = 1 - \exp\left(-\sum_{i=t-a+1}^{t}\lambda_i\right)
$$
The corresponding serosurvey completed at time $t_{sur}$ is informative for the interval $[t_{sur}-a_{max}, t_{sur}]$.
## Model 2. Time-varying FoI - Slow Time-Varying FoI
The *time-varying slow normal model* relies on the following prior distributions for the *FoI* to describe the spread of a given infectious disease within a population over time:
$$
\lambda(t)\sim normal(\lambda(t-1), \sigma) \\
\lambda(t=1) \sim normal(0, 1)
$$
The object `simdata_sw_dec` contains a minimal simulated dataset that emulates a situation where the *FoI* follows a stepwise decreasing tendency (*FoI* panel in Fig. 2). The simulated dataset contains information about 250 samples of individuals between 2 and 47 years old with a number of trials $n=5$. The following code shows how to implement the slow time-varying normal model to this simulated serosurvey:
```{r model_2, include = TRUE, echo = TRUE, results="hide", errors = FALSE, warning = FALSE, message = FALSE, fig.width=4, fig.asp=1.5, fig.align="center", out.width ="50%", fig.keep="none"}
data("simdata_sw_dec")
serodata_sw_dec <- prepare_serodata(simdata_sw_dec)
model_2 <- fit_seromodel(
serodata = serodata_sw_dec,
foi_model = "tv_normal",
iter = 1500
)
plot_seromodel(model_2,
serodata = serodata_sw_dec,
size_text = 6
)
```
```{r model_2_plot, include = TRUE, echo = FALSE, results="hide", errors = FALSE, warning = FALSE, message = FALSE, fig.width=4, fig.asp=1.5, fig.align="center", out.width ="50%", fig.keep="all"}
no_transm <- 0.0000000001
foi_sim_sw_dec <- c(rep(0.2, 25), rep(0.1, 10), rep(no_transm, 17))
model_2_plot <- plot_seromodel(model_2,
serodata = serodata_sw_dec,
foi_sim = foi_sim_sw_dec,
size_text = 6
)
plot(model_2_plot)
```
Figure 2. Slow time-varying serofoi model plot. Simulated (red) vs modelled (blue) *FoI*.
The number of iterations required may depend on the number of years, reflected by the difference between the year of the serosurvey and the maximum age-class sampled.
## Model 3. Time-varying FoI - Fast Epidemic Model
The *time-varying fast epidemic model*, relies on normal prior distributions for the *FoI* in the logarithmic scale, i.e:
$$
\lambda(t)\sim normal(\log(\lambda(t-1)), \sigma) \\
\lambda(t=1) \sim normal(-6, 4)
$$
This is done in order to capture fast changes in the *FoI* trend. Importantly, the standard deviation parameter of this normal distribution of the *FoI* $\lambda(t)$ is set using an upper prior that follows a Cauchy distribution.
In order to test this model we use the minimal simulated dataset contained in the `simdata_large_epi` object. This dataset emulates a hypothetical situation where a three-year epidemic occurs between 2032 and 2035. The simulated serosurvey tests 250 individuals from 0 to 50 years of age in the year 2050. The implementation of the fast epidemic model can be obtained running the following lines of code:
```{r model_3, include = TRUE, echo = TRUE, results="hide", errors = FALSE, warning = FALSE, message = FALSE, fig.width=4, fig.asp=1.5, fig.align="center", out.width ="50%", fig.keep="none"}
data("simdata_large_epi")
serodata_large_epi <- prepare_serodata(simdata_large_epi)
model_3 <- fit_seromodel(
serodata = serodata_large_epi,
foi_model = "tv_normal_log",
iter = 1500
)
model_3_plot <- plot_seromodel(model_3,
serodata = serodata_large_epi,
size_text = 6
)
plot(model_3_plot)
```
```{r model_3_plot, include = TRUE, echo = FALSE, results="hide", errors = FALSE, warning = FALSE, message = FALSE, fig.width=4, fig.asp=1.5, fig.align="center", out.width ="50%", fig.keep="all"}
no_transm <- 0.0000000001
big_outbreak <- 1.5
foi_sim_large_epi <- c(
rep(no_transm, 32),
rep(big_outbreak, 3),
rep(no_transm, 15)
)
model_3_plot <- plot_seromodel(model_3,
serodata = serodata_large_epi,
foi_sim = foi_sim_large_epi,
size_text = 6
)
plot(model_3_plot)
```
Figure 3. *Time-varying fast epidemic serofoi model* plot. Simulated (red) vs modelled (blue) *FoI*.
In Fig 3 we can see that the *fast epidemic serofoi model* is able to identify the large epidemic simulated on the `simdata_large_epi` dataset.
## Models Comparison
The statistical details of the three models are described in Table 2.
::: l-body-outset
| Model Option | Probability of positive case at age $a$ | Prior distribution | Upper priors |
|:----------------|:----------------------------------------|:----------------------------------------------------------|:------------------------|
| `constant` | $\sim binom(n(a,t), P(a,t))$ | $\lambda\sim uniform(0,2)$ | |
| `tv_normal` | $\sim binom(n(a,t), P(a,t))$ | $\lambda\sim normal(\lambda(t-1),\sigma)\\
\lambda(t=1)\sim normal(0,1)$ | $\sigma\sim Cauchy(0,1)$|
| `tv_normal_log` | $\sim binom(n(a,t), P(a,t))$ | $\lambda\sim normal(log(\lambda(t-1)),\sigma)\\
\lambda(t=1)\sim normal(-6,4)$ | $\sigma\sim Cauchy(0,1)$|
:::
Table 2. Statistical characteristics of ***serofoi***'s currently supported models for the FoI ($\lambda$). Here $n$ is the size of an age group $a$ at time-step $t$ and $P$ is its corresponding seroprevalence.
Above we showed that the fast epidemic model (`tv_normal_log`) is able to identify the large epidemic outbreak described by the `simdata_large_epi` dataset, which was simulated according to a step-wise decreasing *FoI* (red line in Fig 3).
Now, we would like to know whether this model actually fits this dataset better than the other available models in ***serofoi***. For this, we also implement both the endemic model (`constant`) and the slow time-varying normal model (`tv_normal`):
```{r model_comparison, include = FALSE, echo = TRUE, results="hide", errors = FALSE, warning = FALSE, message = FALSE }
model_1 <- fit_seromodel(
serodata = serodata_large_epi,
foi_model = "constant",
iter = 800
)
model_1_plot <- plot_seromodel(model_1,
serodata = serodata_large_epi,
size_text = 6
)
model_2 <- fit_seromodel(
serodata = serodata_large_epi,
foi_model = "tv_normal",
iter = 1500
)
model_2_plot <- plot_seromodel(model_2,
serodata = serodata_large_epi,
size_text = 6
)
```
Using the function `cowplot::plot_grid` we can visualise the results of the three models simultaneously:
```{r model_comparison_plot, include = TRUE, echo = TRUE, results="hide", errors = FALSE, warning = FALSE, message = FALSE, fig.width=5, fig.asp=1, fig.align="center", fig.keep="none"}
cowplot::plot_grid(model_1_plot, model_2_plot, model_3_plot,
nrow = 1, ncol = 3, labels = "AUTO"
)
```
```{r model_comparison_plot_, include = TRUE, echo = FALSE, results="hide", errors = FALSE, warning = FALSE, message = FALSE, fig.width=5, fig.asp=1, fig.align="center", fig.keep="all"}
size_text <- 6
model_1_plot <- plot_seromodel(
model_1,
serodata = serodata_large_epi,
foi_sim = foi_sim_large_epi,
size_text = size_text
)
model_2_plot <- plot_seromodel(
model_2,
serodata = serodata_large_epi,
foi_sim = foi_sim_large_epi,
size_text = size_text
)
model_3_plot <- plot_seromodel(
model_3,
serodata = serodata_large_epi,
foi_sim = foi_sim_large_epi,
size_text = size_text
)
cowplot::plot_grid(model_1_plot, model_2_plot, model_3_plot,
nrow = 1, ncol = 3, labels = "AUTO"
)
```
Figure 4. Model comparison between the three serofoi models for a large-epidemic simulated dataset.
A common criterion to decide what model fits the data the best is to choose the one with the larger `elpd`. According to this criterion, in this case the best model is the fast epidemic model, which is the only one that manages to identify the large epidemic (see the second row of panel C in Figure 4).
NOTE: Running the ***serofoi*** models for the first time on your local computer may take a few minutes for the rstan code to compile locally. However, once the initial compilation is complete, there is no further need for local compilation.
## References