plot.hhh4: Plots for Fitted `hhh4`-models

Description

There are six types 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. (2016, 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, 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, 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)
getMaxEV_season(x)
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 = hcl.colors(10),
              labels = FALSE, sp.layout = NULL, ...,
              map = x$stsObj@map, meanHHH = NULL)
plotHHH4_ri(x, component, labels = FALSE, sp.layout = NULL,
            gpar.missing = list(col = "darkgrey", lty = 2, lwd = 2),
            ...)
plotHHH4_neweights(x, plotter = boxplot, ...,
                   exclude = 0, maxlag = Inf)

Arguments

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 al

...

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

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 seasona

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 plotHH

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

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.

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

xlim

numeric vector of length 2 specifying the x-axis range. The default (NULL) is to plot the complete time range.

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 ran

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 us

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

intercept

logical indicating whether to plot seasonality as a multiplicative effect on the respective component (the default intercept=FALSE), or additionally multiplied by the corresponding intercept (intercept=TRUE). The latt

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

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. The default is to use the same scale for the component maps and a separate scale for the map showing the over

col.regions

a vector of colors used to encode the fitted component means (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(gp

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

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 matri

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) spplots 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. 10.1002/bimj.201200037

Meyer, S., Held, L. and H�{oe}hle, M. (2016): Spatio-temporal analysis of epidemic phenomena using the Rpackage surveillance. Journal of Statistical Software. In press. Preprint available at http://arxiv.org/abs/1411.0416

Examples

Run this code

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)

## 'xaxis' option for a nicely formatted time axis
## default tick locations and labels:
plot(measlesFit, units=2, xaxis=list(epochsAsDate=TRUE, line=1))
## an alternative with monthly ticks:
oopts <- surveillance.options(stsTickFactors = c("%m"=0.75, "%Y" = 1.5))
plot(measlesFit, units=2, 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")

## 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 by district averaged over time
if (requireNamespace("gridExtra"))
    plot(measlesFit, type="maps", labels=list(cex=0.6),
         main=c("Total","Endemic","Within district","From other districts"))

## random intercepts of the endemic component
plot(measlesFit, type="ri", component="end", 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"))

## plot completely decomposed mean structure (useless without 'ne' component)
plot(measlesFit, units=bigunits, col=rainbow(measlesFit$nUnit), decompose=TRUE)