In HIDDA/forecasting: Forecasting Based on Surveillance Data

knitr::opts_chunk$set(message = FALSE, warning = FALSE, error = FALSE,
                      fig.align = "center", dev.args = list(pointsize = 10))
options(knitr.duplicate.label = "allow")

options(digits = 4)  # for more compact numerical outputs
library("surveillance")
library("hhh4contacts")

In this vignette, we assess long-term forecasts of the "hhh4" models from vignette("BNV") using:

library("hhh4addon")

This package is developed by Johannes Bracher and available from r packageDescription("hhh4addon", fields = "URL").

Models

We rerun the code from the vignette("BNV") to make the models available:

rmarkdown::render(system.file("doc", "BNV.Rmd", package = "HIDDA.forecasting"),
                  run_pandoc = FALSE, envir = BNV <- new.env(), quiet = TRUE)

(MODELS <- names(BNV$fits))
NMODELS <- length(MODELS)

Predictive distributions

statmoms <- lapply(BNV$fits, stationary_moments, start = 27)

## illustration for age group 00-04
.group <- 1
par(mfrow = n2mfrow(NMODELS), mar = c(5,5,2,1))
invisible(mapply(fanplot_stationary, statmoms, main = names(statmoms),
                 MoreArgs = list(unit = .group, xlab = "Time within season",
                                 mean_col = "white", mean_lty = 1)))

predmoms <- lapply(BNV$fits, predictive_moments, t_condition = 208, lgt = 52,
                   return_Sigma = TRUE)

## illustration for age group 00-04
.group <- 1
par(mfrow = n2mfrow(NMODELS), mar = c(5,5,2,1), las = 1)
invisible(mapply(fanplot_prediction, predmoms, main = names(predmoms),
                 MoreArgs = list(unit = .group, xlab = "Time", ylim = c(0, 55),
                                 mean_col = "white", mean_lty = 1,
                                 l.col = NA, pt.col = 1, pt.cex = 0.4)))

Scaled Dawid-Sebastiani scores of these multivariate forecast distributions:

sapply(predmoms, ds_score_hhh4, detailed = TRUE)

Final size forecasts (aggregate over time)

aggr <- matrix(rep(diag(BNV$NGROUPS), 52), nrow = BNV$NGROUPS,
               dimnames = list(BNV$GROUPS, NULL))
predmoms_aggr <- lapply(predmoms, aggregate_moments, aggregation_matrix = aggr)

par(mfrow = n2mfrow(NMODELS), mar = c(5,5,2,1), las = 2, cex.axis = 0.8)
invisible(mapply(plot_moments_by_unit, predmoms_aggr, main = names(predmoms_aggr),
                 MoreArgs = list(xlab = "Age group", ylim = c(0, 3000),
                                 pt.col = 2, pt.cex = 0.6)))

Scaled Dawid-Sebastiani scores of these multivariate forecast distributions:

sapply(predmoms_aggr, ds_score_hhh4, detailed = TRUE)

HIDDA/forecasting documentation built on Jan. 11, 2024, 1:11 a.m.