hhh4_validation: Predictive Model Assessment for 'hhh4' Models

In jimhester/surveillance: Temporal and Spatio-Temporal Modeling and Monitoring of Epidemic Phenomena

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").
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 the penultimate time point of result$control$subset.
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 one-step-ahead predictions but much faster: "first" will refit the model for the first time point tp[1] only and use this specific fit to calculate all subsequent predictions, whereas "final" will just use result to calculate these. The latter case thus gives nothing else than a subset of result$fitted.values, if the tp's are part of the fitted subset result$control$subset.
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") is to use the parameter estimates from the previous time point, and "final" means to always use the estimates from result as initial values. Alternatively, which.start can be a list of start values as expected by hhh4, which then replace the corresponding estimates from result as initial values. This argument is ignored for “non-rolling” types.
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 verbose-1 if there is more than one time point to predict, otherwise verbose.
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 default setting type="rolling" and which.start="current" (use which.start="final" for this to work).
x	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 score ("dss"), and squared error score ("ses"). The normalized SES ("nses") is also available but it is improper and not computed by default. It is possible to name own scoring rules in which. These must be functions of (x, mu, size), vectorized in all arguments (time x unit matrices) except that size is NULL in case of a Poisson model. See the supplied scoring rules for guidance, e.g., logs or dss.
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 units with individual=FALSE.
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 scoring rule. For the pit-methods, further arguments are passed to pit.default.

Value

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

Author(s)

Michaela Paul and Sebastian Meyer

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

See Also

vignette("hhh4") and hhh4.

Examples

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



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



## plot a (non-randomized) PIT histogram for the predictions
with(pred, pit(x = observed, pdistr = "pnbinom", mu = pred, size = exp(psi)))


## Not run: 
######################################################################
# Do one-step-ahead predictions for the models from the Paul & Held 
# (2011) paper for the influenza data from Bavaria and Baden-Wuerttemberg 
# (this takes some time!)
######################################################################

## see ?hhh4 for a specification of the models

## do 1-step ahead predictions for the last two years

tp <- nrow(fluBYBW)-2*52

val_A0 <- oneStepAhead(res_A0,tp=tp) 
val_B0 <- oneStepAhead(res_B0,tp=tp) 
val_C0 <- oneStepAhead(res_C0,tp=tp) 

val_A1 <- oneStepAhead(res_A1,tp=tp) 
val_B1 <- oneStepAhead(res_B1,tp=tp) 
val_C1 <- oneStepAhead(res_C1,tp=tp) 

val_A2 <- oneStepAhead(res_A2,tp=tp) 
val_B2 <- oneStepAhead(res_B2,tp=tp) 
val_C2 <- oneStepAhead(res_C2,tp=tp) 

val_D <- oneStepAhead(res_D,tp=tp) 


##################################
## compute scores
################################

#scores
vals <- ls(pattern="val_")
nam <- substring(vals,first=5,last=6)

whichScores <- c("logs", "rps", "ses")
scores_i <- list()
meanScores <- NULL
for(i in seq(along.with=vals)){
  sc <- scores(get(vals[i]), which=whichScores, individual=TRUE)
  scores_i[[i]] <- sc
  meanScores <- rbind(meanScores,colMeans(sc, dims=2))
}

names(scores_i) <- nam
rownames(meanScores) <- nam

##comparison with best model B2 

compareWithBest <- function(best, whichModels, nPermut=9999, seed=1234){
  set.seed(seed)
  pVals <- NULL
  for(score in seq_along(whichScores)){
    p <- c()
    for(model in whichModels){
      if(model==best) p <- c(p,NA)
      else p <- c(p,permutationTest(scores_i[[model]][,,score],scores_i[[best]][,,score],
     plot=TRUE,nPermutation=nPermut, verbose=TRUE)$pVal.permut)
    }  
    pVals <- cbind(pVals,p)
  }
  return(pVals)
}

pVals_flu <- compareWithBest(best=6, whichModels=1:10,seed=2059710987)
rownames(pVals_flu) <- nam

## End(Not run)

jimhester/surveillance documentation built on May 19, 2019, 10:33 a.m.