Plot Simulations from "hhh4" Models
hhh4_simulate_plot.RdArrays of simulated counts from simulate.hhh4 can be
visualized as final size boxplots, individual or average time series,
or fan charts (using the fanplot package).
An aggregate-method is also available.
Usage
# S3 method for class 'hhh4sims'
plot(x, ...)
# S3 method for class 'hhh4sims'
aggregate(x, units = TRUE, time = FALSE, ..., drop = FALSE)
as.hhh4simslist(x, ...)
# S3 method for class 'hhh4simslist'
plot(x, type = c("size", "time", "fan"), ...,
groups = NULL, par.settings = list())
# S3 method for class 'hhh4simslist'
aggregate(x, units = TRUE, time = FALSE, ..., drop = FALSE)
plotHHH4sims_size(x, horizontal = TRUE, trafo = NULL, observed = TRUE,
names = base::names(x), ...)
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, ...)
plotHHH4sims_fan(x, which = 1,
fan.args = list(), observed.args = list(), initial.args = list(),
means.args = NULL, key.args = NULL, xlim = NULL, ylim = NULL,
add = FALSE, xaxis = list(), ...)Arguments
- x
an object of class
"hhh4sims"(as resulting from thesimulate-method for"hhh4"models ifsimplify = TRUEwas 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 = TRUEis a convenient shortcut for one plot by unit.- par.settings
a list of graphical parameters for
par. Sensible defaults formfrow,marandlaswill 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.
- names
a character vector of names for
x.- 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
NULLdisables 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, a character vector (providing names for
x), or a list of parameters forlegend.- xlim,ylim
vectors of length 2 determining the axis limits.
- add
a logical indicating if the (mean) simulated time series or the fan chart, respectively, should be added to an existing plot.
- which
a single integer or a character string selecting the model in
xfor which to produce the fan chart. This is only relevant ifxis a"hhh4simslist"of simulations from multiple models. Defaults to the first model.- fan.args
a list of graphical parameters for the
fan, e.g., to employ a differentcolorRampPaletteasfan.col, or to enable contour lines vialn.- observed.args
if a list (of graphical parameters for
lines), the originally observed counts are added to the plot.- means.args
if a list (of graphical parameters for
lines), the point forecasts are added to the plot (by default as a white line within the fan).- key.args
if a list, a color key (in
fan's"boxfan"-style) is added to the fan chart. The list may include positioning parametersstart(the x-position) andylim(the y-range of the color key),spaceto modify the width of the boxfan, andrlabto modify the labels. The color key is disabled by default. An alternative way of labeling the quantiles is via the argumentlninfan.args, see the Examples.- xaxis
if a list of arguments for
addFormattedXAxis, that function is used to draw the time axis, otherwise a default x-axis is drawn.- units
a logical indicating aggregation over units. Can also be a factor (or something convertible to a factor using
as.factor) to aggregate groups of units.- time
a logical indicating if the counts should be summed over the whole simulation period.
- drop
a logical indicating if the unit dimension and the
"hhh4sims"(or"hhh4simslist") class should be dropped after aggregating over (groups of) units.
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
summary(aggregate(salmSims, time = TRUE)) # summary of simulated values
plot(salmSims, type = "size")
## individual and average simulated time series with a confidence interval
plot(salmSims, type = "time", main = "20-weeks-ahead simulation")
## fan chart based on the quantiles of the simulated counts at each time point
## point forecasts are represented by a white line within the fan
if (requireNamespace("fanplot")) {
plot(salmSims, type = "fan", main = "20-weeks-ahead simulation",
fan.args = list(ln = 1:9/10), means.args = list())
}
### 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)
## use a Poisson distribution instead (just for comparison)
measlesFit2 <- update(measlesFit1, family = "Poisson")
## simulate realizations from these models 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
mygroups <- factor(substr(colnames(measlesWeserEms), 4, 4))
apply(aggregate(measlesSims[[1]], time = TRUE, units = mygroups), 1, summary)
plot(measlesSims[[1]], groups = mygroups)
## 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))
## fan charts
if (requireNamespace("fanplot")) {
opar <- par(mfrow = c(2,1))
plot(measlesSims, type = "fan", which = 1, ylim = c(0, 80), main = "NegBin",
key.args = list())
plot(measlesSims, type = "fan", which = 2, ylim = c(0, 80), main = "Poisson")
par(opar)
}