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

Description

There are four 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 and the dominant eigenvalue of the epidemic components.
Plot the complete time course of the"maxEV"(to be used if the neighbourhood component contains time-varying weights).
If the units of the corresponding multivariate"sts"object represent different regions, a map of estimated region-specific intercepts ("ri") of a selected model component can be produced.

Usage

## S3 method for class 'hhh4':
plot(x, type=c("fitted", "season", "maxEV", "ri"), ...)
plotHHH4_fitted(x, units = 1, names = NULL,
                col = c("orange", "blue", "grey85", "black"),
                pch = 19, pt.cex = 0.6,
                par.settings = list(),
                legend = TRUE, legend.args = list(),
                legend.observed = TRUE, ...)
plotHHH4_fitted1(x, unit = 1, main = NULL,
                 col = c("grey30", "grey60", "grey85", "grey0"),
                 pch = 19, pt.cex = 0.6, border = col, 
                 start = x$stsObj@start, end = NULL,
                 xlim = NULL, ylim = NULL, xlab = "", ylab = "No. infected",
                 hide0s = FALSE, meanHHH = NULL)
plotHHH4_season(...,
                components=c("ar", "ne", "end", "maxEV"),
                xlim=NULL, ylim=NULL,
                xlab=NULL, ylab=NULL, main=NULL,
                par.settings = list(), matplot.args = list(),
                legend = NULL, legend.args = list(),
                refline.args = list(), unit = 1)
getMaxEV_season(x, unit = 1)
plotHHH4_maxEV(...,
               matplot.args = list(), refline.args = list(),
               legend.args = list())
getMaxEV(x)
plotHHH4_ri(x, component, sp.layout = NULL,
            gpar.missing = list(col = "darkgrey", lty = 2, lwd = 2),
            ...)

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 the seasonality functions, selection of a unit is only relevant if the model cont

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 4 vector specifying the colors for the spatio-temporal, autoregressive, endemic, and observed elements in this order. For plotHHH4_fitted, a length 3 vector also works, in which case the default "black" will be

pch,pt.cex

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.

start,end

time range to plot specified by vectors of length two in the form c(year,number), see "sts".

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 plots in a list (similar to ylim).

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

vector of components for which to plot estimated seasonality. Additionally, a seasonality plot of the dominant eigenvalue of the epidemic components is available ("maxEV").

matplot.args

list of line style specifications passed to matplot, e.g., lty, lwd, col.

refline.args

list of line style specifications passed to abline to draw the reference line in the plot of the dominant eigenvalue, e.g., lty and col.

component

component for which to plot the estimated region-specific random intercepts. Must partially match one of colnames(ranef(x, tomatrix=TRUE)).

sp.layout

see spplot.

gpar.missing

list of graphical parameters applied to regions with missing random intercepts (i.e., not included in the model). Supplementary regions won't be plotted at all if !is.list(gpar.missing).

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_ri returns the generated spplot, which is a lattice plot of class "trellis".

References

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

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)

## seasonality plot
plot(measlesFit, type="season")

## dominant eigenvalue of the Lambda matrix (cf. Held and Paul, 2012)
getMaxEV_season(measlesFit)$maxEV.const  # here simply exp(AR intercept)
plot(measlesFit, type="maxEV")  # not very exciting here since constant

## random intercepts of the endemic component
plot(measlesFit, type="ri", component="end")

## 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"))