hhh4_validation: Predictive Model Assessment for `hhh4` Models

Description

The function 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.

Usage

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), ...)

Arguments

result

fitted hhh4 model (class "hhh4").

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 t

type

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 types are no true

which.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 defaul

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 s

subset

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

Value

oneStepAhead returns a list (of class "oneStepAhead") with the following components:
predone-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).
observedmatrix with observed counts at the predicted time points. It has the same dimensions and names as pred.
psiin 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.
allConvergedlogical indicating if all successive fits converged.
If keep.estimates=TRUE, there are the following additional elements:
coefficientsmatrix of estimated regression parameters from the successive fits.
Sigma.origmatrix of estimated variance parameters from the successive fits.
logliksmatrix 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 to surveillance 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 reverses 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.

encoding

latin1

References

Czado, C., Gneiting, T. and Held, L. (2009) Predictive model assessment for count data. Biometrics, 65, 1254--1261

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

Examples

Run this code

### 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)

Run the code above in your browser using DataCamp Workspace