twinstim_intensity: Plotting Intensities of Infection over Time or Space

Description

intensityplot method to plot the evolution of the total infection intensity, or its epidemic or endemic components/proportions over time or space (integrated over the other dimension) of fitted twinstim models (or simEpidataCS). The "simEpidataCS"-method is just a wrapper around intensityplot.twinstim by making the "simEpidataCS" object "twinstim"-compatible, i.e. enriching it by the required model components and environment.

The intensity.twinstim auxiliary function returns functions which calculate the endemic or epidemic intensity at a specific time point or location (integrated over the other dimension).

Usage

# S3 method for twinstim
intensityplot(x, which = "epidemic proportion",
    aggregate = c("time", "space"), types = 1:nrow(x$qmatrix),
    tiles, tiles.idcol = NULL, plot = TRUE, add = FALSE,
    tgrid = 101, rug.opts = list(),
    sgrid = 128, polygons.args = list(), points.args = list(),
    cex.fun = sqrt, ...)
# S3 method for simEpidataCS
intensityplot(x, ...)
intensity.twinstim(x,
    aggregate = c("time", "space"), types = 1:nrow(x$qmatrix), 
    tiles, tiles.idcol = NULL)

Value

If plot = FALSE or aggregate = "time", a function is returned, which takes a vector of time points (if aggregate = "time") or a matrix of coordinates (if aggregate = "space"), and returns which on this grid.

intensity.twinstim returns a list containing such functions for the endemic and epidemic intensity (but these are not vectorized).

If plot = TRUE and aggregate = "space", the

trellis.object of the spatial plot is returned.

Arguments

x: an object of class "twinstim" or "simEpidataCS", respectively.
which: "epidemic proportion" (default), "endemic proportion", or "total intensity" (partially matched), or "epidemic intensity" or "endemic intensity". Determines whether to plot the total or component-specific intensities or their proportions as a function of aggregate.
aggregate: One of "time" or "space". The former results in a plot of the evolution of which as a function of time (integrated over the observation region $\bold{W}$), whereas the latter produces a spplot of which over $\bold{W}$ (spanned by tiles). In both cases, which is evaluated on a grid of values, given by tgrid or sgrid, respectively.
types: event types to aggregate. By default, all types of events are aggregated, but one could also be interested in only one specific type or a subset of event types.
tiles: object of class "SpatialPolygons" representing the decomposition of $\bold{W}$ into different regions (as used in the corresponding stgrid of the "epidataCS"). This is only needed for aggregate = "space".
tiles.idcol: either a column index for tiles@data (if tiles is a "SpatialPolygonsDataFrame"), or NULL (default), which refers to the "ID" slot of the polygons, i.e., row.names(tiles). The ID's must correspond to the factor levels of stgrid$tile of the "epidataCS" on which x was fitted.
plot: logical indicating if a plot is desired, which defaults to TRUE. Otherwise, a function will be returned, which takes a vector of time points (if aggregate = "time") or a matrix of coordinates (if aggregate = "space"), and returns which on this grid.
add: logical. If TRUE and aggregate = "time", paths are added to the current plot, using lines. This does not work for aggregate = "space".
tgrid: either a numeric vector of time points when to evaluate which, or a scalar representing the desired number of evaluation points in the observation interval $[t_0, T]$. This argument is unused for aggregate = "space".
rug.opts: if a list, its elements are passed as arguments to the function rug, which will mark the time points of the events if aggregate = "time" (it is unused in the spatial case); otherwise (e.g., NULL), no rug will be produced. By default, the rug argument ticksize is set to 0.02 and quiet is set to TRUE. Note that the argument x of the rug function, which contains the locations for the rug is fixed internally and can not be modified.
sgrid: either an object of class "SpatialPixels" (or coercible to that class) representing the locations where to evaluate which, or a scalar representing the approximate number of points of a grid constructed on the bounding box of tiles. sgrid is internally subsetted to contain only points inside tiles. This argument is unused for aggregate = "time".
polygons.args: if a list, its elements are passed as arguments to sp.polygons, which will add tiles to the plot if aggregate = "space" (it is unused for the temporal plot). The default border colour of the polygons is "darkgrey".
points.args: if a list, its elements are passed as arguments to sp.points, which will add the event locations to the plot if aggregate = "space" (it is unused for the temporal plot). By default, the plot symbol is set to pch=1. The sizes of the points are determined as the product of the argument cex (default: 0.5) of this list and the sizes obtained from the function cex.fun which accounts for multiple events at the same location.
cex.fun: function which takes a vector of counts of events at each unique location and returns a (vector of) cex value(s) for the sizes of the points at the event locations used in points.args. Defaults to the sqrt() function, which for the default circular pch=1 means that the area of each point is proportional to the number of events at its location.
...: further arguments passed to plot or lines (if aggregate = "time"), or to spplot (if aggregate = "space").
For intensityplot.simEpidataCS, arguments passed to intensityplot.twinstim.

Author

Sebastian Meyer

Examples

Run this code

data("imdepi", "imdepifit")

# for the intensityplot we need the model environment, which can be
# easily added by the intelligent update method (no need to refit the model)
imdepifit <- update(imdepifit, model=TRUE)

## path of the total intensity
opar <- par(mfrow=c(2,1))
intensityplot(imdepifit, which="total intensity",
              aggregate="time", tgrid=500)
plot(imdepi, "time", breaks=100)
par(opar)

## time course of the epidemic proportion (default) by event
intensityplot(imdepifit, tgrid=500, types=1)
intensityplot(imdepifit, tgrid=500, types=2,
              add=TRUE, col=2, rug.opts=list(col=2))
legend("topright", legend=levels(imdepi$events$type), lty=1, col=1:2,
       title = "event type")

## endemic and total intensity in one plot
intensityplot(imdepifit, which="total intensity", tgrid=501, lwd=2,
              ylab="intensity")
intensityplot(imdepifit, which="endemic intensity", tgrid=501, lwd=2,
              add=TRUE, col=2, rug.opts=NULL)
text(2500, 0.36, labels="total", col=1, pos=2, font=2)
text(2500, 0.08, labels="endemic", col=2, pos=2, font=2)


## spatial shape of the intensity (aggregated over time)

# need a map of the 'stgrid' tiles, here Germany's districts
load(system.file("shapes", "districtsD.RData", package="surveillance"))

# total intensity (using a rather sparse 'sgrid' for speed)
intensityplot(imdepifit, which="total intensity",
              aggregate="space", tiles=districtsD, sgrid=500,
              col.regions=rev(heat.colors(100)))



if (surveillance.options("allExamples")) {
  # epidemic proportion by type
  maps_epiprop <- lapply(1:2, function (type) {
      intensityplot(imdepifit, which="epidemic", aggregate="space",
                    types=type, tiles=districtsD, sgrid=1000,
                    main=rownames(imdepifit$qmatrix)[type],
                    scales=list(draw=FALSE), at=seq(0,1,by=0.1),
                    col.regions=rev(hcl.colors(10,"YlOrRd")),
                    par.settings=list(par.title.text=list(cex=1)))
  })
  plot(maps_epiprop[[1]], split=c(1,1,2,1), more=TRUE)
  plot(maps_epiprop[[2]], split=c(2,1,2,1))
}