PaperlayerView exact revision on arXiv

Source-linked AI summary

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

Sebastian Meyer, Leonhard Held, Michael Höhle

arXiv:1411.0416v2stat.COcs.CEphysics.data-anstat.AP

TL;DR

Spatio-temporal epidemic data are difficult to analyze because observations are dependent, observationally collected, and often only partially observable. The paper presents surveillance, an open-source R framework implementing regression-oriented models for multiple spatial and temporal data resolutions. It provides visualization, likelihood inference, and simulation while decomposing risk into endemic factors and epidemic dependence, with scope boundaries including age-contact structure and sensitivity to endemic-component misspecification.

  • Problem

    Epidemic surveillance data have self-exciting dependence, are rarely generated by planned experiments, and are often only partially observable, motivating specialized statistical tools.

  • Method

    The paper presents surveillance’s regression-oriented frameworks for continuous point patterns, discrete-unit event histories, and aggregated spatio-temporal counts, with visualization, inference, and simulation tools.

  • Results

    The package implements endemic-epidemic models that separate external endemic risk from dependence on observed past events and supports analyses of infectious-disease surveillance data.

  • Takeaways & Limitations

    surveillance provides an open-source, reproducible framework for analyzing and modeling spatio-temporal epidemic phenomena across individual-level and aggregated surveillance data.

  • Takeaways & Limitations

    The framework does not yet incorporate age-stratified contact structure, and misspecified endemic components can distort model-based reproduction numbers.

Abstract

from arXiv · show

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.

1. Introduction

Epidemic surveillance data challenge classical statistical analysis because observations are dependent, observationally collected, and often only partially observable. The paper presents surveillance as an open-source R framework for monitoring and regression-oriented modeling of spatio-temporal epidemic data.

  • Challenges: Epidemic data are spatio-temporal and self-exciting, with dependent observations, limited observability, and rarely planned experiments.These properties hinder the direct application of classical statistical approaches.
  • The surveillance package: Since 2005, surveillance has provided an open-source R framework for monitoring and modeling epidemic phenomena.The paper focuses on its regression-oriented modeling of spatio-temporal epidemic data.
  • The surveillance package: Regression models in surveillance support assessing environmental, socio-demographic, and control-related influences on endemic and epidemic dynamics.They also support analysis of spatio-temporal event interaction and simulation of epidemic spread from estimated models.
  • Applications: The package’s statistical modeling frameworks have been applied to surveillance data on diseases and health outcomes including influenza, meningococcal disease, measles, rabies, and psychiatric admissions.The cited applications span human, animal, and other surveillance settings.
  • Related software: Compared with related tools, surveillance is presented as uniquely able to estimate regression models for spatio-temporal epidemic data.Other packages address reproduction numbers, phylogenetic transmission, intervention strategies, or temporal epidemic models.
  • Paper organization: The paper reviews three implemented statistical models and illustrates their data handling, visualization, inference, and simulation through infectious-disease surveillance applications.The model-specific sections organize the paper’s methodological and practical examples.

2. Spatio-temporal endemic-epidemic modeling

The paper organizes spatio-temporal epidemic modeling around data resolution, covering continuous point patterns, discrete-unit event histories, and aggregated counts. Its models decompose disease risk into endemic factors and epidemic dependence on observed past events.

  • Model foundations: Epidemic modeling commonly uses stochastic SIR frameworks that represent individuals as susceptible, infectious, or removed.The stochastic model includes a transmission-rate parameter and an assumption about the duration of infectivity.
  • Data types: The paper distinguishes three spatio-temporal data types by spatial and temporal resolution for regression-oriented modeling in surveillance.These include continuous point patterns, discrete-unit event histories, and aggregated count data.
  • Data types: Continuous time-stamped, geo-referenced events are modeled as realizations of self-exciting spatio-temporal point processes.Discrete-unit event histories instead fit a spatial SIR framework represented as a multivariate temporal point process.
  • Endemic-epidemic decomposition: The three model classes are inspired by a Poisson branching process with immigration and additively separate endemic and epidemic risk.The endemic component captures external factors, while the epidemic component is driven by observed past events.
  • Endemic-epidemic decomposition: Endemic risk may vary with seasonality, population density, socio-demographic variables, and vaccination coverage across time or space.Explicit dependence between events is introduced through the epidemic component.
  • Implementation: Each model-specific section combines a theoretical introduction with implementation examples covering data handling, visualization, inference, and simulation.The examples use infectious-disease surveillance data described in the paper’s data overview.

3. Spatio-temporal point pattern of infective events

The surveillance package models individual-level spatio-temporal epidemic events using endemic and self-exciting epidemic components, with tools for inference, visualization, diagnostics, and simulation. The IMD example demonstrates type-specific transmission, spatial interaction assessment, and simulation-based source tracing, while highlighting limitations from tied observations and endemic immigration.

  • Model formulation: The twinstim model represents point-referenced surveillance data through endemic background risk and an observation-driven epidemic component that adds infection pressure from past events.The endemic component captures external factors, while past events trigger offspring during their infectious periods and within spatial interaction ranges.
  • Inference and interpretation: Model-based reproduction numbers are event-specific, depending on event marks, interaction ranges, location, and time, but are adjusted for endemic infections and should be interpreted as effective numbers.Misspecifying the endemic component can distort these estimates; under-reporting and control measures also affect their interpretation.
  • Model formulation: Type-specific twinstim models allow event types to have distinct transmission dynamics while sharing the same endemic pattern.For invasive meningococcal disease, the likelihood separates contributions across event types and uses a transmissibility structure between types.
  • Data preparation: The package prepares IMD case reports as epidataCS objects containing events, the observation region, and a spatio-temporal grid of areal covariates.The example contains 636 individual cases, and projected postcode-centroid coordinates enable Euclidean geometry.
  • Data limitations: Limited sampling resolution can create tied times or locations that conflict with continuous point-process assumptions, so random tie-breaking requires sensitivity analyses.Ties may distort temporal residual diagnostics and can cause power-law spatial kernels to diverge at zero distance.
  • Inference and interpretation: The IMD analysis estimates that serogroup C is less than half as infectious as serogroup B, while ages 3–18 years are estimated to cause twice as many secondary infections as ages 0–2 years.These subject-matter conclusions are derived from rate ratios for endemic and epidemic effects.
  • Model comparison and diagnostics: AIC favors the power-law model over Gaussian and endemic-only alternatives, while the more flexible step function achieves the best AIC but depends strongly on selected knots.The estimated endemic intensity also shows strong seasonality and a slow negative trend, and residual-process diagnostics assess uniformity and serial correlation.
  • Simulation: Simulations from the estimated power-law model retain event-source information, and 21.8% of events are triggered by previous events on average across 30 simulations.Source values distinguish endemic events, prehistory events, and events linked to an infective source.

4. SIR event history of a fixed population

The twinSIR model analyzes individual-level epidemic event histories in a fixed population, combining endemic infection risk with transmission from currently infectious individuals. The Hagelloch measles example illustrates data conversion, covariate construction, constrained inference, model checking, and model selection.

  • Model and data: twinSIR models individual-level surveillance data when the complete SIR event history of a fixed population is assumed known.The Hagelloch example contains sequential event times for 188 affected children and conditions on the index case being initially infectious.
  • Model and data: The infection rate decomposes into an endemic component for external sources and an epidemic component for transmission from currently infectious individuals.The epidemic force can depend on distances, basis functions, and covariates describing interactions between individuals.
  • Inference: The epidemic component gathers transmission terms into a single linear predictor, using non-negative covariates and parameters to preserve a non-negative conditional intensity.This additive construction differs from twinstim’s multiplicative covariate approach and supports constrained penalized likelihood inference.
  • Data preparation: The as.epidata converter standardizes event-history data, checks consistency, and can pre-calculate epidemic terms for twinSIR fitting.Inputs include event times, individual identifiers, coordinates, distance-based functions, interaction covariates, and retained marks.
  • Results: 8.15 times higher transmission was estimated between first- and second-class individuals, while infectious housemates exerted 32.3 times the pressure of infectious children outside the household.The estimated endemic background infection rate without current measles cases was 0.000635.
  • Model assessment: Residual diagnostics revealed deficiencies in the model’s waiting-time description, and the simpler model with only a nothousehold component was preferred.The final infection was necessarily attributed to the endemic component because no infectious children remained in the observed population.

5. Areal time series of counts

The hhh4 framework models infectious-disease counts aggregated by region and period through endemic and epidemic components. The Weser-Ems measles example demonstrates data construction, visualization, spatial coupling, and strong temporal and geographic heterogeneity.

  • Model: The hhh4 model describes infectious-disease counts across units and periods with a conditional negative binomial distribution, allowing Poisson variation as a limiting case.Its variance is µ_it(1+ψ_iµ_it), with shared overdispersion parameters also supported.
  • Model: The model’s mean decomposes into endemic and observation-driven epidemic components, with population-based offsets representing background risk.The endemic component can also account for infections imported from outside the study region or cases not directly linked to observed cases.
  • Spatial dependence: Spatial dependence is represented through transmission weights that describe infection flow between regions and may be informed by movement networks or estimated parametrically.A power-law distance decay based on adjacency order is one suitable choice, with normalized weights distributing prior cases across regions.
  • Inference: Likelihood inference uses analytical scores and Fisher information with a quasi-Newton optimizer, while models with random effects alternate penalized and marginal likelihood maximization.The hhh4 implementation generally converges quickly even with many parameters.
  • Data preparation: The Weser-Ems data are represented by observed-count, population, and neighborhood matrices within an sts object, alongside geographic map information.The example uses weekly counts beginning in 2001, district adjacency orders, and population fractions.
  • Results and visualization: Strong seasonality appeared across the two years, while the districts differed markedly in their reported measles burden.The visualization workflow includes overall and district-level time series, incidence maps, accumulated spatial plots, and animations.

5.3. Modeling and inference

The basic hhh4 model represents measles incidence through endemic, autoregressive, and neighbourhood components, with seasonal variation and population exposure incorporated into the predictors. Inference favors an overdispersed negative binomial model, and most fitted incidence is attributed to within-district autoregression.

  • Model specification: The initial hhh4 model sequentially extends a common-intercept negative binomial model with covariates, transmission weights, and random effects.The model is fitted by penalized maximum likelihood, and the extensions eventually account for unobserved district heterogeneity.
  • Model specification: The endemic predictor includes a trend and seasonal sine-cosine terms, while population fractions enter as multiplicative offsets and epidemic parameters are initially homogeneous.First-order adjacency defines the initial neighbourhood transmission weights.
  • Inference: 0.72 is the estimated epidemic dominant eigenvalue for the basic model.This measure combines the autoregressive and neighbourhood transmission parameters.
  • Inference: 1957 versus 2479 AIC favors the negative binomial basic model over its Poisson update.The estimated overdispersion and its confidence interval also support a negative binomial distribution with overdispersion.
  • Fitted components: The largest fitted-mean contribution comes from within-district autoregression, with little adjacent-district transmission and a small endemic contribution.The comparison is shown for the six districts with more than 20 cases using fitted component plots.

Covariates

The hhh4 framework incorporates regional and temporal covariates into endemic or epidemic components. In the measles analysis, model selection places susceptible proportion in the endemic predictor, where it substantially improves fit and has a positive estimated effect.

  • Covariate setup: Vaccination coverage is used to derive the susceptible proportion, an epidemiologically relevant covariate reflecting the remaining susceptible population.The data include complete vaccination information for all 17 districts in the first complete-data year.
  • Model selection: The analysis compares nine combinations of unchanged, offset, and covariate treatments for the endemic and autoregressive components using AIC.Each combination updates the initial measles model and is evaluated as an hhh4 fit.
  • Model selection: Adding susceptible proportion only to the autoregressive component increases AIC, whereas adding it to the endemic component yields a remarkable improvement.The selected model leaves the autoregressive component unchanged and adds βs log(1 − vi) to the endemic predictor.
  • Estimated effect: 3.29-fold incidence is estimated when a district’s susceptible fraction doubles, with a 95% CI of 2.23–4.86.The estimated exponent is positive and differs from the offset assumption.

Spatial interaction

The spatial-interaction extensions relax direct-adjacency and equal-importation assumptions by incorporating population attraction and distance-dependent transmission weights. Both higher-order transmission and population-based attraction improve model fit, while the power-law model is retained for subsequent extensions.

  • Population attraction: Population-scaled susceptibility addresses commuter-driven spread and attraction toward metropolitan areas.The extension multiplies the neighbourhood transmission parameter φ by a population-based factor.
  • Population attraction: β̂pop = 2.85 (95% CI: 1.83–3.87) and AIC decreases from 1917 to 1887 under the population-attraction model.The result provides strong evidence for an agglomeration effect.
  • Distance-dependent weights: Power-law weights model transmission as a function of adjacency order with normalization across each source district’s neighbourhood.The decay parameter d controls how interaction decreases for higher-order neighbours, and maxlag bounds the spatial interaction order.
  • Distance-dependent weights: 4.10 (95% CI: 2.03–6.17) is the estimated power-law decay parameter, indicating strong decay for higher-order neighbours.A second-order alternative estimates e^ω2 = 0.09 (95% CI: 0.02–0.39).
  • Model comparison: AIC improves when transmission from higher-order neighbours is modeled by either a power law or a second-order specification.The second-order model fits slightly better, but the power-law model is retained because the standalone second-order effect may not be identifiable in complex models and is scientifically implausible.

Random effects

Random effects allow district-specific incidence patterns not explained by observed covariates, including possible under-reporting. They reveal heterogeneous balances between local autoregression and imported transmission, but increase computational cost and remove simple AIC comparisons.

  • Model specification: Random effects are introduced across endemic, autoregressive, and neighbourhood intercepts to represent unobserved district heterogeneity.The framework also supports correlated random effects and conditional autoregressive spatial intercepts.
  • Inference: 0.84 is the epidemic dominant eigenvalue estimated for the random-effects model.The reported fit also includes penalized log-likelihood −869 and marginal log-likelihood −54.2.
  • Estimated heterogeneity: The random-effects maps show pronounced autoregressive heterogeneity, including three western districts associated with large local outbreaks and high case counts.Oldenburg has an estimated autoregressive factor λi = 0.087 and appears to import more cases than explained by population.
  • Estimated heterogeneity: The two districts without reported measles cases have relatively low adjusted endemic incidence and could indicate greater under-reporting.The adjustment accounts for population and susceptible proportion.
  • Limitations: Random effects increase estimation runtime from 0.1 to 4 seconds and invalidate simple AIC-based model comparisons.The authors compare fitted values visually with the initial fit instead.
  • Fitted components: Districts differ in whether fitted cases are mainly attributed to neighbouring transmission or local autoregression.Quantitative model-performance comparisons require more sophisticated techniques than the fitted decompositions.

Predictive model assessment

The paper evaluates competing measles models using proper scoring rules, calibration tests, and PIT histograms for one-week-ahead predictions. Rolling predictions favor the basic model overall, while calibration diagnostics find no evidence of severe miscalibration and suggest a seasonal timing shift.

  • Scoring and prediction: Proper scoring rules compare predictive distributions with later observed counts, including squared-error, logarithmic, ranked-probability, and Dawid–Sebastiani scores.The scores can be averaged across forecasts, and the log-score corresponds to a negative log-likelihood contribution.
  • Scoring and prediction: The assessment compares basic, power-law, and random-effects models on the second quarter of 2002 using one-week-ahead predictions.Rolling predictions refit each model through week t before predicting week t + 1.
  • Scoring and prediction: All scoring rules favor the random-effects model when fitted values from the complete time series are assessed during the test period.This comparison is described as essentially a goodness-of-fit assessment.
  • Scoring and prediction: The basic model gives the best overall mean scores for true rolling one-week-ahead predictions.Differences can be assessed with a paired permutation test or paired t-test.
  • Calibration: There is no clear predictive-performance difference between the basic and random-effects models during the test period.A calibration test for the random-effects model reports z = 0.80671, n = 221, and p-value = 0.4198, providing no evidence of miscalibrated predictions.
  • Calibration: PIT histograms are uniform under calibration, while the observed right-hand decay suggests all models predict higher counts than observed.The authors relate this pattern most likely to the epidemic peak shifting from the second quarter in 2001 to the first quarter in 2002.

Further modeling options

The paper discusses extensions for seasonal epidemic effects, observation intervals, and simulation-based forecasting. It cautions that more complex seasonal models may overfit short measles series and identifies age-structured contact patterns as an important future extension.

  • Seasonality: Seasonal variation can be added to epidemic and neighborhood components using a superposition of harmonic waves.The proposed terms parallel the sinusoidal seasonal structure used for the endemic component.
  • Seasonality: With only two years of measles surveillance and a shifted outbreak start, complex seasonal models are likely to overfit the data.The apparent seasonal shift is observed between 2001 and 2002.
  • Observation intervals: The measles models use a one-week epidemic lag, although measles generation time is approximately 10 days.Some studies therefore aggregate weekly surveillance data into biweekly intervals.
  • Age structure: A useful extension would stratify the epidemic model by age group to represent assortative contact patterns, especially among school children.The paper identifies incorporating this contact structure into hhh4 as current research.
  • Simulation: Simulation from fitted hhh4 models sequentially generates counts using fitted means and previously simulated counts.The procedure requires a fitted model and an initial count vector, and uses negative-binomial or Poisson draws.
  • Simulation: The simulated mean outbreak size is below the 779 observed cases in 2002, although a few large outbreaks are simulated.Figure 21 compares observed weekly cases with the long-term forecast.

6. Conclusion

The paper introduces surveillance as a comprehensive open-source framework for analyzing spatio-temporal surveillance data. It covers individual-level events and aggregated counts, with tools for visualization, likelihood inference, and simulation of endemic–epidemic models.

  • Conclusion: surveillance supports analysis of individual-level event data and aggregated count-data time series.The framework is designed for spatio-temporal surveillance data.
  • Conclusion: The package provides methods for visualization, likelihood inference, and simulation of endemic–epidemic models.The paper presents these capabilities through model-specific applications and examples.
  • Conclusion: Figure 21 presents a simulation-based long-term forecast aggregated over all districts, with weekly simulation means, pointwise 2.5% and 97.5% quantiles, and observed counts.The forecast starts from the last week of 2001.

Involved R packages and versions

The paper records the R package versions and dependencies used to implement its analyses. It identifies CRAN and the surveillance development site as access points for the software.

  • Software versions: The analyses use surveillance 1.10-0 with R version 3.2.2 and knitr for dynamic report generation.The paper also lists first-order dependency versions.
  • Software versions: The listed dependencies include sp, xtable, polyCub, MASS, and spdep.These packages are reported with the versions used in the paper.
  • Software availability: R, surveillance, and the other packages are available from CRAN, while surveillance development is hosted on R-Forge.The paper provides the corresponding project locations.
Loading 1411.0416v2…