Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Use incidence2 + MASS directly instead of i2extras #157

Merged
merged 6 commits into from
Jun 18, 2024
Merged

Use incidence2 + MASS directly instead of i2extras #157

Show file tree
Hide file tree
Changes from 5 commits
Commits
Show all changes
6 commits
Select commit Hold shift + click to select a range
b9e3de2
Update lockfile
Bisaloo Jun 17, 2024
9e6db72
Use MASS+incidence2 directly
Bisaloo Jun 17, 2024
d6e1af8
Rename child document and factor out package name
Bisaloo Jun 17, 2024
3a10213
Add Tim as ctb
Bisaloo Jun 17, 2024
77a9e72
Fix lints
Bisaloo Jun 17, 2024
8a9cfde
Remove paragraph about doubling and halving times
Bisaloo Jun 18, 2024
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
1 change: 1 addition & 0 deletions DESCRIPTION
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -6,6 +6,7 @@ Authors@R: c(
comment = c(ORCID = "0000-0002-4094-1476")),
person("Thibaut", "Jombart", role = "aut"),
person("Carmen", "Tamayo", role = "aut", comment = c(ORCID = "0000-0003-4184-2864")),
person("Tim", "Taylor", role = "ctb", comment = c(ORCID = "0000-0002-8587-7113")),
person("data.org", role = "fnd")
)
Description: A store of outbreak analytics pipelines using different
Expand Down
2 changes: 0 additions & 2 deletions inst/rmarkdown/templates/transmissibility/skeleton/rmdchunks/EpiEstim.Rmd
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,5 +1,3 @@
## _EpiEstim_

```{r lockfile-rt, include = FALSE, message = FALSE}
renv::use("[email protected]")
```
Expand Down
2 changes: 0 additions & 2 deletions inst/rmarkdown/templates/transmissibility/skeleton/rmdchunks/EpiNow2.Rmd
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,5 +1,3 @@
## _EpiNow2_

```{r lockfile-rt, include = FALSE, message = FALSE}
renv::use("[email protected]")
```
Expand Down
104 changes: 57 additions & 47 deletions ...ssibility/skeleton/rmdchunks/i2extras.Rmd → ...nsmissibility/skeleton/rmdchunks/MASS.Rmd
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,24 +1,17 @@
## _i2extras_

```{r lockfile-rt, include = FALSE, message = FALSE}
renv::use(
"[email protected]",
"[email protected]",
"[email protected]" # requires the libsodium system dependency
"[email protected]", # requires the libsodium system dependency
"[email protected]"
)
```

```{r rt-load-pkg}
library(i2extras)
```

### Explanations

The package *i2extras* is used to estimate daily growth rates '_r_', using Poisson
GLMs fitted to epidemic curves. These estimates can be biased by reporting
delays. Here, we exclude the last `r params$incomplete_days` days as incomplete,
and retain a total of `r params$r_estim_window` days to estimate the growth
rate.
We fit negative binomials GLMs to epidemic curves to estimate daily growth rates
'_r_'. These estimates can be biased by reporting delays. Here, we exclude the
last `r params$incomplete_days` days as incomplete, and retain a total of
`r params$r_estim_window` days to estimate the growth rate.

Growth rates can be interpreted as the percent change in daily
incidence. Positive numbers indicate growth, while negative numbers indicate a
Expand Down Expand Up @@ -46,32 +39,62 @@ intervals (CI)
halving time) for each `r group_var`; all results show point estimates and their 95% CI

```{r fig.height = 5 / 3 * n_groups}
# This code uses the {i2extras} function `fit_curve()` to fit a poisson model to
# the epicurves by group_var that represent incidence data for cases included in
# params$data_file
last_counts <- dat_i_day %>%
keep_last(params$r_estim_window + params$incomplete_days)

last_trends <- last_counts %>%
fit_curve(model = "poisson")
# version with negbin model needs more iterations to converge
# fit_curve(model = "negbin", control = glm.control(maxit = 1e3))
plot(last_trends,
angle = 45,
date_format = "%d %b %y",
alpha = 0.8,
n_breaks = 12,
nrow = n_groups,
show_cases = small_counts
)
nest(.key = "data") %>%
mutate(model = lapply(data, function(d) {
MASS::glm.nb(count ~ date_index, d, control = glm.control(maxit = 1e3))
}))

intervals <- last_trends %>%
mutate(result = Map(
function(data, model) {
data %>%
ciTools::add_ci(model, names = c("lower_ci", "upper_ci")) %>%
ciTools::add_pi(model, names = c("lower_pi", "upper_pi")) %>%
as_tibble()
},
data,
model
)) %>%
dplyr::select(any_of(c(group_var, "result"))) %>%
unnest(result)

plot(last_counts, angle = 45) +
geom_line(
aes(date_index, y = pred),
data = intervals,
inherit.aes = FALSE
) +
geom_ribbon(
aes(date_index, ymin = lower_pi, ymax = upper_pi),
alpha = 0.2,
data = intervals,
inherit.aes = FALSE,
fill = "#BBB67E"
)
```

```{r}
# This code uses the {i2extras} function `growth_rate()` to estimate and plot
# growth rates for the last params$r_estim_window + params$incomplete_days days
# included in params$data_file
last_trends %>%
growth_rate() %>%
growth_rates <- last_trends %>%
mutate(
growth = bind_rows(lapply(
model,
function(m) {
ci <- confint(m, 2L, level = 0.95)
tibble(r_lower = ci[1L], r = coef(m)[2L], r_upper = ci[2L])
}
))
) %>%
unpack(growth)

message(colnames(growth_rates))
growth_rates <- growth_rates %>%
dplyr::select(all_of(c(group_var, "r", "r_lower", "r_upper")))

growth_rates %>%
ggplot(aes(y = .data[[group_var]]), fill = custom_grey) +
geom_point(aes(x = r), color = dark_green) +
geom_errorbar(aes(xmin = r_lower, xmax = r_upper), color = dark_green) +
Expand All @@ -91,32 +114,19 @@ last_trends %>%
# This code generates a table where estimated growth rates and doubling/halving
# times are displayed per group_var over the last params$r_estim_window +
# params$incomplete_days days included in params$data_file
last_trends %>%
growth_rate() %>%
dplyr::select(-c(1, 3)) %>%
rename(period = growth_or_decay) %>%
growth_rates %>%
mutate(
"period type" = paste(period, "time"),
r = paste(round(r * 100, 1), "%"),
r_ci = sprintf(
"[%1.1f %% ; %1.1f %% ]",
r_lower * 100,
r_upper * 100
),
time = round(time, 1),
time_ci = sprintf(
"[%.1f ; %.1f]",
time_lower,
time_upper
)
) %>%
dplyr::select(
group_var,
"Daily growth rate" = r,
"Growth rate (95% CI)" = r_ci,
"period type",
"Duration (days)" = time,
"Duration (95% CI)" = time_ci
"Growth rate (95% CI)" = r_ci
) %>%
set_names(toupper) %>%
kbl() %>%
Expand Down
2 changes: 0 additions & 2 deletions inst/rmarkdown/templates/transmissibility/skeleton/rmdchunks/R0.Rmd
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,5 +1,3 @@
## _R0_

```{r lockfile-rt, include = FALSE, message = FALSE}
renv::use("[email protected]")
```
Expand Down
112 changes: 52 additions & 60 deletions inst/rmarkdown/templates/transmissibility/skeleton/skeleton.Rmd
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -41,29 +41,19 @@ params:
rt_estimator:
label: "Which R package to use for Rt estimation"
value: "EpiEstim"
choices: ["EpiEstim", "EpiNow2", "i2extras", "R0"]
choices: ["EpiEstim", "EpiNow2", "MASS", "R0"]
bibliography:
- grateful-refs.bib
---

```{r lockfile, include = FALSE, message = FALSE}
renv::use(
"[email protected]",
"[email protected]",
"[email protected]",
"[email protected]",
"[email protected]",
"[email protected]",
"[email protected]",
"[email protected]",
"[email protected]",
"[email protected]",
"[email protected]",
"[email protected]",
"[email protected]",
"[email protected]",
"[email protected]",
"[email protected]",
"[email protected]",
"[email protected]",
"[email protected]",
"[email protected]",
"[email protected]",
"[email protected]",
"[email protected]",
Expand All @@ -73,101 +63,101 @@ renv::use(
"[email protected]",
"[email protected]",
"[email protected]",
"data.table@1.14.8",
"data.table@1.15.4",
"[email protected]",
"[email protected]",
"distributional@0.3.2",
"distributional@0.4.0",
"[email protected]",
"[email protected]",
"epiverse-trace/epiparameter@328706e", # nolint
"evaluate@0.23",
"evaluate@0.24.0",
"[email protected]",
"[email protected].1",
"fastmap@1.1.1",
"[email protected].2",
"fastmap@1.2.0",
"[email protected]",
"[email protected]",
"[email protected]",
"[email protected].3",
"[email protected].4",
"[email protected]",
"[email protected].0",
"[email protected].1",
"[email protected]",
"[email protected]",
"grates@1.1.0",
"[email protected].4",
"grates@1.2.1",
"[email protected].5",
"[email protected]",
"highr@0.10",
"highr@0.11",
"[email protected]",
"[email protected]",
"[email protected]",
"[email protected]",
"[email protected]",
"[email protected]",
"[email protected]",
"[email protected]",
"[email protected]",
"[email protected]",
"[email protected]",
"[email protected]",
"kableExtra@1.3.4",
"knitr@1.45",
"kableExtra@1.4.0",
"knitr@1.47",
"[email protected]",
"[email protected]",
"lattice@0.21-8",
"lattice@0.22-6",
"[email protected]",
"[email protected].0",
"[email protected].3",
"[email protected]",
"[email protected]",
"[email protected]",
"[email protected]",
"[email protected]",
"[email protected]",
"[email protected]",
"[email protected].0",
"[email protected]163",
"[email protected].1",
"[email protected]164",
"[email protected]",
"[email protected]",
"[email protected]",
"[email protected]",
"[email protected]",
"[email protected]",
"[email protected]",
"[email protected]",
"[email protected]",
"[email protected]",
"[email protected]",
"[email protected]",
"[email protected]",
"[email protected]",
"[email protected]",
"[email protected]",
"[email protected]",
"[email protected]",
"[email protected]",
"[email protected]",
"[email protected]",
"[email protected]",
"[email protected]",
"[email protected]",
"[email protected]",
"[email protected]",
"[email protected]",
"[email protected]",
"[email protected]",
"[email protected]",
"[email protected]",
"[email protected]",
"[email protected]",
"[email protected]",
"[email protected]",
"[email protected]",
"[email protected]",
"[email protected]",
"[email protected]",
"[email protected]",
"[email protected]",
"[email protected]",
"[email protected]",
"[email protected]",
"[email protected].3",
"[email protected].4",
"[email protected]",
"[email protected]",
"[email protected]",
"[email protected]",
"systemfonts@1.0.6",
"systemfonts@1.1.0",
"[email protected]",
"[email protected].0",
"[email protected].0",
"timechange@0.2.0",
"tinytex@0.50",
"[email protected].1",
"[email protected].1",
"timechange@0.3.0",
"tinytex@0.51",
"epiverse-trace/[email protected]", # nolint
"[email protected]",
"[email protected]",
"[email protected]",
"[email protected]",
"[email protected]",
"[email protected]",
"[email protected]",
"writexl@1.4.2",
"xfun@0.42",
"writexl@1.5.0",
"xfun@0.44",
"[email protected]",
"[email protected]",
"[email protected]"
Expand Down Expand Up @@ -479,6 +469,8 @@ dat_i_day <- dat_raw %>%
keep_first(n_distinct(.$date_index) - params$incomplete_days)
```

## _`r params$rt_estimator`_

```{r estimate-rt, child=paste0("rmdchunks/", params$rt_estimator, ".Rmd")}
```

Expand Down
Loading