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

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.


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, Schumacher, and Höhle 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, Held, and Toschke 2008;Paul and Held 2011;Geilhufe, Held, Skrøvseth, Simonsen, and Godtliebsen 2014), meningococcal disease (Paul et al. 2008;Paul and Held 2011;Meyer, Elias, and Höhle 2012), measles (Herzog, Paul, and Held 2011), hand, foot and mouth disease (Chen, Wakefield, Rue, Self, Feng, and Wang 2014), rabies in foxes (Höhle, Paul, and Held 2009), coxiellosis in cows (Schrödle, Held, and Rue 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, Pebesma, and Gómez-Rubio 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 1 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 implementations 2 -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.

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, Höhle, and Hofmann 2005;Held, Hofmann, Höhle, and Schmid 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.
Infective events occur at specific points in continuous time and space, which gives rise to time-space-mark data {(t i , s i , m i ) : i = 1, . . . , n} from a region W observed during a period (0, T ]. Event times and locations can be regarded as a realization of a "self-exciting" spatiotemporal 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) It represents the instantaneous event rate at location s at time point t given all past events, and is often more verbosely denoted by λ * or by explicit conditioning on the "history" H t . Daley and Vere-Jones (2003, Chapter 7) provide a rigorous mathematical definition of this concept. Both the endemic offset ρ [s] [t] and the log-linear predictor ν [s] [t] of endemic covariates are assumed to be piecewise constant on a suitable spatio-temporal grid (e.g., week × district), hence the notation [s] [t] for the region covering 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, where ε j , δ j > 0 are pre-specified temporal and spatial ranges of interaction, possibly infinite. The epidemic log-linear predictor η j of event-specific marks 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+σ) −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 λ(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 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 for infection of a susceptible individual i ∈ S(t). The rate decomposes into a Cox proportional hazards formulation with a log-linear predictor ν i (t) = exp z i (t) β 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 which is represented by a linear combination of basis functions B m , e.g., B-splines (Fahrmeir, Kneib, Lang, and Marx 2013, Section 8.1). The distance-based force of infection is modified additively by a linear predictor of covariates 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 x i (t) α by interchanging the summation order to such that x i (t) comprises all epidemic terms summed over j ∈ I(t). Note that the parameter vector α is constrained to be non-negative to ensure that λ i (t) ≥ 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 and region-specific overdispersion parameters ψ i such that the conditional variance of Y it is µ it (1 + ψ i µ it ). This includes models assuming a common overdispersion ψ i ≡ ψ, as well as the Poisson distribution as the special case ψ i ≡ 0. The endemic component again involves a log-linear predictor (ν 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 ν it , λ it , and φ 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.

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" 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 s i must be defined in a planar coordinate reference system (CRS) to enable Euclidean geometry; see the spTransform-methods in package rgdal (Bivand, Keitt, and Rowlingson 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 ε i (eps.t), the spatial influence radius δ 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 m i of the events, e.g., the age of infected individuals, which may be used as covariates in the epidemic predictor η j influencing infectiousness.
W: a "SpatialPolygons" object representing the observation region 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 precomputes auxiliary variables for model fitting. This includes computing the centered individual influence regions R i = (W ∩b(s i , δ i )) − s i , obtained by intersecting the observation region W with a polygonal approximation of a disc of radius δ i centered at s i . These polygons R i are computed by the function intersectPolyCircle, which wraps functionality of package polyclip (Johnson 2014). Since the 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: 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 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 ρ [s] [t] in the endemic component of Equation 1.

Methods for data handling and visualization
Textual summaries of the data can be obtained by the print and summary methods. Printing "epidataCS" shows the first n (defaults to 6) events: 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).

R> 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 = "") 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")))    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 animatemethod. 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 powerlaw 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.
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

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(ρ [s][t] ) + log(ν [s] [t] ) and log(η 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(η j ) for event-specific infectivity. For instance, epidemic =~log(popdensity) + type corresponds to 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.

surveillance: Spatio-Temporal Analysis of Epidemic Phenomena
The log-likelihood to evaluate is a function of all parameters in the predictors ν [s] [t] and η j and in the interaction functions f and g. Since it has the form it contains a clear computational bottleneck: the two-dimensional integrals R i f ( s ) ds over the polygonal influence regions 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.

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 ρ [s][t] , a (centered) time trend, and a sinusoidal wave of frequency 2π/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 ([s], [t]). Hence the log-likelihood (6) of an endemic-only twinstim simplifies to a sum over all these cells,  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: Coefficients 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) = 1 (0,30] (t), with a Gaussian distance-decay f (x) = exp −x 2 /(2σ 2 ) . 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 ( s ) (and ∂f ( s ) ∂ log σ in the score function) over all domains R i of Equation 6. On Unixalikes, 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.  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 L A T E X 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 L A T E X 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.

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 η j · g(t − t j ) · f ( s − s j ) around its origin. Since this process is independent of the individuals parentage and of other infections, 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, and exponential decay of infectivity (g(t) = e −αt , ε j ≡ ∞). Then, for T → ∞, which corresponds to the basic reproduction number known from the simple SIR model by interpreting |W | as the population size, β as the transmission rate and α as the removal rate. Like in classic epidemic models, the process is sub-critical if µ < 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: Confidence intervals for the estimated reproduction numbersμ 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"). 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 powerlaw 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 + σ) −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)))

Interaction functions
where we no longer apply the upper bound eps.s=200 km for spatial interaction. The powerlaw kernel puts more weight on local transmission and exhibits a heavier tail enabling occasional interaction between events at far distances. 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.

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 AIC imdepifit_endemic 4 19166 imdepifit_Gaussian 9 18967 imdepifit_powerlaw 10 18940 imdepifit_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.: In this example, none of the endemic predictors is removed from the model.

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λ g (t) = 2 k=1 Wλ (s, t, k) ds, 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.   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)
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.

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 (s) = f ( s ), 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 R0method 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> 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").

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 λ(s, t) then becomes λ 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 welldocumented 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.

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: 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 → I event at time tI) one day before the appearance of prodromes, and is removed from the epidemic (I → 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 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"), 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 j∈I(t) B m ( s i − s j ) shall be generated. Here, household (x i,H (t)) and nothousehold (x i,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: 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 λ i (t) of Equation (2) 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:

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(β 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, Welch, and Hunter (2012). Altogether, the CIF for a child i is modeled as where Y i (t) = 1(i ∈ 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 α c1 and α c2 on the force of infection. The model is estimated by maximum likelihood as described in Höhle (2009)   Note that Wald confidence intervals for the epidemic parameters α 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 α c1 and α 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:  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 11a by: R> plot(hagellochFit, which = "epidemic proportion", xlab = "time [days]") 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 11b reveals some deficiencies of the model in describing the waiting times between events, which might be related to the assumption of fixed infection periods.

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

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@RKI 3 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, . . . , 17, t = 1, . . . , 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.

Data structure: "sts"
The epidemic modeling of multivariate count time-series essentially involves two data matrices: the nEpochs×nUnits matrix of observed counts, and a nUnits×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): 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 / 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: 4 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)) 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)
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)

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 districtspecific intercepts α i in the endemic component, we first assume a common intercept α 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.
To account for temporal variation of disease incidence, the endemic log-linear predictor ν t incorporates an overall trend and a sinusoidal wave of frequency ω = 2π/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, λ = exp(α (λ) ) and φ = exp(α (φ) ), 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 = 1(j ∼ i) for the time being, which means that the epidemic can only arrive from directly adjacent districts.

surveillance: Spatio-Temporal Analysis of Epidemic Phenomena
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 ν t , log λ and log φ 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 α (λ) , 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 ψ for all districts. The alternatives are "Poisson" (ψ ≡ 0), and "NegBinM" (district-specific ψ 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) The idx2Exp=1:4 argument yields the estimates for λ, φ, α (ν) and exp(β 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 γ and δ of the sine-cosine terms to the amplitude A and phase shift ϕ of the corresponding sinusoidal wave A sin(ωt + ϕ) in log ν t (Paul et al. 2008). The multiplicative effect of seasonality on ν t is shown in Figure 14  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 ψ = 0. We can underpin this finding by comparing AIC values: R> AIC(measlesFit_basic, update(measlesFit_basic, family = "Poisson")) df AIC measlesFit_basic 7 1957 update(measlesFit_basic, family = "Poisson") 6 2479 The epidemic potential of the process as determined by the parameters λ and φ is best investigated by a combined measure: the dominant eigenvalue (maxEV) of the matrix Λ which has the entries (Λ) ii = λ on the diagonal and (Λ) ij = φw ji for j = 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) 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.

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. 5 This time-constant covariate needs to be transformed to the common matrix structure for incorporation in hhh4: 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  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., β s > 0. Alternatively, this effect could be assumed as an offset, i.e., β s ≡ 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: 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"] The estimated exponentβ 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β s = 3.29 (95% CI: 2.23-4.86).

Spatial interaction
Up to now, the model assumed that the epidemic can only arrive from directly adjacent districts because w ji = 1(j ∼ i), and that all districts have the same coefficient φ 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 commuterdriven spread in our model, we scale the district's susceptibility according to its population fraction by multiplying φ by e βpop i : 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β 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 −d ji , for j = i and w jj = 0, where the decay parameter d is to be estimated. Normalization to w ji / 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 isd = 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 · 1(o ji = 1) + e ω 2 · 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 AIC measlesFit_nepop 9 1887 measlesFit_powerlaw 10 1882 measlesFit_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 16b shows both the power law o −d and the second-order weight eω 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 16a for the power-law model: q q q q q q q q q q q q q q q q q q q q q q qq q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q qq q q q q q q q q q q q q q q q q q q q q q q q q q qq q q q q q q q qq q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q qq q q q q q q q q q q q q q q q q q q q q q q qq q q q q q q q q q q q q q q q q q q qq q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q qq q q q q q q q q q q q q q q q q q q q qq q q q q q q q q q q q q q q q q (a) Normalized weights in the power-law model.  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 k =j o −d jk , 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). 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 underreporting (Bernard, Werber, and Höhle 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 α
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 ( 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))) + } For the autoregressive component in Figure 17a, 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 λ i = exp(α (λ) + α (λ) i ) = 0.087, but it seems to import more cases from other districts than explained by its population (Figure 17b). In Figure 17c, 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 AICbased 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: 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.
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, Gneiting, and Held 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 − µ 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 ses measlesFit_basic 1.09 0.736 5.29 measlesFit_powerlaw 1.10 0.731 5.39 measlesFit_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 ses measlesFit_basic 1.10 0.748 5.40 measlesFit_powerlaw 1.14 0.765 5.87 measlesFit_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 ( 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 ∼ 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: Note that the parameters mu and size (= 1/ψ) 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/ψ), 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 12a). Formal statistical tests for the calibration of count forecasts have been recently proposed by Wei and Held (2014).

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 ω, S s=1 {γ s sin(sωt) + δ s cos(sωt)} , to the autoregressive and/or spatio-temporal log-linear predictors -just like for log ν 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 12a), 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, Bjørnstad, and Grenfell (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.

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 µ 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. By simplify=TRUE, the simulated counts are returned as a 52 × 17 × 100 array instead of a list of 100 "sts" objects. We can, e.g., look at the final size distribution of the simulations:

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.