twinstim_intensity: Plotting Intensities of Infection over Time or Space

Description

intensityplot method to plot the evolution of the total infection intensity, its epidemic proportion or its endemic proportion 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 = c("epidemic proportion", "endemic proportion", "total intensity"),
    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", "endemic proportion", or "total intensity". Partial matching is applied. Determines whether to plot the path of the total intensity or its epidemic or endemic proportions over time or space (which) aggregated over the other dimension and types.
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). By default, the fill colour of the tiles is set to "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 by event
intensityplot(imdepifit, which="epidemic proportion",
              aggregate="time", tgrid=500, types=1)
intensityplot(imdepifit, which="epidemic proportion",
              aggregate="time", tgrid=500, types=2, add=TRUE, col=2)
legend("topright", legend=levels(imdepi$events$type), lty=1, col=1:2,
       title = "event type")

## endemic and total intensity in one plot
intensity_endprop <- intensityplot(imdepifit, which="endemic proportion",
                                   aggregate="time", plot=FALSE)
intensity_total <- intensityplot(imdepifit, which="total intensity",
                                 aggregate="time", tgrid=501, lwd=2)
curve(intensity_endprop(x) * intensity_total(x), add=TRUE, col=2, lwd=2, n=501)
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)
if (surveillance.options("allExamples")) {
  ## load borders of 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)))

  # 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")),
                    colorkey=list(title=list("Epidemic proportion", cex=1)))
  })
  plot(maps_epiprop[[1]], split=c(1,1,2,1), more=TRUE)
  plot(maps_epiprop[[2]], split=c(2,1,2,1))
}