Skip to content

Commit e81a326

Browse files
committed
README rewrite first draft
1 parent abcdd0b commit e81a326

12 files changed

+585
-235
lines changed

DEVELOPMENT.md

+10-4
Original file line numberDiff line numberDiff line change
@@ -32,13 +32,19 @@ Commands for developing the documentation site:
3232
# Basic build and preview
3333
R -e 'pkgdown::clean_site()'
3434
R -e 'devtools::document()'
35-
R -e 'pkgdown::build_site()'
35+
R -e 'pkgdown::build_site(lazy = TRUE, examples = FALSE, devel = TRUE, preview = FALSE)'
3636
```
3737

3838
If you work without R Studio and want to iterate on documentation, you might
39-
find [this
40-
script](https://gist.github.com/gadenbuie/d22e149e65591b91419e41ea5b2e0621)
41-
helpful.
39+
find `Rscript inst/pkgdown-watch.R` helpful to keep a live updating version of the website. Note that you need to have `c("pkgdown", "servr", "devtools", "here", "cli", "fs")` installed.
40+
41+
### Index/homepage
42+
43+
because we are using an `RMD` to make the index, figures are sometimes not added or updated. When in doubt, run `pkgdown::clean_cache()`, `pgkdown::clean_site()`, and delete the following directories/files (path relative to the project directory):
44+
- `readme_files`
45+
- `readme_cache`
46+
- `docs/dev/reference/figures/`
47+
- `man/figures/`
4248

4349
## Versioning
4450

README.Rmd

+241-67
Original file line numberDiff line numberDiff line change
@@ -5,26 +5,95 @@ output: github_document
55
<!-- README.md is generated from README.Rmd. Please edit that file -->
66

77
```{r, include = FALSE}
8-
options(width = 76)
98
knitr::opts_chunk$set(
10-
collapse = TRUE,
11-
comment = "#>",
129
fig.path = "man/figures/README-",
13-
out.width = "100%"
10+
digits = 3,
11+
comment = "#>",
12+
collapse = TRUE,
13+
cache = TRUE,
14+
dev.args = list(bg = "transparent"),
15+
dpi = 300,
16+
cache.lazy = FALSE,
17+
out.width = "90%",
18+
fig.align = "center",
19+
fig.width = 9,
20+
fig.height = 6
21+
)
22+
ggplot2::theme_set(ggplot2::theme_bw())
23+
options(
24+
dplyr.print_min = 6,
25+
dplyr.print_max = 6,
26+
pillar.max_footer_lines = 2,
27+
pillar.min_chars = 15,
28+
stringr.view_n = 6,
29+
pillar.bold = TRUE,
30+
width = 77
1431
)
1532
```
33+
```{r pkgs, include=FALSE, echo=FALSE}
34+
library(epipredict)
35+
library(epidatr)
36+
library(data.table)
37+
library(dplyr)
38+
library(tidyr)
39+
library(ggplot2)
40+
library(magrittr)
41+
library(purrr)
42+
library(scales)
43+
```
44+
45+
```{r coloration, include=FALSE, echo=FALSE}
46+
base <- "#002676"
47+
primary <- "#941120"
48+
secondary <- "#f9c80e"
49+
tertiary <- "#177245"
50+
fourth_colour <- "#A393BF"
51+
fifth_colour <- "#2e8edd"
52+
colvec <- c(base = base, primary = primary, secondary = secondary,
53+
tertiary = tertiary, fourth_colour = fourth_colour,
54+
fifth_colour = fifth_colour)
55+
library(epiprocess)
56+
suppressMessages(library(tidyverse))
57+
theme_update(legend.position = "bottom", legend.title = element_blank())
58+
delphi_pal <- function(n) {
59+
if (n > 6L) warning("Not enough colors in this palette!")
60+
unname(colvec)[1:n]
61+
}
62+
scale_fill_delphi <- function(..., aesthetics = "fill") {
63+
discrete_scale(aesthetics = aesthetics, palette = delphi_pal, ...)
64+
}
65+
scale_color_delphi <- function(..., aesthetics = "color") {
66+
discrete_scale(aesthetics = aesthetics, palette = delphi_pal, ...)
67+
}
68+
scale_colour_delphi <- scale_color_delphi
69+
```
1670

17-
# epipredict
71+
# Epipredict
1872

1973
<!-- badges: start -->
2074
[![R-CMD-check](https://github.com/cmu-delphi/epipredict/actions/workflows/R-CMD-check.yaml/badge.svg)](https://github.com/cmu-delphi/epipredict/actions/workflows/R-CMD-check.yaml)
2175
<!-- badges: end -->
2276

23-
**Note:** This package is currently in development and may not work as expected. Please file bug reports as issues in this repo, and we will do our best to address them quickly.
77+
Epipredict is a framework for building transformation and forecasting pipelines
78+
for epidemiological and other panel time-series datasets.
79+
In addition to tools for building forecasting pipelines, it contains a number of
80+
"canned" forecasters meant to run with little modification as an easy way to get
81+
started forecasting.
82+
83+
It is designed to work well with
84+
[`epiprocess`](https://cmu-delphi.github.io/epiprocess/), a utility for handling
85+
various time series and geographic processing tools in an epidemiological
86+
context.
87+
Both of the packages are meant to work well with the panel data provided by
88+
[`epidatr`](https://cmu-delphi.github.io/epidatr/).
89+
90+
If you are looking for more detail beyond the package documentation, see our
91+
[forecasting book](https://cmu-delphi.github.io/delphi-tooling-book/).
2492

2593
## Installation
2694

27-
To install (unless you're making changes to the package, use the stable version):
95+
To install (unless you're planning on contributing to package development, we
96+
suggest using the stable version):
2897

2998
```r
3099
# Stable version
@@ -33,94 +102,199 @@ pak::pkg_install("cmu-delphi/epipredict@main")
33102
# Dev version
34103
pak::pkg_install("cmu-delphi/epipredict@dev")
35104
```
105+
The documentation for the stable version is at <https://cmu-delphi.github.io/epipredict>, while the development version is at <https://cmu-delphi.github.io/epipredict/dev>.
36106

37-
## Documentation
38-
39-
You can view documentation for the `main` branch at <https://cmu-delphi.github.io/epipredict>.
40107

41-
## Goals for `epipredict`
108+
## Motivating example
42109

43-
**We hope to provide:**
44-
45-
1. A set of basic, easy-to-use forecasters that work out of the box. You should be able to do a reasonably limited amount of customization on them. For the basic forecasters, we currently provide:
46-
* Baseline flatline forecaster
47-
* Autoregressive forecaster
48-
* Autoregressive classifier
49-
* CDC FluSight flatline forecaster
50-
2. A framework for creating custom forecasters out of modular components. There are four types of components:
51-
* Preprocessor: do things to the data before model training
52-
* Trainer: train a model on data, resulting in a fitted model object
53-
* Predictor: make predictions, using a fitted model object
54-
* Postprocessor: do things to the predictions before returning
55-
56-
**Target audiences:**
110+
To demonstrate the kind of forecast epipredict can make, say we're predicting COVID deaths per 100k for each state on
111+
```{r fc_date}
112+
forecast_date <- as.Date("2021-08-01")
113+
```
114+
Below the fold, we construct this dataset as an `epiprocess::epi_df` from JHU data.
57115

58-
* Basic. Has data, calls forecaster with default arguments.
59-
* Intermediate. Wants to examine changes to the arguments, take advantage of
60-
built in flexibility.
61-
* Advanced. Wants to write their own forecasters. Maybe willing to build up
62-
from some components.
116+
<details>
117+
<summary> Creating the dataset using `{epidatr}` and `{epiprocess}` </summary>
118+
This dataset can be found in the package as <TODO DOESN'T EXIST>; we demonstrate some of the typically ubiquitous cleaning operations needed to be able to forecast.
119+
First we pull both jhu-csse cases and deaths from [`{epidatr}` package](https://cmu-delphi.github.io/epidatr/):
120+
```{r case_death}
121+
cases <- pub_covidcast(
122+
source = "jhu-csse",
123+
signals = "confirmed_incidence_prop",
124+
time_type = "day",
125+
geo_type = "state",
126+
time_values = epirange(20200601, 20220101),
127+
geo_values = "*") |>
128+
select(geo_value, time_value, case_rate = value)
63129
64-
The Advanced user should find their task to be relatively easy. Examples of
65-
these tasks are illustrated in the [vignettes and articles](https://cmu-delphi.github.io/epipredict).
130+
deaths <- pub_covidcast(
131+
source = "jhu-csse",
132+
signals = "deaths_incidence_prop",
133+
time_type = "day",
134+
geo_type = "state",
135+
time_values = epirange(20200601, 20220101),
136+
geo_values = "*") |>
137+
select(geo_value, time_value, death_rate = value)
138+
cases_deaths <-
139+
full_join(cases, deaths, by = c("time_value", "geo_value")) |>
140+
as_epi_df(as_of = as.Date("2022-01-01"))
141+
plot_locations <- c("ca", "ma", "ny", "tx")
142+
# plotting the data as it was downloaded
143+
cases_deaths |>
144+
filter(geo_value %in% c("ca", "ma", "ny", "tx")) |>
145+
pivot_longer(cols = c("case_rate", "death_rate"), names_to = "source") |>
146+
ggplot(aes(x = time_value, y = value)) +
147+
geom_line() +
148+
facet_grid(source ~ geo_value, scale = "free") +
149+
scale_x_date(date_breaks = "3 months", date_labels = "%Y %b") +
150+
theme(axis.text.x = element_text(angle = 90, hjust = 1))
151+
```
152+
As with basically any dataset, there is some cleaning that we will need to do to make it actually usable; we'll use some utilities from [`{epiprocess}`](https://cmu-delphi.github.io/epiprocess/) for this.
153+
First, to eliminate some of the noise coming from daily reporting, we do 7 day averaging over a trailing window[^1]:
66154

67-
See also the (in progress) [Forecasting Book](https://cmu-delphi.github.io/delphi-tooling-book/).
155+
[^1]: This makes it so that any given day of the new dataset only depends on the previous week, which means that we avoid leaking future values when making a forecast.
68156

69-
## Intermediate example
157+
```{r smooth}
158+
cases_deaths <-
159+
cases_deaths |>
160+
group_by(geo_value) |>
161+
epi_slide(
162+
cases_7dav = mean(case_rate, na.rm = TRUE),
163+
death_rate_7dav = mean(death_rate, na.rm = TRUE),
164+
.window_size = 7
165+
) |>
166+
ungroup() |>
167+
mutate(case_rate = NULL, death_rate = NULL) |>
168+
rename(case_rate = cases_7dav, death_rate = death_rate_7dav)
169+
```
70170

71-
The package comes with some built-in historical data for illustration, but
72-
up-to-date versions of this could be downloaded with the
73-
[`{epidatr}` package](https://cmu-delphi.github.io/epidatr/)
74-
and processed using
75-
[`{epiprocess}`](https://cmu-delphi.github.io/epiprocess/).[^1]
171+
Then trimming outliers, most especially negative values:
172+
```{r outlier}
173+
cases_deaths <-
174+
cases_deaths |>
175+
group_by(geo_value) |>
176+
mutate(
177+
outlr_death_rate = detect_outlr_rm(time_value, death_rate, detect_negatives = TRUE),
178+
outlr_case_rate = detect_outlr_rm(time_value, case_rate, detect_negatives = TRUE)
179+
) |>
180+
unnest(cols = starts_with("outlr"), names_sep = "_") |>
181+
ungroup() |>
182+
mutate(
183+
death_rate = outlr_death_rate_replacement,
184+
case_rate = outlr_case_rate_replacement) |>
185+
select(geo_value, time_value, case_rate, death_rate)
186+
cases_deaths
187+
```
188+
</details>
76189

77-
[^1]: Other epidemiological signals for non-Covid related illnesses are also
78-
available with [`{epidatr}`](https://github.com/cmu-delphi/epidatr) which
79-
interfaces directly to Delphi's
80-
[Epidata API](https://cmu-delphi.github.io/delphi-epidata/)
190+
After having downloaded and cleaned the data in `cases_deaths`, we plot a subset
191+
of the states, noting the actual forecast date:
81192

82-
```{r epidf, message=FALSE}
83-
library(epipredict)
84-
covid_case_death_rates
193+
<details>
194+
<summary> Plot </summary>
195+
```{r plot_locs}
196+
plot_locations <- c("ca", "ma", "ny", "tx")
197+
forecast_date_label <-
198+
tibble(
199+
geo_value = rep(plot_locations, 2),
200+
source = c(rep("case_rate",4), rep("death_rate", 4)),
201+
dates = rep(forecast_date - 7*2, 2 * length(plot_locations)),
202+
heights = c(rep(150, 4), rep(1.0, 4))
203+
)
204+
processed_data_plot <-
205+
cases_deaths |>
206+
filter(geo_value %in% c("ca", "ma", "ny", "tx")) |>
207+
pivot_longer(cols = c("case_rate", "death_rate"), names_to = "source") |>
208+
ggplot(aes(x = time_value, y = value)) +
209+
geom_line() +
210+
facet_grid(source ~ geo_value, scale = "free") +
211+
geom_vline(aes(xintercept = forecast_date)) +
212+
geom_text(
213+
data = forecast_date_label, aes(x=dates, label = "forecast\ndate", y = heights), size = 3, hjust = "right") +
214+
scale_x_date(date_breaks = "3 months", date_labels = "%Y %b") +
215+
theme(axis.text.x = element_text(angle = 90, hjust = 1))
216+
```
217+
</details>
218+
```{r show-processed-data, warning=FALSE, echo=FALSE}
219+
processed_data_plot
85220
```
86221

87-
To create and train a simple auto-regressive forecaster to predict the death rate two weeks into the future using past (lagged) deaths and cases, we could use the following function.
88-
222+
To make a forecast, we will use a "canned" simple auto-regressive forecaster to predict the death rate four weeks into the future using past (lagged) deaths and cases
89223
```{r make-forecasts, warning=FALSE}
90-
two_week_ahead <- arx_forecaster(
91-
covid_case_death_rates,
224+
four_week_ahead <- arx_forecaster(
225+
cases_deaths |> filter(time_value <= forecast_date),
92226
outcome = "death_rate",
93227
predictors = c("case_rate", "death_rate"),
94228
args_list = arx_args_list(
95229
lags = list(c(0, 1, 2, 3, 7, 14), c(0, 7, 14)),
96-
ahead = 14
230+
ahead = 4 * 7
97231
)
98232
)
99-
two_week_ahead
233+
four_week_ahead
100234
```
101235

102236
In this case, we have used a number of different lags for the case rate, while
103-
only using 3 weekly lags for the death rate (as predictors). The result is both
237+
using zero, one and two weekly lags for the death rate (as predictors). `four_week_ahead` is both
104238
a fitted model object which could be used any time in the future to create
105239
different forecasts, as well as a set of predicted values (and prediction
106-
intervals) for each location 14 days after the last available time value in the
107-
data.
240+
intervals) for each location 28 days after the forecast date.
241+
Plotting the prediction intervals on our subset above[^2]:
242+
243+
[^2]: Alternatively, you could call `auto_plot(four_week_ahead)` to get the full collection of forecasts. This is too busy for the space we have for plotting here.
108244

109-
```{r print-model}
110-
two_week_ahead$epi_workflow
245+
<details>
246+
<summary> Plot </summary>
247+
This is the same kind of plot as `processed_data_plot` above, but with the past data narrowed somewhat
248+
```{r}
249+
narrow_data_plot <-
250+
cases_deaths |>
251+
filter(time_value > "2021-04-01") |>
252+
filter(geo_value %in% c("ca", "ma", "ny", "tx")) |>
253+
pivot_longer(cols = c("case_rate", "death_rate"), names_to = "source") |>
254+
ggplot(aes(x = time_value, y = value)) +
255+
geom_line() +
256+
facet_grid(source ~ geo_value, scale = "free") +
257+
geom_vline(aes(xintercept = forecast_date)) +
258+
geom_text(
259+
data = forecast_date_label, aes(x=dates, label = "forecast\ndate", y = heights), size = 3, hjust = "right") +
260+
scale_x_date(date_breaks = "3 months", date_labels = "%Y %b") +
261+
theme(axis.text.x = element_text(angle = 90, hjust = 1))
111262
```
112263

113-
The fitted model here involved preprocessing the data to appropriately generate
114-
lagged predictors, estimating a linear model with `stats::lm()` and then
115-
postprocessing the results to be meaningful for epidemiological tasks. We can
116-
also examine the predictions.
264+
Putting that together with a plot of the bands, and a plot of the median prediction.
265+
```{r plotting_forecast, warning=FALSE}
266+
epiworkflow <- four_week_ahead$epi_workflow
267+
restricted_predictions <-
268+
four_week_ahead$predictions |>
269+
filter(geo_value %in% plot_locations) |>
270+
rename(time_value = target_date, value = .pred) |>
271+
mutate(source = "death_rate")
272+
forecast_plot <-
273+
narrow_data_plot |>
274+
epipredict:::plot_bands(
275+
restricted_predictions,
276+
fill = "dodgerblue4") +
277+
geom_point(data = restricted_predictions, aes(y = .data$value), color = "orange")
278+
```
279+
</details>
117280

118-
```{r show-preds}
119-
two_week_ahead$predictions
281+
```{r show-single-forecast, warning=FALSE, echo=FALSE}
282+
forecast_plot
120283
```
284+
The orange dot gives the median prediction, while the blue intervals give the 25-75%, 10-90%, and 2.5%-97.5% inter-quantile ranges.
285+
For this particular day and these locations, the forecasts are relatively accurate, with the true data being within the 25-75% interval.
286+
A couple of things to note:
121287

122-
The results above show a distributional forecast produced using data through
123-
the end of 2021 for the 14th of January 2022. A prediction for the death rate
124-
per 100K inhabitants is available for every state (`geo_value`) along with a
125-
90% predictive interval.
288+
1. Our methods are primarily direct forecasters; this means we don't need to
289+
predict 1, 2,..., 27 days ahead to then predict 28 days ahead
290+
2. All of our existing engines are geo-pooled, meaning the training data is
291+
shared across geographies. This has the advantage of increasing the amount of
292+
available training data, with the restriction that the data needs to be on
293+
comparable scales, such as rates.
126294

295+
## Getting Help
296+
If you encounter a bug or have a feature request, feel free to file an [issue on
297+
our github page](https://github.com/cmu-delphi/epipredict/issues).
298+
For other
299+
questions, feel free to contact [Daniel]([email protected]), [David]([email protected]), [Dmitry]([email protected]), or
300+
[Logan]([email protected]), either via email or on the Insightnet slack.

0 commit comments

Comments
 (0)