-
Notifications
You must be signed in to change notification settings - Fork 1
/
README.Rmd
114 lines (98 loc) · 2.85 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
---
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%"
)
```
# epidp
<!-- badges: start -->
[![codecov](https://codecov.io/github/ben18785/epidp/graph/badge.svg?token=STG0INT235)](https://codecov.io/github/ben18785/epidp)
[![R-CMD-check](https://github.com/ben18785/epidp/actions/workflows/R-CMD-check.yaml/badge.svg)](https://github.com/ben18785/epidp/actions/workflows/R-CMD-check.yaml)
<!-- badges: end -->
The goal of `epidp` is to allow covariate information to inform estimates of the time-varying reproduction number, $R_t$.
## Installation
You can install the development version of epidp from [GitHub](https://github.com/) with:
``` r
# install.packages("devtools")
devtools::install_github("ben18785/epidp")
```
## Example
### Step function in $R_t$
We first generate case data assuming a step function for $R_t$.
```{r}
library(epidp)
library(ggplot2)
library(dplyr)
library(magrittr)
library(purrr)
library(tidyr)
rt_fun <- function(t) {
if (t <= 60) {
R <- 2
} else if (t <= 90) {
R <- 0.5
} else {
R <- 1
}
R
}
# simulation parameters
nt <- 200
mean_si <- 6.5
sd_si <- 4.03
i_0 <- 10
# data frame of outputs
epidemic_df <- simulate_renewal_epidemic(rt_fun, nt, mean_si, sd_si, i_0)
# plot
transform_factor <- 300
epidemic_df %>%
select(-c(w_dist, lambda_t)) %>%
mutate(R_t = R_t * transform_factor) %>%
pivot_longer(c(i_t, R_t)) %>%
ggplot(aes(x = t, y = value, colour = name)) +
geom_line() +
scale_y_continuous(
name = "Incidence (i_t)",
sec.axis = sec_axis(~ . / transform_factor, name = "Reproduction Number (R_t)")
) +
labs(x = "Time") +
theme_minimal() +
scale_color_brewer("Series", palette = "Dark2")
```
We now use a Stan version of EpiFilter to estimate the maximum a posteriori estimates
of $R_t$ and overlay these on top of the actual values. Note, these estimates do
not have uncertainty associated with them but the benefit of this is that estimation
is instantaneous. The estimates are close to the actual $R_t$ values after an initial
period when case counts are low.
```{r}
# fit model
fit <- fit_epifilter(
N = length(epidemic_df$i_t),
C = epidemic_df$i_t,
w = epidemic_df$w_dist,
is_sampling = FALSE,
as_vector = FALSE
)
# plot
R <- fit$par$R
epidemic_df %>%
mutate(estimated = R) %>%
rename(true = R_t) %>%
select(t, estimated, true) %>%
pivot_longer(c(estimated, true)) %>%
ggplot(aes(x = t, y = value)) +
geom_line(aes(colour = name)) +
scale_color_brewer("R_t", palette = "Dark2") +
ylab("R_t") +
xlab("Time") +
theme_minimal()
```
## Contributing guidelines
We welcome contributions from collaborators. Before doing so, we ask that you read
our [contributing guidelines](CONTRIBUTING.md) section.