hhh4_simulate_scores {surveillance}R Documentation

Proper Scoring Rules for Simulations from hhh4 Models

Description

Calculate proper scoring rules based on simulated predictive distributions.

Usage

## S3 method for class 'hhh4sims'
scores(x, which = "rps", units = NULL, ..., drop = TRUE)
## S3 method for class 'hhh4simslist'
scores(x, ...)

Arguments

x

an object of class "hhh4sims" (as resulting from the simulate-method for "hhh4" models if simplify = TRUE was set), or an "hhh4simslist", i.e., a list of such simulations potentially obtained from different model fits (using the same simulation period).

which

a character vector indicating which proper scoring rules to compute. By default, only the ranked probability score ("rps") is calculated. Other options include "logs" and "dss".

units

if non-NULL, an integer or character vector indexing the columns of x for which to compute the scores.

drop

a logical indicating if univariate dimensions should be dropped (the default).

...

unused (argument of the generic).

Details

This implementation can only compute univariate scores, i.e., independently for each time point.

The logarithmic score is badly estimated if the domain is large and there are not enough samples to cover the underlying distribution in enough detail (the score becomes infinite when an observed value does not occur in the samples). An alternative is to use kernel density estimation as implemented in logs_sample in the R package scoringRules.

Author(s)

Sebastian Meyer

Examples

data("salmAllOnset")

## fit a hhh4 model to the first 13 years
salmModel <- list(end = list(f = addSeason2formula(~1 + t)),
                  ar = list(f = ~1), family = "NegBin1", subset = 2:678)
salmFit <- hhh4(salmAllOnset, salmModel)

## simulate the next 20 weeks ahead
salmSims <- simulate(salmFit, nsim = 500, seed = 3, subset = 678 + seq_len(20),
                     y.start = observed(salmAllOnset)[678,])

## calculate the RPS at each time point
scores(salmSims, which = "rps")

## produce a PIT histogram based on the empirical distribution function
## of the simulated counts as the forecast distribution at each time point
pit(x = observed(attr(salmSims, "stsObserved")),
    pdistr = apply(salmSims, 1:2, ecdf))

[Package surveillance version 1.16.1 Index]