Ethics statement

The data used in this study were supplied from the UK Health Security Agency (UKHSA) under strict data protection protocols agreed between the University of Warwick and UKHSA. The ethics of the use of these data for these purposes was agreed by UKHSA.

Null model and outbreak detection

A spatio-temporal disease pattern can often be decomposed into two components, endemic and epidemic—although such a distinction can be blurred [34]. The endemic component explains a baseline rate of cases with a persistent yet stable pattern, only reinforced by immigration [35]. The epidemic element can be thought of as realizations of processes with auto-regressive or self-exciting behavior, which arise from a common source, e.g., environmental exposure, increased susceptibility in the population, or invasion of a new pathogen. In other words, the presence of epidemic events is correlated with and reinforces the occurrence of successive events [34].

We argue that the most important characteristic common to all epidemic data is that observations are not independent. Therefore, to state whether cases are observed in atypical frequency, we must compare our observations with those expected if the cases were uncorrelated in both space and time. If the recorded case observations are consistent with such an uncorrelated baseline, no warning should be raised. This suggests that the first step of an outbreak detection system is to develop a naïve model for the observed cases—a null model—which does not consider any spatial or temporal correlations. The second step is then to detect (statistically significant) deviations of the observed pattern from that produced by the null model.

We let X(A) denote the random variable that describes the number of cases detected in a region of space and time () under the null model (X) and define a threshold value θ_A such that (1)

A basic exceedance criterion therefore consists of triggering a warning flag if the number of cases observed over A exceeds a value θ_A. Here, α is the false-positive detection rate, which is the probability of triggering a warning by chance when no outbreak has occurred, and the set represents the entire extent of the spatio-temporal data. The approach described so far is not new [5,36] and can be regarded as the standard in statistical process control [37].

The random neighbourhood covering approach

We assume that the null model is an inhomogeneous Poisson point process with intensity λ(x,y,t), where x and y are spatial coordinates and t is the time. We stress that the rate parameter λ is not necessarily constant in space or time. Indeed, even in the absence of localized outbreaks, we may expect many episodes in highly populated regions and at certain times of the year, as many diseases are known to be influenced by seasonality while weekends or holidays can impact the daily recording practice [10,13] [38]. A Poisson point process across space and time () is defined by the following two properties:

1. X(A₁) and X(A₂) are independent for disjoint subsets A₁ and A₂ of ;
2. For any subset , X(A) has Poisson distribution

(2) where ξ = (x,y,t) and Λ(A) is referred to as the intensity measure [39]. Given the intensity function, it is straightforward to numerically obtain Λ(A), find the αth quantile θ, and state whether the observed number n exceeds this value in the subset A. Those subsets such that the observed event number exceeds θ are tagged as outbreak subsets. For simplicity, it is convenient to treat the time t as an integer variable rounding it up to, e.g., days or weeks, thus is ℝ²×ℕ, its generic subset A is the union ∪_i∈ℕ A_i, A_i = {(x,y,t)∈A|t = i}, and the intensity measure can be easily computed as sum of two-dimensional integrals, i.e., .

The intensity function λ represents the average trend of the endemic infections. We assume that it factorises into a spatial factor λ₀ and a temporal factor λ_time, i.e., λ(x,y,t) = λ₀(x,y) λ_time(t), thus excluding any correlation between the space and time trends. We make the simplifying assumption that λ₀(x,y) is proportional to the population density, reflecting the fact that the number of observed episodes in an area is likely to increase with its population size–although λ₀(x,y) could be matched to historical infection data if these were plentiful. This population-based weighting is arguably the most important factor in the spatial distribution of the episodes, although, across geographical areas, other factors such as the climatic zone or disparities in access to health care could affect the baseline. The temporal factor encodes the disease seasonality, and other long-term trends, and is generally obtained by fitting a periodic function to available data.

The inequality (1) provides a basic exceedance criterion for flagging a subset A (which encompasses a geographical region for a certain time interval) but is far from being ideal. Indeed, a warning is raised with probability α even if the data is actually generated by the baseline model X(A). Further, the aberration can only be tested on the subset A and outbreaks diffused over a larger area are likely to be missed [40]. To address these issues we consider, for each reported case (of given space-time coordinates ), N overlapping subsets , i = 1,2,⋯,N, each of which contains the case (x∈A_i) and apply the criterion of (1) to each one of them in order to test whether there is aberration from the baseline Poisson point process (Fig 1). In the absence of any aberrations and in agreement with property 2, the number of detected cases in any subset A_i has distribution (2). If its intensity measures Λ(A_i) are the same for all the subsets, we expect just a fraction α of them to exceed the threshold (while in the presence of epidemic events the fraction is greater than α). Throughout this manuscript, we stick to the standard false discovery rate α = 0.05 for each single subset. It is worth noting that, by considering the state of multiple subsets, we increase robustness to random fluctuations and make highly unlikely to identify a cluster from random case data. In the null model, we indeed expect just 5% of independent subset to have a warning flag and even with just 20 subsets, the probability of generating a false flag is very small (probability of 19 out of 20 subsets having a flag is 3.6×10⁻²⁴).

For a detected case x, we refer to the proportion w(x) of flagged subsets A∋x as its warning score. By taking this approach, we are able to go beyond the binary warning flag for an outbreak event, and hence present w(x)∈[0, 1] as a continuous warning score specific to a detected case, high values of w corresponding to a high probability that an outbreak is occurring. This suggests raising a warning if w(x) is significantly greater than a chosen threshold, where the statistical significance is tested by means of a binomial test for proportions. Wilson intervals [41] are used as 2.5%-97.5% confidence intervals (CIs) for the warning scores. As N increases, these CIs tighten and the test significancy improves thus allowing practitioners to tune the parameter N and reach their desired significance level.

It is convenient to focus on subsets placed over detected cases, thus testing for aberration only where necessary (Fig 1). This obviously implies that every subset always contains at least one point event and its null-model prediction is shifted to 1+X(A). In space-time, we opt to select cylinders of the same volume ∫_Adξ. In order to allow cylinders of the same volume to enclose outbreaks highly diffused in space or lasting long in time, we randomly generate cylinders of variable radii ϱ and heights h, still subject to the constraint that their volume πϱ²h is constant. The algorithm appears to be robust to the choice of this volume (see simulation experiments); its value is heuristically set to avoid excessively small volumes, expected to contain only one event, and excessively large volumes, which would contain more cases but might not represent well local information.

It is worth noting that a cylinder can accommodate many point events and nearby events are likely to be enclosed by the same cylinder. Consequently, the warning scores are correlated and their values reflect the neighbourhood structure of the detected cases, as demonstrated in the case study of section Results. Also, in simulation experiments performed for validation, clusters are detected as groups of cases with high warning scores. To summarise, for a single infection case x, RaNCover returns a warning score w(x) that measures to what extent its presence alters the local count statistics from those of a Poisson point process of intensity λ. An appropriate baseline intensity needs to be estimated or given as an input.

We now illustrate the method in action by means of synthetic data and baseline. To this end, we assume the baseline intensity λ = λ₀(x,y)λ_time(t) with (3) (4) over a unit square region with parameters (A,B,C,D,E) = (7,0.9,0.5,20,80) where x and y are spatial coordinates, while t is the time index ranging from 1 to 20. This describes a spatial distribution of cases peaked at the centre of the region, modulated by a periodic temporal trend. For each value of t we simulate a spatial Poisson point process with intensity λ, which yields a random pattern of points representing endemic baseline episodes. An outbreak in a delimited polygonal region is simulated by randomly generating additional Poisson point events with intensity 50 λ_time(t) for t = 5,6,7,8,15,16,17,18 and zero elsewhere. Spatial Poisson processes are generated using the thinning method [42]. The resulting point pattern is illustrated in Fig 2A for the first 9 time steps. To identify the events corresponding to the simulated outbreak in this point pattern we apply RaNCover as follows. We draw N = 10,000 random covering cylinders of equal volume, with radii ρ≅(0.10,0.12,0.14,0.17,0.23) and height h = (1,2,3,4,5). The chosen volume is πρ²h≅0.17, which corresponds to about 5 expected baseline events per cylinder. The centre coordinates are randomly sampled with replacement from the episode locations and times, and randomly shifted whilst ensuring that each cylinder contains its assigned episode location. We then flag all the cylinders according to the basic exceedance criterion with α = 0.05. The classification performed by the algorithm is summarised by the warning scores, computed using this set of cylinders. These warning scores indeed appear to capture the synthetic outbreak, as illustrated in Fig 2B.

Performance metrics

Excellent performance is confirmed by the receiver operator characteristic (ROC) curve (panel A in S1 Fig) and the area under the curve (ROC-AUC 0.90, 95%CI 0.85–0.95 by means of bootstrap). The warning scores obtained upon increasing the volume of the cylinders by 40% achieve close performances (ROC-AUC 0.92, 95%CI 0.87–0.96) and are strongly correlated with the first replicate output (Pearson ρ = 0.95, P = 1×10⁻⁸³), thus suggesting that the result is robust to the choice of the cylinder size (panels B and C in S1 Fig).

In a more realistic surveillance scenario, the warning scores are prospectively updated as soon as new cases are detected and included into the dataset, which we take to occur at every point in time. Observing the individual warning scores varying as new observations are included allows us to assess timeliness. This measures the time elapsed between the case detection day and the day at which its warning score stabilises around its final value (corresponding to the retrospectively computed warning score). To simulate this, we iteratively apply the procedure detailed above to all the observations occurred until time τ, τ = 1,2,…,20 and observe how the warning scores of cases detected previous to the update time τ vary with τ, while new cases are included and assessed. Fig 3 shows this process in action, the three panels show the warning scores for outbreak and endemic cases (as defined by the Poisson intensity function) that are detected at times t = 5, t = 6, and t = 7. For this simulation, the values of the warning scores appear to converge after only one or two updates, which suggests that the method can reach excellent timeliness.

Defining iGAS background rate and reporting delay

Our study cohort consists of 10,820 iGAS records comprising 465 different emm subtypes, with the three most prevalent types (viz., 1.0, 89.0, and 12.0) accounting for more than 38% of all records, while other emm types occur one order of magnitude less frequently (Fig 4). A total of 235 different emm subtypes were recorded only once in the cohort.

We estimated the endemic baseline for iGAS as follows. As streptococcal isolates of different emm types are thought to belong to different endemic populations, they cannot form a single outbreak even when they occur in the same geographical location and at the same time [31]. Multiple emm type infections are very uncommon, and our iGAS database does not report presence of patients infected with different emm types. To encode such mutual independence, the baseline intensity λ_M(x,y,t) specific to the emm type M is obtained by multiplying λ(x,y,t) by the relative frequency λ_m of M in the population, computed retrospectively until time τ. In other words, the null-model intensity function is assumed to obey (5)

The temporal and spatial factors (λ_time and λ₀, respectively) remain the same for all emm types. As the incidence of iGAS infection is subject to seasonality [29], λ_time(t) is chosen to be periodic with a period T of one year, viz., (6) where ϑ_τ≔(a,b,c,d) are parameters that must be retrospectively estimated from the data until time τ. To determine these parameters, we use the principle of maximum likelihood estimation (MLE). We denote by t₁,t₂,…,t_n the ordered times of all cases detected until τ, 0<t₁<t₂<⋯<t_n<τ, regardless of their emm types and geographical locations. The parameter τ identifies the time interval used to estimate the temporal trend of the baseline; for iGAS this has to be long enough to capture seasonality, with longer τ arguably producing better estimates. The likelihood function L(ϑ_τ) under the Poisson null model is given by

The value of ϑ_τ that maximises the likelihood was recursively obtained by maximising over a,c,d using the Nelder-Mead method with b held fixed and then finding the global maximum over b with the other parameters held fixed.

The episodes in the database are labelled by patient postcode location. For each postcode, the estimated population size is obtained from the 2011 Census in England and Wales by the Office for National Statistics (ONS), which allows us to tabulate the values of the spatial baseline intensity λ₀(x,y). The corresponding spatial coordinates are taken to be the postcode centroid locations, also retrieved from the ONS. This completes the definition of the endemic intensity function λ(x,y,t).

The available iGAS records also include the week at which the original sample was taken, which, on average, occurred a week earlier than the laboratory emm typing. The time elapsed between the sample date and the emm-type report is referred to as the reporting delay. We generalise RaNCover to include untyped cases by introducing the baseline intensity λ_U≔λ_u(t;τ)λ_time(t)λ₀(x,y) of untyped cases at each update time τ, where λ_u(t;τ) is the average fraction of cases detected at time t and emm typed after time τ. For a generic spatio-temporal subset A, we refer to the number of typed observations, the typed baseline prediction, the number of untyped observations, and the untyped baseline prediction as n_M(A), X_M(A), n_u(A), and X_u(A), respectively. The next step is to define a strategy based on the standard exceedance criterion (1). We flag A if n_M(A) exceeds the threshold value θ_M,A or n_u(A) exceeds the value θ_u,A, i.e., with thresholds defined by Prob(X_M(A)>θ_M,A) = α and Prob(X_u(A)>θ_u,A) = α, respectively. X_M(A) and X_u(A) are assumed to be Poisson point processes of intensities (1−λ_u(t;τ))λ_M(x,y,t) and λ_U(x,y,t;τ), respectively. The warning score of an untyped case x is the fraction w(x) of subsets A∋x flagged accordingly to this criterion.

Simulation study

Benchmarking surveillance systems benefits from realistic simulations of health data [43,44]. As described in Methods, we assume that, in the absence of outbreaks, the detected cases are a realization of a generative Poisson point process with factorised intensity function λ₀(x,y)λ_time(t). Accordingly, we simulated the point process with intensity defined by Eqs (2) and (3) and parameters as estimated in section Methods, thus yielding a synthetic baseline of 4915 endemic cases (distributed over a period of 100 weeks in England) that should not raise new public-health concerns. The temporal intensity function of Eq (6) fitted to all data in the study period is illustrated in Fig 5. We simulated an outbreak by including additional (epidemic) point events localised in 19 selected postcode areas (arbitrarily chosen around the town of St. Albans, north of Greater London) during a chosen time interval (weeks 40 to 59). In order to imitate the initial spread of infections and their subsequent control, the total number of epidemic events per week is initially low, increases and peaks around week 50, and finally decreases until the number of epidemic cases drops to zero (Fig 6A). Aggregating the cases over their geographic locations clearly shows the temporal trend as illustrated in panel A in S2 Fig for the endemic cases and in panel B in S2 Fig for both endemic and epidemic cases. Based on these time series, it is virtually impossible to detect excess cases during the outbreak. Analogously, aggregating over the time yields the geographical distributions of cases, which mainly occur alongside the most populated area (panel C in S2 Fig) and show no appreciable changes if the outbreak in St. Albans is included (panel D in S2 Fig). Yet, whether a single case is part of the outbreak or not is clearly labelled in the synthetic dataset, thus making it ideal for evaluation purposes; this typically is not the case in real data.

The set of these simulated data was used to test whether RaNCover can detect the epidemic components in it. RaNCover for retrospective analysis returns a warning score w(x) for each case x. We set α = 0.05, N = 1,000,000, and cylinder volumes ≅2.605 km²×week, which overall yielded tight CIs (S3 Fig). In Fig 6B, each warning score is plotted vs the case detection time t. In this figure, the markers are coloured in blue if the case was simulated from the baseline and red if simulated from the outbreak. It is possible to appreciate that the baseline events have low warning scores within statistical fluctuations. The bulk of the outbreak events always had warning scores close to one and can be unambiguously discriminated from the baseline based on w(x). Unsurprisingly, few cases at the tails of outbreak can be confused with background, but their warning scores reach relatively large values (w(x)>0.8). The ROC-AUC is 0.99 (95% CI 0.98–1.00), thus suggesting excellent performances. In Fig 6C, we show that the method can also identify the geographic location of the outbreak with precision. The result appears robust to the choice of the cylinder volume, as suggested by the fact that the warning scores calculated from a second replicate where cylinder volume was increased by 40% are strongly correlated with the first replicate (Fig 6D).

As a further test we sequentially calculated the warning scores prospectively, using only the data simulated until the week τ, τ = 40,41,⋯,59. As time progresses and new cases are detected, the warning scores are updated, thus allowing us to observe how quickly the outbreak is detected. The warning-score traces for the epidemic cases are plotted in Fig 6E, where it is possible to appreciate that, for the episode that occurred at t = 40 (which is the first outbreak event and appears isolated at first, see Fig 6A), the warning score is initially very low and increases when new outbreak cases occur (until it converges to a higher stable value). Similarly, the other warning scores become stable after few steps.

For comparison, we performed cluster detection with SaTScan [45] for retrospective space-time analysis using the discrete Poisson model [19]. SaTScan is a mature software widely used in disease surveillance that computes the scan statistics and is written in Java [45]. For our simulation experiment, the cylinder with the optimal likelihood ratio encloses 16 weeks (from 43 to 58) and 33 postcodes (Fig 6C, inset). These SaTScan results are widely consistent with ours; however, it is worth noting that SaTScan took around 3 hours to evaluate this sample dataset on an Intel Core i5 CPU (excluding default Monte Carlo hypothesis testing and not accounting for seasonality, which would further increase the computational burden), thus significantly longer than the 5 minutes circa required by our R [46] implementation of RaNCover. As SaTScan gives cylindrical boundaries of a detected cluster, these do not necessarily separate the outbreak cases from the baseline and, within these boundaries, there are cases correctly labelled only by RaNCover (Fig 6C).

Application to iGAS disease: Cluster analysis and data-point embedding

Warning scores for all episodes of selected emm types are retrospectively generated following the procedure detailed in Methods, with N = 20,000 cylinders per emm type with respect to the baseline intensity of Eq (5). The spatial and temporal factors are obtained from census data and MLE fit, respectively. For each emm type considered, the type-specific factor λ_m is the relative frequencies of cases of that emm type (as illustrated in Fig 4 at each update time). As the average baseline intensity governs the cylinders’ volumes, these are differentially set for each emm type, the rarer the emm type, the larger the volumes.

RaNCover consistently assigns high warning scores to neighbour cases and is thus able to reveal clusters which may belong to outbreaks. As each event is described by its spatial coordinates (e.g., latitude and longitude) and its date, the neighbourhood must be defined in both space and time. We use the t-distributed stochastic neighbour embedding (t-SNE) method to plot these cases (identified by three dimensions, i.e., latitude, longitude, and detection time) in a two-dimensional space in a way that preserves neighbouring identities [47]. In other words, t-SNE is a general statistical-learning method for dimensionality reduction which here is used to pin a two-dimensional point for each point in the three-dimensional space, in such a way that nearby events in space and time appear next to each other in the two-dimensional map. Using a visualisation method based on neighbourhoods rather than on absolute distances is very useful in our context, where the local prevalence of each emm type is important. Indeed, for emm types with low prevalence, recorded episodes that are distant in space and time can be correctly identified as being part of the same outbreak. Conversely, common emm types already established in the host population may show more clusters in the same geographical area and time scale. We illustrate the results for 4 rare emm types (94.0, 108.1, 44.0, and 33.0), and a common emm type that did not appear to cluster (12.0).

Application to iGAS disease: Prospective analysis and emm typing delay

The cluster analysis detailed in the previous subsection is performed retrospectively and is based on cases labelled by the week at which the reference laboratory identified the pathogen’s emm type. As anticipated for synthetic data (Figs 3 and 6E), it is also possible to observe single-case warning scores varying in time as new cases are emm typed and the outbreak progresses. As an example, we focus on the emm type 44.0, which appeared to have caused an outbreak in the East of England (Fig 9) with 23 cases retrospectively flagged at the end of the study period (circled red markers in Fig 9A and 9B). The warning scores computed at each week are illustrated in Fig 12A. As the first case of the cluster is detected, its score is high yet substantially below the threshold of 0.95. As more cases are detected, it approaches the threshold and sharply increases within the first four weeks (Fig 12B), during which nine cases are identified in the region.

It is possible to combine reporting delays in RaNCover. During the interval from the sampling to the emm typing date, we are aware that, at a precise location and time, a case of an emm type of concern might have occurred. Therefore, only incomplete information is available. Even so, it is desirable to account for this, since it may allow earlier outbreak detection. However, it is possible that an untyped case is not actually associated with the clusters of interest, thus including it into a surveillance system may substantially tamper the overall results. It is important to choose a suitable strategy, which we take to flag a subset if the untyped or typed observation exceeds their respective baselines. Prospectively, as the time progresses, new cases are observed and some previously observed cases are emm typed. In such a scenario, at a given point in time, in addition to the typed iGAS cases and their baselines (one for each emm type) one can collect the untyped cases and compute the corresponding baseline. As an example, the number of cases per week and the continuous baselines are plotted against time in Fig 12C for the emm type 44.0—similar snapshots can be obtained for any other point in time. The number of untyped cases is highest in correspondence of the update week and quickly decays to zero as the observation week goes back in time. Using these two baselines and according to the procedure detailed in Methods, we computed the warning scores relative to a potential emm 44.0 outbreak of the untyped cases at each update week. More specifically, we are interested in the early warning scores of those cases which are later found to belong to the 2019 cluster. The early warning scores did not reach the threshold of 0.95, yet they are significantly larger than the baseline warning scores thus potentially providing early warning signals as illustrated in Fig 12D.

Software and reproducibility

Analyses were carried out using R ([46], version 3.4.3). t-distributed stochastic neighbour embedding was performed using the package Rtsne [54] with default parameters. Maps created with Sf [55] using shapefiles from the GADM database (www.gadm.org) and the Ordnance Survey Data Hub (osdatahub.os.uk). All software (also including codes for generating of synthetic data) is archived online at https://github.com/mcavallaro/outbreak-detection.