Citation
- Source:inst/CITATION
+
+ Citation
+ Source:inst/CITATION
+ Held L, Meyer S (2019). -“Forecasting Based on Surveillance Data.” +“Forecasting Based on Surveillance Data.” In Held L, Hens N, O'Neill PD, Wallinga J (eds.), Handbook of Infectious Disease Data Analysis. -Chapman & Hall/CRC. +Chapman \& Hall/CRC.
@InCollection{R:HIDDA.forecasting,
author = {Leonhard Held and Sebastian Meyer},
title = {Forecasting Based on Surveillance Data},
booktitle = {Handbook of Infectious Disease Data Analysis},
editor = {Leonhard Held and Niel Hens and Philip D. O'Neill and Jacco Wallinga},
- publisher = {Chapman & Hall/CRC},
+ publisher = {Chapman \& Hall/CRC},
+ chapter = {25},
year = {2019},
}
- Authors
--
-
-
-
Sebastian Meyer. Author, maintainer. -
-
- -
-
Leonhard Held. Contributor. -
-
-
-HIDDA.forecasting 1.1.1 (2021-03-31)
--
-
-
+ ++
HIDDA.forecasting 1.1.2 (2023-11-29)
+Vignettes have been rebuilt using up-to-date versions of all involved packages in R 4.3.2 (on a new machine). This only gave minor numerical differences in the long-term forecasts of
vignette("CHILI_prophet").
+Accommodate restricted checks without suggested packages.
+
+-HIDDA.forecasting 1.1.1 (2021-03-31)
+-
Vignettes have been rebuilt using up-to-date versions of all involved packages in R 4.0.4, resulting in the following changes:
--
-
vignette("BNV"): age-specific amplitude-shift parameter transformations were wrong insummary(mg_Cpower); bug fixed in surveillance 1.18.0.
+vignette("BNV"): age-specific amplitude-shift parameter transformations were wrong insummary(mg_Cpower); bug fixed in surveillance 1.18.0.surveillance::pit()values for degenerate forecasts with zero probability for observed counts may differ due to a change in surveillance 1.17.1 (they still produce warnings).vignette("CHILI_hhh4"): no rounding ofnfrom 213 to 210 in printedcalibrationTest()(a bug fixed in R 3.6.0).vignette("CHILI_prophet"): minor numerical differences in model fit and predictions due to changes in prophet.
-
-
-+-HIDDA.forecasting 1.1.0 (2019-03-29)
--
-
Use standard PIT for continuous forecasts (
arima,prophet,naive). Differences to the previously used non-randomized PIT histograms for count data are negligible.
+
+- +HIDDA.forecasting 1.1.0 (2019-03-29)
+Use standard PIT for continuous forecasts (
arima,prophet,naive). Differences to the previously used non-randomized PIT histograms for count data are negligible.Add scores for discretized log-normal forecasts, via new function
scores_lnorm_discrete(). These scores are almost identical to the continuous scores, essentially due to the large counts.Vignettes have been rebuilt using up-to-date versions of all involved packages (forecast 8.5, glarma 1.6.0, hhh4contacts 0.13.0, prophet 0.4, scoringRules 0.9.5, surveillance 1.17.0) in R 3.5.3.
-
++HIDDA.forecasting 1.0.0 (2018-09-04)
+This is the version used for the book chapter.
-
The contained vignettes have been built using R 3.5.1 with all dependent packages’ versions as of 25 July 2018 from CRAN. The versions of the main packages were:
--
-
-
+
- forecast 8.4
- glarma 1.6.0 @@ -167,39 +88,31 @@
scoringRules 0.9.4
- surveillance 1.16.2 -
- -
+
Swiss Surveillance Data on Influenza Like Illness, 2000-2016
@@ -135,60 +59,61 @@Swiss Surveillance Data on Influenza Like Illness, 2000-2016
estimated from the Swiss Sentinella Reporting System.data("CHILI")- - -
Format
+data("CHILI")a univariate time series of class zoo,
- where the time index is of class Date
+
Format
+a univariate time series of class zoo,
+ where the time index is of class Date
and always refers to the Tuesday of the notification week
Source
- +Source
The Swiss ILI data has been received on 19 January 2017 by courtesy of:
-Swiss Federal Office of Public Health
- Public Health Directorate
- Communicable Diseases Division
- 3003 Bern
+
Swiss Federal Office of Public Health
+ Public Health Directorate
+ Communicable Diseases Division
+ 3003 Bern
SWITZERLAND
Examples
-+#> Index CHILI -#> Min. :2000-01-04 Min. : 31 -#> 1st Qu.:2004-04-02 1st Qu.: 328 -#> Median :2008-07-01 Median : 973 -#> Mean :2008-07-01 Mean : 4140 -#> 3rd Qu.:2012-09-28 3rd Qu.: 3351 -#> Max. :2016-12-27 Max. :48473
The function dhhh4sims constructs a (non-vectorized)
probability mass function from the result of
-surveillance::simulate.hhh4() (and the corresponding
+surveillance::simulate.hhh4() (and the corresponding
model), as a function of the time point within the simulation period.
The distribution at each time point is obtained as a mixture of
negative binomial (or Poisson) distributions based on the samples from
the previous time point.
dhhh4sims(sims, model)+
dhhh4sims(sims, model)Arguments
+- sims +
a
"hhh4sims"object fromsurveillance::simulate.hhh4().
+
- - model +
the "hhh4" object underlying
sims.
-
Arguments
-| sims | -a |
-
|---|---|
| model | -the "hhh4" object underlying |
-
Value
+Value
+ -a function(x, tp = 1, log = FALSE), which takes a
+
a function(x, tp = 1, log = FALSE), which takes a
vector of model$nUnit counts and calculates the
(log-)probability of observing these counts (given the
model) at the tp'th time point of the simulation
-period (index or character string matching rownames(sims)).
See also
- -logs_hhh4sims() where this function is used.
Author
- +period (index or character string matchingrownames(sims)).
+ See also
+logs_hhh4sims() where this function is used.
Examples
-+#>#>#>CHILI.sts <- sts(observed = CHILI, - epoch = as.integer(index(CHILI)), epochAsDate = TRUE) - -## fit a simple hhh4 model -(f1 <- addSeason2formula(~ 1, period = 365.2425)) -#> ~1 + sin(2 * pi * t/365.2425) + cos(2 * pi * t/365.2425)fit <- hhh4( - stsObj = CHILI.sts, - control = list(ar = list(f = f1), end = list(f = f1), family = "NegBin1") -) - -## simulate the last four weeks (only 200 runs, for speed) -sims <- simulate(fit, nsim = 200, seed = 1, subset = 884:nrow(CHILI.sts), - y.start = observed(CHILI.sts)[883,]) -if (requireNamespace("fanplot")) { - plot(sims, "fan", fan.args = list(ln = c(5,95)/100), - observed.args = list(pch = 19), means.args = list(type = "b")) -} -#>-## derive the weekly forecast distributions -dfun <- dhhh4sims(sims, fit) -dfun(4000, tp = 1) -#> [1] 0.0002194578dfun(4000, tp = 4) -#> [1] 0.0001069622curve(sapply(x, dfun, tp = 4), 0, 30000, type = "h", - main = "4-weeks-ahead forecast", - xlab = "No. infected", ylab = "Probability") --## compare the forecast distributions with the simulated counts -par(mfrow = n2mfrow(nrow(sims))) -for (tp in 1:nrow(sims)) { - MASS::truehist(sims[tp,,], xlab = "counts", ylab = "Probability") - curve(sapply(x, dfun, tp = tp), add = TRUE, lwd = 2) -} --# \dontshow{ -## verify distribution at the first time point (i.e., one-step-ahead NegBin) -stopifnot(identical( - sapply(0:100, dfun, tp = 1), - dnbinom(0:100, - mu = meanHHH(fit$coefficients, terms(fit), subset = 884, total.only = TRUE), - size = sizeHHH(fit$coefficients, terms(fit), subset = 884)) -)) -## check that we have a probability distribution at the last time point -.xgrid <- seq(0, 100000, by = 500) -stopifnot(abs(1 - - integrate(approxfun(.xgrid, sapply(.xgrid, dfun, tp = 4)), 0, 100000)$value -) < 0.001) -# } -
Examples
+library("surveillance")
+#> Loading required package: sp
+#> Loading required package: xtable
+#> This is surveillance 1.22.1; see ‘package?surveillance’ or
+#> https://surveillance.R-Forge.R-project.org/ for an overview.
+CHILI.sts <- sts(observed = CHILI,
+ epoch = as.integer(index(CHILI)), epochAsDate = TRUE)
+
+## fit a simple hhh4 model
+(f1 <- addSeason2formula(~ 1, period = 365.2425))
+#> ~1 + sin(2 * pi * t/365.2425) + cos(2 * pi * t/365.2425)
+fit <- hhh4(
+ stsObj = CHILI.sts,
+ control = list(ar = list(f = f1), end = list(f = f1), family = "NegBin1")
+)
+
+## simulate the last four weeks (only 200 runs, for speed)
+sims <- simulate(fit, nsim = 200, seed = 1, subset = 884:nrow(CHILI.sts),
+ y.start = observed(CHILI.sts)[883,])
+if (requireNamespace("fanplot")) {
+ plot(sims, "fan", fan.args = list(ln = c(5,95)/100),
+ observed.args = list(pch = 19), means.args = list(type = "b"))
+}
+#> Loading required namespace: fanplot
+
+
+## derive the weekly forecast distributions
+dfun <- dhhh4sims(sims, fit)
+dfun(4000, tp = 1)
+#> [1] 0.0002194578
+dfun(4000, tp = 4)
+#> [1] 0.0001069622
+curve(sapply(x, dfun, tp = 4), 0, 30000, type = "h",
+ main = "4-weeks-ahead forecast",
+ xlab = "No. infected", ylab = "Probability")
+
+
+## compare the forecast distributions with the simulated counts
+par(mfrow = n2mfrow(nrow(sims)))
+for (tp in 1:nrow(sims)) {
+ MASS::truehist(sims[tp,,], xlab = "counts", ylab = "Probability")
+ curve(sapply(x, dfun, tp = tp), add = TRUE, lwd = 2)
+}
+
+
+# \dontshow{
+## verify distribution at the first time point (i.e., one-step-ahead NegBin)
+stopifnot(identical(
+ sapply(0:100, dfun, tp = 1),
+ dnbinom(0:100,
+ mu = meanHHH(fit$coefficients, terms(fit), subset = 884, total.only = TRUE),
+ size = sizeHHH(fit$coefficients, terms(fit), subset = 884))
+))
+## check that we have a probability distribution at the last time point
+.xgrid <- seq(0, 100000, by = 500)
+stopifnot(abs(1 -
+ integrate(approxfun(.xgrid, sapply(.xgrid, dfun, tp = 4)), 0, 100000)$value
+) < 0.001)
+# }
+Simulation-Based Forecast Distributions
as a mixture of negative binomial (or Poisson) distributions.dnbmix(means, size = NULL)- -
Arguments
-| means | -a |
-
|---|---|
| size | -the dispersion parameter of the
+
+
+
+ Arguments+
|
-
n.ahead).
- Value
+Value
+ -a function(x, tp = 1, log = FALSE), which takes a vector of
+
a function(x, tp = 1, log = FALSE), which takes a vector of
counts x and calculates the (log-)probabilities of observing
each of these numbers at the tp'th time point of the simulation
period (indexing the rows of means).
See also
- -logs_nbmix() where this function is used.
Author
- +See also
+logs_nbmix() where this function is used.
Examples
-+## a GLARMA example -library("glarma") -y <- as.vector(CHILI) - -## fit a simple NegBin-GLARMA model -X <- t(sapply(2*pi*seq_along(y)/52.1775, - function (x) c(sin = sin(x), cos = cos(x)))) -X <- cbind(intercept = 1, X) -fit <- glarma(y = y[1:883], X = X[1:883,], type = "NegBin", phiLags = 1) - -## simulate the last four weeks (only 500 runs, for speed) -set.seed(1) -means <- replicate(500, { - forecast(fit, n.ahead = 4, newdata = X[884:887,], newoffset = rep(0,4))$mu -}) - -## derive the weekly forecast distributions -dfun <- dnbmix(means, coef(fit, type = "NB")) -dfun(4000, tp = 1) -#> [1] 0.0001458309dfun(4000, tp = 4) -#> [1] 7.550036e-05curve(dfun(x, tp = 4), 0, 30000, type = "h", - main = "4-weeks-ahead forecast", - xlab = "No. infected", ylab = "Probability") --# \dontshow{ -## verify distribution at the first time point (i.e., one-step-ahead NegBin) -stopifnot(all.equal( - dfun(0:100, tp = 1), - dnbinom(0:100, - mu = forecast(fit, n.ahead=1, newdata=X[884,,drop=FALSE])$mu, - size = coef(fit, type = "NB")) -)) -## check that we have a probability distribution at the second time point -.xgrid <- seq(0, 200000, by = 500) -stopifnot(abs(1 - - integrate(approxfun(.xgrid, dfun(.xgrid, tp = 2)), 0, 200000)$value -) < 0.01) -# } -
Examples
+## a GLARMA example
+library("glarma")
+y <- as.vector(CHILI)
+
+## fit a simple NegBin-GLARMA model
+X <- t(sapply(2*pi*seq_along(y)/52.1775,
+ function (x) c(sin = sin(x), cos = cos(x))))
+X <- cbind(intercept = 1, X)
+fit <- glarma(y = y[1:883], X = X[1:883,], type = "NegBin", phiLags = 1)
+
+## simulate the last four weeks (only 500 runs, for speed)
+set.seed(1)
+means <- replicate(500, {
+ forecast(fit, n.ahead = 4, newdata = X[884:887,], newoffset = rep(0,4))$mu
+})
+
+## derive the weekly forecast distributions
+dfun <- dnbmix(means, coef(fit, type = "NB"))
+dfun(4000, tp = 1)
+#> [1] 0.0001458309
+dfun(4000, tp = 4)
+#> [1] 7.550036e-05
+curve(dfun(x, tp = 4), 0, 30000, type = "h",
+ main = "4-weeks-ahead forecast",
+ xlab = "No. infected", ylab = "Probability")
+
+
+# \dontshow{
+## verify distribution at the first time point (i.e., one-step-ahead NegBin)
+stopifnot(all.equal(
+ dfun(0:100, tp = 1),
+ dnbinom(0:100,
+ mu = forecast(fit, n.ahead=1, newdata=X[884,,drop=FALSE])$mu,
+ size = coef(fit, type = "NB"))
+))
+## check that we have a probability distribution at the second time point
+.xgrid <- seq(0, 200000, by = 500)
+stopifnot(abs(1 -
+ integrate(approxfun(.xgrid, dfun(.xgrid, tp = 2)), 0, 200000)$value
+) < 0.01)
+# }
+