# Plots for Fitted `hhh4`

-models

`hhh4_plot.Rd`

There are six `type`

s of plots for fitted `hhh4`

models:

Plot the

`"fitted"`

component means (of selected units) along time along with the observed counts.Plot the estimated

`"season"`

ality of the three components.Plot the time-course of the dominant eigenvalue

`"maxEV"`

.If the units of the corresponding multivariate

`"sts"`

object represent different regions, maps of the fitted mean components averaged over time (`"maps"`

), or a map of estimated region-specific intercepts (`"ri"`

) of a selected model component can be produced.Plot the (estimated) neighbourhood weights (

`"neweights"`

) as a function of neighbourhood order (shortest-path distance between regions), i.e.,`w_ji ~ o_ji`

.

Spatio-temporal `"hhh4"`

models and these plots are illustrated in
Meyer et al. (2017, Section 5), see `vignette("hhh4_spacetime")`

.

## Usage

```
# S3 method for class 'hhh4'
plot(x, type=c("fitted", "season", "maxEV", "maps", "ri", "neweights"), ...)
plotHHH4_fitted(x, units = 1, names = NULL,
col = c("grey85", "blue", "orange"),
pch = 19, pt.cex = 0.6, pt.col = 1,
par.settings = list(),
legend = TRUE, legend.args = list(),
legend.observed = FALSE,
decompose = NULL, total = FALSE, meanHHH = NULL, ...)
plotHHH4_fitted1(x, unit = 1, main = NULL,
col = c("grey85", "blue", "orange"),
pch = 19, pt.cex = 0.6, pt.col = 1, border = col,
start = x$stsObj@start, end = NULL, xaxis = NULL,
xlim = NULL, ylim = NULL, xlab = "", ylab = "No. infected",
hide0s = FALSE, decompose = NULL, total = FALSE, meanHHH = NULL)
plotHHH4_season(..., components = NULL, intercept = FALSE,
xlim = NULL, ylim = NULL,
xlab = NULL, ylab = "", main = NULL,
par.settings = list(), matplot.args = list(),
legend = NULL, legend.args = list(),
refline.args = list(), unit = 1, period = NULL)
getMaxEV_season(x, period = frequency(x$stsObj))
plotHHH4_maxEV(...,
matplot.args = list(), refline.args = list(),
legend.args = list())
getMaxEV(x)
plotHHH4_maps(x, which = c("mean", "endemic", "epi.own", "epi.neighbours"),
prop = FALSE, main = which, zmax = NULL, col.regions = NULL,
labels = FALSE, sp.layout = NULL, ...,
map = x$stsObj@map, meanHHH = NULL)
plotHHH4_ri(x, component, exp = FALSE,
at = list(n = 10), col.regions = cm.colors(100),
colorkey = TRUE, labels = FALSE, sp.layout = NULL,
gpar.missing = list(col = "darkgrey", lty = 2, lwd = 2),
...)
plotHHH4_neweights(x, plotter = boxplot, ...,
exclude = 0, maxlag = Inf)
```

## Arguments

- x
a fitted

`hhh4`

object.- type
type of plot: either

`"fitted"`

component means of selected`units`

along time along with the observed counts, or`"season"`

ality plots of the model components and the epidemic dominant eigenvalue (which may also be plotted along overall time by`type="maxEV"`

, especially if the model contains time-varying neighbourhood weights or unit-specific epidemic effects), or`"maps"`

of the fitted mean components averaged over time, or a map of estimated region-specific random intercepts (`"ri"`

) of a specific model`component`

. The latter two require`x$stsObj`

to contain a map.- ...
For

`plotHHH4_season`

and`plotHHH4_maxEV`

, one or more`hhh4`

-fits, or a single list of these. Otherwise further arguments passed on to other functions.

For the`plot`

-method these go to the specific plot`type`

function.`plotHHH4_fitted`

passes them to`plotHHH4_fitted1`

, which is called sequentially for every unit in`units`

.`plotHHH4_maps`

and`plotHHH4_ri`

pass additional arguments to`spplot`

, and`plotHHH4_neweights`

to the`plotter`

.- units,unit
integer or character vector specifying a single

`unit`

or possibly multiple`units`

to plot. It indexes`colnames(x$stsObj)`

.

In`plotHHH4_fitted`

,`units=NULL`

plots all units.

In the seasonality plot, selection of a unit is only relevant if the model contains unit-specific intercepts or seasonality terms.- names,main
main title(s) for the selected

`unit`

(`s`

) /`components`

. If`NULL`

(default),`plotHHH4_fitted1`

will use the appropriate element of`colnames(x$stsObj)`

, whereas`plotHHH4_season`

uses default titles.- col,border
length 3 vectors specifying the fill and border colors for the endemic, autoregressive, and spatio-temporal component polygons (in this order).

- pch,pt.cex,pt.col
style specifications for the dots drawn to represent the observed counts.

`pch=NA`

can be used to disable these dots.- par.settings
list of graphical parameters for

`par`

. Sensible defaults for`mfrow`

,`mar`

and`las`

will be applied unless overridden or`!is.list(par.settings)`

.- legend
Integer vector specifying in which of the

`length(units)`

frames the legend should be drawn. If a logical vector is supplied,`which(legend)`

determines the frame selection, i.e., the default is to drawn the legend in the first (upper left) frame only, and`legend=FALSE`

results in no legend being drawn.- legend.args
list of arguments for

`legend`

, e.g., to modify the default positioning`list(x="topright", inset=0.02)`

.- legend.observed
logical indicating if the legend should contain a line for the dots corresponding to observed counts.

- decompose
if

`TRUE`

or (a permutation of)`colnames(x$stsObj)`

, the fitted mean will be decomposed into the contributions from each single unit and the endemic part instead of the default endemic + AR + neighbours decomposition.- total
logical indicating if the fitted components should be summed over all units to be compared with the total observed counts at each time point. If

`total=TRUE`

, the`units`

/`unit`

argument is ignored.- start,end
time range to plot specified by vectors of length two in the form

`c(year,number)`

, see`"sts"`

.- xaxis
if this is a list (of arguments for

`addFormattedXAxis`

), the time axis is nicely labelled similar to`stsplot_time`

. Note that in this case or if`xaxis = NA`

, the basic time indexes`1:nrow(x$stsObj)`

will be used as x coordinates, which is different from the long-standing default (`xaxis = NULL`

) with a real time scale.- xlim
numeric vector of length 2 specifying the x-axis range. The default (

`NULL`

) is to plot the complete time range (`type="fitted"`

) or period (`type="season"`

), respectively.- ylim
y-axis range. For

`type="fitted"`

, this defaults to`c(0,max(observed(x$stsObj)[,unit]))`

. For`type="season"`

,`ylim`

must be a list of length`length(components)`

specifying the range for every component plot, or a named list to customize only a subset of these. If only one`ylim`

is specified, it will be recycled for all`components`

plots.- xlab,ylab
axis labels. For

`plotHHH4_season`

,`ylab`

specifies the y-axis labels for all`components`

in a list (similar to`ylim`

). If`NULL`

or incomplete, default mathematical expressions are used. If a single name is supplied such as the default`ylab=""`

(to omit y-axis labels), it is used for all`components`

.- hide0s
logical indicating if dots for zero observed counts should be omitted. Especially useful if there are too many.

- meanHHH
(internal) use different component means than those estimated and available from

`x`

.- components
character vector of component names, i.e., a subset of

`c("ar", "ne", "end")`

, for which to plot the estimated seasonality. If`NULL`

(the default), only components which appear in any of the models in`...`

are plotted.

A seasonality plot of the epidemic dominant eigenvalue is also available by including`"maxEV"`

in`components`

, but it only supports models without epidemic covariates/offsets.- intercept
logical indicating whether to include the global intercept. For

`plotHHH4_season`

, the default (`FALSE`

) means to plot seasonality as a multiplicative effect on the respective component. Multiplication by the intercept only makes sense if there are no further (non-centered) covariates/offsets in the component.- exp
logical indicating whether to

`exp`

-transform the color-key axis labels to show the multiplicative effect of the region-specific random intercept on the respective component. Axis labels are then computed using`log_breaks`

from package scales (if that is available) or`axisTicks`

(as a fallback) respecting the`colorkey$tick.number`

setting (default: 7). The default is`FALSE`

.- at
a numeric vector of breaks for the color levels (see

`levelplot`

), or a list specifying the number of breaks`n`

(default: 10) and their`range`

(default: range of the random effects, extended to be symmetric around 0). In the latter case, breaks are equally spaced (on the original, non-`exp`

scale of the random intercepts). If`exp=TRUE`

, custom breaks (or`range`

) need to be given on the exp-scale.- matplot.args
list of line style specifications passed to

`matplot`

, e.g.,`lty`

,`lwd`

,`col`

.- refline.args
list of line style specifications (e.g.,

`lty`

or`col`

) passed to`abline`

when drawing the reference line (`h=1`

) in plots of seasonal effects (if`intercept=FALSE`

) and of the dominant eigenvalue. The reference line is omitted if`refline.args`

is not a list.- period
a numeric value giving the (longest) period of the harmonic terms in the model. This usually coincides with the

`freq`

of the data (the default), but needs to be adjusted if the model contains harmonics with a longer periodicity.- which
a character vector specifying the components of the mean for which to produce maps. By default, the overall mean and all three components are shown.

- prop
a logical indicating whether the component maps should display proportions of the total mean instead of absolute numbers.

- zmax
a numeric vector of length

`length(which)`

(recycled as necessary) specifying upper limits for the color keys of the maps, using a lower limit of 0. A missing element (`NA`

) means to use a map-specific color key only covering the range of the values in that map (can be useful for`prop = TRUE`

). The default`zmax = NULL`

means to use the same scale for the component maps and a separate scale for the map showing the overall mean.- col.regions
a vector of colors used to encode the fitted component means (see

`levelplot`

). For`plotHHH4_maps`

, the length of this color vector also determines the number of levels, using 10 heat colors by default.- colorkey
a Boolean indicating whether to draw the color key. Alternatively, a list specifying how to draw it, see

`levelplot`

.- map
an object inheriting from

`"SpatialPolygons"`

with`row.names`

covering`colnames(x)`

.- component
component for which to plot the estimated region-specific random intercepts. Must partially match one of

`colnames(ranef(x, tomatrix=TRUE))`

.- labels
determines if and how regions are labeled, see

`layout.labels`

.- sp.layout
optional list of additional layout items, see

`spplot`

.- gpar.missing
list of graphical parameters for

`sp.polygons`

, applied to regions with missing random intercepts, i.e., not included in the model. Such extra regions won't be plotted if`!is.list(gpar.missing)`

.- plotter
the (name of a) function used to produce the plot of weights (a numeric vector) as a function of neighbourhood order (a factor variable). It is called as

`plotter(Weight ~ Distance, ...)`

and defaults to`boxplot`

. A useful alternative is, e.g.,`stripplot`

from package lattice.- exclude
vector of neighbourhood orders to be excluded from plotting (passed to

`factor`

). By default, the neighbourhood weight for order 0 is not shown, which is usually zero anyway.- maxlag
maximum order of neighbourhood to be assumed when computing the

`nbOrder`

matrix. This additional step is necessary iff`neighbourhood(x$stsObj)`

only specifies a binary adjacency matrix.

## Value

`plotHHH4_fitted1`

invisibly returns a matrix of the fitted
component means for the selected `unit`

, and `plotHHH4_fitted`

returns these in a list for all `units`

.`plotHHH4_season`

invisibly returns the plotted y-values, i.e. the
multiplicative seasonality effect within each of `components`

.
Note that this will include the intercept, i.e. the point estimate of
\(exp(intercept + seasonality)\) is plotted and returned.`getMaxEV_season`

returns a list with elements
`"maxEV.season"`

(as plotted by
`plotHHH4_season(..., components="maxEV")`

,
`"maxEV.const"`

and `"Lambda.const"`

(the Lambda matrix and
its dominant eigenvalue if time effects are ignored).`plotHHH4_maxEV`

(invisibly) and `getMaxEV`

return the
dominant eigenvalue of the \(\Lambda_t\) matrix for all time points
\(t\) of `x$stsObj`

.`plotHHH4_maps`

returns a `trellis.object`

if
`length(which) == 1`

(a single `spplot`

), and
otherwise uses `grid.arrange`

from the
gridExtra package to arrange all `length(which)`

`spplot`

s on a single page.
`plotHHH4_ri`

returns the generated `spplot`

, i.e.,
a `trellis.object`

.`plotHHH4_neweights`

eventually calls `plotter`

and
thus returns whatever is returned by that function.

## References

Held, L. and Paul, M. (2012): Modeling seasonality in space-time
infectious disease surveillance data.
*Biometrical Journal*, **54**, 824-843.
doi:10.1002/bimj.201200037

Meyer, S., Held, L. and Höhle, M. (2017):
Spatio-temporal analysis of epidemic phenomena using the R package
surveillance.
*Journal of Statistical Software*, **77** (11), 1-55.
doi:10.18637/jss.v077.i11

## See also

other methods for `hhh4`

fits, e.g., `summary.hhh4`

.

## Examples

```
data("measlesWeserEms")
## fit a simple hhh4 model
measlesModel <- list(
ar = list(f = ~ 1),
end = list(f = addSeason2formula(~0 + ri(type="iid"), S=1, period=52),
offset = population(measlesWeserEms)),
family = "NegBin1"
)
measlesFit <- hhh4(measlesWeserEms, measlesModel)
## fitted values for a single unit
plot(measlesFit, units=2)
## sum fitted components over all units
plot(measlesFit, total=TRUE)
## 'xaxis' option for a nicely formatted time axis
## default tick locations and labels:
plot(measlesFit, total=TRUE, xaxis=list(epochsAsDate=TRUE, line=1))
## an alternative with monthly ticks:
oopts <- surveillance.options(stsTickFactors = c("%m"=0.75, "%Y" = 1.5))
plot(measlesFit, total=TRUE, xaxis=list(epochsAsDate=TRUE,
xaxis.tickFreq=list("%m"=atChange, "%Y"=atChange),
xaxis.labelFreq=list("%Y"=atMedian), xaxis.labelFormat="%Y"))
surveillance.options(oopts)
## plot the multiplicative effect of seasonality
plot(measlesFit, type="season")
## alternative fit with biennial pattern, plotted jointly with original fit
measlesFit2 <- update(measlesFit,
end = list(f = addSeason2formula(~0 + ri(type="iid"), S=2, period=104)))
plotHHH4_season(measlesFit, measlesFit2, components="end", period=104)
## dominant eigenvalue of the Lambda matrix (cf. Held and Paul, 2012)
getMaxEV(measlesFit) # here simply constant and equal to exp(ar.1)
plot(measlesFit, type="maxEV") # not very exciting
## fitted mean components/proportions by district, averaged over time
if (requireNamespace("gridExtra")) {
plot(measlesFit, type="maps", labels=list(cex=0.6),
which=c("endemic", "epi.own"), prop=TRUE, zmax=NA,
main=c("endemic proportion", "autoregressive proportion"))
}
## estimated random intercepts of the endemic component
round(nu0 <- fixef(measlesFit)["end.ri(iid)"], 4) # global intercept
round(ranefs <- ranef(measlesFit, tomatrix = TRUE), 4) # zero-mean deviations
stopifnot(all.equal(
nu0 + ranefs,
ranef(measlesFit, intercept = TRUE) # local intercepts (log-scale)
))
plot(measlesFit, type="ri", component="end",
main="deviations around the endemic intercept (log-scale)")
exp(ranef(measlesFit)) # multiplicative effects, plotted below
plot(measlesFit, type="ri", component="end", exp=TRUE,
main="multiplicative effects",
labels=list(font=3, labels="GEN"))
## neighbourhood weights as a function of neighbourhood order
plot(measlesFit, type="neweights") # boring, model has no "ne" component
## fitted values for the 6 regions with most cases and some customization
bigunits <- tail(names(sort(colSums(observed(measlesWeserEms)))), 6)
plot(measlesFit, units=bigunits,
names=measlesWeserEms@map@data[bigunits,"GEN"],
legend=5, legend.args=list(x="top"), xlab="Time (weekly)",
hide0s=TRUE, ylim=c(0,max(observed(measlesWeserEms)[,bigunits])),
start=c(2002,1), end=c(2002,26), par.settings=list(xaxs="i"))
```