# Spatio-Temporal Analysis of Epidemic Phenomena Using the R Package
surveillance

Sebastian Meyer
University of Zurich Leonhard Held
University of Zurich Michael Höhle
Stockholm University<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>
Sebastian Meyer, Leonhard Held, Michael Höhle Spatio-Temporal Analysis of
Epidemic Phenomena Using the R Package surveillance surveillance: Spatio-
Temporal Analysis of Epidemic Phenomena The availability of geocoded health
data and the inherent temporal structure of communicable diseases have led to
an increased interest in statistical models and software for spatio-temporal
data with epidemic features. The open source R package surveillance can handle
various levels of aggregation at which infective events have been recorded:
individual-level time-stamped geo-referenced data (case reports) in either
continuous space or discrete space, as well as counts aggregated by period and
region. For each of these data types, the surveillance package implements
tools for visualization, likelihoood inference and simulation from recently
developed statistical regression frameworks capturing endemic and epidemic
dynamics. Altogether, this paper is a guide to the spatio-temporal modeling of
epidemic phenomena, exemplified by analyses of public health surveillance data
on measles and invasive meningococcal disease. space-time endemic-epidemic
modeling, infectious disease epidemiology, self-exciting point process,
multivariate time series of counts, branching process with immigration
2014-10-31 Sebastian Meyer
Epidemiology, Biostatistics and Prevention Institute
University of Zurich
Hirschengraben 84
8001 Zurich, Switzerland
E-mail:
URL: http://www.biostat.uzh.ch/aboutus/people/smeyer.html Leonhard Held
Epidemiology, Biostatistics and Prevention Institute
University of Zurich
E-mail:
URL: http://www.biostat.uzh.ch/aboutus/people/held.html Michael Höhle
Department of Mathematics
Stockholm University
E-mail:
URL: http://www.math.su.se/~hoehle

##  1 Introduction

Newly developed statistical methodology often lacks readily available software
to support reproducibility of research and application of the methods in other
settings. For the monitoring and analysis of epidemic phenomena, the open-
source R (R Core Team, 2014) package surveillance fills this gap since 2005.
The package’s functionality can be roughly divided into two categories:
algorithms for prospective aberration detection (as recently reviewed by
Salmon _et al._ , 2014), and statistical modeling frameworks to analyze
surveillance data of epidemic phenomena (the focus of this paper). The spatio-
temporal models implemented in the surveillance package have already been
successfully applied to human influenza (Paul _et al._ , 2008; Paul and Held,
2011; Geilhufe _et al._ , 2014), meningococcal disease (Paul _et al._ , 2008;
Paul and Held, 2011; Meyer _et al._ , 2012), measles (Herzog _et al._ , 2011),
hand, foot and mouth disease (Chen _et al._ , 2014), rabies in foxes (Höhle
_et al._ , 2009), coxiellosis in cows (Schrödle _et al._ , 2012), and the
classical swine fever virus (Höhle, 2009).

During the last few years several R packages for spatio-temporal data have
become available. Basic classes and utility functions are now provided by the
dedicated package spacetime (Pebesma, 2012), which builds upon the quasi
standards sp (Bivand _et al._ , 2013) for spatial data and xts (Ryan and
Ulrich, 2014) for time-indexed data, respectively. The comprehensive and long-
standing package spatstat (Baddeley and Turner, 2005) implements its own
classes and methods specific to spatial and spatio-temporal point patterns,
and offers a wide range of analysis tools including descriptive statistics,
model inference and simulation. For a more complete picture see the CRAN task
view “Handling and Analyzing Spatio-Temporal Data” (Pebesma, 2014). The _R-epi
project_ 111https://sites.google.com/site/therepiproject/ gives an overview of
R packages related to infectious disease epidemiology.

The surveillance package itself deals with three types of spatio-temporal data
distinguished by the available resolution. First, if an entire region is
continuously monitored for infective events (“infections”), which are time-
stamped, geo-referenced, and potentially enriched with further event-specific
data, then a (marked) spatio-temporal point pattern arises. The second data
type we consider comprises the event history of a discrete set of units at
fixed locations followed over time – e.g., farms during livestock epidemics –
while registering when they become susceptible (“at risk”), infected, and
potentially removed (neither at risk nor infectious). Our third data type is
often encountered as a result of privacy protection or reporting regimes, and
is an aggregated version of the individual event data mentioned first: event
counts by region and period.

Section 2 gives an overview of three previously established modeling
frameworks each matching one of the described data types. Sections 3 to 5
illustrate the corresponding implementations222This document is based on
surveillance 1.8-1 in R version 3.1.1 (2014-07-10) using knitr (Xie, 2013b)
for dynamic report generation. surveillance is hosted at
http://surveillance.r-forge.r-project.org/. – including data handling,
visualization, inference, and simulation – by applications to infectious
disease surveillance data. Although all of our case studies originate from
public health surveillance, we stress that our models could also be applied to
other epidemic phenomena such as infectious animal or plant diseases, forest
fires, crimes, or earth quakes. We close with a brief outlook in Section 6.

##  2 Spatio-temporal endemic-epidemic modeling

Epidemic models describe the spread of a communicable disease in a population.
Often, a compartmental view of the population is taken, placing individuals
into one of the three states (S)usceptible, (I)nfectious, or (R)emoved.
Modeling the transitions between these states in a closed population using
deterministic differential equations dates back to the work of Kermack and
McKendrick (1927). This simple SIR model has since been extended in a
multitude of ways, e.g., by addressing additional states or by considering
stochastic versions. Overviews of such SIR modeling can be found in Anderson
and May (1991), Daley and Gani (1999), and Keeling and Rohani (2008).
Statistical inference – i.e., estimation of the model’s parameters from actual
observed data – is, at best, a minor issue in the description of such models.
In contrast, a number of more statistically flavored epidemic models have
emerged recently. This includes, e.g., the TSIR model (Finkenstädt and
Grenfell, 2000), two-component time series models (Held _et al._ , 2005,
2006), and point process models (Lawson and Leimich, 2000; Diggle, 2006).

This article takes a two-component view on epidemic modeling where we use data
as starting point and guidance for what a “useful” epidemic model is. For us,
epidemic data are realizations of a spatio-temporal process with epidemic
(i.e., “self-exciting” or autoregressive) behavior. The surveillance package
provides three regression-oriented modeling frameworks for analyzing such
data, as summarized in Table 1. Although these models are designed for
different data types, they are all inspired by a Poisson branching process
with immigration proposed by Held _et al._ (2005). Its main characteristic is
the additive decomposition of disease risk into “endemic” and “epidemic”
features. The endemic component shall describe the background risk of new
events by external factors (independent of the history of the epidemic), which
in the context of infectious diseases may include seasonality/climate,
population, immigration, socio-demographic variables, and vaccination coverage
(all potentially varying in time and/or space). Explicit dependence between
cases is then introduced through epidemic components driven by the observed
past. In what follows, we give a brief introduction to each of the three
spatio-temporal endemic-epidemic models for a basic understanding, and refer
to the original methodological papers for further details.

|  twinstim (Section 3) |  twinSIR (Section 4) |  hhh4 (Section 5)
---|---|---|---
Data class |  "epidataCS" |  "epidata" |  "sts"
Resolution | individual events in | individual SI[R][S] event | event counts aggregated
continuous space-time | history of fixed population | in space and time
Example | case reports of IMD, | measles outbreak among |  week$\times$district counts of measles,
Germany, 2002–2008 | children in Hagelloch, 1861 | Weser-Ems region, 2001–2002
Model | spatio-temporal | multivariate temporal | multivariate time series
point process | point process | (Poisson or NegBin)
Reference | Meyer _et al._ (2012) | Höhle (2009) | Held and Paul (2012)
Table 1: Spatio-temporal endemic-epidemic models and corresponding data
classes implemented in the package surveillance.

Infective events occur at specific points in continuous time and space, which
gives rise to time-space-mark data
$\\{(t_{i},\bm{s}_{i},\bm{m}_{i}):i=1,\dotsc,n\\}$ from a region $\bm{W}$
observed during a period $(0,T]$. Event times and locations can be regarded as
a realization of a “self-exciting” spatio-temporal point process, for which
Meyer _et al._ (2012) formulated the endemic-epidemic model “twinstim” in
terms of the conditional intensity function (CIF, aka intensity process)

$\lambda(\bm{s},t)=\rho_{[\bm{s}][t]}\,\nu_{[\bm{s}][t]}+\sum_{j\in
I(\bm{s},t)}\eta_{j}\,f(\lVert\bm{s}-\bm{s}_{j}\rVert)\,g(t-t_{j})\>.$ (1)

It represents the instantaneous event rate at location $\bm{s}$ at time point
$t$ given all past events, and is often more verbosely denoted by
$\lambda^{*}$ or by explicit conditioning on the “history” $\mathcal{H}_{t}$.
Daley and Vere-Jones (2003, Chapter 7) provide a rigorous mathematical
definition of this concept. Both the endemic offset $\rho_{[\bm{s}][t]}$ and
the log-linear predictor $\nu_{[\bm{s}][t]}$ of endemic covariates are assumed
to be piecewise constant on a suitable spatio-temporal grid (e.g., week
$\times$ district), hence the notation $[\bm{s}][t]$ for the region covering
$\bm{s}$ during the period which contains $t$. Note that the endemic component
therefore simply equals an inhomogeneous Poisson process for the event counts
by cell of that grid. Furthermore, the observation-driven epidemic component
makes the process “self-exciting” by adding infection pressure from the set of
previous events,

$I(\bm{s},t)=\big{\\{}j:t_{j}<t\>\wedge\>t-t_{j}\leq\varepsilon_{j}\>\wedge\>\lVert\bm{s}-\bm{s}_{j}\rVert\leq\delta_{j}\big{\\}}\>,$

where $\varepsilon_{j},\delta_{j}>0$ are pre-specified temporal and spatial
ranges of interaction, possibly infinite. The epidemic log-linear predictor
$\eta_{j}$ of event-specific marks $\bm{m}_{j}$ weighs the infection pressure,
which is also decreasing over space and time as quantified by parametric
interaction functions $f$ and $g$, respectively. For instance, for the spatial
spread of human infectious diseases, a power-law kernel $f(x)=(x+\sigma)^{-d}$
was found to perform well (Meyer and Held, 2014b), whereas a step function
could be used for the time-course of infectivity (Meyer and Held, 2014a). Data
analysis with this model class is illustrated in Section 3, where we use an
extended version $\lambda(\bm{s},t,k)$ of Equation 1 to account for different
event types with own transmission dynamics (model details in Meyer _et al._ ,
2012).

In the above spatio-temporal point process model, the set of possible event
locations, i.e., the susceptible “population”, consists of the continuous
observation region $\bm{W}$ and is thus infinite. However, if infections can
only occur at specific sites, such as for livestock diseases among farms
(Keeling _et al._ , 2001), an SIR modeling approach can be pursued. In a
similar regression view as above, Höhle (2009) proposed an endemic-epidemic
multivariate temporal point process “twinSIR” modeled by the CIF

$\lambda_{i}(t)\,|\,\\{i\in S(t)\\}=\lambda_{0}(t)\,\nu_{i}(t)+\sum_{j\in
I(t)}\left\\{f(\lVert\bm{s}_{i}-\bm{s}_{j}\rVert)+\bm{w}_{ij}^{\top}\bm{\alpha}^{(w)}\right\\}$
(2)

for infection of a susceptible individual $i\in S(t)$. The rate decomposes
into a Cox proportional hazards formulation with a log-linear predictor
$\nu_{i}(t)=\exp\left(\bm{z}_{i}(t)^{\top}\bm{\beta}\right)$ of covariates to
model infection from external sources, and an additive epidemic component
capturing transmission from the set $I(t)$ of currently infectious
individuals. The force of infection depends on the Euclidean distance to the
infective source through a spatial interaction function

$f(u)=\sum_{m=1}^{M}\alpha^{(f)}_{m}B_{m}(u)\>,$ (3)

which is represented by a linear combination of basis functions $B_{m}$, e.g.,
B-splines (Fahrmeir _et al._ , 2013, Section 8.1). The distance-based force of
infection is modified additively by a linear predictor of covariates
$\bm{w}_{ij}$ describing the interaction of individuals $i$ and $j$ further.
Hence, the whole epidemic component of Equation 2 can be written as a single
linear predictor $\bm{x}_{i}(t)^{\top}\bm{\alpha}$ by interchanging the
summation order to

$\sum_{m=1}^{M}\alpha^{(f)}_{m}\sum_{j\in
I(t)}B_{m}(\lVert\bm{s}_{i}-\bm{s}_{j}\rVert)+\sum_{k=1}^{K}\alpha^{(w)}_{k}\sum_{j\in
I(t)}w_{ijk}=\bm{x}_{i}(t)^{\top}\bm{\alpha}\>,$ (4)

such that $\bm{x}_{i}(t)$ comprises all epidemic terms summed over $j\in
I(t)$. Note that the parameter vector $\bm{\alpha}$ is constrained to be non-
negative to ensure that $\lambda_{i}(t)\geq 0$. Application of this model
class is illustrated in Section 4.

In public health surveillance, the available data is typically limited to
counts $Y_{it}$ of events (cases) within regions and time periods. For such
spatio-temporal count data, a series of papers (Paul _et al._ , 2008; Paul and
Held, 2011; Held and Paul, 2012; Meyer and Held, 2014b) extended the original
endemic-epidemic multivariate time-series model of Held _et al._ (2005),
resulting in the current “hhh4” model class. In its most general formulation,
it assumes that, conditional on past observations, $Y_{it}$ has a negative
binomial distribution with mean

$\mu_{it}=e_{it}\,\nu_{it}+\lambda_{it}\,Y_{i,t-1}+\phi_{it}\sum_{j\neq
i}w_{ji}\,Y_{j,t-1}$ (5)

and region-specific overdispersion parameters $\psi_{i}$ such that the
conditional variance of $Y_{it}$ is $\mu_{it}(1+\psi_{i}\mu_{it})$. This
includes models assuming a common overdispersion $\psi_{i}\equiv\psi$, as well
as the Poisson distribution as the special case $\psi_{i}\equiv 0$. The
endemic component again involves a log-linear predictor ($\nu_{it}$)
proportional to an offset (of expected counts $e_{it}$), whereas the
observation-driven epidemic component here splits up into autoregressive
(local) and neighbourhood effects. For a single time series of counts $Y_{t}$,
hhh4 can be regarded as an extension of glm.nb from package MASS (Venables and
Ripley, 2002) to account for autoregression. With multiple regions, spatio-
temporal dependence is adopted by the third component in Equation 5 with
weights $w_{ji}$ reflecting the flow of infections from region $j$ to region
$i$. These neighbourhood weights can be pre-specified from external data, but
may also be estimated parametrically. Furthermore, possibly correlated random
intercepts in the log-linear predictors $\nu_{it}$, $\lambda_{it}$, and
$\phi_{it}$ are allowed for in addition to offsets, seasonality, and covariate
effects. How to use this model class in practice is illustrated in Section 5.

##  3 twinstim: Geo-referenced events

In this section, we illustrate the application of the endemic-epidemic spatio-
temporal point process intensity model twinstim (1) to case reports of
invasive meningococcal disease (IMD) in Germany, 2002–2008, as analyzed by
Meyer _et al._ (2012). Section 3.1 introduces the basic class "epidataCS" for
such data, Section 3.2 presents the core functionality of fitting and
analyzing twinstims, and Section 3.3 shows how to simulate realizations from a
fitted model.

###  3.1 Data structure: "epidataCS"

####  3.1.1 The "epidataCS" class

The first step toward fitting a twinstim is to turn the relevant data into an
object of the dedicated class "epidataCS". The suffix “CS” indicates that the
underlying point process generating the events is indexed in continuous space.
The class has the following ingredients:

events:

the core data set of point-referenced events as a "SpatialPointsDataFrame"
(package sp). The coordinates $\bm{s}_{i}$ must be defined in a planar
coordinate reference system (CRS) to enable Euclidean geometry; see the
spTransform-methods in package rgdal (Bivand _et al._ , 2014) for how to
project latitude and longitude coordinates. The data slot of the object holds
the following obligatory columns: the event time $t_{i}$, the length of the
infectious period $\varepsilon_{i}$ (eps.t), the spatial influence radius
$\delta_{i}$ (eps.s), and the enclosing tile of the space-time grid stgrid
(see below). For multitype epidemics, an additional type column must specify
the event type $k_{i}$. Further columns may contain additional marks
$\bm{m}_{i}$ of the events, e.g., the age of infected individuals, which may
be used as covariates in the epidemic predictor $\eta_{j}$ influencing
infectiousness.

W:

a "SpatialPolygons" object representing the observation region $\bm{W}$
containing the events. It must use the same CRS as events and its area must be
consistent with the sum of the tiles’ areas of the stgrid described next.

stgrid:

a data frame of covariates on a full space-time grid for the endemic component
of a twinstim (1), e.g., socio-demographic or ecological information by week
and district. Time intervals are given by a start (and stop) column, regions
are indexed by a factor variable tile, and a column area gives the area of the
respective region in a unit compatible with the CRS used for events and W.
Note that a geographic representation of the regions in stgrid is not required
for model estimation, and is thus not part of the "epidataCS" class. It is,
however, useful to have a polygon map of stgrid since it is necessary for
simulation from the estimated conditional intensity model.

qmatrix (for multitype epidemics):

a square indicator matrix for possible transmission between event types. The
default identity matrix corresponds to an independent spread of all types as,
e.g., assumed for the two meningococcal finetypes in the example below.

The converter function as.epidataCS checks the consistency of the supplied
data and pre-computes auxiliary variables for model fitting. This includes
computing the centered individual influence regions
$\bm{R}_{i}=(\bm{W}\cap\bar{b}(\bm{s}_{i},\delta_{i}))-\bm{s}_{i}$, obtained
by intersecting the observation region $\bm{W}$ with a polygonal approximation
of a disc of radius $\delta_{i}$ centered at $\bm{s}_{i}$. These polygons
$\bm{R}_{i}$ are computed by the function intersectPolyCircle, which wraps
functionality of package polyclip (Johnson, 2014). Since the $\bm{R}_{i}$ are
involved as integration domains in the point process likelihood, the more
detailed the polygons of W represent the observation region the longer it will
take to fit a twinstim. It is thus advisable to sacrifice some shape details
for speed by reducing the polygon complexity, e.g., by applying one of the
simplification methods available at MapShaper.org (Harrower and Bloch, 2006).
Within R, there is the function thinnedSpatialPoly in package maptools (Bivand
and Lewin-Koh, 2014), which implements the Douglas and Peucker (1973)
reduction method, and spatstat’s simplify.owin procedure.

Generation of an "epidataCS" object is illustrated below for the IMD data,
which is already contained as data("imdepi") in our package. The main input
component is the "SpatialPointsDataFrame" of individual point-referenced cases
of IMD:

R> summary(events)

Object of class SpatialPointsDataFrameCoordinates: min maxx40394665y27103525Is
projected: TRUEproj4string :[+init=epsg:3035 +units=km +proj=laea +lat_0=52
+lon_0=10 +x_0=4321000 +y_0=3210000+ellps=GRS80 +no_defs]Number of points:
636Data attributes: time tile type eps.t eps.s sex agegrpMin. : 0 05354 : 34
B:336 Min. :30 Min. :200 female:292 [0,3) :1941st Qu.: 539 05370 : 27 C:300
1st Qu.:30 1st Qu.:200 male :339 [3,19) :279Median :1155 11000 : 27 Median :30
Median :200 NA’s : 5 [19,Inf):162Mean :1193 05358 : 13 Mean :30 Mean :200 NA’s
: 13rd Qu.:1808 05162 : 12 3rd Qu.:30 3rd Qu.:200Max. :2543 05382 : 12 Max.
:30 Max. :200 (Other):511 Event coordinates are available as residence
postcode centroids projected in the _European Terrestrial Reference System
1989_ in kilometer units. Time is measured in days and the tile column
contains the official key of the enclosing administrative district (matching
stgrid$tile$). There are two types of events labeled as "B" and "C", which
refer to the serogroups of the two meningococcal finetypes _B:P1.7-2,4:F1-5_
and _C:P1.5,2:F3-3_ contained in the data. The infectious period is assumed to
last up to eps.t = 30 days and spatial interaction is limited to a eps.s = 200
km radius for all events. The latter has numerical advantages for a Gaussian
interaction function $f$ with a relatively small standard deviation. For a
power-law kernel, however, this restriction will be dropped to enable
occasional long-range transmission. The last two data attributes displayed in
the above event summary are covariates from the case reports: the gender and
age group of the patient.

As observation window W, we use a simplified polygon representation of
Germany’s boundary:

R> load(system.file("shapes", "districtsD.RData", package = "surveillance"))
This file contains both the "SpatialPolygonsDataFrame" districtsD describing
Germany’s 413 administrative districts as at January 1, 2009, as well as their
union stateD obtained by the call rgeos::gUnaryUnion(districtsD) (Bivand and
Rundel, 2014).

The stgrid of covariates for the endemic model component is a simplified
version of what has been used by Meyer _et al._ (2012):

start stop tile area popdensity1 0 31 01001 56.4 1557.12 0 31 01002 118.7
1996.63 0 31 01003 214.2 987.64 ... ... ... ... ...34690 2526 2557 16075
1148.5 79.234691 2526 2557 16076 843.5 133.634692 2526 2557 16077 569.1 181.5
To reduce computational cost and package size, we did not include any time-
dependent external covariates such as (lagged) weekly influenza counts by
district, and we use a temporal resolution of monthly (instead of weekly)
intervals to model endemic time trends and seasonality. Thus, the district-
specific population density (popdensity) is the only grid-based covariate and
will be included as the offset $\rho_{[\bm{s}][t]}$ in the endemic component
of Equation 1.

Finally, assuming a qmatrix without transmission across types, we construct an
object of class "epidataCS" as follows:

R> imdepi <\- as.epidataCS(events = events, W = stateD, stgrid = stgrid,\+
qmatrix = diag(2), nCircle2Poly = 16, verbose = FALSE)

####  3.1.2 Methods for data handling and visualization

Display | Subset | Extract | Modify | Convert
---|---|---|---|---
print |  [ |  nobs |  update |  as.epidata
summary |  head |  marks |  untie |  epidataCS2sts
plot |  tail |  |  |
animate |  subset |  |  |
as.stepfun |  |  |  |
Table 2: Generic and non-generic functions applicable to "epidataCS" objects.

Textual summaries of the data can be obtained by the print and summary
methods. Printing "epidataCS" shows the first n (defaults to 6) events:

R> imdepi

Observation period: 0 - 2557Observation window (bounding box): [4031, 4672] x
[2684, 3550]Spatio-temporal grid (not shown): 84 time blocks x 413 tilesTypes
of events: "B" "C"Overall number of events: 636 coordinates time tile type
eps.t eps.s sex agegrp BLOCK start popdensity1 (4110, 3200) 0.212 05554 B 30
200 male [3,19) 1 0 2612 (4120, 3080) 0.712 05382 C 30 200 male [3,19) 1 0
5193 (4410, 2920) 5.591 09574 B 30 200 female [19,Inf) 1 0 2094 (4200, 2880)
7.117 08212 B 30 200 female [3,19) 1 0 16665 (4130, 3220) 22.060 05554 C 30
200 male [3,19) 1 0 2616 (4090, 3180) 24.954 05170 C 30 200 male [3,19) 1 0
455[....] Note that during conversion to "epidataCS", the columns BLOCK (time
interval index), start and popdensity have been copied from stgrid to the
events data frame for the respective times and locations. The event marks
(including time and location) can be extracted by marks(imdepi), and these are
summarized when printing the summary(imdepi).

Visualizations are provided by the as.stepfun, plot and animate methods. The
step function shows the time course of the number of infectives, which is
driven by the event times and lengths of the infectious periods. For the IMD
data we obtain Figure 1 by:

R> plot(as.stepfun(imdepi), xlim = summary(imdepi)$timeRange, xaxs = "i",\+
xlab = "Time [days]", ylab = "Current number of infectives", main = "")

Figure 1: Time course of the number of infectives assuming infectious periods
of 30 days.

The plot-method offers aggregation over time or space, e.g. (Figure 2):

R> plot(imdepi, "time", breaks = c(0, unique(imdepi$stgrid$stop)),\+ col =
c("indianred", "darkblue"), ylim = c(0, 20))R> plot(imdepi, "space", lwd =
2,\+ points.args = list(pch = c(1, 19), col = c("indianred", "darkblue")))

(a)

(b)

Figure 2: Occurrence of the two finetypes in the temporal and spatial
dimensions.

The time series plot shows the monthly aggregated number of cases by finetype
in a stacked histogram as well as each type’s cumulative number over time. The
spatial plot shows the observation window W with the locations of all cases
(by type), where the areas of the points are proportional to the number of
cases at the respective location.

The above static plots do not capture the space-time dynamics of epidemic
spread. An animation may provide additional insight and can be produced by the
corresponding animate-method. For instance, to look at the first year of the
B-type in a weekly sequence of snapshots in a web browser (using facilities of
the animation package (Xie, 2013a)):

R> library("animation")R> saveHTML(animate(subset(imdepi, type=="B"),
interval=c(0,365), time.spacing=7),\+ htmlfile = "imdepi_animation.html",\+
nmax = Inf, interval = 0.2, loop = FALSE,\+ title = "Animation of the first
year of type B events") Selecting events from "epidataCS" is enabled by the
[-method (used by the subset-method), which returns a new "epidataCS" object
containing only the selected subset of events.

A limited data sampling resolution may lead to tied event times or locations,
which are in conflict with a continuous spatio-temporal point process model.
For instance, a temporal residual analysis would suggest model deficiencies
(Meyer _et al._ , 2012, Figure 4), and a power-law kernel for spatial
interaction may diverge if there are events with zero distance to potential
source events (Meyer and Held, 2014b). The function untie breaks ties by
random shifts. This has already been applied to the event _times_ in imdepi by
subtracting a U(0,1)-distributed random number from the original dates. To
break tied locations, we shift all locations by a random vector with length up
to half the observed minimum spatial separation.

R> eventDists <\- dist(coordinates(imdepi$events))R> (minsep <\-
min(eventDists[eventDists > 0]))

[1] 1.17

R> set.seed(321)R> imdepi_untied <\- untie(imdepi, amount = list(s =
minsep/2)) Such random tie-breaking requires sensitivity analyses as discussed
by Meyer and Held (2014b).

The update-method is useful to change the values of the maximum interaction
ranges eps.t and eps.s, since it takes care of the necessary updates of the
hidden auxiliary variables in an "epidataCS" object. For an unbounded
interaction radius:

R> imdepi_untied_infeps <\- update(imdepi_untied, eps.s = Inf)

Last but not least, "epidataCS" can be converted to the other classes
"epidata" (Section 4) and "sts" (Section 5) by aggregation. The method
as.epidata.epidataCS aggregates events by region (tile), and the function
epidataCS2sts yields counts by region and time interval. The data could then,
e.g., be analyzed by the multivariate time series model presented in Section
5. We can also use visualization tools of the "sts" class, e.g., to produce
Figure 3:

R> imdepi_sts <\- epidataCS2sts(imdepi, freq=12, start=c(2002,1),
tiles=districtsD)R> plot(imdepi_sts, type = observed~time, legend.opts =
NULL)R> plot(imdepi_sts, type = observed~unit, population =
districtsD$POPULATION/100000)

(a)

(b)

Figure 3: IMD cases (joint types) aggregated as an "sts" object by month and
district.

###  3.2 Modeling and inference

Having prepared the data as an object of class "epidataCS", the function
twinstim can be used to perform likelihood inference for conditional intensity
models matching Equation 1. The main arguments for twinstim are the formulae
of the endemic and epidemic predictors,
$\log(\rho_{[\bm{s}][t]})+\log(\nu_{[\bm{s}][t]})$ and $\log(\eta_{j})$, and
the spatial and temporal interaction functions siaf ($f$) and tiaf ($g$),
respectively. The endemic formula incorporates covariates from stgrid, and
usually contains at least the population density offset, i.e., endemic =
offset(log(popdensity)). There can be additional effects of time, which are
functions of the variable start from stgrid, or effects of, e.g., socio-
demographic and ecological variables. The second, epidemic formula specifies
$\log(\eta_{j})$ for event-specific infectivity. For instance, epidemic =
log(popdensity) + type corresponds to
$\eta_{j}=\rho_{[\bm{s}_{j}][t_{j}]}^{\gamma_{\rho}}\exp(\gamma_{0}+\gamma_{C}\mathbbm{1}(k_{j}=C))$,
which models different infectivity of the event types, and scales with
population density (a grid-based covariate) to reflect higher contact rates
and thus infectivity in more densly populated regions. For the interaction
functions, several alternatives as listed in Table 3 are predefined. They are
applicable out-of-the-box and illustrated later as part of the case study.

Spatial (siaf.*) |  Temporal (tiaf.*)
---|---
constant |  constant
gaussian |  exponential
powerlaw |  step
powerlawL |
step |
student |
Table 3: Predefined spatial and temporal interaction functions.

The log-likelihood to evaluate is a function of all parameters in the
predictors $\nu_{[\bm{s}][t]}$ and $\eta_{j}$ and in the interaction functions
$f$ and $g$. Since it has the form

$\left[\sum_{i=1}^{n}\log\lambda(\bm{s}_{i},t_{i})\right]-\int_{0}^{T}\int_{\bm{W}}\lambda(\bm{s},t)\,\mathrm{d}\bm{s}\,\mathrm{d}t\>,$
(6)

it contains a clear computational bottleneck: the two-dimensional integrals
$\int_{\bm{R}_{i}}f(\lVert\bm{s}\rVert)\,\mathrm{d}\bm{s}$ over the polygonal
influence regions $\bm{R}_{i}$. Calculation is trivial for (piecewise)
constant $f$, but otherwise requires numerical integration. For this purpose,
the R package polyCub (Meyer, 2014) offers cubature methods over polygonal
domains as described in Meyer and Held (2014c, Section 2). One of these
methods also takes analytical advantage of the assumed isotropy in such a way
that numerical integration remains in only one dimension.

####  3.2.1 Basic example

To illustrate statistical inference with twinstim, we will estimate several
models for the simplified and “untied” IMD data presented in Section 3.1. In
the endemic component, we include the population density as a multiplicative
offset $\rho_{[\bm{s}][t]}$, a (centered) time trend, and a sinusoidal wave of
frequency $2\pi/365$ to capture seasonality (cf. Held and Paul, 2012, Section
2.2):

R> (endemic <\- addSeason2formula(~offset(log(popdensity)) + I(start/365 -
3.5),\+ period = 365, timevar = "start"))

~offset(log(popdensity)) + I(start/365 - 3.5) + sin(2 * pi * start/365) +
cos(2 * pi * start/365) Because of the aforementioned integrations in the log-
likelihood (6), it is advisable to first fit an endemic-only model to obtain
reasonable start values for more complex epidemic models:

R> imdepifit_endemic <\- twinstim(endemic = endemic, epidemic = ~0,\+ data =
imdepi_untied, subset = !is.na(agegrp)) Note that this twinstim _without_ an
epidemic component can actually be represented by a Poisson regression model
for aggregated counts, which provides a nice link to the count data model hhh4
(5) illustrated in Section 5. To see this, recall that the endemic component
of a twinstim (1) is piecewise constant on the spatio-temporal grid with cells
$([\bm{s}],[t])$. Hence the log-likelihood (6) of an endemic-only twinstim
simplifies to a sum over all these cells,

$\sum_{[\bm{s}],[t]}\left\\{Y_{[\bm{s}][t]}\log\lambda([\bm{s}],[t])-\lvert[\bm{s}]\rvert\,\lvert[t]\rvert\,\lambda([\bm{s}],[t])\right\\}\>,$

where $Y_{[\bm{s}][t]}$ is the aggregated number of events observed in cell
$([\bm{s}],[t])$, and $\lvert[\bm{s}]\rvert$ and $\lvert[t]\rvert$ denote area
and length of a cell, respectively. Except for an additive constant, the above
log-likelihood is equivalently obtained from the Poisson model
$Y_{[\bm{s}][t]}\sim\operatorname{Po}(\lvert[\bm{s}]\rvert\,\lvert[t]\rvert\,\lambda([\bm{s}],[t]))$.
This relation offers a means of code validation by passing via the established
glm function to fit an endemic-only twinstim model to "epidataCS". This is
provided by the wrapper function glm_epidataCS:

R> imdepifit_glm <\- glm_epidataCS(formula = endemic,\+ data =
subset(imdepi_untied, !is.na(agegrp)))

Display | Extract | Modify | Other
---|---|---|---
print |  nobs |  update |  simulate
summary |  vcov |  add1 |
xtable |  logLik |  drop1 |
plot |  extractAIC |  stepComponent |
intensityplot |  profile |  |
iafplot |  residuals |  |
checkResidualProcess |  terms |  |
|  R0 |  |
Table 4: Generic and non-generic functions applicable to "twinstim" objects.
Note that there is no need for specific coef or confint methods, since the
respective default methods from package stats apply outright.

Many of the standard functions to access model fits in R are also implemented
for "twinstim" fits (Table 4 gives an overview). For instance we can produce a
summary as usual:

R> summary(imdepifit_endemic)

Call:twinstim(endemic = endemic, epidemic = ~0, data = imdepi_untied, subset =
!is.na(agegrp))Coefficients of the endemic component: Estimate Std. Error z
value Pr(>|z|)h.(Intercept) -20.3683 0.0419 -486.24 < 2e-16 ***h.I(start/365 -
3.5) -0.0444 0.0200 -2.22 0.027 *h.sin(2 * pi * start/365) 0.2733 0.0576 4.75
2.0e-06 ***h.cos(2 * pi * start/365) 0.3509 0.0581 6.04 1.5e-09 ***\---Signif.
codes: 0 ’***’ 0.001 ’**’ 0.01 ’*’ 0.05 ’.’ 0.1 ’ ’ 1No epidemic
component.AIC: 19166Log-likelihood: -9579Number of log-likelihood evaluations:
5Number of score function evaluations: 5Runtime: 0.167 seconds

Because of the aforementioned equivalence of the endemic component with a
Poisson regression model, the coefficients can be interpreted as log rate
ratios in the usual way. For instance, the endemic rate is estimated to
decrease by 1-exp(coef(imdepifit_endemic)[2]) $=$ 4.3% per year. Coefficient
correlations can be retrieved by the argument correlation=TRUE in the summary
call (just like for summary.glm), but may also be extracted via the standard
cov2cor(vcov(imdepifit_endemic)).

We now update the endemic model to take additional spatio-temporal dependence
between events into account. Infectivity shall depend on the meningococcal
finetype and the age group of the patient, and is assumed to be constant over
time (default), $g(t)=\mathbbm{1}_{(0,30]}(t)$, with a Gaussian distance-decay
$f(x)=\exp\left\\{-x^{2}/(2\sigma^{2})\right\\}$. This model was originally
selected by Meyer _et al._ (2012) and can be fitted as follows:

R> imdepifit_Gaussian <\- update(imdepifit_endemic,\+ epidemic = ~type +
agegrp, siaf = siaf.gaussian(),\+ control.siaf = list(F = list(adapt = 0.25),
Deriv = list(nGQ = 13)),\+ start = list(epidemic = c("(Intercept)"=-12.5,
"typeC"=-1, "siaf.1"=2.8)),\+ model = TRUE, cores = 2 * (.Platform$OS.type ==
"unix")) To reduce the runtime of this example, we specified suitable start
values for some parameters (others start at 0) and set control.siaf with a
rather low number of nodes for the cubature of $f(\lVert\bm{s}\rVert)$ (and
$\frac{\partial f(\lVert\bm{s}\rVert)}{\partial\log\sigma}$ in the score
function) over all domains $\bm{R}_{i}$ of Equation 6. On Unix-alikes, these
numerical integrations can be performed in parallel via the “multicore”
functions (mclapply et al.) from the base package parallel, here using cores=2
CPU’s. For later generation of an intensityplot, the model environment is
retained with the result.

| RR | 95% CI | p-value
---|---|---|---
h.I(start/365 - 3.5) | 0.956 | 0.91–1.00 | 0.04
h.sin(2 * pi * start/365) | 1.243 | 1.09–1.41 | 0.0008
h.cos(2 * pi * start/365) | 1.375 | 1.21–1.56 |  $<$0.0001
e.typeC | 0.401 | 0.24–0.68 | 0.0006
e.agegrp[3,19) | 1.943 | 1.03–3.66 | 0.04
e.agegrp[19,Inf) | 0.748 | 0.30–1.84 | 0.53
Table 5: Estimated rate ratios (RR) and associated Wald confidence intervals
(CI) for endemic (h.) and epidemic (e.) terms. This table was generated by
xtable(imdepifit_Gaussian).

Table 5 shows a LaTeX version of the model summary as generated by its xtable
method (Dahl, 2014), which provides rate ratios for the endemic and epidemic
effects. There is an alternative toLatex method that simply translates the
summary table of coefficients to LaTeX without transformation. On the subject
matter level, we can conclude from Table 5 that the meningococcal finetype of
serogroup C is less than half as infectious as the B-type, and that patients
in the age group 3 to 18 years cause approximately twice as many secondary
infections as infants aged 0 to 2 years.

####  3.2.2 The reproduction number

Similar to the simple SIR model (see, e.g., Keeling and Rohani, 2008, Section
2.1), the point process model (1) features a reproduction number derived from
its branching process interpretation. As soon as an individual becomes
infected, it spreads the disease according to an inhomogeneous Poisson process
with rate $\eta_{j}\cdot g(t-t_{j})\cdot f(\lVert\bm{s}-\bm{s_{j}}\rVert)$
around its origin. Since this process is independent of the individuals
parentage and of other infections,

$\mu_{j}=\eta_{j}\cdot\left[\int_{0}^{\min(T-t_{j},\varepsilon_{j})}g(t)\,dt\right]\cdot\left[\int_{\bm{R}_{j}}f(\bm{s})\,d\bm{s}\right]$
(7)

is the expected number of secondary infections triggered by the infective $j$
within its observed range of interaction. Equation 7 can also be motivated by
looking at a spatio-temporal version of the simple SIR model wrapped into the
"twinstim" class (1): no endemic component, homogeneous force of infection
($\eta_{j}\equiv\beta$), homogeneous mixing in space ($f(\bm{s})\equiv 1$,
$\delta_{j}\equiv\infty$), and exponential decay of infectivity
($g(t)=e^{-\alpha t}$, $\varepsilon_{j}\equiv\infty$). Then, for
$T\rightarrow\infty$,

$\mu=\beta\cdot\left[\int_{0}^{\infty}e^{-\alpha
t}\,dt\right]\cdot\left[\int_{\bm{W}-\bm{s_{j}}}1\,d\bm{s}\right]=\beta\cdot\lvert\bm{W}\rvert/\alpha\>,$

which corresponds to the basic reproduction number known from the simple SIR
model by interpreting $\lvert\bm{W}\rvert$ as the population size, $\beta$ as
the transmission rate and $\alpha$ as the removal rate. Like in classic
epidemic models, the process is sub-critical if $\mu<1$ holds, which means
that its eventual extinction is almost sure. However, it is crucial to
understand that in a full model with an endemic component, new infections may
always occur via “immigration”. Hence, reproduction numbers in "twinstim" are
adjusted for infections occurring independently of previous infections.
Furthermore, under-reporting and implemented control measures imply that the
estimates are to be thought of as effective reproduction numbers.

From fitted "twinstim" objects, we can extract the estimated event-specific
reproduction numbers via the R0 method. For instance in our data example,
_type_ -specific summaries of the expected numbers of secondary infections are
calculated by:

R> R0_events <\- R0(imdepifit_Gaussian)R> tapply(R0_events,
marks(imdepi_untied)[names(R0_events),"type"], summary)

$B Min. 1st Qu. Median Mean 3rd Qu. Max. 0.050 0.138 0.184 0.218 0.320 0.358$C
Min. 1st Qu. Median Mean 3rd Qu. Max. 0.0304 0.0552 0.0737 0.0962 0.1430
0.1430 Confidence intervals for the estimated reproduction numbers
$\hat{\mu}_{j}$ can be obtained by simulation, i.e., evaluating Equation 7 for
samples from the asymptotic multivariate normal distribution of the parameter
estimates (Mandel, 2013). For this purpose, the R0-method takes an argument
newcoef, which is exemplified in help("R0").

####  3.2.3 Interaction functions

Figure 4: Estimates of spatial interaction $f(x)$ scaled by $e^{\gamma_{0}}$.
The standard deviation of the bivariate Gaussian kernel was estimated as
$\hat{\sigma}=$ 16.30 (95% CI: 13.89–19.13), and the estimates for the power
law are $\hat{\sigma}=$ 4.48 (95% CI: 1.76–11.41) and $\hat{d}=$ 2.45 (95% CI:
1.79–3.36).

Figure 4 shows several estimated spatial interaction functions, which can be
plotted by, e.g., plot(imdepifit_Gaussian, which="siaf"). Meyer and Held
(2014b) found that a power-law decay of spatial interaction is more
appropriate than a Gaussian kernel to describe the spread of infectious
diseases among humans. To use the power-law kernel $f(x)=(x+\sigma)^{-d}$ the
Gaussian model can be updated by

R> imdepifit_powerlaw <\- update(imdepifit_Gaussian, data =
imdepi_untied_infeps,\+ siaf = siaf.powerlaw(), control.siaf = NULL,\+ start =
list(epidemic = c("(Intercept)"=-6.2, "siaf.1"=1.5, "siaf.2"=0.9))) where we
no longer apply the upper bound eps.s=200 km for spatial interaction. The
power-law kernel puts more weight on local transmission and exhibits a heavier
tail enabling occasional interaction between events at far distances. Table 3
also lists the step function kernel as an alternative, which is particularly
useful for two reasons: first, it is a more flexible approach since it
estimates interaction between the given knots without assuming a functional
form, and second, the spatial integrals with this kernel are fast to compute
analytically. Therefore, siaf.step can be used for a quick inspection of
spatial interaction. We update the model to use four steps at log-equidistant
knots up to an interaction range of 100 km:

R> imdepifit_step4 <\- update(imdepifit_Gaussian, data =
imdepi_untied_infeps,\+ siaf = siaf.step(exp((1:4)*log(100)/5), maxRange=100),
control.siaf = NULL,\+ start = c("e.(Intercept)" = -10, setNames(-2:-5,
paste0("e.siaf.", 1:4)))) Figure 4 suggests that the estimated step function
is in line with the power law.

For the temporal interaction function $g(t)$, model updates and plots are
equally possible in principle, e.g., update(imdepifit_Gaussian,
tiaf=tiaf.exponential()). However, the events in the IMD data are too rare to
infer the time-course of infectivity with confidence. More sophisticated
B-spline kernels of higher order have been investigated by Meyer and Held
(2014a) for spatial as well as temporal interaction, but they require far more
computational cost and are currently not part of the package.

####  3.2.4 Model selection

Model performance can be easily compared by Akaike’s Information Criterion
(AIC):

R> AIC(imdepifit_endemic, imdepifit_Gaussian, imdepifit_powerlaw,
imdepifit_step4)

df AICimdepifit_endemic 4 19166imdepifit_Gaussian 9 18967imdepifit_powerlaw 10
18940imdepifit_step4 12 18933 This suggests superiority of the power-law vs.
the Gaussian model and the endemic-only model. The more flexible step function
yields the best AIC value but its shape strongly depends on the chosen knots
and is not guaranteed to be monotonically decreasing.

The convenience function stepComponent – a wrapper around the step function in
the base package stats – can be used to perform AIC-based stepwise model
selection, e.g.:

R> imdepifit_endemic_sel <\- stepComponent(imdepifit_endemic, component =
"endemic") In this example, none of the endemic predictors is removed from the
model.

####  3.2.5 Diagnostic plots

There are two other types of plots implemented for "twinstim" objects:
intensityplot for fitted intensity curves, and checkResidualProcess to
investigate the goodness of fit. Figure 5 shows the temporal intensityplot
$\hat{\lambda}_{g}(t)=\sum_{k=1}^{2}\int_{\bm{W}}\hat{\lambda}(\bm{s},t,k)\,\mathrm{d}\bm{s}$,
i.e., the fitted ground intensity aggregated over space and both types, which
can be produced by the following code:

R> intensityplot(imdepifit_powerlaw, aggregate = "time", types = 1:2, tgrid =
501) Note that this represents a realization of a stochastic process, since it
depends on the occurred cases. The estimated endemic intensity was also added
to the plot and exhibits strong seasonality and a slow negative trend. The
proportion of the endemic intensity is rather constant along time since there
have not been any major outbreaks. This proportion can be visualized by
specifying which="endemic proportion" in the above call. Setting
aggregate="space" produces spatial intensityplots as in Figure 6, which shows
the accummulated epidemic proportion by type.

Figure 5: Fitted intensity process aggregated over space and both types.

(a)

(b)

Figure 6: Epidemic proportion of the fitted intensity process accumulated over
time by type.

The function checkResidualProcess transforms the temporal “residual process”
in such a way that it exhibits a uniform distribution and lacks serial
correlation if the fitted model describes the true CIF well (Ogata, 1988,
Section 3.3). These properties can be checked graphically as in Figure 7
produced by:

R> checkResidualProcess(imdepifit_powerlaw)

Figure 7: Inspecting the transformed temporal residual process $u_{i}$. The
left plot shows no evidence for deviation from uniformity and the right plot
suggests absence of serial correlation.

The left plot shows the empirical cumulative distribution function of the
transformed residuals with a 95% confidence band obtained by inverting the
corresponding Kolmogorov-Smirnov test. The right-hand plot suggests absence of
serial correlation of the transformed residuals.

###  3.3 Simulation

Simulation from the fitted model is another useful tool to investigate the
goodness of fit. To identify regions with unexpected IMD dynamics, Meyer _et
al._ (2012) compared the observed numbers of cases by district to the
respective 2.5% and 97.5% quantiles of 100 simulations from the selected
model. Simulations also allow us to investigate the stochastic volatility of
the endemic-epidemic process, to obtain probabilistic forecasts, and to
perform parametric bootstrap of the spatio-temporal point pattern.

The simulation algorithm we apply is described in Meyer _et al._ (2012,
Section 4). It requires a geographic representation of the tiles of stgrid,
and functionality for sampling locations according to the selected spatial
kernel $f_{2}(\bm{s})=f(\lVert\bm{s}\rVert)$, which is implemented for the
ones listed in Table 3. The following code simulates a marked spatio-temporal
point pattern for the first 90 days of our estimated "twinstim":

R> imdepisim <\- simulate(imdepifit_powerlaw, nsim = 1, seed = 1, t0 = 0, T =
90,\+ data = imdepi_untied_infeps, tiles = districtsD) The simulated object
belongs to the "simEpidataCS" class (extending "epidataCS"), which carries
some additional components related to the generating model. This enables an
R0-method as well as intensityplots for simulated data. All methods for
"epidataCS" are of course applicable, e.g., for visualization. Event marks are
per default sampled from their respective empirical distribution in the
original data, but a customized generator can be supplied as argument rmarks.
Because we have started simulation at time t0=0, no events from data have been
used as the prehistory, i.e., the first simulated event is necessarily driven
by the endemic model component. The source of each simulated event can
actually be seen from imdepisim$events$source, which is 0 for endemic events
and refers to the infective event by the row index otherwise. In the above
example, we obtain:

R> table(imdepisim$events$source)

0 3 11 18 2127 1 1 2 1 Thus, out of 32 simulated events, 5 have been triggered
by previous events.

In case of multiple simulations (nsim > 1), the result is by default
simplified in that only the events instead of a full "epidataCS" object are
retained from every run to save memory and computation time. All other
components, which do not vary between simulations, e.g., the stgrid, are only
stored from the first run. There is a [[-method for such "simEpidataCSlist"s
in order to extract single simulations as full "simEpidataCS" objects also
from the simplified structure.

For more detailed simulation examples, we refer to help("simulate.twinstim").

##  4 twinSIR: SIR event history of a fixed population

When the units where the events can occur are a known discrete set, the point
process is no longer indexed in a continuous spatial domain, but is of
multivariate nature. The conditional intensity function $\lambda(\bm{s},t)$
then becomes $\lambda_{i}(t)$, which characterizes the instantaneous risk of
infection of unit $i$ at time $t$ given the history of the epidemic. For
twinSIR models based on the CIF (2), the individual units are assumed to
belong to a fixed population for which the complete SIR event history is
known. As a case study, we will use a particularly well-documented measles
outbreak among children of the isolated German village Hagelloch in the year
1861, which has previously been analyzed by, e.g., Neal and Roberts (2004).
Section 4.1 introduces the basic class "epidata" for the associated event
history data, and Section 4.2 presents the core functionality of fitting and
analyzing twinSIR point process models. Due to the many similarities with the
twinstim framework covered in Section 3, we condense the present twinSIR
treatment accordingly.

###  4.1 Data structure: "epidata"

New SIR-type event data typically arrive in the form of a simple data frame
with one row per individual and the time points of the sequential events of
the individual as columns. For the 1861 Hagelloch measles epidemic, such a
data set of the 188 affected children is contained in the surveillance
package:

R> data("hagelloch")R> head(hagelloch.df, n = 4)

PN NAME FN HN AGE SEX PRO ERU CL DEAD IFTO SI1 1 Mueller 41 61 7 female
1861-11-21 1861-11-25 1st class <NA> 45 102 2 Mueller 41 61 6 female
1861-11-23 1861-11-27 1st class <NA> 45 123 3 Mueller 41 61 4 female
1861-11-28 1861-12-02 preschool <NA> 172 94 4 Seibold 61 62 13 male 1861-11-27
1861-11-28 2nd class <NA> 180 10 C PR CA NI GE TD TM x.loc y.loc tPRO tERU
tDEAD tR tI1 no complicatons 4 4 3 1 NA NA 142 100 22.7 26.2 NA 29.2 21.72 no
complicatons 4 4 3 1 3 40.3 142 100 24.2 28.8 NA 31.8 23.23 no complicatons 4
4 3 2 1 40.5 142 100 29.6 33.7 NA 36.7 28.64 no complicatons 1 1 1 1 3 40.7
165 102 28.1 29.0 NA 32.0 27.1 The help("hagelloch") contains a description of
all columns. Here we concentrate on the event columns PRO (appearance of
prodromes), ERU (eruption), and DEAD (day of death if during the outbreak). We
take the day on which the index case developed first symptoms, 30 October 1861
(min(hagelloch.df$PRO)$), as the start of the epidemic, i.e., we condition on
this case being initially infectious. As for twinstim, the property of point
processes that concurrent events have zero probability requires special
treatment. Ties are due to the interval censoring of the data to a daily basis
– we broke these ties by adding random jitter to the event times within the
given days. The resulting columns tPRO, tERU, and tDEAD are relative to the
defined start time. Following Neal and Roberts (2004), we assume that each
child becomes infectious (S $\rightarrow$ I event at time tI) one day before
the appearance of prodromes, and is removed from the epidemic (I $\rightarrow$
R event at time tR) three days after the appearance of rash or at the time of
death, whichever comes first.

For further processing of the data, we convert hagelloch.df to the
standardized "epidata" structure for twinSIR. This is done by the converter
function as.epidata, which also checks consistency and optionally pre-
calculates epidemic terms $\bm{x}_{i}(t)$ of Equation 4 to be incorporated in
a twinSIR model. The following call generates the hagelloch "epidata" object:

R> hagelloch <\- as.epidata(\+ hagelloch.df, t0 = 0, tI.col = "tI", tR.col =
"tR",\+ id.col = "PN", coords.cols = c("x.loc", "y.loc"),\+ f = list(\+
household = function(u) u == 0,\+ nothousehold = function(u) u > 0\+ ),\+ w =
list(\+ c1 = function (CL.i, CL.j) CL.i == "1st class" & CL.j == CL.i,\+ c2 =
function (CL.i, CL.j) CL.i == "2nd class" & CL.j == CL.i\+ ),\+ keep.cols =
c("SEX", "AGE", "CL")\+ ) The coordinates (x.loc, y.loc) correspond to the
location of the household the child lives in and are measured in meters. Note
that "twinSIR" allows tied locations of individuals and assumes the relevant
spatial location to be fixed during the entire observation period. The
argument f lists distance-dependent basis functions $B_{m}$ for which the
epidemic terms $\sum_{j\in I(t)}B_{m}(\lVert\bm{s}_{i}-\bm{s}_{j}\rVert)$
shall be generated. Here, household ($x_{i,H}(t)$) and nothousehold
($x_{i,\bar{H}}(t)$) count for each child the number of currently infective
children in its household and outside its household, respectively. Similar to
Neal and Roberts (2004), we also calculate the covariate-based epidemic terms
c1 ($x_{i,c1}(t)$) and c2 ($x_{i,c2}(t)$) counting the number of currently
infective classmates. Note from the corresponding definitions of $w_{ij1}$ and
$w_{ij2}$ in w that c1 is always zero for children of the second class and c2
is always zero for children of the first class. For pre-school children, both
variables equal zero over the whole period. By the last argument keep.cols, we
choose to only keep the covariates SEX, AGE, and school CLass from
hagelloch.df. The first few rows of the generated "epidata" object are shown
below:

R> head(hagelloch, n = 4)

BLOCK id start stop atRiskY event Revent x.loc y.loc SEX AGE CL1 1 1 0 1.14 1
0 0 142 100 female 7 1st class2 1 2 0 1.14 1 0 0 142 100 female 6 1st class3 1
3 0 1.14 1 0 0 142 100 female 4 preschool4 1 4 0 1.14 1 0 0 165 102 male 13
2nd class household nothousehold c1 c21 0 1 0 02 0 1 0 03 0 1 0 04 0 1 0 1

The "epidata" structure inherits from counting processes as implemented by the
class "Surv" of package survival (Therneau, 2014) and also used in, e.g., the
timereg package (Scheike and Zhang, 2011). Specifically, the observation
period is splitted up into consecutive time intervals (start; stop] of
constant conditional intensities. As the CIF $\lambda_{i}(t)$ of Equation (2)
only changes at time points, where the set of infectious individuals $I(t)$ or
some endemic covariate in $\nu_{i}(t)$ change, those occurrences define the
break points of the time intervals. Altogether, the hagelloch event history
consists of 375 time BLOCKs of 188 rows, where each row describes the state of
the individual id during the corresponding time interval. The susceptibility
status and the I- and R-events are captured by the columns atRiskY, event and
Revent, respectively. The atRiskY column indicates if the individual is at
risk of becoming infected in the current interval. The event columns indicate,
which individual was infected or removed at the stop time. Note that at most
one entry in the event and Revent columns is 1, all others are 0.

Display | Subset | Modify
---|---|---
print |  [ |  update
summary |  |
plot |  |
animate |  |
stateplot |  |
Table 6: Generic and non-generic functions applicable to "epidata" objects.

Apart from being the input format for twinSIR models, "epidata" has several
associated methods (Table 6), which are similar in spirit to the methods
described for the "epidataCS" class. For example, Figure 8 illustrates the
course of the Hagelloch measles epidemic by counting processes for the number
of susceptible, infectious and recovered (or “removed”) children,
respectively:

R> plot(hagelloch, xlab = "Time [days]")

Figure 8: Evolution of the 1861 Hagelloch measles epidemic in terms of the
numbers of susceptible, infectious, and recovered children. The bottom rug
marks the infection times tI.

Figure 9 shows the locations of the households, where the size of each dot is
proportional to the number of children in the household:

R> hagelloch_coords <\- summary(hagelloch)$coordinatesR>
plot(hagelloch_coords, xlab = "x [m]", ylab = "y [m]",\+ pch = 15, asp = 1,
cex = sqrt(multiplicity(hagelloch_coords)))R> legend(x = "topleft", legend =
c(1,4,8), pt.cex = sqrt(c(1,4,8)), pch = 15,\+ title = "Household size")

Figure 9: Spatial locations of the Hagelloch households. The size of each dot
is proportional to the number of children in the household.

###  4.2 Modeling and inference

####  4.2.1 Basic example

To illustrate the flexibility of twinSIR we will analyze the Hagelloch data
using class room and household indicators similar to Neal and Roberts (2004).
We include an additional endemic background rate $\exp(\beta_{0})$, which
allows for multiple outbreaks triggered by external sources. Consequently, we
do not need to ignore the child that got infected about one month after the
end of the main epidemic (see the last event mark in Figure 8), as, e.g., done
in a thorough network-based analysis of the Hagelloch data by Groendyke _et
al._ (2012). Altogether, the CIF for a child $i$ is modeled as

$\lambda_{i}(t)=Y_{i}(t)\cdot\left[\exp(\beta_{0})+\alpha_{H}x_{i,H}(t)+\alpha_{c1}x_{i,c1}(t)+\alpha_{c2}x_{i,c2}(t)+\alpha_{\bar{H}}x_{i,\bar{H}}(t)\right]\>,$
(8)

where $Y_{i}(t)=\mathbbm{1}(i\in S(t))$ is the at-risk indicator. By counting
the number of infectious classmates separately for both school classes as
described in the previous section, we allow for class-specific effects
$\alpha_{c1}$ and $\alpha_{c2}$ on the force of infection. The model is
estimated by maximum likelihood as described in Höhle (2009) using the
following call:

R> hagellochFit <\- twinSIR(~ household + c1 + c2 + nothousehold, data =
hagelloch)R> summary(hagellochFit)

Call:twinSIR(formula = ~household + c1 + c2 + nothousehold, data =
hagelloch)Coefficients: Estimate Std. Error z value Pr(>|z|)household 0.026868
0.006113 4.39 1.1e-05 ***c1 0.023892 0.005026 4.75 2.0e-06 ***c2 0.002932
0.000755 3.88 0.0001 ***nothousehold 0.000831 0.000142 5.87 4.3e-09
***cox(logbaseline) -7.362644 0.887989 -8.29 < 2e-16 ***\---Signif. codes: 0
’***’ 0.001 ’**’ 0.01 ’*’ 0.05 ’.’ 0.1 ’ ’ 1Total number of infections:
187One-sided AIC: 1245(simulated penalty weights)Log-likelihood: -619Number of
log-likelihood evaluations: 119

The results show, e.g., a 0.0239 / 0.0029 $=$ 8.15 times higher transmission
between individuals in the 1st class than in the 2nd class. Furthermore, an
infectious housemate adds 0.0269 / 0.0008 $=$ 32.3 times as much infection
pressure as infectious children outside the household. The endemic background
rate of infection in a population with no current measles cases is estimated
to be $\exp(\hat{\beta}_{0})=\exp(-7.36)=0.000635$. An associated Wald
confidence interval (CI) based on the asymptotic normality of the maximum
likelihood estimator (MLE) can be obtained by exp-transforming the confint for
$\beta_{0}$:

R> exp(confint(hagellochFit, parm = "cox(logbaseline)"))

2.5 % 97.5 %cox(logbaseline) 0.000111 0.00362

Note that Wald confidence intervals for the epidemic parameters $\bm{\alpha}$
are to be treated carefully, because their construction does not take the
restricted parameter space into account. For more adequate statistical
inference, the behavior of the log-likelihood near the MLE can be investigated
using the profile-method for "twinSIR" objects. For instance, to evaluate the
normalized profile log-likelihood of $\alpha_{c1}$ and $\alpha_{c2}$ on an
equidistant grid of 25 points within the corresponding 95% Wald CIs, we do:

R> prof <\- profile(hagellochFit,\+ list(c(match("c1",
names(coef(hagellochFit))), NA, NA, 25),\+ c(match("c2",
names(coef(hagellochFit))), NA, NA, 25))) The profiling result contains 95%
highest likelihood based CIs for the parameters, as well as the Wald CIs for
comparison:

R> prof$ci.hl

idx hl.low hl.up wald.low wald.up mlec1 2 0.01522 0.03497 0.01404 0.03374
0.02389c2 3 0.00158 0.00454 0.00145 0.00441 0.00293 An advantage of the
profile likelihood CIs for the $\bm{\alpha}$’s is that the parameter
constraints are more naturally taken into account. The entire functional form
of the normalized profile log-likelihood on the requested grid as stored in
prof$lp$ can be visualized by:

R> plot(prof)

Figure 10: Normalized log-likelihood for $\alpha_{c1}$ and $\alpha_{c2}$ when
fitting the twinSIR model formulated in Equation (8) to the Hagelloch data.

####  4.2.2 Model Diagnostics

Display | Extract | Other
---|---|---
print |  vcov |  simulate
summary |  logLik |
plot |  AIC |
intensityplot |  extractAIC |
checkResidualProcess |  profile |
|  residuals |
Table 7: Generic and non-generic functions for "twinSIR". There are no
specific coef or confint methods, since the respective default methods from
package stats apply outright.

Table 7 lists all methods for the "twinSIR" class. For example, to investigate
how the CIF decomposes into endemic and epidemic intensity over time, we
produce Figure LABEL:fig:hagellochFit_plot1 by:

R> plot(hagellochFit, which = "epidemic proportion", xlab = "time [days]")

(a)

(b)

Figure 11: Diagnostic plots for the twinSIR model formulated in Equation 8.

Note that the last infection was necessarily caused by the endemic component
since there were no more infectious children in the observed population which
could have triggered the new case. We can also inspect temporal Cox-Snell-like
residuals of the fitted point process using the function checkResidualProcess
as for the spatio-temporal point process models "twinstim" in Section 3.2. The
resulting Figure LABEL:fig:hagellochFit_plot2 reveals some deficiencies of the
model in describing the waiting times between events, which might be related
to the assumption of fixed infection periods.

Finally, "twinSIR"’s AIC-method computes the one-sided AIC (Hughes and King,
2003) as described in Höhle (2009), which can be used for model selection. For
instance, we may consider a more flexible model for local spread using a step
function for the distance kernel $f(u)$ in Equation 3. An updated model with
$B_{1}=I_{(0;100)}(u)$, $B_{2}=I_{[100;200)}(u)$, $B_{3}=I_{[200;\infty)}(u)$
can be fitted as follows:

R> knots <\- c(100, 200)R> fstep <\- list(B1 = function(D) D > 0 & D <
knots[1],\+ B2 = function(D) D >= knots[1] & D < knots[2],\+ B3 = function(D)
D >= knots[2])R> hagellochFit_fstep <\- twinSIR(~ household + c1 + c2 + B1 +
B2 + B3,\+ data = update(hagelloch, f = fstep)) According to

R> AIC(hagellochFit, hagellochFit_fstep)

df AIChagellochFit 5 1245hagellochFit_fstep 7 1246 the simpler model with just
a nothousehold component instead of the more flexible distance-based step
function is preferred.

###  4.3 Simulation

Simulation from fitted twinSIR models is possible as described in detail in
Höhle (2009, Section 4). The implementation is made available by an
appropriate simulate-method for class "twinSIR". Because both the algorithm
and the call are similar to the invocation on "twinstim" objects already
described in Section 3.3, we skip the illustration here and refer to
help("simulate.twinSIR").

##  5 hhh4: Multivariate time series of counts

In public health surveillance, the most common type of spatio-temporal data
arise as part of routine case reporting to public health authorities. The
Robert Koch Institute (RKI) in Germany for example maintains a database of
cases of notifiable diseases under the German _Protection against Infection
Act_ (“Infektionsschutzgesetz”). Geographically and temporally aggregated data
can be retrieved through online queries via the
SurvStat@RKI333https://survstat.rki.de application, which we used to obtain
counts of measles by week and district in the Weser-Ems region of Lower
Saxony, Germany, during the years 2001 and 2002 (as of Annual Report 2005).
These spatio-temporal count data constitute the response $Y_{it}$,
$i=1,\dotsc,17$, $t=1,\dotsc,104$, for our illustration of the endemic-
epidemic multivariate time-series model “hhh4” (5). In Section 5.1, we briefly
present the associated data class "sts" (“surveillance time series”). In
Section 5.2, a simple model for the measles data based on the original
analysis by Held _et al._ (2005) is introduced, which is then sequentially
improved by suitable model extensions. The final Section 5.3 illustrates
simulation from fitted "hhh4" models.

###  5.1 Data structure: "sts"

The epidemic modeling of multivariate count time-series essentially involves
two data matrices: the nEpochs$\times$nUnits matrix of observed counts, and a
nUnits$\times$nUnits matrix quantifying the interaction of the units. In this
paper, we focus on the spatio-temporal setting, in which case the “units” are
geographical districts whose adjacency matrix provides the simplest measure of
spatial interaction. A "SpatialPolygons" representation of the districts
enables nice visualizations and can be used to derive the districts’ adjacency
matrix automatically using the function poly2adjmat, which wraps around
functionality of package spdep (Bivand, 2014):

R> weserems_adjmat <\- poly2adjmat(map) Given this binary matrix, the function
nbOrder computes the distance matrix to enable modeling of spatial interaction
as a function of neighbourhood order (Meyer and Held, 2014b):

R> weserems_nbOrder <\- nbOrder(weserems_adjmat, maxlag = 10)

In infectious disease epidemiology, it is natural to accompany the counts with
district-specific population numbers in order to report disease incidence or
to account for population size in statistical models. The population data must
obey the same matrix structure as the observed counts, and for numerical
reasons, we incorporate population fractions $e_{i}=n_{i}/\sum_{i}n_{i}$
rather than absolute numbers $n_{i}$. Constructing an object of the S4-class
"sts" from the aforementioned ingredients is as simple as:444The
data(”measlesWeserEms”) is part of the surveillance package and constitutes a
corrected version of data(”measles.weser”) used by Held _et al._ (2005).
Differences are documented on the associated help page.

R> measlesWeserEms <\- new("sts", start = c(2001, 1), freq = 52, epoch =
1:nrow(counts),\+ observed = counts, neighbourhood = weserems_nbOrder,\+ map =
map, populationFrac = populationFrac)

Here, start and freq have the same meaning as for classical time-series
objects of class "ts", i.e., (year, sample number) of the first observation
and the number of observations per year. See Höhle and Mazick (2010) and
Salmon _et al._ (2014) for detailed descriptions of the "sts" class, also used
for the prospective aberration detection facilities of the surveillance
package.

We can visualize such aggregated surveillance data in four ways: individual
time series, overall time series, map of accumulated counts by district, or
animated maps. For instance, Figure 12 has been produced by the following
code:

R> plot(measlesWeserEms, type = observed ~ time, legend.opts = NULL)R>
plot(measlesWeserEms, type = observed ~ unit, labels = list(font = 2))

(a)

(b)

Figure 12: Counts of measles in the Weser-Ems region.

The default plot type is observed time | unit in order to plot the individual time series by district. The output is shown in Figure 13 excluding the districts 03401 (SK Delmenhorst) and 03405 (SK Wilhelmshaven) without any reported cases:

R> measlesWeserEms15 <\- measlesWeserEms[, colSums(observed(measlesWeserEms))
> 0]R> plot(measlesWeserEms15)

Figure 13: Individual time series of the 15 affected districts.

An animation of the data can be easily produced as well. We recommend to use
converters of the animation package, e.g., to watch the series of plots in a
web browser. The following code will generate weekly disease maps during the
year 2001 with the respective total number of cases shown in a legend and an
evolving time series plot at the bottom:

R> library("animation")R> saveHTML(animate(measlesWeserEms, tps = 1:52,
total.args = list()),\+ title = "Evolution of the measles epidemic in the
Weser-Ems region",\+ ani.width = 500, ani.height = 600)

###  5.2 Modeling and inference

We start by modeling the measles counts in the Weser-Ems region by a slightly
simplified version of the original negative binomial model by Held _et al._
(2005). Instead of district-specific intercepts $\alpha_{i}$ in the endemic
component, we first assume a common intercept $\alpha$ in order to not be
forced to restrict the model to the 15 affected districts. After the
estimation and illustration of this basic model, we will discuss the following
sequential extensions: covariates (district-specific vaccination coverage),
neighbourhood weights, and random effects to eventually account for unobserved
heterogeneity of the districts.

####  5.2.1 Basic model

Our initial model has the following mean structure:

$\displaystyle\mu_{it}$
$\displaystyle=e_{i}\,\nu_{t}+\lambda\,Y_{i,t-1}+\phi\sum_{j\neq
i}w_{ji}Y_{j,t-1}\>,$ (9) $\displaystyle\log(\nu_{t})$
$\displaystyle=\alpha^{(\nu)}+\beta_{t}t+\gamma\sin(\omega
t)+\delta\cos(\omega t)\>.$ (10)

To account for temporal variation of disease incidence, the endemic log-linear
predictor $\nu_{t}$ incorporates an overall trend and a sinusoidal wave of
frequency $\omega=2\pi/52$. As a basic district-specific measure of disease
incidence, the population fraction $e_{i}$ is included as a multiplicative
offset. The effects of the epidemic components,
$\lambda=\exp(\alpha^{(\lambda)})$ and $\phi=\exp(\alpha^{(\phi)})$, are
assumed homogeneous across districts and constant over time. Note that the
superscripts in brackets are used to distinguish the component-specific
intercept parameters. Furthermore, we define $w_{ji}=\mathbbm{1}(j\sim i)$ for
the time being, which means that the epidemic can only arrive from directly
adjacent districts.

This "hhh4" model transforms into the following list of control arguments:

R> measlesModel_basic <\- list(\+ end = list(f = addSeason2formula(~1 + t,
period = measlesWeserEms@freq),\+ offset = population(measlesWeserEms)),\+ ar
= list(f = ~1),\+ ne = list(f = ~1, weights = neighbourhood(measlesWeserEms)
== 1),\+ family = "NegBin1")

The formulae of the three predictors $\log\nu_{t}$, $\log\lambda$ and
$\log\phi$ are specified as element f of the end, ar, and ne lists,
respectively. For the endemic formula we use the convenient function
addSeason2formula to generate the sine-cosine terms, and we take the
multiplicative population offset $e_{i}$ from the measlesWeserEms object. The
autoregressive part only consists of the intercept $\alpha^{(\lambda)}$,
whereas the spatio-temporal component additional specifies the matrix of
neighbourhood weights $(w_{ji})$ to use – here an indication of first-order
neighbourhood. By family = "NegBin1" we choose a negative binomial model with
a common overdispersion parameter $\psi$ for all districts. The alternatives
are "Poisson" ($\psi\equiv 0$), and "NegBinM" (district-specific $\psi_{i}$).
The complete list of control arguments is then fed into the hhh4 function
together with the data to estimate the model, a summary of which is printed
below.

R> measlesFit_basic <\- hhh4(stsObj = measlesWeserEms, control =
measlesModel_basic)R> summary(measlesFit_basic, idx2Exp = 1:4, amplitudeShift
= TRUE, maxEV = TRUE)

Call:hhh4(stsObj = measlesWeserEms, control = measlesModel_basic)Coefficients:
Estimate Std. Errorexp(ar.1) 0.64540 0.07927exp(ne.1) 0.01581
0.00420exp(end.1) 1.08025 0.27884exp(end.t) 1.00119 0.00426end.A(2 * pi *
t/52) 1.16423 0.19212end.s(2 * pi * t/52) -0.63436 0.13350overdisp 2.01384
0.28544Epidemic dominant eigenvalue: 0.72Log-likelihood: -972AIC: 1957BIC:
1996Number of units: 17Number of time points: 103 The idx2Exp=1:4 argument
yields the estimates for $\lambda$, $\phi$, $\alpha^{(\nu)}$ and
$\exp(\beta_{t})$ instead of their respective internal log-values. For
instance, exp(end.t) represents the seasonality-adjusted factor by which the
basic endemic incidence increases per week. The amplitudeShift argument
transforms the internal coefficients $\gamma$ and $\delta$ of the sine-cosine
terms to the amplitude $A$ and phase shift $\varphi$ of the corresponding
sinusoidal wave $A\sin(\omega t+\varphi)$ in $\log\nu_{t}$ (Paul _et al._ ,
2008). The multiplicative effect of seasonality on $\nu_{t}$ is shown in
Figure 14 produced by:

R> plot(measlesFit_basic, type = "season", components = "end", main = "")

Figure 14: Estimated multiplicative effect of seasonality on the endemic mean.

The overdisp parameter and its 95% confidence interval obtained by

R> confint(measlesFit_basic, parm = "overdisp")

2.5 % 97.5 %overdisp 1.45 2.57 suggest that a negative binomial distribution
with overdispersion is more adequate than a Poisson model corresponding to
$\psi=0$. We can underpin this finding by comparing AIC values:

R> AIC(measlesFit_basic, update(measlesFit_basic, family = "Poisson"))

df AICmeaslesFit_basic 7 1957update(measlesFit_basic, family = "Poisson") 6
2479

The epidemic potential of the process as determined by the parameters
$\lambda$ and $\phi$ is best investigated by a combined measure: the dominant
eigenvalue (maxEV) of the matrix $\bm{\Lambda}$ which has the entries
$(\Lambda)_{ii}=\lambda$ on the diagonal and $(\Lambda)_{ij}=\phi w_{ji}$ for
$j\neq i$ (Paul _et al._ , 2008). If the dominant eigenvalue is smaller than
unity, it can be interpreted as the epidemic proportion of total disease
incidence. In the above model, the estimate is 72%. A more intuitive way of
judging the relative importance of the three model components is offered by
Figure 15. It shows the fitted component means along with the observed time
series for the six districts with more than 20 cases, and has been produced
by:

R> districts2plot <\- names(which(colSums(observed(measlesWeserEms)) > 20))R>
plot(measlesFit_basic, type = "fitted", units = districts2plot, hide0s = TRUE)

Figure 15: Fitted components in the initial model measlesFit_basic for the six
districts with more than 20 cases. Dots are only drawn for positive weekly
counts.

Other plot types and methods for fitted "hhh4" models as listed in Table 8
will be applied in the course of the following model extensions.

Display | Extract | Modify | Other
---|---|---|---
print |  nobs |  update |  predict
summary |  coef |  |  simulate
plot |  fixef |  |  oneStepAhead
|  ranef |  |
|  vcov |  |
|  confint |  |
|  logLik |  |
|  residuals |  |
|  terms |  |
|  formula |  |
Table 8: Functions applicable to "hhh4" objects.

####  5.2.2 Covariates

The "hhh4" model framework allows for covariate effects on the endemic or
epidemic contributions to disease incidence. Covariates may vary over both
regions and time and thus obey the same matrix structure as the observed
counts. For infectious disease models, the regional vaccination coverage is an
important example of such a covariate, since it reflects the (remaining)
susceptible population. In a thorough analysis of measles occurrence in the
German federal states, Herzog _et al._ (2011) found vaccination coverage to be
associated with outbreak size. We follow their approach of using the district-
specific proportion $1-v_{i}$ of unvaccinated children just starting school as
a proxy for the susceptible population. As $v_{i}$ we use the proportion of
children vaccinated with at least one dose among the ones presenting their
vaccination card at school entry in district $i$ in the year 2004.555This is
the first year with complete data for all 17 districts, available from the
public health department of Lower Saxony at
http://www.nlga.niedersachsen.de/portal/live.php?navigation_id=27093 (in
German). This time-constant covariate needs to be transformed to the common
matrix structure for incorporation in hhh4:

R> Sprop <\- matrix(1 - measlesWeserEms<EMAIL_ADDRESS>byrow = TRUE,\+
nrow = nrow(measlesWeserEms), ncol = ncol(measlesWeserEms))R>
summary(Sprop[1,])

Min. 1st Qu. Median Mean 3rd Qu. Max. 0.0306 0.0481 0.0581 0.0675 0.0830
0.1400

There are several ways to account for the susceptible proportion in our model,
among which the simplest is to update the endemic population offset $e_{i}$ by
multiplication with $(1-v_{i})$. Herzog _et al._ (2011) found that the
susceptible proportion is best added as a covariate in the autoregressive
component in the form

$\lambda_{i}\,Y_{i,t-1}=\exp\big{(}\alpha^{(\lambda)}+\beta_{s}\log(1-v_{i})\big{)}\,Y_{i,t-1}=\exp\big{(}\alpha^{(\lambda)}\big{)}\,(1-v_{i})^{\beta_{s}}\,Y_{i,t-1}$

according to the mass action principle (Keeling and Rohani, 2008). A higher
proportion of susceptibles in district $i$ is expected to boost the generation
of new infections, i.e., $\beta_{s}>0$. Alternatively, this effect could be
assumed as an offset, i.e., $\beta_{s}\equiv 1$. To choose between endemic
and/or autoregressive effects, and multiplicative offset vs. covariate
modeling, we perform AIC-based model selection. First, set up a grid of all
combinations of envisaged extensions for the endemic and autoregressive
components:

R> Soptions <\- c("unchanged", "Soffset", "Scovar")R> SmodelGrid <\-
expand.grid(end = Soptions, ar = Soptions)R> row.names(SmodelGrid) <\-
do.call("paste", c(SmodelGrid, list(sep = "|"))) Then update the initial model
measlesFit_basic according to each row of SmodelGrid:

R> measlesFits_vacc <\- apply(X = SmodelGrid, MARGIN = 1, FUN = function
(options) {\+ updatecomp <\- function (comp, option)\+ switch(option,
"unchanged" = list(),\+ "Soffset" = list(offset = comp$offset * Sprop),\+
"Scovar" = list(f = update(comp$f, ~. + log(Sprop))))\+
update(measlesFit_basic,\+ end = updatecomp(measlesFit_basic$control$end,
options[1]),\+ ar = updatecomp(measlesFit_basic$control$ar, options[2]),\+
data = list(Sprop = Sprop))\+ }) The resulting object measlesFits_vacc is a
list of 9 "hhh4" fits, which are named according to the corresponding Soptions
used for the endemic and autoregressive component. We construct a call of the
function AIC taking all list elements as arguments:

R> AIC.list <\- function (object, ..., k = 2)\+ {\+ if
(is.null(names(object))) stop("the list of models must have names")\+
eval(as.call(c(lapply(c("AIC", names(object)), as.name), list(k = k))),\+
envir = object)\+ }R> aics_vacc <\- AIC(measlesFits_vacc)R>
aics_vacc[order(aics_vacc[, "AIC"]), ]

df AIC`Scovar|unchanged` 8 1917`Scovar|Scovar` 9 1919`Soffset|unchanged` 7
1922`Soffset|Scovar` 8 1924`Scovar|Soffset` 8 1934`Soffset|Soffset` 7
1937unchanged|unchanged 7 1957`unchanged|Scovar` 8 1959`unchanged|Soffset` 7
1967

Hence, AIC increases if the susceptible proportion is only added to the
autoregressive component, but we see a remarkable improvement when adding it
to the endemic component. The best model is obtained by leaving the
autoregressive component unchanged ($\lambda$) and adding the term
$\beta_{s}\log(1-v_{i})$ to the endemic predictor in Equation 10.

R> measlesFit_vacc <\- measlesFits_vacc[["Scovar|unchanged"]]R>
coef(measlesFit_vacc, se = TRUE)["end.log(Sprop)", ]

Estimate Std. Error 1.718 0.288 The estimated exponent $\hat{\beta}_{s}$ is
both clearly positive and different from the offset assumption. In other
words, if a district’s fraction of susceptibles is doubled, the endemic
measles incidence is estimated to multiply by $2^{\hat{\beta}_{s}}=$ 3.29 (95%
CI: 2.23–4.86).

####  5.2.3 Spatial interaction

Up to now, the model assumed that the epidemic can only arrive from directly
adjacent districts because $w_{ji}=\mathbbm{1}(j\sim i)$, and that all
districts have the same coefficient $\phi$ for importing cases from
neighbouring regions. Given the ability of humans to travel further and
preferrably to metropolitan areas, both assumptions seem overly simplistic.
First, to reflect commuter-driven spread in our model, we scale the district’s
susceptibility according to its population fraction by multiplying $\phi$ by
$e_{i}^{\beta_{pop}}$:

R> measlesFit_nepop <\- update(measlesFit_vacc,\+ ne = list(f = ~log(pop)),\+
data = list(pop = population(measlesWeserEms))) As in an analysis of influenza
in Southern Germany’s districts (Meyer and Held, 2014b), we find strong
evidence for such an agglomeration effect: the estimated exponent is
$\hat{\beta}_{pop}=$ 2.85 (95% CI: 1.83–3.87) and AIC decreases from 1917 to
1887. Geilhufe _et al._ (2014) additionally incorporated data on incoming road
and air traffic to model immigration more explicitly.

To account for long-range transmission of cases, Meyer and Held (2014b)
proposed the estimation of the weights $w_{ji}$ as a function of the order of
neighbourhood $o_{ji}$ between districts. Specifically, a power-law model was
suggested, which assumes the form $w_{ji}=o_{ji}^{-d}$, for $j\neq i$ and
$w_{jj}=0$, where the decay parameter $d$ is to be estimated. Normalization to
$w_{ji}/\sum_{k}w_{jk}$ is recommended and applied by default when specifying
W_powerlaw as neighbourhood weights:

R> measlesFit_powerlaw <\- update(measlesFit_nepop,\+ ne = list(weights =
W_powerlaw(maxlag = 5))) The argument maxlag sets an upper bound for spatial
interaction in terms of neighbourhood order. Here we set no limit since
max(neighbourhood(measlesWeserEms)) is 5. The resulting parameter estimate is
$\hat{d}=$ 4.10 (95% CI: 2.03–6.17), which represents a strong decay of
spatial interaction for higher-order neighbours. As an alternative to the
parametric power law, order-specific weights up to maxlag can be modeled by
using W_np instead of W_powerlaw. For instance, W_np(maxlag=2) corresponds to
a second-order model, i.e.,
$w_{ji}=1\cdot\mathbbm{1}(o_{ji}=1)+e^{\omega_{2}}\cdot\mathbbm{1}(o_{ji}=2)$,
which is also row-normalized by default:

R> measlesFit_np2 <\- update(measlesFit_nepop,\+ ne = list(weights =
W_np(maxlag = 2)))R> AIC(measlesFit_nepop, measlesFit_powerlaw,
measlesFit_np2)

df AICmeaslesFit_nepop 9 1887measlesFit_powerlaw 10 1882measlesFit_np2 10 1881
In spite of the second-order model resulting in a slightly better fit, we will
use the power-law model as a basis for further model extensions since the
stand-alone second-order effect is not always identifiable in more complex
models and is scientifically implausible.

Figure LABEL:fig:measlesFit_neweights2 shows both the power law $o^{-\hat{d}}$
and the second-order weight $e^{\hat{\omega}_{2}}=$ 0.09 (95% CI: 0.02–0.39).
Alternatively, the plot type="neweights" for "hhh4" objects can produce a
stripplot (Sarkar, 2008) of $w_{ji}$ against $o_{ji}$ as shown in Figure
LABEL:fig:measlesFit_neweights1 for the power-law model:

R> library("lattice")R> plot(measlesFit_powerlaw, type = "neweights", plotter
= stripplot,\+ panel = function (...) {panel.stripplot(...);
panel.average(...)},\+ jitter.data = TRUE, xlab = expression(o[ji]), ylab =
expression(w[ji]))

(a)

(b)

Figure 16: Estimated weights as a function of neighbourhood order.

Note that only horizontal jitter is added in this case. Because of
normalization, the weight $w_{ji}$ for transmission from district $j$ to
district $i$ is determined not only by the districts’ neighbourhood $o_{ji}$
but also by the total amount of neighbourhood of district $j$ in the form of
$\sum_{k\neq j}o_{jk}^{-d}$, which causes some variation of the weights for a
specific neighbourhood order.

The neighbourhood-based approach to modeling spatial interaction is especially
attractive if movement network data are not available. Otherwise, a matrix
quantifying the epidemiological coupling between different regions (Keeling
and Rohani, 2008, Chapter 7) could enter as weights in the neighbourhood
component (Paul _et al._ , 2008; Geilhufe _et al._ , 2014).

####  5.2.4 Random effects

Paul and Held (2011) introduced random effects for "hhh4" models, which are
useful if the districts exhibit heterogeneous incidence levels not explained
by observed covariates. For infectious disease surveillance data, a typical
example of unobserved heterogeneity is under-reporting (Bernard _et al._ ,
2014). Our measles data from the Weser-Ems region even contain two districts
without any reported cases within the two years, while the district with the
smallest population (03402, SK Emden) had the second-largest number of cases
reported (Figure 13). Hence, allowing for district-specific intercepts in the
endemic or epidemic components is expected to improve the model fit. To
estimate independent random effects $\alpha_{i}^{(\nu)}\stackrel{{\scriptstyle
iid}}{{\sim}}\operatorname{N}(\alpha^{(\nu)},\sigma_{\nu}^{2})$,
$\alpha_{i}^{(\lambda)}\stackrel{{\scriptstyle
iid}}{{\sim}}\operatorname{N}(\alpha^{(\lambda)},\sigma_{\lambda}^{2})$, and
$\alpha_{i}^{(\phi)}\stackrel{{\scriptstyle
iid}}{{\sim}}\operatorname{N}(\alpha^{(\phi)},\sigma_{\phi}^{2})$ in all three
components, we update the corresponding formulae as follows:

R> measlesFit_ri <\- update(measlesFit_powerlaw,\+ end = list(f =
update(formula(measlesFit_powerlaw)$end, ~ . + ri() - 1)),\+ ar = list(f =
update(formula(measlesFit_powerlaw)$ar, ~ . + ri() - 1)),\+ ne = list(f =
update(formula(measlesFit_powerlaw)$ne, ~ . + ri() - 1)))

R> summary(measlesFit_ri, amplitudeShift = TRUE, maxEV = TRUE)

Call:hhh4(stsObj = object$stsObj, control = control)Random effects: Var
Corrar.ri(iid) 1.076ne.ri(iid) 1.294 0end.ri(iid) 1.312 0 0Fixed effects:
Estimate Std. Errorar.ri(iid) -1.61389 0.38197ne.log(pop) 3.42406
1.07722ne.ri(iid) 6.62429 2.81553end.t 0.00578 0.00480end.A(2 * pi * t/52)
1.20359 0.20149end.s(2 * pi * t/52) -0.47916 0.14205end.log(Sprop) 1.79350
0.69159end.ri(iid) 4.42260 1.94605neweights.d 3.60640 0.77602overdisp 0.97723
0.15132Epidemic dominant eigenvalue: 0.84Penalized log-likelihood:
-869Marginal log-likelihood: -54.2Number of units: 17Number of time points:
103 The summary now contains an extra section with the estimated variance
components $\sigma_{\lambda}^{2}$, $\sigma_{\phi}^{2}$, and $\sigma_{\nu}^{2}$
of the random effects. We did not assume correlation between the three
intercepts, but this is possible by specifying ri(corr="all") in the component
formulae. The implementation also supports a conditional autoregressive
formulation (Besag _et al._ , 1991) for spatially correlated intercepts by
using ri(type="car").

The estimated district-specific intercepts of the three components are not
shown in the above coefficient summary. Their values can be extracted by the
ranef-method,

R> head(ranef(measlesFit_ri, tomatrix = TRUE))

ar.ri(iid) ne.ri(iid) end.ri(iid)03401 0.000 -0.0567 -1.00403402 1.223 0.0431
1.52603403 -0.827 1.5588 -0.62003404 -0.377 -0.5013 -1.11903405 0.000 -0.0911
-1.28203451 -0.524 -0.1135 0.496 and they can also be visualized in a map by
the plot type="ri":

R> for (comp in c("ar", "ne", "end")) {\+ print(plot(measlesFit_ri, type =
"ri", component = comp, xlim = c(6.6,8.8),\+ labels = list(cex = 0.6), at =
seq(-1.7, 1.7, length.out = 15)))\+ }

(a)

(b)

(c)

Figure 17: Estimated district-specific random effects in measlesFit_ri.

For the autoregressive component in Figure LABEL:fig:measlesFit_ri_map1, we
see a pronounced heterogeneity between the three western districts in blue and
the remaining districts. These three districts have been affected by large
local outbreaks and are also the ones with the highest overall numbers of
cases. In contrast, the city of Oldenburg (03403) is estimated with a
relatively low autoregressive factor
$\lambda_{i}=\exp(\alpha^{(\lambda)}+\alpha_{i}^{(\lambda)})=0.087$, but it
seems to import more cases from other districts than explained by its
population (Figure LABEL:fig:measlesFit_ri_map2). In Figure
LABEL:fig:measlesFit_ri_map3, the two districts without any reported measles
cases (03401 and 03405) appear in dark pink, which means that they exhibit a
relatively low endemic incidence after adjusting for the population and
susceptible proportion. Such districts could be suspected of a larger amount
of under-reporting.

Note that the extra flexiblility of the random effects model comes at a price.
First, the estimation runtime increases considerably from 0.2 seconds for the
previous power-law model measlesFit_powerlaw to 4.9 seconds with additional
random effects. Furthermore, we no longer obtain AIC values in the model
summary, since random effects invalidate simple AIC-based model comparisons
(Greven and Kneib, 2010). Of course we can plot the fitted values and visually
compare their quality with the initial fit shown in Figure 15:

R> plot(measlesFit_ri, type = "fitted", units = districts2plot, hide0s = TRUE)

Figure 18: Fitted components in the random effects model measlesFit_ri for the
six districts with more than 20 cases. Compare to Figure 15.

For some of these districts, a great amount of cases is now explained via
transmission from neighbouring regions while others are mainly influenced by
the local autoregression. However, for quantitative model comparisons we have
to resort to more sophisticated techniques presented in the next section.

####  5.2.5 Predictive model assessment

Paul and Held (2011) suggest to evaluate one-step-ahead forecasts from
competing models by proper scoring rules for count data (Czado _et al._ ,
2009). These scores measure the discrepancy between the predictive
distribution $P$ from a fitted model and the later observed value $y$. A well-
known example is the squared error score (“ses”) $(y-\mu_{P})^{2}$, which is
usually averaged over a suitable set of forecasts to obtain the mean squared
error. More elaborate scoring rules such as the logarithmic score (“logs”) or
the ranked probability score (“rps”) take into account the whole predictive
distribution to assess calibration and sharpness simultaneously – see the
recent review by Gneiting and Katzfuss (2014). Lower scores correspond to
better predictions.

In the "hhh4" framework, predictive model assessment is made available by the
functions oneStepAhead and scores. We will use the second quarter of 2002 as
the test period, and compare the basic model, the power-law model, and the
random effects model. First, we use the "final" fits on the complete time
series to compute the predictions, which then simply correspond to the fitted
values during the test period:

R> tp <\- c(65, 77)R> models2compare <\- paste0("measlesFit_", c("basic",
"powerlaw", "ri"))R> measlesPreds1 <\- lapply(mget(models2compare),
oneStepAhead, tp=tp, type="final")

R> SCORES <\- c("logs", "rps", "ses")R> measlesScores1 <\-
lapply(measlesPreds1, scores, which = SCORES, individual = TRUE) Note that in
this case, the log-score for a model’s prediction in district $i$ in week $t$
equals the associated negative log-likelihood contribution. Comparing the mean
scores from different models is thus essentially a goodness-of-fit assessment:

R> t(sapply(measlesScores1, colMeans, dims = 2))

logs rps sesmeaslesFit_basic 1.09 0.736 5.29measlesFit_powerlaw 1.10 0.731
5.39measlesFit_ri 1.01 0.638 4.82 All scoring rules claim that the random
effects model gives the best fit during the second quarter of 2002. Now we
turn to true one-week-ahead predictions of type="rolling", which means that we
always refit the model up to week $t$ to get predictions for the next week
$t+1$:

R> measlesPreds2 <\- lapply(mget(models2compare), oneStepAhead, tp = tp,\+
type = "rolling", which.start = "final",\+ cores = 2 * (.Platform$OS.type ==
"unix"))R> measlesScores2 <\- lapply(measlesPreds2, scores, which = SCORES,
individual = TRUE)R> t(sapply(measlesScores2, colMeans, dims = 2))

logs rps sesmeaslesFit_basic 1.10 0.748 5.40measlesFit_powerlaw 1.14 0.765
5.87measlesFit_ri 1.11 0.763 7.08 Thus, the most parsimonious initial model
measlesFit_basic gives the best one-week-ahead predictions in terms of overall
mean scores. Statistical significance of the differences in mean scores can be
investigated by a permutation test for paired data or a paired t-test, both of
which are computed by the permutationTest function:

R> sapply(SCORES, function (score) {\+
permutationTest(measlesScores2$measlesFit_ri[,,score],\+
measlesScores2$measlesFit_basic[,,score])\+ })

logs rps sesdiffObs 0.00782 0.0154 1.68pVal.permut 0.863 0.723 0.18pVal.t
0.854 0.717 0.171 Hence, there is no clear evidence for a difference between
the basic and the random effects model with regard to predictive performance
during the test period. In longer surveillance time series such assessments
might reveal more obvious differences.

Czado _et al._ (2009) also describe an informal graphical approach to assess
calibration of the predictions: (non-randomized) probability integral
transform (PIT) histograms for count data. Under the hypothesis of
calibration, i.e., $y_{it}\sim P_{it}$ for all predictive distributions
$P_{it}$ in the test period, the PIT histogram is uniform. Underdispersed
predictive distributions lead to U-shaped histograms, and biased predictions
cause skewness. Such non-randomized PIT histograms for count data are
available by calling the function pit:

R> for (m in models2compare) {\+ with(measlesPreds2[[m]],\+ pit(x = observed,
pdistr = "pnbinom", mu = pred, size = exp(psi),\+ plot = list(ylim =
c(0,1.25), main = m)))\+ }

Figure 19: Non-randomized PIT histograms of competing models to check
calibration of the one-week-ahead predictions during the second quarter of
2002.

Note that the parameters mu and size ($=1/\psi$) of the predictive
distribution are different for each predicted time point. The psi component of
the predictions measlesPreds2[[m]] contains the estimated overdispersion
parameter on the internal scale $\log(1/\psi)$, thus we have to specify
size=exp(psi) to obtain the parametrization expected by the pnbinom function.
Figure 19 shows the PIT histograms of the competing models generated by the
above code. In this aggregate view of the predictions over all districts and
weeks of the test period, the predictive performance is indeed comparable
between the models, and there is no evidence of badly dispersed predictions.
However, the right-hand decay in all histograms suggests that all models tend
to predict higher counts than observed. This is most likely related to the
seasonal shift between the years 2001 and 2002. In 2001, the peak of the
epidemic was in the second quarter, while it already occured in the first
quarter in 2002 (cp. Figure LABEL:fig:measlesWeserEms1). Formal statistical
tests for the calibration of count forecasts have been recently proposed by
Wei and Held (2014).

####  5.2.6 Further modeling options

In the previous sections we extended our model for measles in the Weser-Ems
region with respect to spatial variation of the counts and their interaction.
Temporal variation was only accounted for in the endemic component, which
included a long-term trend and a sinusoidal wave on the log-scale. Held and
Paul (2012) suggest to also allow seasonal variation of the epidemic force by
adding a superposition of $S$ harmonic waves of fundamental frequency
$\omega$,

$\sum_{s=1}^{S}\left\\{\gamma_{s}\sin(s\omega t)+\delta_{s}\cos(s\omega
t)\right\\}\>,$

to the autoregressive and/or spatio-temporal log-linear predictors – just like
for $\log\nu_{t}$ in Equation 10 with $S=1$. However, given only two years of
measles surveillance and the apparent shift of seasonality with regard to the
start of the outbreak in 2002 compared to 2001 (Figure
LABEL:fig:measlesWeserEms1), more complex seasonal models are likely to
overfit the data. Concerning the coding in R, adding sine-cosine terms in the
epidemic components poses no difficulties by again using the convenience
function addSeason2formula. Updating a previous model for different numbers of
harmonics is even simpler, since the update-method has a corresponding
argument S. The plots of type="season" and type="maxEV" for "hhh4" fits can
visualize the estimated component seasonality.

All of our models for the measles surveillance data incorporated an epidemic
effect of the counts from the local district and its neighbours. Without
further notice, we thereby assumed a lag equal to the observation interval of
one week. However, the generation time of measles is around 10 days (Anderson
and May, 1991), which is why some studies, e.g., Finkenstädt _et al._ (2002)
or Herzog _et al._ (2011), aggregate their weekly measles surveillance data
into biweekly intervals. Fine and Clarkson (1982) used weekly counts in their
analysis and report that biweekly aggregation would have little effect on the
results. We can also perform such a sensitivity analysis by running the whole
code of the current section based on aggregate(measlesWeserEms, nfreq = 26).
Doing so, the parameter estimates of the various models retain their order of
magnitude and conclusions remain the same. However, with the number of time
points halved, the complex random effects model would no longer be
identifiable for all subsets of the data when calculating one-week-ahead
predictions during the test period.

###  5.3 Simulation

Simulation from fitted "hhh4" models is enabled by an associated simulate-
method. Compared to the point process models of Sections 3 and 4, simulation
is less complex since it essentially consists of sequential calls of rnbinom
(or rpois) from the base package stats. At each time point $t$, the mean
$\mu_{it}$ is determined by plugging in the parameter estimates and the counts
$Y_{i,t-1}$ simulated at the previous time point. In addition to a model fit,
we thus need to pass an initial vector of counts y.start to the simulate-
method. As an example, we simulate 100 realizations of the evolution of
measles during the year 2002 based on the fitted random effects model and the
counts of the last week of the year 2001 in the 17 districts:

R> (y.start <\- observed(measlesWeserEms)[52,])

03401 03402 03403 03404 03405 03451 03452 03453 03454 03455 03456 03457 03458
03459 0 0 0 0 0 0 0 0 0 0 0 25 0 003460 03461 03462 0 0 0

R> measlesSim <\- simulate(measlesFit_ri, nsim = 100, seed = 1, subset =
53:104,\+ y.start = y.start, simplify = TRUE) By simplify=TRUE, the simulated
counts are returned as a $52\times 17\times 100$ array instead of a list of
100 "sts" objects. We can, e.g., look at the final size distribution of the
simulations:

R> summary(colSums(measlesSim, dims = 2))

Min. 1st Qu. Median Mean 3rd Qu. Max. 223 326 424 550 582 3970 A few large
outbreaks have been simulated, but the mean size is below the observed number
of sum(observed(measlesWeserEms)[53:104,]) $=779$ cases in the year 2002. We
refer to help("simulate.hhh4") and Meyer and Held (2014b, Section 4.3) for
further examples.

##  6 Outlook

In the present work we have introduced the surveillance package as a
comprehensive framework for the analysis of spatio-temporal surveillance data
covering individual-level event data as well as aggregated count data time
series. Altogether, the package offers a multitude of methods for
visualization, likelihood inference and simulation of endemic-epidemic models.
Additional functionality beyond the illustrations used in Sections 3 to 5 can
be found in the package’s help pages. The extra vignette("hhh4") contains
examples of modeling univariate and bivariate count time series using hhh4. We
hope that the package can help to serve an increased need in analyzing spatio-
temporal epidemic data using statistical models. From a methodological point
of view, future work will focus on using the endemic-epidemic model classes as
a framework for model-based cluster detection in spatio-temporal data.

## Acknowledgements

The implementation of the hhh4 model is mainly due to Michaela Paul, to whom
we are thankful for all methodological advances and code contributions in the
past years. We also acknowledge all other code contributors in the long
history of the surveillance package (in alphabetical order): Thais Correa,
Mathias Hofmann, Christian Lang, Juliane Manitz, Andrea Riebler, Daniel
Sabanés Bové, Maëlle Salmon, Dirk Schumacher, Stefan Steiner, Mikko Virtanen,
Valentin Wimmer. Many have also helped us by investigating the package and
giving feedback: Doris Altmann, Johannes Dreesman, Johannes Elias, Marc
Geilhufe, Kurt Hornik, Mayeul Kauffmann, Marcos Prates, Brian D. Ripley, Barry
Rowlingson, Christopher W. Ryan, Klaus Stark, Yann Le Strat, André Michael
Toschke, Wei Wei, Achim Zeileis.

Financial support by the Munich Center of Health Sciences (2007–2010) and the
Swiss National Science Foundation (since 2007) is gratefully acknowledged.

## References

  * Anderson and May (1991) Anderson RM, May RM (1991).  _Infectious Diseases of Humans: Dynamics and Control_.  Oxford University Press.
  * Baddeley and Turner (2005) Baddeley A, Turner R (2005).  “spatstat: An R Package for Analyzing Spatial Point Patterns.”  _Journal of Statistical Software_ , 12(6), 1–42.  URL http://www.jstatsoft.org/v12/i06/.
  * Bernard _et al._ (2014) Bernard H, Werber D, Höhle M (2014).  “Estimating the Under-Reporting of Norovirus Illness in Germany Utilizing Enhanced Awareness of Diarrhoea During a Large Outbreak of Shiga Toxin-Producing E. Coli O104:H4 in 2011 – A Time Series Analysis.”  _BMC Infectious Diseases_ , 14(1), 116.  10.1186/1471-2334-14-116.
  * Besag _et al._ (1991) Besag J, York J, Mollié A (1991).  “Bayesian Image-Restoration, with Two Applications in Spatial Statistics.”  _Annals of the Institute of Statistical Mathematics_ , 43(1), 1–20.
  * Bivand (2014) Bivand R (2014).  _spdep: Spatial dependence: weighting schemes, statistics and models_.  R package version 0.5-77, URL http://CRAN.R-project.org/package=spdep.
  * Bivand _et al._ (2014) Bivand R, Keitt T, Rowlingson B (2014).  _rgdal: Bindings for the Geospatial Data Abstraction Library_.  R package version 0.9-1, URL http://CRAN.R-project.org/package=rgdal.
  * Bivand and Lewin-Koh (2014) Bivand R, Lewin-Koh N (2014).  _maptools: Tools for reading and handling spatial objects_.  R package version 0.8-30, URL http://CRAN.R-project.org/package=maptools.
  * Bivand and Rundel (2014) Bivand R, Rundel C (2014).  _rgeos: Interface to Geometry Engine - Open Source (GEOS)_.  R package version 0.3-8, URL http://CRAN.R-project.org/package=rgeos.
  * Bivand _et al._ (2013) Bivand RS, Pebesma E, Gómez-Rubio V (2013).  _Applied Spatial Data Analysis with R_ , volume 10 of _Use R!_ 2nd edition. Springer, New York.  ISBN 1-4614-7617-8.  URL http://www.asdar-book.org.
  * Chen _et al._ (2014) Chen CX, Wakefield JC, Rue H, Self S, Feng ZJ, Wang Y (2014).  “Bayesian Penalized Spline Models for the Analysis of Spatial-Temporal Count Data.”  Submitted to Statistics in Medicine.
  * Czado _et al._ (2009) Czado C, Gneiting T, Held L (2009).  “Predictive Model Assessment for Count Data.”  _Biometrics_ , 65(4), 1254–1261.  10.1111/j.1541-0420.2009.01191.x.
  * Dahl (2014) Dahl DB (2014).  _xtable: Export tables to LaTeX or HTML_.  R package version 1.7-4, URL http://CRAN.R-project.org/package=xtable.
  * Daley and Gani (1999) Daley DJ, Gani J (1999).  _Epidemic Modelling: An Introduction_.  Cambridge University Press.
  * Daley and Vere-Jones (2003) Daley DJ, Vere-Jones D (2003).  _An Introduction to the Theory of Point Processes_ , volume I: Elementary Theory and Methods of _Probability and its Applications_.  2nd edition. Springer-Verlag, New York.  ISBN 0-387-95541-0.
  * Diggle (2006) Diggle PJ (2006).  “Spatio-Temporal Point Processes, Partial Likelihood, Foot and Mouth Disease.”  _Statistical Methods in Medical Research_ , 15(4), 325–336.  10.1191/0962280206sm454oa.
  * Douglas and Peucker (1973) Douglas DH, Peucker TK (1973).  “Algorithms for the Reduction of the Number of Points Required to Represent a Digitized Line or its Caricature.”  _Cartographica: The International Journal for Geographic Information and Geovisualization_ , 10(2), 112–122.  10.3138/FM57-6770-U75U-7727.
  * Fahrmeir _et al._ (2013) Fahrmeir L, Kneib T, Lang S, Marx B (2013).  _Regression: Models, Methods and Applications_.  Springer-Verlag.  ISBN 978-3-642-34332-2.  10.1007/978-3-642-34333-9.
  * Fine and Clarkson (1982) Fine PEM, Clarkson JA (1982).  “Measles in England and Wales—I: An Analysis of Factors Underlying Seasonal Patterns.”  _International Journal of Epidemiology_ , 11(1), 5–14.  10.1093/ije/11.1.5.
  * Finkenstädt _et al._ (2002) Finkenstädt BF, Bjørnstad ON, Grenfell BT (2002).  “A Stochastic Model for Extinction and Recurrence of Epidemics: Estimation and Inference for Measles Outbreaks.”  _Biostatistics_ , 3(4), 493–510.  10.1093/biostatistics/3.4.493.
  * Finkenstädt and Grenfell (2000) Finkenstädt BF, Grenfell BT (2000).  “Time Series Modelling of Childhood Diseases: A Dynamical Systems Approach.”  _Journal of the Royal Statistical Society C_ , 49(2), 187–205.  10.1111/1467-9876.00187.
  * Geilhufe _et al._ (2014) Geilhufe M, Held L, Skrøvseth SO, Simonsen GS, Godtliebsen F (2014).  “Power Law Approximations of Movement Network Data for Modeling Infectious Disease Spread.”  _Biometrical Journal_ , 56(3), 363–382.  10.1002/bimj.201200262.
  * Gneiting and Katzfuss (2014) Gneiting T, Katzfuss M (2014).  “Probabilistic Forecasting.”  _Annual Review of Statistics and Its Application_ , 1(1), 125–151.  10.1146/annurev-statistics-062713-085831.
  * Greven and Kneib (2010) Greven S, Kneib T (2010).  “On the Behaviour of Marginal and Conditional AIC in Linear Mixed Models.”  _Biometrika_ , 97(4), 773–789.  10.1093/biomet/asq042.
  * Groendyke _et al._ (2012) Groendyke C, Welch D, Hunter DR (2012).  “A Network-Based Analysis of the 1861 Hagelloch Measles Data.”  _Biometrics_ , 68(3), 755–765.  10.1111/j.1541-0420.2012.01748.x.
  * Harrower and Bloch (2006) Harrower M, Bloch M (2006).  “Mapshaper.org: A Map Generalization Web Service.”  _IEEE Computer Graphics and Applications_ , 26(4), 22–27.  10.1109/MCG.2006.85.
  * Held _et al._ (2006) Held L, Hofmann M, Höhle M, Schmid V (2006).  “A Two-Component Model for Counts of Infectious Diseases.”  _Biostatistics_ , 7(3), 422–437.  10.1093/biostatistics/kxj016.
  * Held _et al._ (2005) Held L, Höhle M, Hofmann M (2005).  “A Statistical Framework for the Analysis of Multivariate Infectious Disease Surveillance Counts.”  _Statistical Modelling_ , 5(3), 187–199.  10.1191/1471082X05st098oa.
  * Held and Paul (2012) Held L, Paul M (2012).  “Modeling Seasonality in Space-Time Infectious Disease Surveillance Data.”  _Biometrical Journal_ , 54(6), 824–843.  10.1002/bimj.201200037.
  * Herzog _et al._ (2011) Herzog SA, Paul M, Held L (2011).  “Heterogeneity in Vaccination Coverage Explains the Size and Occurrence of Measles Epidemics in German Surveillance Data.”  _Epidemiology and Infection_ , 139(04), 505–515.  10.1017/S0950268810001664.
  * Hughes and King (2003) Hughes AW, King ML (2003).  “Model Selection Using AIC in the Presence of One-Sided Information.”  _Journal of Statistical Planning and Inference_ , 115(2), 397–411.  10.1016/S0378-3758(02)00159-3.
  * Höhle (2009) Höhle M (2009).  “Additive-Multiplicative Regression Models for Spatio-Temporal Epidemics.”  _Biometrical Journal_ , 51(6), 961–978.  10.1002/bimj.200900050.
  * Höhle and Mazick (2010) Höhle M, Mazick A (2010).  “Aberration Detection in R Illustrated by Danish Mortality Monitoring.”  In TA Kass-Hout, X Zhang (eds.), _Biosurveillance: Methods and Case Studies_ , pp. 215–238. Chapman & Hall/CRC.
  * Höhle _et al._ (2009) Höhle M, Paul M, Held L (2009).  “Statistical Approaches to the Monitoring and Surveillance of Infectious Diseases for Veterinary Public Health.”  _Preventive Veterinary Medicine_ , 91(1), 2–10.  10.1016/j.prevetmed.2009.05.017.
  * Johnson (2014) Johnson A (2014).  _polyclip: Polygon Clipping_.  R package version 1.3-0, ported to R by Adrian Baddeley and Brian Ripley, URL http://CRAN.R-project.org/package=polyclip.
  * Keeling and Rohani (2008) Keeling MJ, Rohani P (2008).  _Modeling Infectious Diseases in Humans and Animals_.  Princeton University Press.  URL http://www.modelinginfectiousdiseases.org/.
  * Keeling _et al._ (2001) Keeling MJ, Woolhouse MEJ, Shaw DJ, Matthews L, Chase-Topping M, Haydon DT, Cornell SJ, Kappey J, Wilesmith J, Grenfell BT (2001).  “Dynamics of the 2001 UK Foot and Mouth Epidemic: Stochastic Dispersal in a Heterogeneous Landscape.”  _Science_ , 294(5543), 813–817.  10.1126/science.1065973.
  * Kermack and McKendrick (1927) Kermack WO, McKendrick AG (1927).  “A Contribution to the Mathematical Theory of Epidemics.”  _Proceedings of the Royal Society of London A_ , 115(772), 700–721.  10.1098/rspa.1927.0118.
  * Lawson and Leimich (2000) Lawson AB, Leimich P (2000).  “Approaches to the Space-Time Modelling of Infectious Disease Behaviour.”  _IMA Journal of Mathematics Applied in Medicine and Biology_ , 17(1), 1–13.
  * Mandel (2013) Mandel M (2013).  “Simulation-Based Confidence Intervals for Functions with Complicated Derivatives.”  _The American Statistician_ , 67(2), 76–81.  10.1080/00031305.2013.783880.
  * Meyer (2014) Meyer S (2014).  _polyCub: Cubature over Polygonal Domains_.  R package version 0.5-1, URL http://CRAN.R-project.org/package=polyCub.
  * Meyer _et al._ (2012) Meyer S, Elias J, Höhle M (2012).  “A Space-Time Conditional Intensity Model for Invasive Meningococcal Disease Occurrence.”  _Biometrics_ , 68(2), 607–616.  10.1111/j.1541-0420.2011.01684.x.
  * Meyer and Held (2014a) Meyer S, Held L (2014a).  “Flexible Estimation of Spatio-Temporal Interaction in a Point Process Model for Infectious Disease Spread.”  In T Kneib, F Sobotka, J Fahrenholz, H Irmer (eds.), _Proceedings of the 29th International Workshop on Statistical Modelling_ , volume 1. Statistical Modelling Society, Göttingen, Germany.  10.5167/uzh-99895.
  * Meyer and Held (2014b) Meyer S, Held L (2014b).  “Power-Law Models for Infectious Disease Spread.”  _The Annals of Applied Statistics_ , 8(3), 1612–1639.  10.1214/14-AOAS743.
  * Meyer and Held (2014c) Meyer S, Held L (2014c).  “Supplement B of ‘Power-Law Models for Infectious Disease Spread’.”  10.1214/14-AOAS743SUPPB.  http://www.biostat.uzh.ch/research/manuscripts/powerlaw.html.
  * Neal and Roberts (2004) Neal PJ, Roberts GO (2004).  “Statistical Inference and Model Selection for the 1861 Hagelloch Measles Epidemic.”  _Biostatistics_ , 5(2), 249–261.  10.1093/biostatistics/5.2.249.
  * Ogata (1988) Ogata Y (1988).  “Statistical Models for Earthquake Occurrences and Residual Analysis for Point Processes.”  _Journal of the American Statistical Association_ , 83(401), 9–27.  URL http://www.jstor.org/stable/2288914.
  * Paul and Held (2011) Paul M, Held L (2011).  “Predictive Assessment of a Non-Linear Random Effects Model for Multivariate Time Series of Infectious Disease Counts.”  _Statistics in Medicine_ , 30(10), 1118–1136.  10.1002/sim.4177.
  * Paul _et al._ (2008) Paul M, Held L, Toschke A (2008).  “Multivariate Modelling of Infectious Disease Surveillance Data.”  _Statistics in Medicine_ , 27(29), 6250–6267.  10.1002/sim.3440.
  * Pebesma (2012) Pebesma E (2012).  “spacetime: Spatio-Temporal Data in R.”  _Journal of Statistical Software_ , 51(7), 1–30.  URL http://www.jstatsoft.org/v51/i07/.
  * Pebesma (2014) Pebesma E (2014).  “CRAN Task View: Handling and Analyzing Spatio-Temporal Data.”  http://CRAN.R-project.org/web/views/SpatioTemporal.html.  Version 2014-01-12.
  * R Core Team (2014) R Core Team (2014).  _R: A Language and Environment for Statistical Computing_.  R Foundation for Statistical Computing, Vienna, Austria.  URL http://www.R-project.org/.
  * Ryan and Ulrich (2014) Ryan JA, Ulrich JM (2014).  _xts: eXtensible Time Series_.  R package version 0.9-7, URL http://CRAN.R-project.org/package=xts.
  * Salmon _et al._ (2014) Salmon M, Schumacher D, Höhle M (2014).  “Monitoring Count Time Series in R: Aberration Detection in Public Health Surveillance.”  Submitted to the _Journal of Statistical Software_.
  * Sarkar (2008) Sarkar D (2008).  _Lattice: Multivariate Data Visualization with R_.  Springer, New York.  ISBN 978-0-387-75968-5, URL http://lmdvr.r-forge.r-project.org.
  * Scheike and Zhang (2011) Scheike TH, Zhang MJ (2011).  “Analyzing Competing Risk Data Using the R timereg package.”  _Journal of Statistical Software_ , 38(2), 1–15.  URL http://www.jstatsoft.org/v38/i02.
  * Schrödle _et al._ (2012) Schrödle B, Held L, Rue H (2012).  “Assessing the Impact of a Movement Network on the Spatiotemporal Spread of Infectious Diseases.”  _Biometrics_ , 68(3), 736–744.  10.1111/j.1541-0420.2011.01717.x.
  * Therneau (2014) Therneau TM (2014).  _A Package for Survival Analysis in S_.  R package version 2.37-7, URL http://CRAN.R-project.org/package=survival.
  * Venables and Ripley (2002) Venables WN, Ripley BD (2002).  _Modern Applied Statistics with S_.  Fourth edition. Springer, New York.  ISBN 0-387-95457-0, URL http://www.stats.ox.ac.uk/pub/MASS4.
  * Wei and Held (2014) Wei W, Held L (2014).  “Calibration Tests for Count Data.”  _TEST_.  10.1007/s11749-014-0380-8.  In press.
  * Xie (2013a) Xie Y (2013a).  “animation: An R Package for Creating Animations and Demonstrating Statistical Methods.”  _Journal of Statistical Software_ , 53(1), 1–27.  URL http://www.jstatsoft.org/v53/i01/.
  * Xie (2013b) Xie Y (2013b).  _Dynamic Documents with R and knitr_.  The R Series. Chapman and Hall/CRC, Boca Raton, Florida.  ISBN 978-1482203530.  URL http://yihui.name/knitr/.