### univariate salmonella agona count time series
data("salmonella.agona")
## convert from old "disProg" to new "sts" class
salmonella <- disProg2sts(salmonella.agona)
## generate formula for temporal and seasonal trends
f.end <- addSeason2formula(~1 + t, S=1, period=52)
model <- list(ar = list(f = ~1), end = list(f = f.end), family = "NegBin1")
## fit the model
result <- hhh4(salmonella, model)
## do sequential one-step-ahead predictions for the last 5 weeks
pred <- oneStepAhead(result, nrow(salmonella)-5, type="rolling",
which.start="final", verbose=FALSE)
pred
quantile(pred)
confint(pred)
## simple plot of the 80% one-week-ahead prediction interval
## and point forecasts
if (requireNamespace("fanplot"))
plot(pred, probs = c(.1,.9), means.args = list())
# \dontshow{
## test equivalence of parallelized version
if (.Platform$OS.type == "unix" && isTRUE(parallel::detectCores() > 1))
stopifnot(identical(pred,
oneStepAhead(result, nrow(salmonella)-5, type="rolling",
which.start="final", verbose=FALSE, cores=2)))
# }
## note: oneStepAhead(..., type="final") just means fitted values
stopifnot(identical(
unname(oneStepAhead(result, nrow(salmonella)-5, type="final")$pred),
unname(tail(fitted(result), 5))))
## compute scores of the one-step-ahead predictions
(sc <- scores(pred))
## the above uses the scores-method for "oneStepAhead" predictions,
## which is a simple wrapper around the default method:
scores(x = pred$observed, mu = pred$pred, size = exp(pred$psi))
## scores with respect to the fitted values are similar
(scFitted <- scores(result, subset = nrow(salmonella)-(4:0)))
# \dontshow{
## test that scFitted is equivalent to scores(oneStepAhead(..., type = "final"))
stopifnot(all.equal(
scFitted,
scores(oneStepAhead(result, nrow(salmonella)-5, type="final")),
check.attributes = FALSE))
# }
## test if the one-step-ahead predictions are calibrated
calibrationTest(pred) # p = 0.8746
## the above uses the calibrationTest-method for "oneStepAhead" predictions,
## which is a simple wrapper around the default method:
calibrationTest(x = pred$observed, mu = pred$pred, size = exp(pred$psi))
## we can also test calibration of the fitted values
## using the calibrationTest-method for "hhh4" fits
calibrationTest(result, subset = nrow(salmonella)-(4:0))
## plot a (non-randomized) PIT histogram for the predictions
pit(pred)
## the above uses the pit-method for "oneStepAhead" predictions,
## which is a simple wrapper around the default method:
pit(x = pred$observed, pdistr = "pnbinom", mu = pred$pred, size = exp(pred$psi))
### multivariate measles count time series
## (omitting oneStepAhead forecasts here to keep runtime low)
data("measlesWeserEms")
## simple hhh4 model with random effects in the endemic component
measlesModel <- list(
end = list(f = addSeason2formula(~0 + ri(type="iid"))),
ar = list(f = ~1),
family = "NegBin1")
measlesFit <- hhh4(measlesWeserEms, control = measlesModel)
## assess overall (in-sample) calibration of the model, i.e.,
## if the observed counts are from the fitted NegBin distribution
calibrationTest(measlesFit) # default is DSS (not suitable for low counts)
calibrationTest(measlesFit, which = "logs") # p = 0.7238
## to assess calibration in the second year for a specific district
calibrationTest(measlesFit, subset = 53:104, units = "03452", which = "rps")
pit(measlesFit, subset = 53:104, units = "03452")
### For a more sophisticated multivariate analysis of
### areal time series of influenza counts - data("fluBYBW") -
### see the (computer-intensive) demo("fluBYBW") script:
demoscript <- system.file("demo", "fluBYBW.R", package = "surveillance")
#file.show(demoscript)
Run the code above in your browser using DataLab