`hhh4`

Models`oneStepAhead`

computes successive one-step-ahead
predictions for a (random effects) HHH model fitted by `hhh4`

. The function `scores`

computes a number of (strictly) proper
scoring rules based on such one-step-ahead predictions
(or, for the `hhh4`

-method, directly based on the fitted values).

See Paul and Held (2011) for further details.

There are also `calibrationTest`

and `pit`

methods for `oneStepAhead`

predictions.

```
oneStepAhead(result, tp, type = c("rolling", "first", "final"),
which.start = c("current", "final"),
keep.estimates = FALSE, verbose = TRUE, cores = 1)
```scores(x, ...)
## S3 method for class 'oneStepAhead':
scores(x, which = c("logs", "rps", "dss", "ses"),
units = NULL, sign = FALSE, individual = FALSE, reverse = TRUE, ...)
## S3 method for class 'hhh4':
scores(x, which = c("logs", "rps", "dss", "ses"),
subset = x$control$subset, units = seq_len(x$nUnit), sign = FALSE, ...)

## S3 method for class 'oneStepAhead':
calibrationTest(x, ...)

## S3 method for class 'oneStepAhead':
pit(x, ...)
## S3 method for class 'hhh4':
pit(x, subset = x$control$subset, units = seq_len(x$nUnit), ...)

result

fitted

`hhh4`

model (class `"hhh4"`

).tp

numeric vector of length 1 or 2.
One-step-ahead predictions are computed from time points

`tp[1]`

, ..., `tp[2]`

(yielding predictions for time
points `tp[1]+1`

, ...), where
`tp[2]`

defaults to ttype

The default

`"rolling"`

procedure sequentially
refits the model up to each time point in `tp`

and computes
the one-step-ahead predictions for the respective next time point.
The alternative `type`

s are no truewhich.start

Which initial parameter values should be used when successively
refitting the model to subsets of the data (up to time point

`tp[1]`

, up to `tp[1]+1`

, ...) if `type="rolling"`

?
Default (`"current"`

) keep.estimates

logical indicating if parameter estimates and log-likelihoods from
the successive fits should be returned.

verbose

non-negative integer (usually in the range

`0:3`

) specifying
the amount of tracing information to output.
During `hhh4`

model updates, the following verbosity is used:
`0`

if `cores > 1`

, otherwise <cores

the number of cores to use when computing
the predictions for the set of time points

`tp`

in parallel
(with `mclapply`

).
Note that parallelization is not possible in the defaulx

an object of class

`"oneStepAhead"`

or `"hhh4"`

.which

character vector determining which scores to compute.
The package surveillance implements the following proper
scoring rules: logarithmic score (

`"logs"`

), ranked probability
score (`"rps"`

), Dawid-Sebastiani ssubset

subset of time points for which to calculate the scores
(or the PIT histogram).
Defaults to the subset used for fitting the model.

units

integer or character vector indexing the units for which
the scores (or the PIT histogram) should be computed.
By default (

`NULL`

) all units are considered.sign

logical indicating if the function should also return

`sign(x-mu)`

, i.e., the sign of the difference between
the observed counts and corresponding predictions.
This does not really make sense when averaging over multiple
`u`

individual

logical indicating if the individual scores of the

`units`

should be returned. By default (`FALSE`

), the
individual scores are averaged over all `units`

.reverse

logical indicating if the rows (time points) should be
reversed in the result (historical default for the

`oneStepAhead`

-method)....

Unused by the

`scores`

-methods.
For the `calibrationTest`

-method, further arguments are passed
to `calibrationTest.default`

, e.g., `which`

to
select a `oneStepAhead`

returns a list (of class`"oneStepAhead"`

) with the following components:pred one-step-ahead predictions in a matrix, where each row corresponds to one of the time points requested via the argument `tp`

, and which has`ncol(result$stsObj)`

unit-specific columns. The rownames indicate the predicted time points and the column names are identical to`colnames(result$stsObj)`

.observed matrix with observed counts at the predicted time points. It has the same dimensions and names as `pred`

.psi in case of a negative-binomial model, a matrix of the estimated overdispersion parameter(s) at each time point on the internal -log-scale (1 column if `"NegBin1"`

,`ncol(observed)`

columns if`"NegBinM"`

or shared overdispersion). For a`"Poisson"`

model, this component is`NULL`

.allConverged logical indicating if all successive fits converged. - If
`keep.estimates=TRUE`

, there are the following additional elements: coefficients matrix of estimated regression parameters from the successive fits. Sigma.orig matrix of estimated variance parameters from the successive fits. logliks matrix with columns `"loglikelihood"`

and`"margll"`

with their obvious meanings.- The function
`scores`

computes the scoring rules specified in the argument`which`

. If multiple`units`

are selected and`individual=TRUE`

, the result is an array of dimensions`c(nrow(pred),length(units),5+sign)`

(up tosurveillance 1.8-0, the first two dimensions were collapsed to give a matrix). Otherwise, the result is a matrix with`nrow(pred)`

rows and`5+sign`

columns. If there is only one predicted time point, the first dimension is dropped in both cases. CAVE: For historical reasons,`scores.oneStepAhead`

by default`reverse`

s the order of the rows (time points)! The`hhh4`

-method of`scores`

does not reverse rows.The

`calibrationTest`

- and`pit`

-methods are just convenient wrappers around the respective default methods. See also`calibrationTest.hhh4`

.

Paul, M. and Held, L. (2011) Predictive assessment of a non-linear
random effects model for multivariate time series of infectious
disease counts. Statistics in Medicine, **30**, 1118--1136

`vignette("hhh4")`

and `vignette("hhh4_spacetime")`

### univariate salmonella agona data data("salmonella.agona") ## convert to sts class salmonella <- disProg2sts(salmonella.agona) ## generate formula for temporal and seasonal trends f.end <- addSeason2formula(f = ~1 + t, S=1, period=52) model1 <- list(ar = list(f=~1), end = list(f=f.end), family = "NegBin1") ## fit the model result <- hhh4(salmonella, model1) ## do sequential one-step-ahead predictions for the last 5 weeks pred <- oneStepAhead(result, nrow(salmonella)-5, type="rolling", which.start="final", verbose=FALSE) ## 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))) ## oneStepAhead(..., type="final") just means fitted values stopifnot(identical( unname(oneStepAhead(result, nrow(salmonella)-5, type="final", verbose=FALSE)$pred), unname(tail(fitted(result), 5)))) ## compute scores (sc <- scores(pred, reverse = FALSE)) ## scores with respect to the fitted values are similar (scFitted <- scores(result, subset = nrow(salmonella)-c(4:0))) ## test that scFitted is equivalent to scores(oneStepAhead(..., type = "final")) stopifnot(all.equal( scFitted, scores(oneStepAhead(result, nrow(salmonella)-5, type="final", verbose=FALSE), reverse = FALSE), check.attributes = FALSE)) ## plot a (non-randomized) PIT histogram for the predictions with(pred, pit(x = observed, pdistr = "pnbinom", mu = pred, size = exp(psi))) ### 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(file.path("demo", "fluBYBW.R"), package = "surveillance") demoscript #file.show(demoscript)