hhh4_simulate_plot: Summarize Simulations from '"hhh4"' Models

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

Description

Arrays of simulated counts from simulate.hhh4 can be visualized in various levels of aggregation: final size, time series. Furthermore, proper scoring rules can be calculated based on the simulated predictive distributions. Be aware, though, that the current implementation can only compute univariate scores, i.e., it treats the predictions at the various time points as independent.

Usage

## S3 method for class 'hhh4sims'
plot(x, ...)

as.hhh4simslist(x, ...)
## S3 method for class 'hhh4simslist'
plot(x, type = c("size", "time"), ...,
     groups = NULL, par.settings = list())

plotHHH4sims_size(x, horizontal = TRUE, trafo = NULL, observed = TRUE, ...)

plotHHH4sims_time(x, average = mean, individual = length(x) == 1,
    conf.level = if (individual) 0.95 else NULL,
    matplot.args = list(), initial.args = list(), legend = length(x) > 1,
    xlim = NULL, ylim = NULL, add = FALSE, ...)

## 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).
type	a character string indicating the summary plot to produce.
...	further arguments passed to methods.
groups	an optional factor to produce stratified plots by groups of units. The special setting groups = TRUE is a convenient shortcut for one plot by unit.
par.settings	a list of graphical parameters for par. Sensible defaults for mfrow, mar and las will be applied unless overridden or !is.list(par.settings).
horizontal	a logical indicating if the boxplots of the final size distributions should be horizontal (the default).
trafo	an optional transformation function from the scales package, e.g., sqrt_trans.
observed	a logical indicating if a line and axis value for the observed size of the epidemic should be added to the plot. Alternatively, a list with graphical parameters can be specified to modify the default values.
average	scalar-valued function to apply to the simulated counts at each time point.
individual	a logical indicating if the individual simulations should be shown as well.
conf.level	a scalar in (0,1), which determines the level of the pointwise quantiles obtained from the simulated counts at each time point. A value of NULL disables the confidence interval.
matplot.args	a list of graphical parameters for matlines.
initial.args	if a list (of graphical parameters for lines), a bar for the initial number of cases is added to the plot.
legend	a logical or a list of parameters for legend.
xlim,ylim	vectors of length 2 determining the axis limits.
add	a logical indicating if the (mean) simulated time series should be added to an existing plot.
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 to compute the scores for that subset only.
drop	a logical indicating if univariate dimensions should be dropped (the default).

Author(s)

Sebastian Meyer

Examples

### univariate example
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 = 300, seed = 3, subset = 678 + seq_len(20),
                     y.start = observed(salmAllOnset)[678,])

## compare final size distribution to observed value
plot(salmSims)

## simulated time series
plot(salmSims, type = "time", main = "2-weeks-ahead simulation")


### multivariate example
data("measlesWeserEms")

## fit a hhh4 model to the first year
measlesModel <- list(
    end = list(f = addSeason2formula(~1), offset = population(measlesWeserEms)),
    ar = list(f = ~1),
    ne = list(f = ~1 + log(pop),
        weights = W_powerlaw(maxlag = 5, normalize = TRUE)),
    family = "NegBin1", subset = 2:52,
    data = list(pop = population(measlesWeserEms)))
measlesFit1 <- hhh4(measlesWeserEms, control = measlesModel)
measlesFit2 <- update(measlesFit1, family = "Poisson")

## simulate realizations from this model during the second year
measlesSims <- lapply(X = list(NegBin = measlesFit1, Poisson = measlesFit2),
                      FUN = simulate, nsim = 50, seed = 1, subset = 53:104,
                      y.start = observed(measlesWeserEms)[52,])

## final size of the first model
plot(measlesSims[[1]])

## stratified by groups of districts
plot(measlesSims[[1]], groups = factor(substr(colnames(measlesWeserEms), 4, 4)))

## a class and plot-method for a list of simulations from different models
measlesSims <- as.hhh4simslist(measlesSims)
plot(measlesSims)

## simulated time series
plot(measlesSims, type = "time", individual = TRUE, ylim = c(0, 80))

## compare proper scoring rules for a specific subset of the regions
## (CAVE: these are univariate scores for each time point and region,
##        which do not account for dependence over time)
measlesScores5 <- scores(measlesSims, which = "rps",
                         units = substr(colnames(measlesWeserEms), 4, 4) == "5")
sapply(measlesScores5, mean)

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