diff --git a/.Rbuildignore b/.Rbuildignore index 8f951b3..368b0a4 100644 --- a/.Rbuildignore +++ b/.Rbuildignore @@ -8,3 +8,6 @@ ^.cache$ ^inst/pkg-structure$ ^_pkgdown.yml$ +^\.github$ +^_pkgdown\.yml$ +^pkgdown$ diff --git a/.github/.gitignore b/.github/.gitignore new file mode 100644 index 0000000..2d19fc7 --- /dev/null +++ b/.github/.gitignore @@ -0,0 +1 @@ +*.html diff --git a/.github/workflows/check-full.yaml b/.github/workflows/check-full.yaml new file mode 100644 index 0000000..7475bf8 --- /dev/null +++ b/.github/workflows/check-full.yaml @@ -0,0 +1,85 @@ +on: + push: + branches: + - master + pull_request: + branches: + - master + +name: R-CMD-check + +jobs: + R-CMD-check: + runs-on: ${{ matrix.config.os }} + + name: ${{ matrix.config.os }} (${{ matrix.config.r }}) + + strategy: + fail-fast: false + matrix: + config: + - {os: macOS-latest, r: 'devel'} + - {os: macOS-latest, r: '4.0'} + - {os: windows-latest, r: '4.0'} + - {os: ubuntu-16.04, r: '4.0', rspm: "https://demo.rstudiopm.com/all/__linux__/xenial/latest"} + - {os: ubuntu-16.04, r: '3.6', rspm: "https://demo.rstudiopm.com/all/__linux__/xenial/latest"} + - {os: ubuntu-16.04, r: '3.5', rspm: "https://demo.rstudiopm.com/all/__linux__/xenial/latest"} + + env: + R_REMOTES_NO_ERRORS_FROM_WARNINGS: true + RSPM: ${{ matrix.config.rspm }} + + steps: + - uses: actions/checkout@v2 + + - uses: r-lib/actions/setup-r@master + with: + r-version: ${{ matrix.config.r }} + + - uses: r-lib/actions/setup-pandoc@master + + - name: Query dependencies + run: | + install.packages('remotes') + saveRDS(remotes::dev_package_deps(dependencies = TRUE), ".github/depends.Rds", version = 2) + shell: Rscript {0} + + - name: Cache R packages + if: runner.os != 'Windows' + uses: actions/cache@v1 + with: + path: ${{ env.R_LIBS_USER }} + key: ${{ runner.os }}-r-${{ matrix.config.r }}-1-${{ hashFiles('.github/depends.Rds') }} + restore-keys: ${{ runner.os }}-r-${{ matrix.config.r }}-1- + + - name: Install system dependencies + if: runner.os == 'Linux' + env: + RHUB_PLATFORM: linux-x86_64-ubuntu-gcc + run: | + sudo add-apt-repository -y "ppa:marutter/rrutter3.5" + sudo add-apt-repository -y "ppa:marutter/c2d4u3.5" + sudo apt update + sudo apt install r-cran-rstan + Rscript -e "remotes::install_github('r-hub/sysreqs')" + sysreqs=$(Rscript -e "cat(sysreqs::sysreq_commands('DESCRIPTION'))") + sudo -s eval "$sysreqs" + + - name: Install dependencies + run: | + remotes::install_deps(dependencies = TRUE) + remotes::install_cran("rcmdcheck") + shell: Rscript {0} + + - name: Check + env: + _R_CHECK_CRAN_INCOMING_: false + run: rcmdcheck::rcmdcheck(args = c("--no-manual", "--as-cran"), error_on = "warning", check_dir = "check") + shell: Rscript {0} + + - name: Upload check results + if: failure() + uses: actions/upload-artifact@master + with: + name: ${{ runner.os }}-r${{ matrix.config.r }}-results + path: check diff --git a/.github/workflows/pkgdown.yaml b/.github/workflows/pkgdown.yaml new file mode 100644 index 0000000..5e9a811 --- /dev/null +++ b/.github/workflows/pkgdown.yaml @@ -0,0 +1,44 @@ +on: + push: + branches: master + +name: pkgdown + +jobs: + pkgdown: + runs-on: macOS-latest + env: + GITHUB_PAT: ${{ secrets.GITHUB_TOKEN }} + steps: + - uses: actions/checkout@v2 + + - uses: r-lib/actions/setup-r@master + + - uses: r-lib/actions/setup-pandoc@master + + - name: Query dependencies + run: | + install.packages('remotes') + saveRDS(remotes::dev_package_deps(dependencies = TRUE), ".github/depends.Rds", version = 2) + shell: Rscript {0} + + - name: Cache R packages + uses: actions/cache@v1 + with: + path: ${{ env.R_LIBS_USER }} + key: macOS-r-4.0-1-${{ hashFiles('.github/depends.Rds') }} + restore-keys: macOS-r-4.0-1- + + - name: Install dependencies + run: | + install.packages("remotes") + remotes::install_deps(dependencies = TRUE) + remotes::install_dev("pkgdown") + shell: Rscript {0} + + - name: Install package + run: R CMD INSTALL . + + - name: Deploy package + run: pkgdown::deploy_to_branch(new_process = FALSE) + shell: Rscript {0} diff --git a/.gitignore b/.gitignore index 35448b7..2f6dba8 100644 --- a/.gitignore +++ b/.gitignore @@ -2,3 +2,5 @@ .RData .Rproj.user inst/doc +.DS_Store +docs diff --git a/DESCRIPTION b/DESCRIPTION index 0bd1d54..1937c38 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -7,14 +7,14 @@ Authors@R: c( role = c("aut", "cre"), email = "contact@samabbott.co.uk", comment = c(ORCID = "0000-0001-8057-8037")), - person("Joel", "Hellewell", - email = "joel.hellewell@lshtm.ac.uk", - role = c("aut"), - comment = c(ORCID = "0000-0003-2683-0849")), person(given = "Nikos", family = "Bosse", role = c("aut"), email = "nikosbosse@gmail.com"), + person("Joel", "Hellewell", + email = "joel.hellewell@lshtm.ac.uk", + role = c("aut"), + comment = c(ORCID = "0000-0003-2683-0849")), person("Katharine", "Sherratt", email = "katharine.sherratt@lshtm.ac.uk", role = c("aut")), @@ -24,24 +24,38 @@ Authors@R: c( person("Robin", "Thompson", email = "robin.thompson@lshtm.ac.uk", role = c("aut")), + person("Aurelien", "Chateigner", + email = "aurelien.chateigner@inrae.fr", + role = c("aut")), + person("Sylvain", "Mareschal", + email = "mareschal@ovsa.fr", + role = c("aut")), + person("Andrea", "Rau", + email = "andrea.rau@inrae.fr", + role = c("aut")), + person("Nathalie", "Vialaneix", + email = "nathalie.vialaneix@inrae.fr", + role = c("aut")), + person(given = "Michael", + family = "DeWitt", + role = c("aut"), + email = "me.dewitt.jr@gmail.com", + comment = c(ORCID = "0000-0001-8940-1967")), person("Sebastian", "Funk", email = "sebastian.funk@lshtm.ac.uk", role = c("aut"))) -Description: To forecast the time-varying reproduction number and using this to forecast reported case counts. Includes +Description: To forecast the time-varying reproduction number and use this to forecast reported case counts. Includes tools to evaluate a range of models across samples and time series using proper scoring rules. License: MIT + file LICENSE Encoding: UTF-8 LazyData: true Roxygen: list(markdown = TRUE) -RoxygenNote: 7.0.2 +RoxygenNote: 7.1.0 Depends: R (>= 3.3.0) Imports: scoringutils, scoringRules, - bsts, - fabletools, - tsibble, purrr, furrr, magrittr, @@ -60,7 +74,14 @@ Remotes: Suggests: knitr, rmarkdown, + fabletools, + tsibble, + forecastHybrid, fable, + bsts, feasts, - future.apply + future.apply, + brms, + tidybayes, + testthat VignetteBuilder: knitr diff --git a/NAMESPACE b/NAMESPACE index 6258deb..2324183 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -1,12 +1,14 @@ # Generated by roxygen2: do not edit by hand export("%>%") +export(brms_model) export(bsts_model) export(compare_models) export(compare_timeseries) export(draw_from_si_prob) export(evaluate_model) export(fable_model) +export(forecastHybrid_model) export(forecast_cases) export(forecast_rt) export(iterative_case_forecast) @@ -23,10 +25,10 @@ export(summarise_forecast) export(summarise_scores) importFrom(HDInterval,hdi) importFrom(R.utils,withTimeout) -importFrom(bsts,bsts) -importFrom(bsts,predict.bsts) importFrom(cowplot,theme_cowplot) +importFrom(data.table,`:=`) importFrom(dplyr,arrange) +importFrom(dplyr,bind_cols) importFrom(dplyr,bind_rows) importFrom(dplyr,filter) importFrom(dplyr,group_by) @@ -40,8 +42,6 @@ importFrom(dplyr,select_if) importFrom(dplyr,slice) importFrom(dplyr,summarise) importFrom(dplyr,ungroup) -importFrom(fabletools,forecast) -importFrom(fabletools,model) importFrom(furrr,future_map) importFrom(furrr,future_pmap) importFrom(ggplot2,aes) @@ -67,9 +67,14 @@ importFrom(scoringutils,bias) importFrom(scoringutils,interval_score) importFrom(scoringutils,pit) importFrom(scoringutils,sharpness) +importFrom(stats,median) +importFrom(stats,quantile) +importFrom(stats,rpois) +importFrom(stats,rt) +importFrom(stats,sd) +importFrom(stats,setNames) importFrom(tibble,tibble) importFrom(tidyr,expand_grid) importFrom(tidyr,gather) importFrom(tidyr,spread) importFrom(tidyr,unnest) -importFrom(tsibble,tsibble) diff --git a/R/brms_model.R b/R/brms_model.R new file mode 100644 index 0000000..9b9691a --- /dev/null +++ b/R/brms_model.R @@ -0,0 +1,76 @@ +#' brms model wrapper +#' +#' Allows users to specify a model using the [brms::bf()] wrapper from `brms` +#' Note that `brms` and `tidybayes` must both be installed for this +#' model wrapper to be functional. +#' +#' @param y Numeric vector of time points to forecast +#' @param samples Numeric, number of samples to take. +#' @param model A `brms` model wrapped in the [brms::bf()] function +#' @param horizon Numeric, the time horizon over which to predict. +#' @param n_cores Numeric, the number of cores to use, default of 1 +#' @param n_chains Numeric, the number of chains to use, default of 4 +#' @param n_iter Numeric, the number of iterations in the sampler to use, +#' default of 4000 +#' @param ... additional arguments passed to `brms` (e.g. priors or family) +#' @return A dataframe of predictions (with columns representing the time horizon and rows representing samples). +#' @export +#' @importFrom data.table `:=` +#' @examples \dontrun{ +#' +#' ## Used on its own +#' ## Note: More iterations and chains should be used +#' library(brms) +#'brms_model(y = EpiSoon::example_obs_rts[1:10, ]$rt, +#' model = brms::bf(y ~ gp(time)), +#' samples = 10, horizon = 7, n_iter = 40, n_chains = 1, refresh =0) +#' +#' ## Used for forecasting +#' ## Note that the timeout parameter has been increased to allow +#' ## for the time for the code to be compiled +#' ## Note: More iterations and chains should be used +#' +#' forecast_rt(EpiSoon::example_obs_rts[1:10, ], +#' model = function(...){ +#' brms_model(model = brms::bf(y ~ gp(time)), n_iter = 40, n_chains = 1, ...)}, +#' horizon = 7, samples = 10, timeout = 300) +#'} + +brms_model <- function(y = NULL, samples = NULL, + horizon = NULL, model = NULL, n_cores = 1, + n_chains = 4, n_iter = 2000, ...) { + + check_suggests("brms") + check_suggests("tidybayes") + check_suggests("tsibble") + + ## Make input numeric into correct tsibble format + timeseries <- tsibble::tsibble(y = y, time = 1:length(y), index = time) + + ## Fit the model + fit <- brms::brm(formula = model, data = timeseries, + chains = n_chains, iter = n_iter, + cores = n_cores, ...) + + # Create Prediction Data Frame + dat_new <- data.frame(time = (length(y)+1):(length(y)+horizon)) + + prediction <- tidybayes::add_fitted_draws(newdata = dat_new, + model = fit, n = samples) + + prediction <- dplyr::mutate(prediction, row_id = dplyr::row_number()) + + prediction <- dplyr::ungroup(prediction) + + prediction <- data.table::as.data.table(prediction) + + prediction <- data.table::dcast(prediction, row_id~time, value.var = ".value") + + prediction <- prediction[,row_id:=NULL] + + ## Extract samples and tidy format + samples <- as.data.frame(prediction) + + return(samples) + +} diff --git a/R/bsts_model.R b/R/bsts_model.R index c837253..af4e619 100644 --- a/R/bsts_model.R +++ b/R/bsts_model.R @@ -6,8 +6,9 @@ #' @param horizon Numeric, the time horizon over which to predict. #' @return A dataframe of predictions (with columns representing the time horizon and rows representing samples). #' @export -#' @importFrom bsts bsts predict.bsts -#' @examples +#' @examples \dontrun{ +#' +#' library(bsts) #' #' ## Used on its own #' bsts_model(y = EpiSoon::example_obs_rts[1:10, ]$rt, @@ -21,10 +22,14 @@ #' function(ss, y){ #' bsts::AddAr(ss, y = y, lags = 3)}, ...)}, #' horizon = 7, samples = 10) +#'} bsts_model <- function(y = NULL, samples = NULL, horizon = NULL, model = NULL) { + check_suggests("bsts") + + model <- model(list(), y) ## Fit the model @@ -44,9 +49,4 @@ bsts_model <- function(y = NULL, samples = NULL, return(samples) - - - - - } diff --git a/R/fable_model.R b/R/fable_model.R index 9d0d15a..a47650d 100644 --- a/R/fable_model.R +++ b/R/fable_model.R @@ -1,8 +1,9 @@ #' Fable model wrapper #' #' -#' @description Provides an interface for models from the `fable` package. Not the `ARIMA` model -#' requires the `feast` package. If `future` is being used `fable` will require `future.apply` in +#' @description Provides an interface for models from the `fable` package. +#' Note the `feasts::ARIMA` model requires the `feast` package. If `future` +#' is being used `fable` will require `future.apply` in #' order to not silently fail. #' #' @param model A `fable` model object. For models that use a formula interface time @@ -11,11 +12,9 @@ #' @return A dataframe of predictions (with columns representing the #' time horizon and rows representing samples). #' @export -#' @importFrom tsibble tsibble -#' @importFrom fabletools model forecast #' @importFrom purrr map #' @importFrom dplyr bind_rows -#' @examples +#' @examples \dontrun{ #' ## Used on its own #' fable_model(y = EpiSoon::example_obs_rts[1:10, ]$rt, #' model = fable::ARIMA(y ~ time), @@ -26,9 +25,15 @@ #' model = function(...){ #' fable_model(model = fable::ARIMA(y ~ time), ...)}, #' horizon = 7, samples = 10) +#'} fable_model <- function(y = NULL, samples = NULL, horizon = NULL, model = NULL) { +check_suggests("tsibble") +check_suggests("fable") +check_suggests("fabletools") +check_suggests("feasts") +check_suggests("future.apply") ## Make input numeric into correct tsibble format timeseries <- tsibble::tsibble(y = y, time = 1:length(y), index = time) diff --git a/R/forecastHybrid_model.R b/R/forecastHybrid_model.R new file mode 100644 index 0000000..740fa93 --- /dev/null +++ b/R/forecastHybrid_model.R @@ -0,0 +1,89 @@ +#' forecastHybrid model wrapper +#' +#' Allows users to forecast using ensembles from the `forecastHybrid` package. Note that +#' whilst weighted ensembles can be created this is not advised when samples > 1 as currently +#' samples are derived assuming a normal distribution using the upper and lower confidence intervals of the ensemble. +#' These confidence intervals are themselves either based on the unweighted mean of the ensembled +#' models or the maximum/minimum from the candiate models. Note that `forecastHybrid` must be installed for this +#' model wrapper to be functional. +#' @param y Numeric vector of time points to forecast +#' @param samples Numeric, number of samples to take. +#' @param horizon Numeric, the time horizon over which to predict. +#' @param model_params List of parameters to pass to `forecastHybrid::hybridModel`. +#' @param forecast_params List of parameters to pass to `forecastHybrid:::forecast.hybridModel`. +#' @return A dataframe of predictions (with columns representing the time horizon and rows representing samples). +#' @export +#' @importFrom purrr map2 +#' @importFrom dplyr bind_cols +#' @examples \dontrun{ +#' +#' library(forecastHybrid) +#' +#' ## Used on its own +#' forecastHybrid_model(y = EpiSoon::example_obs_rts$rt, +#' samples = 10, horizon = 7) +#' +#' +#'## Used with non-default arguments +#'## Note that with the current sampling from maximal confidence intervals model +#'## Weighting using cross-validation will only have an impact when 1 sample is used. +#'forecastHybrid_model(y = EpiSoon::example_obs_rts$rt, +#' samples = 1, horizon = 7, +#' model_params = list(cvHorizon = 7, windowSize = 7, +#' rolling = TRUE, models = "zeta")) +#' +#' +#' ## Used for forecasting +#' forecast_rt(EpiSoon::example_obs_rts, +#' model = EpiSoon::forecastHybrid_model, +#' horizon = 7, samples = 1) +#' +#'## Used for forcasting with non-default arguments +#'forecast_rt(EpiSoon::example_obs_rts, +#' model = function(...){EpiSoon::forecastHybrid_model( +#' model_params = list(models = "zte"), +#' forecast_params = list(PI.combination = "mean"), ...) +#' }, +#' horizon = 7, samples = 10) +#'} +forecastHybrid_model <- function(y = NULL, samples = NULL, + horizon = NULL, model_params = NULL, + forecast_params = NULL) { + + + check_suggests("forecastHybrid") + + + ## Fit the model + fitted_model <- suppressMessages( + suppressWarnings( + do.call(forecastHybrid::hybridModel, c(list(y = y, parallel = FALSE, + num.cores = 1), model_params)) + ) + ) + + ## Predict using the model + prediction <- do.call(forecastHybrid:::forecast.hybridModel, + c(list(object = fitted_model, h = horizon), + forecast_params)) + + ## Extract samples and tidy format + sample_from_model <- prediction + + if (samples == 1) { + sample_from_model <- data.frame(t(as.data.frame(sample_from_model$mean))) + rownames(sample_from_model) <- NULL + }else{ + mean <- as.numeric(prediction$mean) + upper <- prediction$upper[, ncol(prediction$upper)] + lower <- prediction$lower[, ncol(prediction$lower)] + sd <- (upper - lower) / 3.92 + sample_from_model <- purrr::map2(mean, sd, + ~ rnorm(samples, mean = .x, sd = .y)) + + sample_from_model <- dplyr::bind_cols(sample_from_model) + } + + return(sample_from_model) + +} diff --git a/R/forecast_rt.R b/R/forecast_rt.R index 255f346..641300d 100644 --- a/R/forecast_rt.R +++ b/R/forecast_rt.R @@ -25,7 +25,7 @@ #' horizon = 7, samples = 10) forecast_rt <- function(rts, model = model, horizon = 7, samples = 1000, - bound_rt = TRUE, timeout = 30) { + bound_rt = TRUE, timeout = 100) { ## Set up for model fitting y <- rts$rt diff --git a/R/utils.R b/R/utils.R new file mode 100644 index 0000000..b362855 --- /dev/null +++ b/R/utils.R @@ -0,0 +1,51 @@ +#' @importFrom stats median quantile rpois rt sd setNames + +globalVariables( + c("timeseries", + "bottom", + "cases", + "calibration", + "ci", + "crps", + "dss", + "forecast_date", + "hdi_50", + "hdi_90", + "horizon", + "index", + "infectiousness", + "iqr", + "logs", + "lower", + "model_name", + "obs_sample", + "row_id", + "score", + "time", + "top", + "upper", + "time", + "rts", + "value", + "y", + "model") +) + +check_suggests <- function(pkg_name, dev_message = NULL){ + if (!requireNamespace(pkg_name, quietly = TRUE)) { + msg <- sprintf("This function requires `%s` to work.", + pkg_name) + + if (!is.null(dev_message)) { + msg <- paste(msg, dev_message) + } else{ + msg <- paste(msg, "Please install it.\n") + } + stop(msg, call. = FALSE) + } +} + + + + + diff --git a/README.Rmd b/README.Rmd index b070d91..8d25e67 100644 --- a/README.Rmd +++ b/README.Rmd @@ -13,11 +13,10 @@ knitr::opts_chunk$set( # EpiSoon [![DOI](https://zenodo.org/badge/248311916.svg)](https://zenodo.org/badge/latestdoi/248311916) +[![R build status](https://github.com/epiforecasts/EpiSoon/workflows/R-CMD-check/badge.svg)](https://github.com/epiforecasts/EpiSoon) [![Build Status](https://travis-ci.com/epiforecasts/EpiSoon.svg?branch=master)](https://travis-ci.com/epiforecasts/EpiSoon) -*Warning: This package is a work in progress and is currently developed solely with the COVID-19 outbreak in mind. Breaking changes may occur and the authors cannot guarantee support.* - -**Aim:** To forecast the time-varying reproduction number and using this to forecast reported case counts. +**Aim:** To forecast the time-varying reproduction number and use this to forecast reported case counts. ## Installation diff --git a/README.md b/README.md index 98dd09d..a29bad5 100644 --- a/README.md +++ b/README.md @@ -1,14 +1,13 @@ # EpiSoon +[![DOI](https://zenodo.org/badge/248311916.svg)](https://zenodo.org/badge/latestdoi/248311916) +[![R build +status](https://github.com/epiforecasts/EpiSoon/workflows/R-CMD-check/badge.svg)](https://github.com/epiforecasts/EpiSoon) [![Build Status](https://travis-ci.com/epiforecasts/EpiSoon.svg?branch=master)](https://travis-ci.com/epiforecasts/EpiSoon) -*Warning: This package is a work in progress and is currently developed -solely with the COVID-19 outbreak in mind. Breaking changes may occur -and the authors cannot guarantee support.* - -**Aim:** To forecast the time-varying reproduction number and using this +**Aim:** To forecast the time-varying reproduction number and use this to forecast reported case counts. ## Installation @@ -84,70 +83,72 @@ forecasts <- EpiSoon::compare_timeseries(obs_rts, obs_cases, models, forecasts #> $forecast_rts -#> # A tibble: 791 x 12 -#> timeseries model forecast_date date horizon median mean sd bottom -#> -#> 1 Region 1 AR 3 2020-03-04 2020-03-05 1 2.12 2.00 0.599 0.386 -#> 2 Region 1 AR 3 2020-03-04 2020-03-06 2 1.92 1.90 0.503 0.733 -#> 3 Region 1 AR 3 2020-03-04 2020-03-07 3 1.74 1.67 0.722 0 -#> 4 Region 1 AR 3 2020-03-04 2020-03-08 4 1.63 1.57 0.735 0 -#> 5 Region 1 AR 3 2020-03-04 2020-03-09 5 1.49 1.59 0.620 0.871 -#> 6 Region 1 AR 3 2020-03-04 2020-03-10 6 1.44 1.42 0.881 0 -#> 7 Region 1 AR 3 2020-03-04 2020-03-11 7 1.41 1.38 0.879 0 -#> 8 Region 1 AR 3 2020-03-05 2020-03-06 1 2.20 2.20 0.00600 2.19 -#> 9 Region 1 AR 3 2020-03-05 2020-03-07 2 2.13 2.13 0.0145 2.12 -#> 10 Region 1 AR 3 2020-03-05 2020-03-08 3 2.07 2.07 0.0137 2.05 -#> # … with 781 more rows, and 3 more variables: lower , upper , +#> # A tibble: 798 x 12 +#> timeseries model forecast_date date horizon median mean sd bottom +#> +#> 1 Region 1 AR 3 2020-03-04 2020-03-05 1 0 1.08 2.15 0 +#> 2 Region 1 AR 3 2020-03-04 2020-03-06 2 0.203 0.955 1.59 0 +#> 3 Region 1 AR 3 2020-03-04 2020-03-07 3 0.232 0.790 1.26 0 +#> 4 Region 1 AR 3 2020-03-04 2020-03-08 4 0.259 1.70 2.69 0 +#> 5 Region 1 AR 3 2020-03-04 2020-03-09 5 0.238 1.40 2.35 0 +#> 6 Region 1 AR 3 2020-03-04 2020-03-10 6 0.862 0.830 0.880 0 +#> 7 Region 1 AR 3 2020-03-04 2020-03-11 7 1.22 1.92 2.95 0 +#> 8 Region 1 AR 3 2020-03-05 2020-03-06 1 0.378 1.18 1.88 0 +#> 9 Region 1 AR 3 2020-03-05 2020-03-07 2 0.536 1.87 3.64 0 +#> 10 Region 1 AR 3 2020-03-05 2020-03-08 3 0.758 1.87 3.44 0 +#> # … with 788 more rows, and 3 more variables: lower , upper , #> # top #> #> $rt_scores -#> # A tibble: 629 x 10 -#> timeseries model forecast_date date horizon dss crps logs bias -#> -#> 1 Region 1 AR 3 2020-03-04 2020-03-05 1 -0.914 0.0985 -0.312 0.4 -#> 2 Region 1 AR 3 2020-03-04 2020-03-06 2 -1.11 0.148 0.148 0.4 -#> 3 Region 1 AR 3 2020-03-04 2020-03-07 3 -0.265 0.223 0.558 0.4 -#> 4 Region 1 AR 3 2020-03-04 2020-03-08 4 -0.143 0.259 0.736 0.4 -#> 5 Region 1 AR 3 2020-03-04 2020-03-09 5 -0.413 0.297 0.822 0.4 -#> 6 Region 1 AR 3 2020-03-04 2020-03-10 6 0.172 0.340 1.08 0.4 -#> 7 Region 1 AR 3 2020-03-04 2020-03-11 7 0.217 0.347 1.09 0.4 -#> 8 Region 1 AR 3 2020-03-05 2020-03-06 1 -6.97 0.00727 -3.21 1 -#> 9 Region 1 AR 3 2020-03-05 2020-03-07 2 -7.38 0.0100 -2.70 0.100 -#> 10 Region 1 AR 3 2020-03-05 2020-03-08 3 -1.37 0.0278 -0.281 0 -#> # … with 619 more rows, and 1 more variable: sharpness +#> # A tibble: 630 x 14 +#> timeseries model forecast_date date horizon dss crps logs bias +#> +#> 1 Region 1 AR 3 2020-03-04 2020-03-05 1 1.76 1.42 13.8 -0.6 +#> 2 Region 1 AR 3 2020-03-04 2020-03-06 2 1.49 1.14 2.48 -0.8 +#> 3 Region 1 AR 3 2020-03-04 2020-03-07 3 1.65 1.15 2.49 -0.8 +#> 4 Region 1 AR 3 2020-03-04 2020-03-08 4 1.90 0.914 2.04 -0.400 +#> 5 Region 1 AR 3 2020-03-04 2020-03-09 5 1.70 0.999 2.45 -0.400 +#> 6 Region 1 AR 3 2020-03-04 2020-03-10 6 1.70 0.851 1.80 -0.8 +#> 7 Region 1 AR 3 2020-03-04 2020-03-11 7 2.06 0.564 1.62 -0.400 +#> 8 Region 1 AR 3 2020-03-05 2020-03-06 1 1.48 1.18 6.29 -0.6 +#> 9 Region 1 AR 3 2020-03-05 2020-03-07 2 2.48 1.00 2.29 -0.6 +#> 10 Region 1 AR 3 2020-03-05 2020-03-08 3 2.37 0.779 1.75 -0.6 +#> # … with 620 more rows, and 5 more variables: sharpness , +#> # calibration , median , iqr , ci #> #> $forecast_cases -#> # A tibble: 629 x 12 +#> # A tibble: 630 x 12 #> timeseries model forecast_date date horizon median mean sd bottom #> -#> 1 Region 1 AR 3 2020-03-04 2020-03-05 1 57 52.3 15.4 11 -#> 2 Region 1 AR 3 2020-03-04 2020-03-06 2 58 57.9 17.1 22 -#> 3 Region 1 AR 3 2020-03-04 2020-03-07 3 70 63.7 28.5 0 -#> 4 Region 1 AR 3 2020-03-04 2020-03-08 4 62.5 68.6 36.6 0 -#> 5 Region 1 AR 3 2020-03-04 2020-03-09 5 64.5 74.8 41.3 28 -#> 6 Region 1 AR 3 2020-03-04 2020-03-10 6 71 80.2 58.4 0 -#> 7 Region 1 AR 3 2020-03-04 2020-03-11 7 81.5 91.8 69.5 0 -#> 8 Region 1 AR 3 2020-03-05 2020-03-06 1 65 66.6 7.63 57 -#> 9 Region 1 AR 3 2020-03-05 2020-03-07 2 83 78.3 9.52 59 -#> 10 Region 1 AR 3 2020-03-05 2020-03-08 3 89.5 90.1 8.24 76 -#> # … with 619 more rows, and 3 more variables: lower , upper , +#> 1 Region 1 AR 3 2020-03-04 2020-03-05 1 0 32.3 66.2 0 +#> 2 Region 1 AR 3 2020-03-04 2020-03-06 2 7.5 36.3 62.3 0 +#> 3 Region 1 AR 3 2020-03-04 2020-03-07 3 6 33.9 70.3 0 +#> 4 Region 1 AR 3 2020-03-04 2020-03-08 4 6.5 60.3 129. 0 +#> 5 Region 1 AR 3 2020-03-04 2020-03-09 5 4.5 67.7 170. 0 +#> 6 Region 1 AR 3 2020-03-04 2020-03-10 6 16 38.7 50.0 0 +#> 7 Region 1 AR 3 2020-03-04 2020-03-11 7 38.5 53.9 64.1 0 +#> 8 Region 1 AR 3 2020-03-05 2020-03-06 1 13 46.4 74.3 0 +#> 9 Region 1 AR 3 2020-03-05 2020-03-07 2 17.5 88.6 178. 0 +#> 10 Region 1 AR 3 2020-03-05 2020-03-08 3 30.5 67.7 107. 0 +#> # … with 620 more rows, and 3 more variables: lower , upper , #> # top #> #> $case_scores -#> # A tibble: 629 x 11 +#> # A tibble: 630 x 15 #> timeseries model sample forecast_date date horizon dss crps logs #> -#> 1 Region 1 AR 3 1 2020-03-04 2020-03-05 1 5.90 5.07 3.78 -#> 2 Region 1 AR 3 1 2020-03-04 2020-03-06 2 6.44 8.47 4.09 -#> 3 Region 1 AR 3 1 2020-03-04 2020-03-07 3 7.40 12.8 4.40 -#> 4 Region 1 AR 3 1 2020-03-04 2020-03-08 4 8.02 18.8 4.88 -#> 5 Region 1 AR 3 1 2020-03-04 2020-03-09 5 8.44 23.6 5.07 -#> 6 Region 1 AR 3 1 2020-03-04 2020-03-10 6 9.25 34.7 5.45 -#> 7 Region 1 AR 3 1 2020-03-04 2020-03-11 7 9.68 41.5 5.59 -#> 8 Region 1 AR 3 1 2020-03-05 2020-03-06 1 4.74 4.38 3.57 -#> 9 Region 1 AR 3 1 2020-03-05 2020-03-07 2 5.55 4.89 3.42 -#> 10 Region 1 AR 3 1 2020-03-05 2020-03-08 3 6.43 7.49 3.76 -#> # … with 619 more rows, and 2 more variables: bias , sharpness +#> 1 Region 1 AR 3 1 2020-03-04 2020-03-05 1 8.52 40.8 35.8 +#> 2 Region 1 AR 3 1 2020-03-04 2020-03-06 2 8.54 36.3 5.79 +#> 3 Region 1 AR 3 1 2020-03-04 2020-03-07 3 9.06 55.9 8.89 +#> 4 Region 1 AR 3 1 2020-03-04 2020-03-08 4 9.74 58.9 6.74 +#> 5 Region 1 AR 3 1 2020-03-04 2020-03-09 5 10.3 77.9 10.3 +#> 6 Region 1 AR 3 1 2020-03-04 2020-03-10 6 12.4 78.2 6.29 +#> 7 Region 1 AR 3 1 2020-03-04 2020-03-11 7 11.7 91.1 7.27 +#> 8 Region 1 AR 3 1 2020-03-05 2020-03-06 1 8.65 38.1 7.27 +#> 9 Region 1 AR 3 1 2020-03-05 2020-03-07 2 10.3 43.1 6.13 +#> 10 Region 1 AR 3 1 2020-03-05 2020-03-08 3 9.36 38.6 5.66 +#> # … with 620 more rows, and 6 more variables: bias , sharpness , +#> # calibration , median , iqr , ci ``` - Plot an evaluation of Rt forecasts using iterative @@ -182,24 +183,20 @@ EpiSoon::plot_forecast_evaluation(forecasts$forecast_cases, obs_cases, c(7)) + ``` r EpiSoon::summarise_scores(forecasts$case_scores) -#> # A tibble: 15 x 9 -#> score model bottom lower median mean upper top sd -#> -#> 1 bias AR 3 0 0 0.100 0.256 0.4 1 0.354 -#> 2 bias ARIMA 0 0 0.0500 0.236 0.4 1 0.342 -#> 3 bias Semi-local linear t… 0 0 0.100 0.293 0.6 1 0.354 -#> 4 crps AR 3 4.08 10.9 19.7 35.1 37.3 159. 39.9 -#> 5 crps ARIMA 3.06 9.12 21.5 35.2 39.4 157. 39.7 -#> 6 crps Semi-local linear t… 2.95 9.01 20.0 32.8 35.7 170. 39.9 -#> 7 dss AR 3 5.02 6.89 9.48 15.8 16.6 69.1 17.3 -#> 8 dss ARIMA 4.57 6.57 9.57 16.9 17.0 78.5 19.9 -#> 9 dss Semi-local linear t… 4.53 6.70 8.57 13.3 13.0 64.5 15.1 -#> 10 logs AR 3 3.52 4.52 5.22 14.4 9.26 62.7 44.5 -#> 11 logs ARIMA 3.35 4.24 5.25 12.4 9.29 76.9 19.1 -#> 12 logs Semi-local linear t… 3.20 4.28 5.06 9.59 6.67 60.1 15.7 -#> 13 sharpness AR 3 4.60 10.4 15.6 18.8 23.7 51.0 13.1 -#> 14 sharpness ARIMA 4.45 9.64 14.1 15.6 19.8 33.2 7.78 -#> 15 sharpness Semi-local linear t… 5.36 11.1 16.3 20.4 24.5 55.3 15.1 +#> # A tibble: 27 x 9 +#> score model bottom lower median mean upper top sd +#> +#> 1 bias AR 3 -1.00e+0 -8.75e-1 -8.00e-1 -6.94e-1 -6.00e-1 -0.200 2.65e-1 +#> 2 bias ARIMA -1.00e+0 2.00e-1 8.00e-1 5.33e-1 1.00e+0 1 5.70e-1 +#> 3 bias Semi-l… -1.00e+0 3.00e-1 8.00e-1 5.37e-1 1.00e+0 1 5.63e-1 +#> 4 calib… AR 3 8.57e-5 8.57e-5 8.57e-5 2.83e-2 2.00e-3 0.504 1.26e-1 +#> 5 calib… ARIMA 8.57e-5 8.57e-5 8.57e-5 3.59e-2 1.50e-4 0.509 1.15e-1 +#> 6 calib… Semi-l… 8.57e-5 8.57e-5 8.57e-5 4.30e-2 1.50e-4 0.587 1.28e-1 +#> 7 ci AR 3 1.41e+2 3.76e+2 6.91e+2 1.60e+3 1.63e+3 8488. 2.28e+3 +#> 8 ci ARIMA 2.28e+1 4.62e+1 8.40e+1 1.04e+3 1.57e+3 6341. 1.72e+3 +#> 9 ci Semi-l… 2.51e+1 5.67e+1 1.26e+2 1.05e+3 1.62e+3 5913. 1.75e+3 +#> 10 crps AR 3 3.63e+1 9.19e+1 1.40e+2 1.58e+2 2.00e+2 429. 9.54e+1 +#> # … with 17 more rows ``` ## Docker diff --git a/_pkgdown.yml b/_pkgdown.yml index a8c0a72..7feaacb 100644 --- a/_pkgdown.yml +++ b/_pkgdown.yml @@ -43,6 +43,8 @@ reference: contents: - bsts_model - fable_model + - brms_model + - forecastHybrid_model - title: Predicting cases from reproduction numbers desc: Functions to for mapping reproduction number estimates to forecast cases contents: diff --git a/docs/404.html b/docs/404.html index 34ba4ac..0963de0 100644 --- a/docs/404.html +++ b/docs/404.html @@ -10,23 +10,27 @@ - + - + - + + + + + - - + + - + - - + + @@ -53,7 +57,7 @@ - +
+ +
@@ -127,7 +137,7 @@

Page not found (404)

-

Site built with pkgdown 1.4.1.

+

Site built with pkgdown 1.5.1.

diff --git a/docs/LICENSE-text.html b/docs/LICENSE-text.html index abced4e..f1d2678 100644 --- a/docs/LICENSE-text.html +++ b/docs/LICENSE-text.html @@ -10,23 +10,27 @@ - + - + - + + + + + - - + + - + - - + + @@ -53,7 +57,7 @@ - +
+ +
@@ -129,7 +139,7 @@

License

-

Site built with pkgdown 1.4.1.

+

Site built with pkgdown 1.5.1.

diff --git a/docs/LICENSE.html b/docs/LICENSE.html index bcb316e..7ab2c6c 100644 --- a/docs/LICENSE.html +++ b/docs/LICENSE.html @@ -10,23 +10,27 @@ - + - + - + + + + + - - + + - + - - + + @@ -53,7 +57,7 @@ - +
+ +
@@ -133,7 +143,7 @@

MIT License

-

Site built with pkgdown 1.4.1.

+

Site built with pkgdown 1.5.1.

diff --git a/docs/articles/index.html b/docs/articles/index.html index 0d4e376..a0bcdb8 100644 --- a/docs/articles/index.html +++ b/docs/articles/index.html @@ -10,23 +10,27 @@ - + - + - + + + + + - - + + - + - - + + @@ -53,7 +57,7 @@ - +
@@ -131,7 +136,7 @@

All vignettes

-

Site built with pkgdown 1.4.1.

+

Site built with pkgdown 1.5.1.

diff --git a/docs/articles/introduction.html b/docs/articles/introduction.html index 611f050..d7b4ceb 100644 --- a/docs/articles/introduction.html +++ b/docs/articles/introduction.html @@ -6,19 +6,19 @@ Getting started • EpiSoon - - - - + + + + + - - + - +
-

Getting started

+

Getting started

@@ -91,11 +91,15 @@

  • Load the package (bsts for models, ggplot2 for plotting, and cowplot for theming)
-
library(EpiSoon)
-library(bsts)
-library(fable)
-library(cowplot)
-library(dplyr)
+
library(EpiSoon)
+library(bsts)
+#> Warning: package 'BoomSpikeSlab' was built under R version 3.6.2
+#> Warning: package 'Boom' was built under R version 3.6.2
+#> Warning: package 'MASS' was built under R version 3.6.2
+library(fable)
+#> Warning: package 'fable' was built under R version 3.6.2
+library(cowplot)
+library(dplyr)

@@ -103,110 +107,116 @@

  • We use an example dataframe built into the package but this could be replaced with your own data.
- +
EpiSoon::example_obs_rts
+#>          rt       date
+#> 1  2.490547 2020-03-01
+#> 2  2.442588 2020-03-02
+#> 3  2.402473 2020-03-03
+#> 4  2.335572 2020-03-04
+#> 5  2.266551 2020-03-05
+#> 6  2.192293 2020-03-06
+#> 7  2.146429 2020-03-07
+#> 8  2.104371 2020-03-08
+#> 9  2.059281 2020-03-09
+#> 10 2.027134 2020-03-10
+#> 11 2.014678 2020-03-11
+#> 12 1.998946 2020-03-12
+#> 13 1.968350 2020-03-13
+#> 14 1.947376 2020-03-14
+#> 15 1.906984 2020-03-15
+#> 16 1.812842 2020-03-16
+#> 17 1.718532 2020-03-17
+#> 18 1.665646 2020-03-18
+#> 19 1.639927 2020-03-19
+#> 20 1.633795 2020-03-20
+#> 21 1.682025 2020-03-21
+#> 22 1.561653 2020-03-22
  • Fit a bsts model and produce a Rt forecast. Any appropriately wrapped model can be used (see bsts_model and fable_model for an examples).
- +
rt_forecast <- forecast_rt(EpiSoon::example_obs_rts[1:10, ],
+                       model = function(...){EpiSoon::bsts_model(model = function(ss, y){bsts::AddAutoAr(ss, y = y, lags = 10)}, ...)},
+                       horizon = 21, samples = 10)
+
+rt_forecast
+#> # A tibble: 210 x 4
+#>    sample date          rt horizon
+#>     <int> <date>     <dbl>   <int>
+#>  1      1 2020-03-11  1.99       1
+#>  2      2 2020-03-11  1.99       1
+#>  3      3 2020-03-11  1.97       1
+#>  4      4 2020-03-11  1.95       1
+#>  5      5 2020-03-11  1.97       1
+#>  6      6 2020-03-11  1.96       1
+#>  7      7 2020-03-11  1.98       1
+#>  8      8 2020-03-11  2.06       1
+#>  9      9 2020-03-11  2.05       1
+#> 10     10 2020-03-11  2.00       1
+#> # … with 200 more rows
  • Score the forecast
- +
rt_scores <- score_forecast(rt_forecast, EpiSoon::example_obs_rts)
+
+rt_scores
+#> # A tibble: 12 x 11
+#>    date       horizon   dss   crps   logs   bias sharpness calibration median
+#>    <date>       <int> <dbl>  <dbl>  <dbl>  <dbl>     <dbl>       <dbl>  <dbl>
+#>  1 2020-03-11       1 -6.29 0.0200 -1.46  -0.6      0.0290       0.878 0.0573
+#>  2 2020-03-12       2 -5.34 0.0179 -1.87  -0.200    0.0516       0.878 0.0468
+#>  3 2020-03-13       3 -4.69 0.0249 -1.63  -0.400    0.0601       0.878 0.0628
+#>  4 2020-03-14       4 -4.38 0.0341 -1.33  -0.8      0.0556       0.878 0.0881
+#>  5 2020-03-15       5 -4.19 0.0325 -1.08  -0.200    0.118        0.878 0.0452
+#>  6 2020-03-16       6 -3.89 0.0336 -0.997  0.200    0.125        0.878 0.0424
+#>  7 2020-03-17       7 -3.64 0.0306 -1.36   0.200    0.0718       0.878 0.0482
+#>  8 2020-03-18       8 -3.38 0.0507 -0.863  0.8      0.0995       0.878 0.109 
+#>  9 2020-03-19       9 -3.10 0.0587 -0.846  0.6      0.0938       0.878 0.151 
+#> 10 2020-03-20      10 -3.25 0.0451 -0.802  0.4      0.172        0.878 0.0858
+#> 11 2020-03-21      11 -3.18 0.0405 -0.965 -0.200    0.0966       0.878 0.0900
+#> 12 2020-03-22      12 -3.02 0.0428 -1.02   0.6      0.124        0.878 0.0752
+#> # … with 2 more variables: iqr <dbl>, ci <dbl>
  • Summarise the forecast scores
- +
summarise_scores(rt_scores)
+#> Warning: attributes are not identical across measure variables;
+#> they will be dropped
+#> # A tibble: 9 x 8
+#>   score        bottom   lower  median    mean   upper     top     sd
+#> * <chr>         <dbl>   <dbl>   <dbl>   <dbl>   <dbl>   <dbl>  <dbl>
+#> 1 bias        -0.745  -0.250   0       0.0333  0.45    0.745  0.510 
+#> 2 calibration  0.878   0.878   0.878   0.878   0.878   0.878  0     
+#> 3 ci           0.128   0.278   0.491   0.441   0.586   0.708  0.202 
+#> 4 crps         0.0185  0.0292  0.0338  0.0359  0.0434  0.0565 0.0122
+#> 5 dss         -6.03   -4.46   -3.77   -4.03   -3.24   -3.04   1.01  
+#> 6 iqr          0.0623  0.0942  0.150   0.146   0.171   0.264  0.0642
+#> 7 logs        -1.80   -1.38   -1.05   -1.19   -0.939  -0.814  0.342 
+#> 8 median       0.0432  0.0478  0.0690  0.0752  0.0886  0.140  0.0322
+#> 9 sharpness    0.0353  0.0590  0.0952  0.0915  0.119   0.159  0.0401
  • Summarise the forecast
- +
summarised_rt_forecast <- summarise_forecast(rt_forecast)
+
+summarised_rt_forecast
+#> # A tibble: 21 x 9
+#>    date       horizon median  mean     sd bottom lower upper   top
+#>    <date>       <int>  <dbl> <dbl>  <dbl>  <dbl> <dbl> <dbl> <dbl>
+#>  1 2020-03-11       1   1.99  1.99 0.0369   1.95  1.96  1.99  2.06
+#>  2 2020-03-12       2   1.98  1.97 0.0642   1.86  1.96  2.02  2.06
+#>  3 2020-03-13       3   1.94  1.92 0.0851   1.74  1.92  2.00  2.03
+#>  4 2020-03-14       4   1.90  1.88 0.0893   1.70  1.87  1.94  2.03
+#>  5 2020-03-15       5   1.88  1.85 0.116    1.63  1.87  2.00  2.00
+#>  6 2020-03-16       6   1.83  1.82 0.150    1.51  1.75  1.89  2.05
+#>  7 2020-03-17       7   1.74  1.76 0.165    1.46  1.70  1.80  2.04
+#>  8 2020-03-18       8   1.72  1.75 0.166    1.43  1.67  1.81  2.00
+#>  9 2020-03-19       9   1.72  1.75 0.183    1.43  1.64  1.74  2.08
+#> 10 2020-03-20      10   1.68  1.69 0.197    1.35  1.54  1.69  2.03
+#> # … with 11 more rows
  • Plot the forecast against observed data
-
plot_forecast(summarised_rt_forecast, EpiSoon::example_obs_rts)
+
plot_forecast(summarised_rt_forecast, EpiSoon::example_obs_rts)

@@ -215,106 +225,112 @@

  • Forecasting cases requires the observed cases on which the observed Rt estimates were based
- +
EpiSoon::example_obs_cases
+#> # A tibble: 63 x 2
+#>    cases date      
+#>    <dbl> <date>    
+#>  1     1 2020-01-20
+#>  2     0 2020-01-21
+#>  3     1 2020-01-22
+#>  4     0 2020-01-23
+#>  5     0 2020-01-24
+#>  6     0 2020-01-25
+#>  7     1 2020-01-26
+#>  8     0 2020-01-27
+#>  9     0 2020-01-28
+#> 10     0 2020-01-29
+#> # … with 53 more rows
  • It also requires an assumption to be made about the serial interval (defined using probability distribution).
- +
EpiSoon::example_serial_interval
+#>         1    2    3    4    5    6    7    8    9   10   11   12   14 
+#> 0.00 0.03 0.25 0.17 0.09 0.15 0.13 0.05 0.05 0.03 0.02 0.01 0.01 0.01
  • Forecast cases (using the case data on which the observed Rt estimates were based)
- +
case_forecast <- forecast_cases(EpiSoon::example_obs_cases, rt_forecast,
+                                serial_interval = EpiSoon::example_serial_interval)
+
+case_forecast
+#> # A tibble: 210 x 4
+#>    sample date       cases horizon
+#>  *  <dbl> <date>     <int>   <int>
+#>  1      1 2020-03-11   171       1
+#>  2      1 2020-03-12   210       2
+#>  3      1 2020-03-13   238       3
+#>  4      1 2020-03-14   281       4
+#>  5      1 2020-03-15   323       5
+#>  6      1 2020-03-16   352       6
+#>  7      1 2020-03-17   414       7
+#>  8      1 2020-03-18   445       8
+#>  9      1 2020-03-19   526       9
+#> 10      1 2020-03-20   539      10
+#> # … with 200 more rows
  • Score the cases forecast
- +
case_scores <- score_case_forecast(case_forecast, EpiSoon::example_obs_cases)
+
+case_scores
+#> # A tibble: 12 x 11
+#>    date       horizon   dss   crps  logs  bias sharpness calibration median
+#>    <date>       <int> <dbl>  <dbl> <dbl> <dbl>     <dbl>       <dbl>  <dbl>
+#>  1 2020-03-11       1  5.41   3.61  3.69   0.3      13.3   0.0000500      5
+#>  2 2020-03-12       2  5.56   5.97  4.09   0.4      13.3   0.0000500     21
+#>  3 2020-03-13       3  7.49  18.1   5.25   0.6      20.0   0.0000500     58
+#>  4 2020-03-14       4  6.58  10.3   4.57   0.5      22.2   0.0000500     36
+#>  5 2020-03-15       5  7.85  22.4   5.27   0.4      28.2   0.0000500     94
+#>  6 2020-03-16       6 10.4   58.9   6.15   0.8      56.3   0.0000500    176
+#>  7 2020-03-17       7 10.7   65.1   5.98   1        70.4   0.0000500    205
+#>  8 2020-03-18       8 10.7   75.6   6.26   0.8      86.0   0.0000500    209
+#>  9 2020-03-19       9 10.7   78.1   6.32   0.8     102.    0.0000500    250
+#> 10 2020-03-20      10 10.8   77.3   6.63   0.8      66.7   0.0000500    206
+#> 11 2020-03-21      11 10.8   41.8   5.75   0.4      96.4   0.0000500     75
+#> 12 2020-03-22      12 13.2  243.    7.47   0.8     139.    0.0000500    679
+#> # … with 2 more variables: iqr <dbl>, ci <dbl>
  • Summarise the cases scores
- +
summarise_scores(case_scores)
+#> Warning: attributes are not identical across measure variables;
+#> they will be dropped
+#> # A tibble: 9 x 8
+#>   score        bottom     lower     median      mean     upper       top      sd
+#> * <chr>         <dbl>     <dbl>      <dbl>     <dbl>     <dbl>     <dbl>   <dbl>
+#> 1 bias        3.28e-1   4.00e-1    7.00e-1   6.33e-1   8.00e-1   9.45e-1   0.227
+#> 2 calibrat…   5.00e-5   5.00e-5    5.00e-5   5.00e-5   5.00e-5   5.00e-5   0    
+#> 3 ci          4.41e+1   7.53e+1    2.34e+2   4.29e+2   5.44e+2   1.71e+3 573.   
+#> 4 crps        4.26e+0   1.61e+1    5.04e+1   5.84e+1   7.60e+1   1.98e+2  65.0  
+#> 5 dss         5.45e+0   7.26e+0    1.05e+1   9.18e+0   1.07e+1   1.25e+1   2.50 
+#> 6 iqr         1.56e+1   4.97e+1    1.92e+2   2.51e+2   3.20e+2   9.02e+2 300.   
+#> 7 logs        3.80e+0   5.08e+0    5.87e+0   5.62e+0   6.27e+0   7.24e+0   1.10 
+#> 8 median      9.40e+0   5.25e+1    1.35e+2   1.68e+2   2.07e+2   5.61e+2 182.   
+#> 9 sharpness   1.33e+1   2.17e+1    6.15e+1   5.96e+1   8.86e+1   1.29e+2  41.2
  • Summarise the cases forecast
- +
summarised_case_forecast <- summarise_case_forecast(case_forecast)
+
+summarised_case_forecast
+#> # A tibble: 21 x 9
+#>    date       horizon median  mean    sd bottom lower upper   top
+#>    <date>       <int>  <dbl> <dbl> <dbl>  <dbl> <dbl> <dbl> <dbl>
+#>  1 2020-03-11       1   170.  168.  15.7    141   165   181   186
+#>  2 2020-03-12       2   204.  200.  15.6    170   202   217   217
+#>  3 2020-03-13       3   237   234.  26.5    190   222   249   269
+#>  4 2020-03-14       4   269   264.  24.0    220   251   281   292
+#>  5 2020-03-15       5   320   300.  42.4    211   318   348   348
+#>  6 2020-03-16       6   354   349.  55.6    239   308   358   424
+#>  7 2020-03-17       7   398.  399.  73.0    298   334   414   553
+#>  8 2020-03-18       8   448.  461.  95.4    300   407   507   627
+#>  9 2020-03-19       9   524.  526. 122.     321   418   549   750
+#> 10 2020-03-20      10   557   586. 159.     335   508   598   915
+#> # … with 11 more rows
  • Plot the forecast against observed case data
-
plot_forecast(summarised_case_forecast, EpiSoon::example_obs_cases)
+
plot_forecast(summarised_case_forecast, EpiSoon::example_obs_cases)

@@ -323,49 +339,49 @@

  • To explore the quality of a models forecast it can help to iteratively forecast from each available data point. This is supported in EpiSoon using the following:
- +
it_rt_forecast <- iterative_rt_forecast(EpiSoon::example_obs_rts,
+                                        model = function(...){
+                                          EpiSoon::bsts_model(model = function(ss, y){bsts::AddAutoAr(ss, y = y, lags = 10)}, ...)
+                                          },
+                                        horizon = 7, samples = 10, min_points = 4)
+
+it_rt_forecast
+#> # A tibble: 1,260 x 5
+#>    forecast_date sample date          rt horizon
+#>  * <chr>          <int> <date>     <dbl>   <int>
+#>  1 2020-03-05         1 2020-03-06  2.31       1
+#>  2 2020-03-05         2 2020-03-06  2.16       1
+#>  3 2020-03-05         3 2020-03-06  2.22       1
+#>  4 2020-03-05         4 2020-03-06  2.22       1
+#>  5 2020-03-05         5 2020-03-06  2.33       1
+#>  6 2020-03-05         6 2020-03-06  2.25       1
+#>  7 2020-03-05         7 2020-03-06  2.16       1
+#>  8 2020-03-05         8 2020-03-06  2.23       1
+#>  9 2020-03-05         9 2020-03-06  2.21       1
+#> 10 2020-03-05        10 2020-03-06  2.17       1
+#> # … with 1,250 more rows
  • We can then iteratively forecast cases using the following:
- +
it_cases_forecast <- iterative_case_forecast(it_fit_samples = it_rt_forecast,
+                                             cases = EpiSoon::example_obs_cases,
+                                             serial_interval = EpiSoon::example_serial_interval)
+
+it_cases_forecast
+#> # A tibble: 1,260 x 5
+#>    forecast_date sample date       cases horizon
+#>  * <chr>          <dbl> <date>     <int>   <int>
+#>  1 2020-03-05         1 2020-03-06    79       1
+#>  2 2020-03-05         1 2020-03-07   101       2
+#>  3 2020-03-05         1 2020-03-08   108       3
+#>  4 2020-03-05         1 2020-03-09   142       4
+#>  5 2020-03-05         1 2020-03-10   173       5
+#>  6 2020-03-05         1 2020-03-11   207       6
+#>  7 2020-03-05         1 2020-03-12   257       7
+#>  8 2020-03-05         2 2020-03-06    81       1
+#>  9 2020-03-05         2 2020-03-07    78       2
+#> 10 2020-03-05         2 2020-03-08   104       3
+#> # … with 1,250 more rows
  • All functionality shown above is also supported for iterative forecasting.
@@ -374,93 +390,93 @@

Evaluate a model

In real world use we are likely to want to evaluate a model by iteratively forecasting Rts and cases, summarising these forecasts, scoring them and then returning them in a sensible format. These steps are all contained in the evaluate_model function.

-
model_eval <- evaluate_model(EpiSoon::example_obs_rts,
-                             EpiSoon::example_obs_cases,
-                             model = function(...){
-                                          EpiSoon::bsts_model(model = function(ss, y){bsts::AddAutoAr(ss, y = y, lags = 10)}, ...)
-                                          },
-                             horizon = 21, samples = 10,
-                             serial_interval = EpiSoon::example_serial_interval)
-
-model_eval
-#> $forecast_rts
-#> # A tibble: 399 x 10
-#>    forecast_date date       horizon median  mean    sd bottom lower upper   top
-#>    <chr>         <date>       <int>  <dbl> <dbl> <dbl>  <dbl> <dbl> <dbl> <dbl>
-#>  1 2020-03-04    2020-03-05       1   2.18 2.21  0.114 2.07    2.07  2.21  2.40
-#>  2 2020-03-04    2020-03-06       2   2.10 1.99  0.287 1.49    2.02  2.30  2.30
-#>  3 2020-03-04    2020-03-07       3   1.92 1.76  0.493 0.923   1.92  2.31  2.31
-#>  4 2020-03-04    2020-03-08       4   1.74 1.48  0.756 0.0800  1.69  2.33  2.33
-#>  5 2020-03-04    2020-03-09       5   1.63 1.30  0.809 0       1.60  2.24  2.24
-#>  6 2020-03-04    2020-03-10       6   1.55 1.19  0.821 0       1.51  2.23  2.23
-#>  7 2020-03-04    2020-03-11       7   1.41 1.05  0.869 0       1.35  2.19  2.19
-#>  8 2020-03-04    2020-03-12       8   1.29 0.980 0.880 0       1.16  2.15  2.15
-#>  9 2020-03-04    2020-03-13       9   1.39 0.985 0.875 0       1.28  2.10  2.10
-#> 10 2020-03-04    2020-03-14      10   1.31 0.924 0.823 0       1.30  2.05  2.05
-#> # … with 389 more rows
-#> 
-#> $rt_scores
-#> # A tibble: 171 x 12
-#>    forecast_date date       horizon     dss   crps   logs   bias sharpness
-#>    <chr>         <date>       <int>   <dbl>  <dbl>  <dbl>  <dbl>     <dbl>
-#>  1 2020-03-04    2020-03-05       1 -4.14   0.0475 -0.818 -0.400     0.127
-#>  2 2020-03-04    2020-03-06       2 -2.05   0.0886 -0.213 -0.400     0.232
-#>  3 2020-03-04    2020-03-07       3 -0.828  0.190   0.286 -0.6       0.452
-#>  4 2020-03-04    2020-03-08       4  0.0948 0.287   0.661 -0.6       0.539
-#>  5 2020-03-04    2020-03-09       5  0.441  0.373   0.879 -0.8       0.582
-#>  6 2020-03-04    2020-03-10       6  0.643  0.442   0.986 -0.8       0.754
-#>  7 2020-03-04    2020-03-11       7  0.994  0.554   1.12  -0.8       0.898
-#>  8 2020-03-04    2020-03-12       8  1.13   0.595   1.17  -0.8       1.00 
-#>  9 2020-03-04    2020-03-13       9  1.03   0.564   1.12  -0.8       0.855
-#> 10 2020-03-04    2020-03-14      10  1.22   0.625   1.18  -0.8       0.780
-#> # … with 161 more rows, and 4 more variables: calibration <dbl>, median <dbl>,
-#> #   iqr <dbl>, ci <dbl>
-#> 
-#> $forecast_cases
-#> # A tibble: 171 x 10
-#>    forecast_date date       horizon median  mean     sd bottom lower upper   top
-#>    <chr>         <date>       <int>  <dbl> <dbl>  <dbl>  <dbl> <dbl> <dbl> <dbl>
-#>  1 2020-03-04    2020-03-05       1   65.5  65.7   3.37     60    64    68    71
-#>  2 2020-03-04    2020-03-06       2   81    76.4  12.1      59    79    94    94
-#>  3 2020-03-04    2020-03-07       3   77    78.1  20.6      47    65    94   108
-#>  4 2020-03-04    2020-03-08       4   90.5  81.7  45.3       5    82   135   135
-#>  5 2020-03-04    2020-03-09       5  103    84.1  54.6       0    73   129   152
-#>  6 2020-03-04    2020-03-10       6  104.   95.6  71.9       0   101   198   198
-#>  7 2020-03-04    2020-03-11       7   99    87    77.0       0    82   195   195
-#>  8 2020-03-04    2020-03-12       8  101   102.   99.3       0     0   119   256
-#>  9 2020-03-04    2020-03-13       9  141   124.  123.        0     0   165   356
-#> 10 2020-03-04    2020-03-14      10  124   121.  124.        0     0   135   364
-#> # … with 161 more rows
-#> 
-#> $case_scores
-#> # A tibble: 171 x 13
-#>    sample forecast_date date       horizon   dss  crps  logs   bias sharpness
-#>    <chr>  <chr>         <date>       <int> <dbl> <dbl> <dbl>  <dbl>     <dbl>
-#>  1 1      2020-03-04    2020-03-05       1  3.04  1.69  2.53  0.6        3.71
-#>  2 1      2020-03-04    2020-03-06       2  4.97  4.82  3.86  0.200     13.3 
-#>  3 1      2020-03-04    2020-03-07       3  6.20  7.35  4.28 -0.400     21.5 
-#>  4 1      2020-03-04    2020-03-08       4  7.74 13.4   4.93 -0.200     51.1 
-#>  5 1      2020-03-04    2020-03-09       5  8.27 15.4   5.04 -0.300     41.5 
-#>  6 1      2020-03-04    2020-03-10       6  8.90 27.3   5.48 -0.200     92.7 
-#>  7 1      2020-03-04    2020-03-11       7  9.78 45.1   5.63 -0.8      113.  
-#>  8 1      2020-03-04    2020-03-12       8 10.0  52.9   5.93 -0.400    147.  
-#>  9 1      2020-03-04    2020-03-13       9 10.0  51.8   6.08 -0.6      162.  
-#> 10 1      2020-03-04    2020-03-14      10 10.8  88.4   6.39 -0.8      171.  
-#> # … with 161 more rows, and 4 more variables: calibration <dbl>, median <dbl>,
-#> #   iqr <dbl>, ci <dbl>
+
model_eval <- evaluate_model(EpiSoon::example_obs_rts,
+                             EpiSoon::example_obs_cases,
+                             model = function(...){
+                                          EpiSoon::bsts_model(model = function(ss, y){bsts::AddAutoAr(ss, y = y, lags = 10)}, ...)
+                                          },
+                             horizon = 21, samples = 10,
+                             serial_interval = EpiSoon::example_serial_interval)
+
+model_eval
+#> $forecast_rts
+#> # A tibble: 399 x 10
+#>    forecast_date date       horizon median  mean    sd bottom lower upper   top
+#>  * <chr>         <date>       <int>  <dbl> <dbl> <dbl>  <dbl> <dbl> <dbl> <dbl>
+#>  1 2020-03-04    2020-03-05       1   2.29  2.29 0.164  2.01   2.25  2.35  2.66
+#>  2 2020-03-04    2020-03-06       2   2.20  2.24 0.404  1.82   2.12  2.34  3.27
+#>  3 2020-03-04    2020-03-07       3   2.10  2.10 0.413  1.49   2.03  2.29  2.95
+#>  4 2020-03-04    2020-03-08       4   1.98  1.95 0.387  1.13   1.89  2.11  2.49
+#>  5 2020-03-04    2020-03-09       5   1.81  1.80 0.456  0.868  1.58  2.05  2.41
+#>  6 2020-03-04    2020-03-10       6   1.81  1.67 0.504  0.678  1.81  2.23  2.23
+#>  7 2020-03-04    2020-03-11       7   1.83  1.61 0.614  0.520  1.82  2.56  2.56
+#>  8 2020-03-04    2020-03-12       8   1.78  1.52 0.709  0.392  1.19  1.95  2.75
+#>  9 2020-03-04    2020-03-13       9   1.59  1.32 0.819  0.122  1.04  1.98  2.65
+#> 10 2020-03-04    2020-03-14      10   1.52  1.27 0.864  0      1.00  2.06  2.73
+#> # … with 389 more rows
+#> 
+#> $rt_scores
+#> # A tibble: 171 x 12
+#>    forecast_date date       horizon    dss   crps   logs   bias sharpness
+#>    <chr>         <date>       <int>  <dbl>  <dbl>  <dbl>  <dbl>     <dbl>
+#>  1 2020-03-04    2020-03-05       1 -3.70  0.0276 -1.33   0.200    0.0839
+#>  2 2020-03-04    2020-03-06       2 -1.90  0.0638 -0.422  0        0.169 
+#>  3 2020-03-04    2020-03-07       3 -1.86  0.0695 -0.563 -0.400    0.198 
+#>  4 2020-03-04    2020-03-08       4 -1.82  0.0862 -0.310 -0.400    0.169 
+#>  5 2020-03-04    2020-03-09       5 -1.32  0.146   0.239 -0.6      0.352 
+#>  6 2020-03-04    2020-03-10       6 -0.921 0.170   0.370 -0.400    0.483 
+#>  7 2020-03-04    2020-03-11       7 -0.595 0.214   0.631 -0.6      0.688 
+#>  8 2020-03-04    2020-03-12       8 -0.281 0.264   0.801 -0.8      0.563 
+#>  9 2020-03-04    2020-03-13       9  0.185 0.351   0.981 -0.6      0.701 
+#> 10 2020-03-04    2020-03-14      10  0.276 0.394   1.05  -0.6      0.785 
+#> # … with 161 more rows, and 4 more variables: calibration <dbl>, median <dbl>,
+#> #   iqr <dbl>, ci <dbl>
+#> 
+#> $forecast_cases
+#> # A tibble: 171 x 10
+#>    forecast_date date       horizon median  mean    sd bottom lower upper   top
+#>  * <chr>         <date>       <int>  <dbl> <dbl> <dbl>  <dbl> <dbl> <dbl> <dbl>
+#>  1 2020-03-04    2020-03-05       1   67.5  69.3  10.4     58    58    69    93
+#>  2 2020-03-04    2020-03-06       2   81.5  83.7  13.5     71    76    86   118
+#>  3 2020-03-04    2020-03-07       3  100.  104.   27.1     70    82   105   167
+#>  4 2020-03-04    2020-03-08       4  109   104    27.1     47    91   116   151
+#>  5 2020-03-04    2020-03-09       5  133   130.   44.3     55   114   167   206
+#>  6 2020-03-04    2020-03-10       6  143   139.   54.1     37   104   173   213
+#>  7 2020-03-04    2020-03-11       7  164.  153.   73.5     38    86   177   260
+#>  8 2020-03-04    2020-03-12       8  185   169.  105.      11    83   220   370
+#>  9 2020-03-04    2020-03-13       9  182.  177.  139.      14    14   194   444
+#> 10 2020-03-04    2020-03-14      10  214.  209   196.       0    91   303   662
+#> # … with 161 more rows
+#> 
+#> $case_scores
+#> # A tibble: 171 x 13
+#>    sample forecast_date date       horizon   dss  crps  logs   bias sharpness
+#>  * <chr>  <chr>         <date>       <int> <dbl> <dbl> <dbl>  <dbl>     <dbl>
+#>  1 1      2020-03-04    2020-03-05       1  4.99  3.05  3.29  0.200      9.64
+#>  2 1      2020-03-04    2020-03-06       2  5.80  5.13  3.41  0.6        7.41
+#>  3 1      2020-03-04    2020-03-07       3  6.86  8.09  4.14  0.4       17.0 
+#>  4 1      2020-03-04    2020-03-08       4  6.50  6.    4.09  0.200     19.3 
+#>  5 1      2020-03-04    2020-03-09       5  7.59 12.8   4.92  0.200     45.2 
+#>  6 1      2020-03-04    2020-03-10       6  7.88 13.2   5.04  0         51.1 
+#>  7 1      2020-03-04    2020-03-11       7  8.53 17.7   5.42 -0.200     97.9 
+#>  8 1      2020-03-04    2020-03-12       8  9.26 23.9   5.68  0        102.  
+#>  9 1      2020-03-04    2020-03-13       9  9.82 39.3   6.10 -0.200    153.  
+#> 10 1      2020-03-04    2020-03-14      10 10.5  46.1   6.22 -0.200    157.  
+#> # … with 161 more rows, and 4 more variables: calibration <dbl>, median <dbl>,
+#> #   iqr <dbl>, ci <dbl>
  • All functionality outlined above can be applied to this output but a special plotting function (plot_forecast_evaluation) is also provided. First evaluate the Rt forecast against observed values.
- +
plot_forecast_evaluation(model_eval$forecast_rts,
+                         EpiSoon::example_obs_rts,
+                         horizon_to_plot = 7)

  • Then evaluate forecast cases against observed values.
- +
plot_forecast_evaluation(model_eval$forecast_cases,
+                         EpiSoon::example_obs_cases,
+                         horizon_to_plot = 7)

@@ -475,22 +491,11 @@

- @@ -501,7 +506,7 @@

-

Site built with pkgdown 1.4.1.

+

Site built with pkgdown 1.5.1.

diff --git a/docs/articles/introduction_files/figure-html/unnamed-chunk-15-1.png b/docs/articles/introduction_files/figure-html/unnamed-chunk-15-1.png index 0ea5977..be5dac0 100644 Binary files a/docs/articles/introduction_files/figure-html/unnamed-chunk-15-1.png and b/docs/articles/introduction_files/figure-html/unnamed-chunk-15-1.png differ diff --git a/docs/articles/introduction_files/figure-html/unnamed-chunk-19-1.png b/docs/articles/introduction_files/figure-html/unnamed-chunk-19-1.png index 2712736..303c73a 100644 Binary files a/docs/articles/introduction_files/figure-html/unnamed-chunk-19-1.png and b/docs/articles/introduction_files/figure-html/unnamed-chunk-19-1.png differ diff --git a/docs/articles/introduction_files/figure-html/unnamed-chunk-20-1.png b/docs/articles/introduction_files/figure-html/unnamed-chunk-20-1.png index 3914863..f3c0fa7 100644 Binary files a/docs/articles/introduction_files/figure-html/unnamed-chunk-20-1.png and b/docs/articles/introduction_files/figure-html/unnamed-chunk-20-1.png differ diff --git a/docs/articles/introduction_files/figure-html/unnamed-chunk-8-1.png b/docs/articles/introduction_files/figure-html/unnamed-chunk-8-1.png index b63b044..147ab7d 100644 Binary files a/docs/articles/introduction_files/figure-html/unnamed-chunk-8-1.png and b/docs/articles/introduction_files/figure-html/unnamed-chunk-8-1.png differ diff --git a/docs/authors.html b/docs/authors.html index ccb4fd3..92cf86b 100644 --- a/docs/authors.html +++ b/docs/authors.html @@ -10,23 +10,27 @@ - + - + - + + + + + - - + + - + - - + + @@ -53,7 +57,7 @@ - +
-

Site built with pkgdown 1.4.1.

+

Site built with pkgdown 1.5.1.

diff --git a/docs/bootstrap-toc.css b/docs/bootstrap-toc.css new file mode 100644 index 0000000..5a85941 --- /dev/null +++ b/docs/bootstrap-toc.css @@ -0,0 +1,60 @@ +/*! + * Bootstrap Table of Contents v0.4.1 (http://afeld.github.io/bootstrap-toc/) + * Copyright 2015 Aidan Feldman + * Licensed under MIT (https://github.com/afeld/bootstrap-toc/blob/gh-pages/LICENSE.md) */ + +/* modified from https://github.com/twbs/bootstrap/blob/94b4076dd2efba9af71f0b18d4ee4b163aa9e0dd/docs/assets/css/src/docs.css#L548-L601 */ + +/* All levels of nav */ +nav[data-toggle='toc'] .nav > li > a { + display: block; + padding: 4px 20px; + font-size: 13px; + font-weight: 500; + color: #767676; +} +nav[data-toggle='toc'] .nav > li > a:hover, +nav[data-toggle='toc'] .nav > li > a:focus { + padding-left: 19px; + color: #563d7c; + text-decoration: none; + background-color: transparent; + border-left: 1px solid #563d7c; +} +nav[data-toggle='toc'] .nav > .active > a, +nav[data-toggle='toc'] .nav > .active:hover > a, +nav[data-toggle='toc'] .nav > .active:focus > a { + padding-left: 18px; + font-weight: bold; + color: #563d7c; + background-color: transparent; + border-left: 2px solid #563d7c; +} + +/* Nav: second level (shown on .active) */ +nav[data-toggle='toc'] .nav .nav { + display: none; /* Hide by default, but at >768px, show it */ + padding-bottom: 10px; +} +nav[data-toggle='toc'] .nav .nav > li > a { + padding-top: 1px; + padding-bottom: 1px; + padding-left: 30px; + font-size: 12px; + font-weight: normal; +} +nav[data-toggle='toc'] .nav .nav > li > a:hover, +nav[data-toggle='toc'] .nav .nav > li > a:focus { + padding-left: 29px; +} +nav[data-toggle='toc'] .nav .nav > .active > a, +nav[data-toggle='toc'] .nav .nav > .active:hover > a, +nav[data-toggle='toc'] .nav .nav > .active:focus > a { + padding-left: 28px; + font-weight: 500; +} + +/* from https://github.com/twbs/bootstrap/blob/e38f066d8c203c3e032da0ff23cd2d6098ee2dd6/docs/assets/css/src/docs.css#L631-L634 */ +nav[data-toggle='toc'] .nav > .active > ul { + display: block; +} diff --git a/docs/bootstrap-toc.js b/docs/bootstrap-toc.js new file mode 100644 index 0000000..1cdd573 --- /dev/null +++ b/docs/bootstrap-toc.js @@ -0,0 +1,159 @@ +/*! + * Bootstrap Table of Contents v0.4.1 (http://afeld.github.io/bootstrap-toc/) + * Copyright 2015 Aidan Feldman + * Licensed under MIT (https://github.com/afeld/bootstrap-toc/blob/gh-pages/LICENSE.md) */ +(function() { + 'use strict'; + + window.Toc = { + helpers: { + // return all matching elements in the set, or their descendants + findOrFilter: function($el, selector) { + // http://danielnouri.org/notes/2011/03/14/a-jquery-find-that-also-finds-the-root-element/ + // http://stackoverflow.com/a/12731439/358804 + var $descendants = $el.find(selector); + return $el.filter(selector).add($descendants).filter(':not([data-toc-skip])'); + }, + + generateUniqueIdBase: function(el) { + var text = $(el).text(); + var anchor = text.trim().toLowerCase().replace(/[^A-Za-z0-9]+/g, '-'); + return anchor || el.tagName.toLowerCase(); + }, + + generateUniqueId: function(el) { + var anchorBase = this.generateUniqueIdBase(el); + for (var i = 0; ; i++) { + var anchor = anchorBase; + if (i > 0) { + // add suffix + anchor += '-' + i; + } + // check if ID already exists + if (!document.getElementById(anchor)) { + return anchor; + } + } + }, + + generateAnchor: function(el) { + if (el.id) { + return el.id; + } else { + var anchor = this.generateUniqueId(el); + el.id = anchor; + return anchor; + } + }, + + createNavList: function() { + return $('
    '); + }, + + createChildNavList: function($parent) { + var $childList = this.createNavList(); + $parent.append($childList); + return $childList; + }, + + generateNavEl: function(anchor, text) { + var $a = $(''); + $a.attr('href', '#' + anchor); + $a.text(text); + var $li = $('
  • '); + $li.append($a); + return $li; + }, + + generateNavItem: function(headingEl) { + var anchor = this.generateAnchor(headingEl); + var $heading = $(headingEl); + var text = $heading.data('toc-text') || $heading.text(); + return this.generateNavEl(anchor, text); + }, + + // Find the first heading level (`

    `, then `

    `, etc.) that has more than one element. Defaults to 1 (for `

    `). + getTopLevel: function($scope) { + for (var i = 1; i <= 6; i++) { + var $headings = this.findOrFilter($scope, 'h' + i); + if ($headings.length > 1) { + return i; + } + } + + return 1; + }, + + // returns the elements for the top level, and the next below it + getHeadings: function($scope, topLevel) { + var topSelector = 'h' + topLevel; + + var secondaryLevel = topLevel + 1; + var secondarySelector = 'h' + secondaryLevel; + + return this.findOrFilter($scope, topSelector + ',' + secondarySelector); + }, + + getNavLevel: function(el) { + return parseInt(el.tagName.charAt(1), 10); + }, + + populateNav: function($topContext, topLevel, $headings) { + var $context = $topContext; + var $prevNav; + + var helpers = this; + $headings.each(function(i, el) { + var $newNav = helpers.generateNavItem(el); + var navLevel = helpers.getNavLevel(el); + + // determine the proper $context + if (navLevel === topLevel) { + // use top level + $context = $topContext; + } else if ($prevNav && $context === $topContext) { + // create a new level of the tree and switch to it + $context = helpers.createChildNavList($prevNav); + } // else use the current $context + + $context.append($newNav); + + $prevNav = $newNav; + }); + }, + + parseOps: function(arg) { + var opts; + if (arg.jquery) { + opts = { + $nav: arg + }; + } else { + opts = arg; + } + opts.$scope = opts.$scope || $(document.body); + return opts; + } + }, + + // accepts a jQuery object, or an options object + init: function(opts) { + opts = this.helpers.parseOps(opts); + + // ensure that the data attribute is in place for styling + opts.$nav.attr('data-toggle', 'toc'); + + var $topContext = this.helpers.createChildNavList(opts.$nav); + var topLevel = this.helpers.getTopLevel(opts.$scope); + var $headings = this.helpers.getHeadings(opts.$scope, topLevel); + this.helpers.populateNav($topContext, topLevel, $headings); + } + }; + + $(function() { + $('nav[data-toggle="toc"]').each(function(i, el) { + var $nav = $(el); + Toc.init($nav); + }); + }); +})(); diff --git a/docs/index.html b/docs/index.html index c0b9daf..b3c2419 100644 --- a/docs/index.html +++ b/docs/index.html @@ -6,20 +6,20 @@ Forecast Cases Using Reproduction Numbers • EpiSoon - - - - + + + + + - - +
    -
    -

    Site built with pkgdown 1.4.1.

    +

    Site built with pkgdown 1.5.1.

    diff --git a/docs/reference/summarise_forecast.html b/docs/reference/summarise_forecast.html index b0a0443..a55c40e 100644 --- a/docs/reference/summarise_forecast.html +++ b/docs/reference/summarise_forecast.html @@ -10,23 +10,27 @@ - + - + - + + + + + - - + + - + - - + + @@ -37,7 +41,6 @@ - @@ -55,7 +58,7 @@ - +
    - @@ -174,7 +173,7 @@

    Contents

    -

    Site built with pkgdown 1.4.1.

    +

    Site built with pkgdown 1.5.1.

    diff --git a/docs/reference/summarise_scores.html b/docs/reference/summarise_scores.html index e57b369..754b9ed 100644 --- a/docs/reference/summarise_scores.html +++ b/docs/reference/summarise_scores.html @@ -10,23 +10,27 @@ - + - + - + + + + + - - + + - + - - + + @@ -37,7 +41,6 @@ - @@ -55,7 +58,7 @@ - +
    - @@ -236,7 +246,7 @@

    Contents

    -

    Site built with pkgdown 1.4.1.

    +

    Site built with pkgdown 1.5.1.

    diff --git a/docs/sitemap.xml b/docs/sitemap.xml index 8a4efed..fe93046 100644 --- a/docs/sitemap.xml +++ b/docs/sitemap.xml @@ -3,6 +3,9 @@ epiforecasts.io/EpiSoon//index.html + + epiforecasts.io/EpiSoon//reference/brms_model.html + epiforecasts.io/EpiSoon//reference/bsts_model.html diff --git a/man/brms_model.Rd b/man/brms_model.Rd new file mode 100644 index 0000000..d00f8ae --- /dev/null +++ b/man/brms_model.Rd @@ -0,0 +1,64 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/brms_model.R +\name{brms_model} +\alias{brms_model} +\title{brms model wrapper} +\usage{ +brms_model( + y = NULL, + samples = NULL, + horizon = NULL, + model = NULL, + n_cores = 1, + n_chains = 4, + n_iter = 2000, + ... +) +} +\arguments{ +\item{y}{Numeric vector of time points to forecast} + +\item{samples}{Numeric, number of samples to take.} + +\item{horizon}{Numeric, the time horizon over which to predict.} + +\item{model}{A \code{brms} model wrapped in the \code{\link[brms:bf]{brms::bf()}} function} + +\item{n_cores}{Numeric, the number of cores to use, default of 1} + +\item{n_chains}{Numeric, the number of chains to use, default of 4} + +\item{n_iter}{Numeric, the number of iterations in the sampler to use, +default of 4000} + +\item{...}{additional arguments passed to \code{brms} (e.g. priors or family)} +} +\value{ +A dataframe of predictions (with columns representing the time horizon and rows representing samples). +} +\description{ +Allows users to specify a model using the \code{\link[brms:bf]{brms::bf()}} wrapper from \code{brms} +Note that \code{brms} and \code{tidybayes} must both be installed for this +model wrapper to be functional. +} +\examples{ +\dontrun{ + +## Used on its own +## Note: More iterations and chains should be used +library(brms) +brms_model(y = EpiSoon::example_obs_rts[1:10, ]$rt, + model = brms::bf(y ~ gp(time)), + samples = 10, horizon = 7, n_iter = 40, n_chains = 1, refresh =0) + +## Used for forecasting +## Note that the timeout parameter has been increased to allow +## for the time for the code to be compiled +## Note: More iterations and chains should be used + +forecast_rt(EpiSoon::example_obs_rts[1:10, ], + model = function(...){ + brms_model(model = brms::bf(y ~ gp(time)), n_iter = 40, n_chains = 1, ...)}, + horizon = 7, samples = 10, timeout = 300) +} +} diff --git a/man/bsts_model.Rd b/man/bsts_model.Rd index ac9afff..fb00b78 100644 --- a/man/bsts_model.Rd +++ b/man/bsts_model.Rd @@ -22,6 +22,9 @@ A dataframe of predictions (with columns representing the time horizon and rows BSTS model wrapper } \examples{ +\dontrun{ + +library(bsts) ## Used on its own bsts_model(y = EpiSoon::example_obs_rts[1:10, ]$rt, @@ -36,3 +39,4 @@ bsts_model(y = EpiSoon::example_obs_rts[1:10, ]$rt, bsts::AddAr(ss, y = y, lags = 3)}, ...)}, horizon = 7, samples = 10) } +} diff --git a/man/example_obs_cases.Rd b/man/example_obs_cases.Rd index 21aa9e5..4438045 100644 --- a/man/example_obs_cases.Rd +++ b/man/example_obs_cases.Rd @@ -4,7 +4,9 @@ \name{example_obs_cases} \alias{example_obs_cases} \title{Example Observed Cases} -\format{A data frame containing cases reported on each date.} +\format{ +A data frame containing cases reported on each date. +} \usage{ example_obs_cases } diff --git a/man/example_obs_rts.Rd b/man/example_obs_rts.Rd index 5d508be..cad3deb 100644 --- a/man/example_obs_rts.Rd +++ b/man/example_obs_rts.Rd @@ -4,7 +4,9 @@ \name{example_obs_rts} \alias{example_obs_rts} \title{Example Observed Rts} -\format{A data frame containing Rts estimated for each date.} +\format{ +A data frame containing Rts estimated for each date. +} \usage{ example_obs_rts } diff --git a/man/example_serial_interval.Rd b/man/example_serial_interval.Rd index 1fee3c2..9117483 100644 --- a/man/example_serial_interval.Rd +++ b/man/example_serial_interval.Rd @@ -4,7 +4,9 @@ \name{example_serial_interval} \alias{example_serial_interval} \title{Example Serial Interval} -\format{A vector giviing the probability for each day} +\format{ +A vector giviing the probability for each day +} \usage{ example_serial_interval } diff --git a/man/fable_model.Rd b/man/fable_model.Rd index f82e50a..002f807 100644 --- a/man/fable_model.Rd +++ b/man/fable_model.Rd @@ -21,11 +21,13 @@ A dataframe of predictions (with columns representing the time horizon and rows representing samples). } \description{ -Provides an interface for models from the \code{fable} package. Not the \code{ARIMA} model -requires the \code{feast} package. If \code{future} is being used \code{fable} will require \code{future.apply} in +Provides an interface for models from the \code{fable} package. +Note the \code{feasts::ARIMA} model requires the \code{feast} package. If \code{future} +is being used \code{fable} will require \code{future.apply} in order to not silently fail. } \examples{ +\dontrun{ ## Used on its own fable_model(y = EpiSoon::example_obs_rts[1:10, ]$rt, model = fable::ARIMA(y ~ time), @@ -37,3 +39,4 @@ forecast_rt(EpiSoon::example_obs_rts[1:10, ], fable_model(model = fable::ARIMA(y ~ time), ...)}, horizon = 7, samples = 10) } +} diff --git a/man/figures/unnamed-chunk-7-1.png b/man/figures/unnamed-chunk-7-1.png index a6ba98c..4727b68 100644 Binary files a/man/figures/unnamed-chunk-7-1.png and b/man/figures/unnamed-chunk-7-1.png differ diff --git a/man/figures/unnamed-chunk-8-1.png b/man/figures/unnamed-chunk-8-1.png index ef6ef88..0fcd2e3 100644 Binary files a/man/figures/unnamed-chunk-8-1.png and b/man/figures/unnamed-chunk-8-1.png differ diff --git a/man/forecastHybrid_model.Rd b/man/forecastHybrid_model.Rd new file mode 100644 index 0000000..aea16a6 --- /dev/null +++ b/man/forecastHybrid_model.Rd @@ -0,0 +1,69 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/forecastHybrid_model.R +\name{forecastHybrid_model} +\alias{forecastHybrid_model} +\title{forecastHybrid model wrapper} +\usage{ +forecastHybrid_model( + y = NULL, + samples = NULL, + horizon = NULL, + model_params = NULL, + forecast_params = NULL +) +} +\arguments{ +\item{y}{Numeric vector of time points to forecast} + +\item{samples}{Numeric, number of samples to take.} + +\item{horizon}{Numeric, the time horizon over which to predict.} + +\item{model_params}{List of parameters to pass to \code{forecastHybrid::hybridModel}.} + +\item{forecast_params}{List of parameters to pass to \code{forecastHybrid:::forecast.hybridModel}.} +} +\value{ +A dataframe of predictions (with columns representing the time horizon and rows representing samples). +} +\description{ +Allows users to forecast using ensembles from the \code{forecastHybrid} package. Note that +whilst weighted ensembles can be created this is not advised when samples > 1 as currently +samples are derived assuming a normal distribution using the upper and lower confidence intervals of the ensemble. +These confidence intervals are themselves either based on the unweighted mean of the ensembled +models or the maximum/minimum from the candiate models. Note that \code{forecastHybrid} must be installed for this +model wrapper to be functional. +} +\examples{ +\dontrun{ + +library(forecastHybrid) + +## Used on its own +forecastHybrid_model(y = EpiSoon::example_obs_rts$rt, + samples = 10, horizon = 7) + + +## Used with non-default arguments +## Note that with the current sampling from maximal confidence intervals model +## Weighting using cross-validation will only have an impact when 1 sample is used. +forecastHybrid_model(y = EpiSoon::example_obs_rts$rt, + samples = 1, horizon = 7, + model_params = list(cvHorizon = 7, windowSize = 7, + rolling = TRUE, models = "zeta")) + + +## Used for forecasting + forecast_rt(EpiSoon::example_obs_rts, + model = EpiSoon::forecastHybrid_model, + horizon = 7, samples = 1) + +## Used for forcasting with non-default arguments +forecast_rt(EpiSoon::example_obs_rts, + model = function(...){EpiSoon::forecastHybrid_model( + model_params = list(models = "zte"), + forecast_params = list(PI.combination = "mean"), ...) + }, + horizon = 7, samples = 10) +} +} diff --git a/man/forecast_rt.Rd b/man/forecast_rt.Rd index 88627de..6ea3cc9 100644 --- a/man/forecast_rt.Rd +++ b/man/forecast_rt.Rd @@ -10,7 +10,7 @@ forecast_rt( horizon = 7, samples = 1000, bound_rt = TRUE, - timeout = 30 + timeout = 100 ) } \arguments{ diff --git a/tests/testthat.R b/tests/testthat.R new file mode 100644 index 0000000..efda7e7 --- /dev/null +++ b/tests/testthat.R @@ -0,0 +1,4 @@ +library("testthat") +library("EpiSoon") + +test_check("EpiSoon") diff --git a/tests/testthat/test_compare_models.R b/tests/testthat/test_compare_models.R new file mode 100644 index 0000000..cc8fec2 --- /dev/null +++ b/tests/testthat/test_compare_models.R @@ -0,0 +1,54 @@ +##------------------------------------------------------------ +context("Test of the function 'compare_models'.... testing outputs") + +models <- list("AR 3" = function(...) { + EpiSoon::bsts_model(model = function(ss, y){ + bsts::AddAr(ss, y = y, lags = 2) + }, ...)}, + "Semi-local linear trend" = function(...) { + EpiSoon::bsts_model(model = function(ss, y){ + bsts::AddSemilocalLinearTrend(ss, y = y) + }, ...)} +) + +set.seed(1234) + +out <- compare_models(EpiSoon::example_obs_rts, + EpiSoon::example_obs_cases, + models, + horizon = 7, samples = 10, + serial_interval = EpiSoon::example_serial_interval) + +test_that("Outputs have proper lengths and names", { + expect_length(out, 4) + + expect_named(out, c("forecast_rts", "rt_scores", "forecast_cases", "case_scores")) + + # 'return_raw' version + out_raw <- compare_models(EpiSoon::example_obs_rts, + EpiSoon::example_obs_cases, + models, + horizon = 7, samples = 10, + serial_interval = EpiSoon::example_serial_interval, + return_raw = TRUE) + + expect_length(out_raw, 6) + + expect_named(out_raw, c("forecast_rts", "rt_scores", "forecast_cases", + "case_scores", "raw_rt_forecast", "raw_case_forecast")) + +}) + +test_that("Outputs return results for all models", { + expect_identical(names(models), unique(out$forecast_rts$model)) + expect_equal(sum(is.na(out$forecast_rts)), 0) + + expect_identical(names(models), unique(out$rt_scores$model)) + expect_equal(sum(is.na(out$rt_scores)), 0) + + expect_identical(names(models), unique(out$forecast_cases$model)) + expect_equal(sum(is.na(out$forecast_cases)), 0) + + expect_identical(names(models), unique(out$case_scores$model)) + expect_equal(sum(is.na(out$case_scores)), 0) +}) diff --git a/tests/testthat/test_compare_timeseries.R b/tests/testthat/test_compare_timeseries.R new file mode 100644 index 0000000..d2fb192 --- /dev/null +++ b/tests/testthat/test_compare_timeseries.R @@ -0,0 +1,58 @@ +##------------------------------------------------------------ +context("Test of the function 'compare_timeseries'.... testing outputs") + +obs_rts <- EpiSoon::example_obs_rts %>% + dplyr::mutate(timeseries = "Region 1") %>% + dplyr::bind_rows(EpiSoon::example_obs_rts %>% + dplyr::mutate(timeseries = "Region 2")) + +obs_cases <- EpiSoon::example_obs_cases %>% + dplyr::mutate(timeseries = "Region 1") %>% + dplyr::bind_rows(EpiSoon::example_obs_cases %>% + dplyr::mutate(timeseries = "Region 2")) + +models <- list("AR 3" = function(...){ + EpiSoon::bsts_model(model = function(ss, y){ + bsts::AddAr(ss, y = y, lags = 3) + }, ...)}, + "Semi-local linear trend" = function(...) { + EpiSoon::bsts_model(model = function(ss, y){ + bsts::AddSemilocalLinearTrend(ss, y = y) + }, ...)} +) + +out <- compare_timeseries(obs_rts, obs_cases, models, + horizon = 7, samples = 10, + serial_interval = EpiSoon::example_serial_interval) + +test_that("Outputs have proper lengths and names", { + expect_length(out, 4) + + expect_named(out, c("forecast_rts", "rt_scores", "forecast_cases", "case_scores")) + + # 'return_raw' version + out_raw <- compare_timeseries(obs_rts, obs_cases, models, + horizon = 7, samples = 10, + serial_interval = EpiSoon::example_serial_interval, + return_raw = TRUE) + + expect_length(out_raw, 6) + + expect_named(out_raw, c("forecast_rts", "rt_scores", "forecast_cases", + "case_scores", "raw_rt_forecast", "raw_case_forecast")) + +}) + +test_that("Outputs return results for all models", { + expect_identical(names(models), unique(out$forecast_rts$model)) + expect_equal(sum(is.na(out$forecast_rts)), 0) + + expect_identical(names(models), unique(out$rt_scores$model)) + expect_equal(sum(is.na(out$rt_scores)), 0) + + expect_identical(names(models), unique(out$forecast_cases$model)) + # expect_equal(sum(is.na(out$forecast_cases)), 0) + + expect_identical(names(models), unique(out$case_scores$model)) + # expect_equal(sum(is.na(out$case_scores)), 0) +}) diff --git a/tests/testthat/test_draw_from_si_prob.R b/tests/testthat/test_draw_from_si_prob.R new file mode 100644 index 0000000..f6afe05 --- /dev/null +++ b/tests/testthat/test_draw_from_si_prob.R @@ -0,0 +1,17 @@ +# draw_from_si_prob +# This function returns result for requested days. It simply calls for +# serial_interval[days_ago] but checking first the 'days_ago' are in the +# interval. Otherwise, it returns 0. +# Things to test: +## In tests +# - that the function return the proper result +# - that the result is of the proper class + +exampleResult <- draw_from_si_prob(c(1, 2, 4, 10), + EpiSoon::example_serial_interval) +# expectedResult <- EpiSoon::example_serial_interval[c(1, 2, 4, 10)] +test_that("The expected draw is obtained", { + # expect_equal(exampleResult, expectedResult) + expect_true(is.numeric(exampleResult)) + expect_length(exampleResult, 4) +}) diff --git a/tests/testthat/test_evaluate_model.R b/tests/testthat/test_evaluate_model.R new file mode 100644 index 0000000..06ae80d --- /dev/null +++ b/tests/testthat/test_evaluate_model.R @@ -0,0 +1,60 @@ +##------------------------------------------------------------ +context("Test of the function 'evaluate_model'.... testing inputs") + +sampled_obs <- EpiSoon::example_obs_rts %>% + dplyr::mutate(sample = 1) %>% + dplyr::bind_rows(EpiSoon::example_obs_rts %>% + dplyr::mutate(sample = 2)) + +sampled_cases <- EpiSoon::example_obs_cases %>% + dplyr::mutate(sample = 1) %>% + dplyr::bind_rows(EpiSoon::example_obs_cases %>% + dplyr::mutate(sample = 2)) + +a_model <- function(...) { + EpiSoon::bsts_model(model = function(ss, y) { + bsts::AddSemilocalLinearTrend(ss, y = y) + }, ...) +} + +test_that("Inputs with unequal lengths return error", { + sampled_cases_sub <- sampled_cases %>% filter(sample != 2) + + expect_error( + evaluate_model(sampled_obs, + sampled_cases_sub, + model = a_model, + horizon = 7, samples = 10, + serial_interval = EpiSoon::example_serial_interval), + "Must have the same number of Rt and case samples." + ) +}) + +context("Test of the function 'evaluate_model'.... testing outputs") + +test_that("Outputs have proper lenghts and names", { + out <- evaluate_model(sampled_obs, + sampled_cases, + model = a_model, + horizon = 7, samples = 10, + serial_interval = EpiSoon::example_serial_interval) + + expect_length(out, 4) + + expect_named(out, c("forecast_rts", "rt_scores", "forecast_cases", "case_scores")) + + # raw outputs added + out <- evaluate_model(sampled_obs, + sampled_cases, + model = a_model, + horizon = 7, samples = 10, + serial_interval = EpiSoon::example_serial_interval, + return_raw = TRUE) + + expect_length(out, 6) + + expect_named(out, + c("forecast_rts", "rt_scores", "forecast_cases", "case_scores", + "raw_rt_forecast", "raw_case_forecast")) + +}) diff --git a/tests/testthat/test_forecast_cases.R b/tests/testthat/test_forecast_cases.R new file mode 100644 index 0000000..b2719bb --- /dev/null +++ b/tests/testthat/test_forecast_cases.R @@ -0,0 +1,152 @@ + +context("Test of forecast_cases()") + +# Fit and forecast to use (provided example) +forecast <- forecast_rt( + EpiSoon::example_obs_rts[1:10,], + model = function(...) { + EpiSoon::bsts_model( + model = function(ss, y){ + bsts::AddAutoAr(ss, y=y, lags=10) + }, + ... + ) + }, + horizon = 7, + samples = 10 +) + + + +test_that("forecast_cases() output is of expected format", { + # Provided example + out <- forecast_cases( + EpiSoon::example_obs_cases, + fit_samples = forecast, + serial_interval = EpiSoon::example_serial_interval + ) + + expect_s3_class(out, c("tbl_df", "tbl", "data.frame")) + expect_named(out, c("sample", "date", "cases", "horizon")) + expect_equal(nrow(forecast), nrow(out)) + + # FIXME : recommended to preserve original order + forecast <- forecast[ order(forecast$sample, forecast$date) ,] + expect_equal(forecast$sample, out$sample) + expect_equal(forecast$date, out$date) + expect_equal(forecast$horizon, out$horizon) + + # FIXME : recommended to preserve column classes (sample casted from integer to numeric) +}) + +test_that("forecast_cases() handles missing arguments as expected", { + + expect_error( + forecast_cases(), + "is missing" + ) + + # FIXME : recommended to print a more explicit message + expect_error( + forecast_cases( + fit_samples = forecast, + serial_interval = EpiSoon::example_serial_interval + ) + ) + + # FIXME : recommended to print a more explicit message + expect_error( + forecast_cases( + cases = EpiSoon::example_obs_cases, + serial_interval = EpiSoon::example_serial_interval + ) + ) + + # FIXME : recommended to print a more explicit message + expect_error( + forecast_cases( + cases = EpiSoon::example_obs_cases, + fit_samples = forecast + ) + ) + + expect_identical( + { + set.seed(42) + forecast_cases( + EpiSoon::example_obs_cases, + fit_samples = forecast, + serial_interval = EpiSoon::example_serial_interval + ) + },{ + set.seed(42) + forecast_cases( + EpiSoon::example_obs_cases, + fit_samples = forecast, + serial_interval = EpiSoon::example_serial_interval, + forecast_date = forecast$date[1] + ) + } + ) + + expect_identical( + { + set.seed(42) + forecast_cases( + EpiSoon::example_obs_cases, + fit_samples = forecast, + serial_interval = EpiSoon::example_serial_interval + ) + },{ + set.seed(42) + forecast_cases( + EpiSoon::example_obs_cases, + fit_samples = forecast, + serial_interval = EpiSoon::example_serial_interval, + horizon = max(forecast$horizon) + ) + } + ) + + expect_identical( + { + set.seed(42) + forecast_cases( + EpiSoon::example_obs_cases, + fit_samples = forecast, + serial_interval = EpiSoon::example_serial_interval + ) + },{ + set.seed(42) + forecast_cases( + EpiSoon::example_obs_cases, + fit_samples = forecast, + serial_interval = EpiSoon::example_serial_interval, + rdist = rpois + ) + } + ) + +}) + +# FIXME : 'horizon' appears to have no effect + +test_that("forecast_cases() can handle custom sampling functions", { + + expect <- function(FUN) { + expect_silent( + forecast_cases( + EpiSoon::example_obs_cases, + fit_samples = forecast, + serial_interval = EpiSoon::example_serial_interval, + rdist = FUN + ) + ) + } + + expect(rpois) + expect(rnorm) + expect(function(n,mean) { sample(-5:5+mean, n, replace=TRUE) }) + +}) + diff --git a/tests/testthat/test_forecast_rt.R b/tests/testthat/test_forecast_rt.R new file mode 100644 index 0000000..ca678ba --- /dev/null +++ b/tests/testthat/test_forecast_rt.R @@ -0,0 +1,40 @@ +##------------------------------------------------------------ +context("Test of the function 'forecast_rt'.... testing outputs") + +mod_test <- function(...) { + EpiSoon::bsts_model(model = function(ss, y){bsts::AddAutoAr(ss, y = y, lags = 10)}, ...) +} +rts_test <- EpiSoon::example_obs_rts[1:10, ] +horizon_test <- 7 +samples_test <- 10 +run_test <- forecast_rt(rts = rts_test, + model = mod_test, + horizon = horizon_test, + samples = samples_test, + bound_rt = TRUE, + timeout = 30) +test_that("Output has proper length, names, class", { + expect_s3_class(run_test, c("tbl_df", "tbl", "data.frame")) + expect_named(run_test) + expect_equal(nrow(run_test), horizon_test*samples_test) + expect_length(run_test, 4) + expect_named(run_test, c("sample", "date", "rt", "horizon")) +}) +test_that("bound_rt = TRUE has appropriate bounding behavior", { + expect_gte(min(run_test$rt), 0) +}) + +##------------------------------------------------------------ +context("Test of the function 'forecast_rt'.... testing inputs") +misspecified_model <- function(...) { + lm(y~x) +} + +test_that("Outputs returns error for inappropriately specified model", { + expect_error(forecast_rt(rts = rts_test, + model = misspecified_model, + horizon = horizon_test, + samples = samples_test, + bound_rt = TRUE, + timeout = 30)) +}) diff --git a/tests/testthat/test_iterative_case_forecast.R b/tests/testthat/test_iterative_case_forecast.R new file mode 100644 index 0000000..d166aa5 --- /dev/null +++ b/tests/testthat/test_iterative_case_forecast.R @@ -0,0 +1,44 @@ + +context("Test of iterative_case_forecast()") + +# Iterative RT forecast to use (provided example) +it_forecast <- iterative_rt_forecast( + EpiSoon::example_obs_rts, + model = function(...) { + EpiSoon::bsts_model( + model = function(ss, y) { bsts::AddSemilocalLinearTrend(ss, y = y) }, + ... + ) + }, + horizon = 7, + samples = 10 +) + + + +test_that("iterative_case_forecast() output is of expected format", { + # Provided example + out <- iterative_case_forecast( + it_fit_samples = it_forecast, + cases = EpiSoon::example_obs_cases, + serial_interval = EpiSoon::example_serial_interval + ) + + expect_s3_class(out, c("tbl_df", "tbl", "data.frame")) + expect_named(out, c("forecast_date", "sample", "date", "cases", "horizon")) + expect_equal(nrow(it_forecast), nrow(out)) + + # FIXME : recommended to preserve original order + sorted_forecast <- it_forecast[ order(it_forecast$forecast_date, it_forecast$sample, it_forecast$date) ,] + expect_equal(sorted_forecast$sample, out$sample) + expect_equal(sorted_forecast$date, out$date) + expect_equal(sorted_forecast$horizon, out$horizon) + expect_equal(sorted_forecast$forecast_date, out$forecast_date) + + # FIXME : recommended to preserve column classes (sample casted from integer to numeric) +}) + +# FIXME : no 'horizon' argument to pass to forecast_cases() ? + +# 'rdist' tested in forecast_cases() (passed untouched) + diff --git a/tests/testthat/test_iterative_rt_forecast.R b/tests/testthat/test_iterative_rt_forecast.R new file mode 100644 index 0000000..96b8142 --- /dev/null +++ b/tests/testthat/test_iterative_rt_forecast.R @@ -0,0 +1,32 @@ +##------------------------------------------------------------ +context("Test of the function 'iterative_rt_forecast'.... testing outputs") + +rts_test <- EpiSoon::example_obs_rts +mod_test <- function(...) { + EpiSoon::bsts_model(model = function(ss, y) { + bsts::AddSemilocalLinearTrend(ss, y = y) + }, ...) +} +horizon_test <- 7 +samples_test <- 10 +min_points_test <- 4 + +run_test <- iterative_rt_forecast(rts_test, + model = mod_test, + horizon = horizon_test, + samples = samples_test, + min_points = min_points_test) + +test_that("Output has proper length, names, class", { + expect_s3_class(run_test, c("tbl_df", "tbl", "data.frame")) + expect_named(run_test) + expect_equal(nrow(run_test), (nrow(rts_test)-min_points_test)*horizon_test*samples_test) + expect_length(run_test, 5) + expect_named(run_test, c("forecast_date", "sample", "date", "rt", "horizon")) +}) +test_that("bound_rt = TRUE has appropriate bounding behavior", { + expect_gte(min(run_test$rt), 0) +}) +test_that("Earliest forecast date is `min_points`+1 days after the first rts input date", { + expect_equal(min(rts_test$date)+min_points_test+1, min(run_test$date)) +}) diff --git a/tests/testthat/test_plot_forecast.R b/tests/testthat/test_plot_forecast.R new file mode 100644 index 0000000..5c5d1f6 --- /dev/null +++ b/tests/testthat/test_plot_forecast.R @@ -0,0 +1,28 @@ +##------------------------------------------------------------ +context("Test of the function 'plot_forecast'.... testing output") + +samples <- forecast_rt(EpiSoon::example_obs_rts[1:10, ], + model = function(...) { + EpiSoon::bsts_model(model = function(ss, y) { + bsts::AddSemilocalLinearTrend(ss, y = y) + }, ...)}, + horizon = 21, samples = 10) + +test_that("Plot is obtained for summarized forecast", { + summarised_forecast <- summarise_forecast(samples) + p <- plot_forecast(summarised_forecast, EpiSoon::example_obs_rts) + + expect_s3_class(p, c("gg", "ggplot")) + expect_silent(plot(p)) +}) + +test_that("Plot is obtained for summarized cases", { + pred_cases <- forecast_cases(EpiSoon::example_obs_cases, samples, + serial_interval = EpiSoon::example_serial_interval) + + summarised_case_forecast <- summarise_case_forecast(pred_cases) + p <- plot_forecast(summarised_case_forecast, EpiSoon::example_obs_cases) + + expect_s3_class(p, c("gg", "ggplot")) + expect_silent(plot(p)) +}) diff --git a/tests/testthat/test_plot_forecast_evaluation.R b/tests/testthat/test_plot_forecast_evaluation.R new file mode 100644 index 0000000..fde6511 --- /dev/null +++ b/tests/testthat/test_plot_forecast_evaluation.R @@ -0,0 +1,20 @@ +##------------------------------------------------------------ +context("Test of the function 'plot_forecast_evaluation'.... testing output") + +forecast_eval <- evaluate_model(EpiSoon::example_obs_rts, + EpiSoon::example_obs_cases, + model = function(...) { + EpiSoon::bsts_model(model = function(ss, y) { + bsts::AddSemilocalLinearTrend(ss, y = y) + }, ...)}, + serial_interval = EpiSoon::example_serial_interval, + horizon = 7, samples = 10) + +test_that("Plot is obtained for outputs of models", { + p <- plot_forecast_evaluation(forecast_eval$forecast_rts, + EpiSoon::example_obs_rts, + horizon_to_plot = 7) + + expect_s3_class(p, c("gg", "ggplot")) + expect_silent(plot(p)) +}) diff --git a/tests/testthat/test_predict_cases.R b/tests/testthat/test_predict_cases.R new file mode 100644 index 0000000..aac7408 --- /dev/null +++ b/tests/testthat/test_predict_cases.R @@ -0,0 +1,41 @@ +# predict_cases +# This function predicts the cases from one sample of forecast Rt. It needs a +# data frame with two columns (date and cases), another data frame with the +# columns rt and date (rt is numerical and date is a date format), a vector of +# numerical values that describe the probability distribution of the interval. +# The other arguments does not seem to be mandatory (to check) +# Things to test: +## In tests +# - that the function returns the proper result +# - that the result class is correct + +# set.seed(1234) + +forecast <- forecast_rt(EpiSoon::example_obs_rts[1:10, ], model = function(...){ + EpiSoon::bsts_model(model = function(ss, y){ + bsts::AddSemilocalLinearTrend(ss, y = y) + }, ...) + }, horizon = 7, samples = 1) + +predictedCases <- predict_cases( + cases = EpiSoon::example_obs_cases, + rts = forecast, + forecast_date = as.Date("2020-03-10"), + serial_interval = EpiSoon::example_serial_interval +) + +# expectedTable <- tibble::as_tibble( +# data.frame( +# date = as.Date(c("2020-03-11", "2020-03-12", "2020-03-13", "2020-03-14", +# "2020-03-15", "2020-03-16", "2020-03-17")), +# cases = c(120, 204, 132, 131, 219, 151, 193) +# ) +# ) + +test_that("The expected Rt sample forecasts predict cases are obtained", { + expect_s3_class(predictedCases, "data.frame") + expect_named(predictedCases, c("date", "cases")) + expect_length(predictedCases$date, 7) # must be equal to horizon x length(forecast_date) + # expect_equal(predictedCases$date, expectedTable$date) + # expect_gte(cor(predictedCases$cases, expectedTable$cases), .9) +}) diff --git a/tests/testthat/test_predict_current_cases.R b/tests/testthat/test_predict_current_cases.R new file mode 100644 index 0000000..aad390a --- /dev/null +++ b/tests/testthat/test_predict_current_cases.R @@ -0,0 +1,30 @@ +# predict_current_cases +# same remark as for the tests of 'test_predict_cases' +predictedCurrentCases <- predict_current_cases( + cases = EpiSoon::example_obs_cases, + rts = EpiSoon::example_obs_rts, + serial_interval = EpiSoon::example_serial_interval +) + +# expectedCurrentCases <- data.frame( +# rt = EpiSoon::example_obs_rts$rt, +# date = EpiSoon::example_obs_rts$date, +# infectiousness = c(14.145, 14.495, 14.925, 15.435, 16.635, 18.705, 22.805, +# 27.555, 41.465, 60.965, 74.735, 104.655, 151.405, 158.155, +# 158.155, 160.445, 162.955, 165.845, 172.575, 184.005, +# 209.880, 223.155), +# cases = c(30, 53, 36, 37, 33, 39, 54, 63, 90, 119, 152, 201, 302, 337, 318, +# 282, 320, 296, 298, 311, 359, 349) +# ) + +test_that("The expected Rts based on observed data predict cases are obtained", { + expect_s3_class(predictedCurrentCases, "data.frame") + expect_named(predictedCurrentCases, + c("rt", "date", "infectiousness", "cases")) + expect_length(predictedCurrentCases$date, nrow(EpiSoon::example_obs_rts)) + # expect_equal(predictedCurrentCases$rt, expectedCurrentCases$rt) + # expect_equal(predictedCurrentCases$date, expectedCurrentCases$date) + # expect_equal(predictedCurrentCases$infectiousness, + # expectedCurrentCases$infectiousness) + # expect_gte(cor(predictedCurrentCases$cases, expectedCurrentCases$cases), .9) +}) diff --git a/tests/testthat/test_score_case_forecast.R b/tests/testthat/test_score_case_forecast.R new file mode 100644 index 0000000..04c0677 --- /dev/null +++ b/tests/testthat/test_score_case_forecast.R @@ -0,0 +1,30 @@ +##------------------------------------------------------------ +context("Test of the function 'score_case_forecast'.... testing outputs") + +rts_test <- EpiSoon::example_obs_rts[1:10, ] +mod_test <- function(...) { + EpiSoon::bsts_model(model = function(ss, y){ + bsts::AddSemilocalLinearTrend(ss, y = y) + }, ...) +} +samples_test <- forecast_rt(rts = rts_test, + model = mod_test, + horizon = 7, + samples = 10) + +case_test <- EpiSoon::example_obs_cases +serial_interval_test <- EpiSoon::example_serial_interval +pred_cases <- forecast_cases(cases = case_test, + fit_samples = samples_test, + serial_interval = serial_interval_test) + +run_test <- score_case_forecast(pred_cases, case_test) + +test_that("Output has proper length, names, class", { + expect_s3_class(run_test, c("tbl_df", "tbl", "data.frame")) + expect_named(run_test) + expect_equal(nrow(run_test), length(unique(pred_cases$date))) + expect_named(run_test, c("date", "horizon", "dss", "crps", "logs", "bias", + "sharpness", "calibration", "median", "iqr", "ci")) + expect_equal(nrow(run_test), 7) +}) diff --git a/tests/testthat/test_score_forecast.R b/tests/testthat/test_score_forecast.R new file mode 100644 index 0000000..6bce938 --- /dev/null +++ b/tests/testthat/test_score_forecast.R @@ -0,0 +1,25 @@ +##------------------------------------------------------------ +context("Test of the function 'score_forecast'.... testing outputs") + +## Fit a model (using a subset of observations) +samples <- forecast_rt(EpiSoon::example_obs_rts[1:10, ], + model = function(...) { + EpiSoon::bsts_model(model = function(ss, y){ + bsts::AddSemilocalLinearTrend(ss, y = y) + }, ...) + }, + horizon = 7, samples = 10) + +## Score the model fit (with observations during the time horizon of the forecast) +obs_test <- EpiSoon::example_obs_rts +run_test <- score_forecast(samples, + observations = obs_test) + +test_that("Output has proper length, names, class", { + expect_s3_class(run_test, c("tbl_df", "tbl", "data.frame")) + expect_named(run_test) + expect_equal(nrow(run_test), length(intersect(samples$date, obs_test$date))) + expect_named(run_test, c("date", "horizon", "dss", "crps", "logs", "bias", + "sharpness", "calibration", "median", "iqr", "ci")) +}) + diff --git a/tests/testthat/test_summarise_case_forecast.R b/tests/testthat/test_summarise_case_forecast.R new file mode 100644 index 0000000..6bef8cc --- /dev/null +++ b/tests/testthat/test_summarise_case_forecast.R @@ -0,0 +1,77 @@ + +context("Test of summarise_forecast() and summarise_case_forecast()") + +# FIXME : recommended to define classes and corresponding summary() methods to guarantee input class and comply to R standards + +# horizon = 7 (provided example) +forecast_h7 <- forecast_rt( + example_obs_rts, + model = function(...) { + EpiSoon::bsts_model( + model = function(ss, y) { bsts::AddSemilocalLinearTrend(ss, y = y) }, + ... + ) + }, + horizon = 7, + samples = 10 +) +case_forecast_h7 <- forecast_cases( + EpiSoon::example_obs_cases, + forecast_h7, + EpiSoon::example_serial_interval +) +summary_forecast_h7 <- summarise_forecast(forecast_h7) +summary_case_forecast_h7 <- summarise_case_forecast(case_forecast_h7) + +# horizon = 1 +forecast_h1 <- forecast_rt( + example_obs_rts, + model = function(...) { + EpiSoon::bsts_model( + model = function(ss, y) { bsts::AddSemilocalLinearTrend(ss, y = y) }, + ... + ) + }, + horizon = 1, + samples = 10 +) +case_forecast_h1 <- forecast_cases( + EpiSoon::example_obs_cases, + forecast_h1, + EpiSoon::example_serial_interval +) +summary_forecast_h1 <- summarise_forecast(forecast_h1) +summary_case_forecast_h1 <- summarise_case_forecast(case_forecast_h1) + + + +test_that("summarise_forecast() output is of expected format", { + expect_s3_class(summary_forecast_h7, c("tbl_df", "tbl", "data.frame")) + expect_named(summary_forecast_h7, c("date", "horizon", "median", "mean", "sd", "bottom", "lower", "upper", "top")) +}) + +test_that("summarise_case_forecast() output is of expected format", { + expect_s3_class(summary_case_forecast_h7, c("tbl_df", "tbl", "data.frame")) + expect_named(summary_case_forecast_h7, c("date", "horizon", "median", "mean", "sd", "bottom", "lower", "upper", "top")) +}) + +test_that("summarise_forecast() rows vary with 'horizon'", { + # Ascending horizon + expect_equal(summary_forecast_h1$horizon, 1) + expect_equal(summary_forecast_h7$horizon, 1:7) + + # Ascending date + expect_equal(summary_forecast_h1$date, sort(unique(forecast_h1$date))) + expect_equal(summary_forecast_h7$date, sort(unique(forecast_h7$date))) +}) + +test_that("summarise_case_forecast() rows vary with 'horizon'", { + # Ascending horizon + expect_equal(summary_case_forecast_h1$horizon, 1) + expect_equal(summary_case_forecast_h7$horizon, 1:7) + + # Ascending date + expect_equal(summary_case_forecast_h1$date, sort(unique(forecast_h1$date))) + expect_equal(summary_case_forecast_h7$date, sort(unique(forecast_h7$date))) +}) + diff --git a/tests/testthat/test_summarise_scores.R b/tests/testthat/test_summarise_scores.R new file mode 100644 index 0000000..d2e5665 --- /dev/null +++ b/tests/testthat/test_summarise_scores.R @@ -0,0 +1,146 @@ +##------------------------------------------------------------ +context("Test of the function 'summarise_scores'.... testing outputs") + +rts <- EpiSoon::example_obs_rts %>% + dplyr::mutate(timeseries = "Region 1") %>% + dplyr::bind_rows(EpiSoon::example_obs_rts %>% + dplyr::mutate(timeseries = "Region 2")) + +cases <- EpiSoon::example_obs_cases %>% + dplyr::mutate(timeseries = "Region 1") %>% + dplyr::bind_rows(EpiSoon::example_obs_cases %>% + dplyr::mutate(timeseries = "Region 2")) + +models <- list("AR 3" = function(...) { + EpiSoon::bsts_model(model = function(ss, y) { + bsts::AddAr(ss, y = y, lags = 3) + }, ...)}, + "Semi-local linear trend" = function(...) { + EpiSoon::bsts_model(model = function(ss, y) { + bsts::AddSemilocalLinearTrend(ss, y = y) + }, ...)} +) + +out <- compare_timeseries(rts, cases, models, horizon = 7, samples = 10, + serial_interval = EpiSoon::example_serial_interval) +exp_scores <- setdiff(names(out$rt_scores), + c("timeseries", "model", "forecast_date", "date", + "horizon")) + +test_that("Default outputs have proper lengths and names", { + summary_out <- summarise_scores(out$rt_scores) + + # test that the summary for all scores are obtained for every model + expect_equal(nrow(summary_out), length(models) * length(exp_scores)) + + expect_named(summary_out, c("score", "model", "bottom", "lower", "median", + "mean", "upper", "top", "sd")) + + expect_s3_class(summary_out, c("tbl_df", "tbl", "data.frame")) +}) + +test_that("Summary over variables has proper lengths and names", { + ## horizon + summary_horizon <- summarise_scores(out$rt_scores, "horizon") + + # summary for all scores are obtained for every horizon and every model + expect_equal(nrow(summary_horizon), + length(models) * length(exp_scores) * max(out$rt_scores$horizon)) + + expect_named(summary_horizon, c("horizon", "score", "model", "bottom", "lower", + "median", "mean", "upper", "top", "sd")) + + ## timeseries + summary_timeseries <- summarise_scores(out$rt_scores, "timeseries") + + # summary for all scores are obtained for every region and every model + nb_timeseries <- length(unique(out$rt_scores$timeseries)) + expect_equal(nrow(summary_timeseries), + length(models) * length(exp_scores) * nb_timeseries) + + expect_named(summary_timeseries, c("timeseries", "score", "model", "bottom", + "lower", "median", "mean", "upper", "top", + "sd")) + + ## timeseries AND horizon + summary_variables <- summarise_scores(out$rt_scores, c("timeseries", "horizon")) + + # summary for all scores are obtained for every horizon, every region, every model + nb_timeseries <- length(unique(out$rt_scores$timeseries)) + expect_equal(nrow(summary_variables), + length(models) * length(exp_scores) * + max(out$rt_scores$horizon) * nb_timeseries) + + expect_named(summary_variables, c("timeseries", "horizon", "score", "model", + "bottom", "lower", "median", "mean", "upper", + "top", "sd")) +}) + +## Instead summarise across region and summarise case scores +test_that("Specific score summary has proper lengths and names", { + ## one score + summary_logs <- summarise_scores(out$case_scores, sel_scores = "logs") + + # summary for all scores are obtained for every horizon and every model + expect_equal(nrow(summary_logs), length(models)) + + ## more than one score + summary_scores <- summarise_scores(out$case_scores, + sel_scores = c("logs", "dss")) + + # summary for all scores are obtained for every horizon and every model + expect_equal(nrow(summary_scores), length(models) * 2) +}) + +test_that("Summary are handled properly when there is just one model or no field 'model'", { + ## one model + one_model <- out$rt_scores %>% filter(model == "Semi-local linear trend") + + summary_one_model <- summarise_scores(one_model) + + # test that the summary for all scores are obtained for every model + expect_equal(nrow(summary_one_model), length(exp_scores)) + + ## no field 'model' + out$rt_scores$model <- NULL + + summary_no_model <- summarise_scores(out$rt_scores) + + # test that the summary for all scores are obtained for every model + expect_equal(nrow(summary_no_model), length(exp_scores)) + + expect_named(summary_no_model, c("score", "bottom", "lower", "median", "mean", + "upper", "top", "sd")) +}) + +test_that("Summary works for outputs of 'evaluate_model'", { + sampled_obs <- EpiSoon::example_obs_rts %>% + dplyr::mutate(sample = 1) %>% + dplyr::bind_rows(EpiSoon::example_obs_rts %>% + dplyr::mutate(sample = 2)) + sampled_cases <- EpiSoon::example_obs_cases %>% + dplyr::mutate(sample = 1) %>% + dplyr::bind_rows(EpiSoon::example_obs_cases %>% + dplyr::mutate(sample = 2)) + a_model <- function(...) { + EpiSoon::bsts_model(model = function(ss, y) { + bsts::AddSemilocalLinearTrend(ss, y = y) + }, ...) + } + + out <- evaluate_model(sampled_obs, sampled_cases, model = a_model, + horizon = 7, samples = 10, + serial_interval = EpiSoon::example_serial_interval) + + summary_models <- summarise_scores(out$rt_scores) + + # test that the summary for all scores are obtained for every model + exp_scores <- setdiff(names(out$rt_scores), + c("sample", "model", "forecast_date", "date", "horizon")) + expect_equal(nrow(summary_models), length(exp_scores)) + + expect_named(summary_models, c("score", "bottom", "lower", "median", "mean", + "upper", "top", "sd")) + + expect_s3_class(summary_models, c("tbl_df", "tbl", "data.frame")) +})