-
Notifications
You must be signed in to change notification settings - Fork 3
/
README.Rmd
186 lines (125 loc) · 8.21 KB
/
README.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
---
output: github_document
---
<!-- README.md is generated from README.Rmd. Please edit that file -->
```{r, include = FALSE}
knitr::opts_chunk$set(
collapse = TRUE,
comment = "#>",
fig.path = "man/figures/README-",
out.width = "100%"
)
```
# Forecast case and sequence notifications using variant of concern strain dynamics
[](https://github.com/epiforecasts/forecast.vocs/actions/workflows/R-CMD-as-cran-check.yaml) [](https://app.codecov.io/gh/epiforecasts/forecast.vocs)
[](https://epiforecasts.r-universe.dev/) [](https://github.com/epiforecasts/forecast.vocs/blob/master/LICENSE.md/) [](https://github.com/epiforecasts/forecast.vocs/graphs/contributors)
[](https://zenodo.org/badge/latestdoi/383161374)
Contains models and tools to produce short-term forecasts for both case and sequence notifications assuming circulation of either one or two variants. Tools are also provided to allow the evaluation of the use of sequence data for short-term forecasts in both real-world settings and in user generated scenarios.
## Installation
### Installing the package
Install the stable development version of the package with:
```{r, eval = FALSE}
install.packages("forecast.vocs", repos = "https://epiforecasts.r-universe.dev")
```
Install the unstable development from GitHub using the following,
```{r, eval = FALSE}
remotes::install_github("epiforecasts/forecast.vocs", dependencies = TRUE)
```
### Installing CmdStan
If you don't already have CmdStan installed then, in addition to installing `forecast.vocs`, it is also necessary to install CmdStan using CmdStanR's
`install_cmdstan()` function to enable model fitting in `forecast.vocs`. A suitable C++ toolchain is also required. Instructions are provided in the [_Getting started with
CmdStanR_](https://mc-stan.org/cmdstanr/articles/cmdstanr.html) vignette. See the [CmdStanR documentation](https://mc-stan.org/cmdstanr/) for further details and support.
```{r, eval = FALSE}
cmdstanr::install_cmdstan()
```
## Quick start
This quick start uses data from Germany that includes COVID-19 notifications and sequences with sequences either being positive or negative for the Delta variant. It shows how to produce forecasts for both a one and two strain model for the 19th of June 2021 when the latest available data estimated that approximately 7% of COVID-19 were positive for the Delta variant. Note that estimated growth rates and reproduction numbers shown here have been rescaled using an assumed generation time of 5.5 days and a weakly informative prior centred around Delta being 50% more transmissible than non-Delta cases has been used.
```{r}
library(forecast.vocs)
options(mc.cores = 4)
obs <- filter_by_availability(
germany_covid19_delta_obs,
date = as.Date("2021-06-19")
)
obs
current_obs <- latest_obs(germany_covid19_delta_obs)
```
### Forecast all-in-one
Run a forecast for both one and two strain models (or optionally just one of these) using the `forecast()` function. This provides a wrapper around other package tooling to initialise, fit, and summarise forecasts. Multiple forecasts can be performed efficiently across dates and scenarios using `forecast_across_dates()` and `forecast_across_scenarios()`.
```{r, message = FALSE}
forecasts <- forecast(obs,
strains = c(1, 2), voc_scale = c(0.4, 0.2),
voc_label = "Delta", scale_r = 5.5 / 7,
adapt_delta = 0.99, max_treedepth = 15,
refresh = 0, show_messages = FALSE
)
forecasts
```
Summarised posterior and forecast estimates can be directly extracted from the output using the included `summary()` method. Alternatively, the complete
summarised posterior can be extracted using the following (which also has its own `summary()` and `plot()` options to make exploring the results easier).
See the documentation for further details.
```{r, eval = FALSE}
summary(forecasts, target = "posterior", type = "all")
```
Plot the posterior prediction for cases from the single strain model ("Overall"), from the two strain model ("Combined"), and the unobserved estimates for each strain.
```{r}
plot(forecasts, current_obs)
```
Plot the posterior prediction for the fraction of cases that have the Delta variant from the two strain model.
```{r}
plot(forecasts, current_obs, type = "voc_frac", voc_label = "Delta variant")
```
Plot the posterior estimate for the growth rate over the mean of the generation time for COVID-19 cases (here assumed to be 5.5 days).
```{r}
plot(forecasts, type = "growth")
```
Plot the posterior estimate for the effective reproduction number of Delta and non-Delta cases.
```{r}
plot(forecasts, type = "rt")
```
Plot the estimated transmission advantage for Delta vs non-Delta over time. Note that as here a model with a correlated variant growth rates has been used meaning that the difference between variants may vary over time (potentially suitable where imports or present or a variant has greater immune escape than the baseline). If it is likely that the growth advantage is constant over time consider the fixed scaling model. If the variants are likely to be circulating in independent populations consider the independent variant model. See the documentation and model definition vignette for details.
```{r}
plot(forecasts, type = "voc_advantage")
```
Alternatively a list of all plots with sensible defaults can be produced using the following.
```{r, eval = FALSE}
plot(forecasts, current_obs, type = "all", voc_label = "Delta variant")
```
A potentially useful optional integration with [`scoringutils`](https://epiforecasts.io/scoringutils/) package is also provided to streamline scoring forecasts using proper scoring rules. See the documentation of the `scoringutils` package for more on this but as a simple example below we calculate average scores both the single and two strain models for this single Forecast
```{r}
library(scoringutils)
library(knitr)
summary(forecasts, target = "forecast", type = "cases") |>
fv_score_forecast(current_obs) |>
summarise_scores(by = "strains") |>
kable()
```
### Step by step forecast
Rather than using the all-in-one `forecast()` function individual package functions can be used to produce a forecast as follows.
```{r, eval = FALSE}
dt <- fv_as_data_list(obs, horizon = 4)
model <- fv_model(strains = 2)
inits <- fv_inits(dt, strains = 2)
fit <- fv_sample(
data = dt, model = model, init = inits,
voc_scale = c(0.4, 0.2),
adapt_delta = 0.99, max_treedepth = 15,
refresh = 0, show_messages = FALSE
)
posterior <- fv_tidy_posterior(fit, voc_label = "Delta", scale_r = 5.5 / 7)
```
As for forecasts produced with the `forecast()` function summary plots and estimates can then be produced using `summary()`, and `plot()` methods. For example case forecasts can be returned using the following,
```{r, eval = FALSE}
summary(posterior, type = "cases", forecast = TRUE)
```
## Citation
If using `forecast.vocs` in your work please consider citing it using the following,
```{r, echo = FALSE}
citation("forecast.vocs")
```
## How to make a bug report or feature request
Please briefly describe your problem and what output you expect in an [issue](https://github.com/epiforecasts/forecast.vocs/issues). If you have a question, please don't open an issue. Instead, ask on our [Q and A page](https://github.com/epiforecasts/forecast.vocs/discussions/categories/q-a).
## Contributing
We welcome contributions and new contributors! We particularly appreciate help on priority problems in the [issues](https://github.com/epiforecasts/forecast.vocs/issues). Please check and add to the issues, and/or add a [pull request](https://github.com/epiforecasts/forecast.vocs/pulls).
## Code of Conduct
Please note that the `forecast.vocs` project is released with a [Contributor Code of Conduct](epiforecasts.io/forecast.vocs/CODE_OF_CONDUCT.html). By contributing to this project, you agree to abide by its terms.