Epidemiologists responsible for routine public health surveillance seek early indications of possible disease outbreaks. A growing collection of data types has become available to aid this surveillance, but with the analysis of these data comes the need to reduce the search possibilities. With cluster detection methodology, analysis of data sources can help by indicating the likely location and extent of a potential outbreak. Spatial and spatiotemporal scan statistics have the additional advantages of controlling the number of alarms resulting from multiple testing and accounting for usual clustering patterns, as found in the data [1].

Some natural clustering is expected, as in urban areas with high population density. If cluster detection methods do not account for this heterogeneity, customarily high concentrations will be signaled as anomalous. Only departures from the usual distribution should be flagged. In the implementation of the widely used SaTScan software [2], a detection statistic is applied to all candidate clusters, and the cluster with the highest statistical value is tested for significance. The statistic is a function of observed and expected values both inside and outside the candidate cluster. Thus, an accurate estimate of the expected spatial distribution is necessary for determination of relevant clusters with high sensitivity at reasonable false alarm rates. This requirement is consistent whether the detection approach uses scan statistics or another methodology.

A natural choice for the expected spatial distribution is the census population by sub-region of the geographic area under surveillance. In the United States, for example, population counts stratified by 3-digit zip code, 5-digit zip code, or census tract are freely downloadable [3]. This choice is appropriate for some applications of scan statistics to chronic disease surveillance where the data give true population rates relative to census totals and long-term effects are of interest. One problem with using these counts is that they represent the population measured at the most recent census, taken only every 10 years, so they do not include more recent population shifts. A more serious problem is that some biosurveillance data sources are not representative of the entire population, so they might produce expectations that are biased relative to the data of interest. This bias typically occurs because the data represent only one segment of the population that is not evenly dispersed over the region. Examples are data covering only military personnel and their dependents, as in this study, participants in a health maintenance organization or other medical plan, or customers of a pharmacy chain. Organizations furnishing such data typically do not or cannot furnish details of their spatial coverage. Therefore, system developers have used past data to form spatial expectations.

If substantial (typically multi-year) quality historical data are available, detailed modeling of counts at the sub-region level could be attempted to calculate expected distributed counts. Kleinman et al [4] obtained significantly reduced false alarm rates with this approach. However, health monitors typically need to apply data with little data history, and longer data history is likely to be unrepresentative of current spatial patterns because of changes in population behavior, in participation of data providers, and in data collection systems. This article focuses on data-based methods of computing expected count distributions assuming that only 1–8 weeks of recent data are available.

For clarity, we distinguish two separate problems in applying data with spatial information to detection of disease clusters. The first problem is how to aggregate the data into sub-regions. Ozonoff et al and Olson et al. [5,6] demonstrated that the level of aggregation can strongly affect the power of scan statistics for cluster detection. The second problem is that given the total number of cases occurring in a monitored region and the chosen partition of sub-regions, what is the expected distribution of cases among these sub-regions if no public health event of interest is underway? Reasonable agreement with this distribution might be seen as the null hypothesis that is rejected if a case cluster is found at a level of significance above a preset threshold. Practical data considerations often drive the answer to the first problem. Privacy laws typically preclude the use of patient addresses or other exact location fields, so zip codes, census tracts, and locations of clinics specified in data records have been used. The second problem is the subject of this article.

The choice of an effective method for estimating the case distribution is a function of the monitored data source. Even within the same data source, for example, estimation of the current distribution for a broad, seasonal syndromic filtering of clinical records is likely different from the same estimation for a rare sub-syndrome occurring evenly throughout the year. For the current study, we used data from BioSense, a US national biosurveillance system operated by the Centers for Disease Control and Prevention (CDC) [7]. This system uses computer databases for automated disease surveillance. The source of data for the study was the BioSense collection of US Department of Defense (DoD) military outpatient clinic data records. Since 2004, these data records have been sent daily to CDC through the TriCare Management Activity, which manages clinical patient data for all US military treatment facilities. The patient base includes personnel on active duty, retirees, and dependents. Patient record fields provide key information including the visit date, age, gender, residence and treatment clinic zip codes, and for medical classification, a code from the International Classification of Diseases, Ninth Revision (commonly called ICD-9), used by hospital and physicians' office personnel to report billing information to insurance companies [8].

The study dataset included 3 years (1/1/2004–12/31/2006) of DoD patient records from military outpatient clinics in the state of Texas, which has a large, widespread military population. These records were classified according to their ICD-9 code fields into 11 standard BioSense syndrome groups [9]. The respiratory and rash syndromes have been chosen for this study to represent high-count, seasonal data behavior and sparse, nonseasonal behavior, respectively. The record fields allowing spatial aggregation were the zip code of the clinic where the patient was seen and the zip code of patient residence. The patient residence zip codes were distributed nationwide, mostly outside of Texas. Figure 1 illustrates the number of records with residence zip codes inside and outside Texas for each DOD facility in Texas.

One obvious difference between the geographic patient distribution for the study dataset and the general census distribution results from the large fraction of non-Texas residents among these patients who live near the military bases, with the result that zip codes near those bases are overrepresented relative to census figures. However, for reasons of both patient privacy and military operational security, the spatial distribution of people eligible for treatment at base clinics is unavailable for analysis. Because of commuting by Texas residents, it cannot be assumed from Figure 1 that nearly all patient residences are close to treatment facilities. Many commuters are eligible for treatment at military clinics throughout the continental United States, and in western states such as Texas, commuting distances of over 50 miles are common. For these reasons, data-based spatial estimation methods are required. The Methods section discusses the data preparation and the estimation methods tested for cluster detection.

Using location fields in patient records to seek infectious disease clusters, the underlying assumption is that the field value approximates the location of exposure to the disease-causing agent. The location of exposure is difficult to determine in individual-based studies and much more so in automated, population-based systems. The records in the study dataset contain fields for zip codes of the treatment facility and of the patient's home. In the DoD population whose records compose this dataset, most of the active duty personnel work on military bases that include treatment facilities. Therefore, a working hypothesis is that the facility zip code should be used to detect an outbreak in which exposures occur in the work environment; whereas, the home zip code should be used to detect an outbreak resulting from residence-based exposures. However, the patient home zip codes in the study data often refer to the home town of the patient, not to the zip code of military duty at the time of treatment. In such cases, the facility zip code is a better choice for the current residence. There is no practical way to know which zip code to use for the residence-based analysis, so in view of the number of patients who commute, we have used the following "100-mile rule" on patient zip codes to prepare the data before the analysis.

▪ First, we included all records that have either a facility zip code or a home zip code inside Texas. We also included records from four facilities outside of Texas but near the Texas-Oklahoma and Texas-Louisiana borders because of the large numbers of Texas residents commuting across the border and using those facilities.

▪ Second, we calculated for each record the distance between the home and facility zip codes.

▪ For the residence-based analysis, we used the patient home zip code if this distance was less than 100 miles. Otherwise, we used the facility zip code.

▪ If the home zip code field was empty, the facility zip code was used as the residence zip code. If the facility zip code was missing or outside Texas, the home zip code was used.

We used this reasonably simple rule because a more detailed assignment would require information that was available neither for the study nor for routine monitoring. In a normal situation most people would probably not travel more than a couple of hours to a facility unless there were no other choices.

Restriction of the study data set to those records classified in syndrome groups yielded 1,233 residence zip codes and 32 facility zip codes. The residence zip codes were well distributed in the eastern half of the state with high concentrations around the cities, but much sparser to the west. Most of the 32 facility zip codes were to the east with a high concentration in the San Antonio area. The west was covered only by two facilities at El Paso in the western tip. The mean daily record counts per residence zip code were 0.29 and 2.67 for rash and respiratory syndromes, respectively. The corresponding mean daily record counts per facility zip code were 4 and 38. Thus, respiratory counts exceeded rash counts by a factor of >9. Although rash counts have a relatively smooth pattern with a slightly increasing trend during the 3-year period, respiratory counts have a strong seasonal pattern (Figure 2). Both syndromes have a strong day-of-week effect (Figure 3).

We used two types of estimation approaches to calculate expected counts: a baseline-mean approach and a space-time permutation approach. Figure 4 illustrates these two approaches. In each half of this figure, columns geo-1, geo-2 etc. are used to indicate subregions. Time is shown vertically, with a set of rows for each block of days. The "start date" and "end date" give the study date limits including multiple runs. For a given run in the series, "baseline-begin" and "baseline-end" show the data dates bounding the baseline, followed by a 2-day buffer and then by the dates whose data are being tested for clusters. After each run, the arrows indicate that these blocks move forward by a day.

We implemented the algorithms described in the previous sections into C and MatLab computer programs both for research purposes and for possible use in BioSense. Development of codes giving identical random draws in both software environments allowed perfect agreement of results despite the diverse coding methods. These methods were also developed to facilitate comparison of the spatial estimation methods, to allow inspection of all variables at all stages, and to incorporate two of the following features: (1) computation of the statistical significance of maximum likelihood clusters using the extreme value (or Gumbel) distribution [12], and (2) for the capability to inject realistic signals on the background data for detection performance testing. The C program runs in PC and/or Unix platforms and is callable from SAS code. This code does not have all of the features of SaTScan. Readers interested in obtaining the code should contact the corresponding author.

Our implementation searches cylinders with circular bases. The cylinder height is the number of days in the current test period, so that for a 7-day cylinder, a candidate cluster might have from 1 to 7 days of data. The analysis that follows focuses on comparing spatial estimation methods and is restricted to 1-day cylinders. For the spatial search, any point grid of possible cluster centers could be used, and our study used the set of centroids of zip codes in the data. For analyses based on facility and presumed residence zip codes, the numbers of zip codes were 32 and 1,233, respectively. We defined the radius of a candidate cluster as the maximum distance from the center to another centroid in the cluster, and we allowed radii as large as 100 km for the residence zip code runs and 300 km for the clinic zip runs, based on possible care seeking behavior.

We applied the Poisson Generalized Likelihood Ratio, GLR [11]. Presume that on test day t, there are Obs cases observed within a candidate cluster, Exp cases expected (calculated by one of the methods of this paper or otherwise modelled) in the cluster, and Tot is the total number of cases in the entire surveillance region since the beginning of the baseline. Our implementation followed the common practice of restricting attention to clusters for which more cases are observed than expected. When this condition fails, we set the test statistic to zero, and when Obs > Exp, the test statistic is GLR, or the logarithm of:

This statistic gives a measure of the likelihood of the observed counts inside and outside the cluster given the expected distribution and the total number of observed cases. Candidate clusters within the search space are ranked in descending order by their GLR values. The next step is to decide which if any of these clusters should be flagged as statistically significant.

The LLR values used as the test statistic do not comply with a known distribution in general, so thresholds for significance are usually determined empirically. The empirical approach is to randomly generate a large set of simulated data distributions under the null hypothesis of the expected spatial distribution for the test period. For each simulated distribution k, apply the same search procedure used for the observed data to obtain the maximum test statistics max(LLR(k)) over all candidate clusters. Derive the LLR threshold for significance from the set of max(LLR(k)) and apply it to decide if one or more of the clusters found using the observed data are significant. Specific steps include:

Reliable quantitative measurement of the sensitivity and specificity of a cluster detection method requires multiple signals distributed in time over the sub-regions of interest. Documentation of authentic outbreaks sufficient to provide enough detail for this purpose is very rare. Previous power studies [15,16] aimed at chronic disease applications focused on increased risk in a few sub-regions. However, infectious disease surveillance is concerned with detection of transient events in which cases are variably distributed in the region of interest. For a controlled, quantifiable detection performance capability, we therefore added to the software an injection component that allows the user the capability to: (1) specify the injected signal in detail, including the number of cases to inject, the epidemic curve of these injected cases by day, and the geometry of distribution of these cases; (2) produce the injected cases for each day using random draws and add them to the background data; (3) run the scan statistics for a set of days including the days of injected cases to determine whether clusters are detected at critical significance levels for the injected signals; and, (4) conduct sets of repeated trials of these tests at different start dates to allow computation of detection sensitivity and timeliness. The stochastic generation of signals and changing data background over time provides realistic variability for these trials.

The signal generation method is as follows:

▪ For each trial, a signal start date is selected, and we generate a stochastic epidemic curve using random draws from a lognormal distribution [17,18].

▪ Parameters of this distribution are supplied by the user to control the median peak day and time spread of the injected cases.

▪ These random draws are tabulated to give the total number of injected cases on each day. For each day with injected cases, these cases must be assigned to sub-regions, or zip codes in this study.

▪ Three geometries have been implemented to reflect hypotheses of spread of an aerosol pathogen. The implemented geometries include radial, wedge-shaped, and hourglass-shaped spread at rates that the user may specify.

▪ For the results given here, we used the radial spread, in which the drop-off of cases from the centre zip code is proportional to e^-d/k, for a dispersion factor k that determines the radial decay rate.

In the current study, total signals consisted of either 50 or 100 cases added to the original respiratory counts for either home or clinic zip codes. For epicurve generation, we used lognormal parameters ζ = 1.3, σ = .4. These parameters gave region-wide time series epidemic curves of 8–13 days in duration. The radial spatial distribution was calculated using distance from the signal centre in kilometres, and for our Texas DoD dataset, we used a dispersion factor of k = 18 for a plausible spread of cases. We modified this spatial distribution procedure to avoid producing obvious signals by injecting large numbers of counts in zip codes where there were few or no baseline counts. For this reason, the maximum allowed signal for each zip code for each day was at least 2, but no more than the maximum number of background cases in that zip code for the case day of week. We ran the distribution algorithm imposing this criterion recursively until all cases were distributed. For the repeated trials to measure detection sensitivity, we chose a start day for the first signal, and in each subsequent trial, the start date was advanced by 8 days to sample the background counts at different seasons, avoid day-of-week bias, and reduce correlation among successive trials. We generated these trials over approximately 3 years of background data to obtain 127 trials. For comparison of detection capability using residence zip codes and facility zip codes, we first produced injected signals at the resolution of the 1,233 residence zip codes and ran the scan statistics at that resolution. We then distributed these cases to the 32 facilities using distributions based on the distance to each facility, for comparison runs at the coarser resolution.

For the sequential runs without injects to assess background cluster rates, the update procedure is depicted in Figure 4. Our primary purpose was to compare spatial estimation methods, so the test period for all runs was a single day with no look-back, the purely spatial case. For dates 1/1/2004 to 12/31/2006, we started the first run at the end of the 2-day buffer after the first baseline period. The maximum spatial cluster sizes were 300 km and 50 km for DoD facility zip scanning and patient zip scanning, respectively.

We applied the four estimation methods to scan all test dates and all coordinates in this study dataset. Because they are to be used for daily surveillance, and disease outbreaks are expected to be rare, a primary concern for using these methods is the actual cluster determination rate, or how often the method will indicate a need for investigation. If alarms occur too often, then the alarms corresponding to true health events might be ignored.

We ran the purely spatial cluster detection algorithm for the 3 years of data with no correction for temporal multiple testing to get the relative cluster determination rates of the four methods. In the absence of information about any true outbreaks during the data period, we could not refer to significant clusters as true or false alarms, so as a surrogate for a false alarm rate, we measured the background alarm rate as the number of significant clusters per unit time for the p-value 0.003.

Table 1 shows the number of significant clusters found per 100 weekdays and per 100 weekends, counting as significant those clusters with at least 2 sub-regions per cluster and with Gumbel p-value ≤ 0.003. Figure 5 shows the background cluster rates using the facility-level zip codes for (a) the rash syndrome and (b) the respiratory syndrome. Analogous plots for runs using residence-level zip codes are given in Figures 6(a) and 6(b). Key observations are:

▪ From Table 1, in general, the rash syndrome yield lower cluster rates than the respiratory syndrome, probably because of the seasonality and high mean and overdispersion of the respiratory counts.

▪ For the rash syndrome, the W2 estimation method with the 28-day baseline has a clear advantage over the unstratified averaging of the C2 method at the facility level, but not at the residence level, at which the day-of-week effect is less clear. The plotted levels of cluster alerts suggest that the full respiratory syndrome is too noisy for this dataset in the sense that its high mean and overdispersed variance cause much more significant but random clustering than would be expected for a sparser outcome variable, or especially for a dataset in which spatial distribution is similar to that of census data. Replacement by a more specific sub-syndrome such as influenza-like illness (ILI) should be considered; even though ILI counts are also seasonal, counts are much smaller, and their spatial distribution is likely to be fairly stable for most of the year.

▪ For the respiratory syndrome, management of day-of-week (DOW) effects is a primary concern, and the W2 and DOW methods clearly yield fewer cluster alerts than methods that ignore these effects. See Figures 5(b) and 6(b).

▪ The longest baseline of 56 days gives reduced background clustering only in certain circumstances, such as when the DOW method is used at the residence zip level. The finer temporal DOW stratification, combined with the finer spatial stratification, appears to require a longer baseline. For coarser aggregation, the longer baseline can actually increase the number of cluster alerts, likely because of seasonal effects.

▪ Extending the temporal stratification from the W2 to the DOW method never yields a clear overall advantage in reduced alerting.

▪ The W2 method with a 28-day baseline is a safe choice for controlling the number of cluster alerts overall and gives as clear advantage in some situations. However, for a sparse outcome variable such as the rash visit counts and a fine spatial subdivision of records, a longer baseline might reduce the alerting.

As described in the Methods section, we randomly generated 127 visit count signals with a region-wide lognormal epidemic curve and a radically decaying spatial distribution. Separate signal sets were generated for totals of 50 and 100 outbreak cases. For each set, we added these signals to the authentic respiratory syndrome background data for a series of detection trials.

Figure 7 shows the sensitivity versus p-value threshold for the C2 and W2 detection methods associated with no injects (background) and 100 injects. When no injects were used, C2 had a higher background cluster rates than W2 for a given p-value. A background alert rate of 0.02 (indicated by the dashed line in Figure 7) corresponded to a p-value of about 10^-4 for W2 and 10^-6 for C2, apparently because neglecting day-of-week effects produces spurious cluster alerts when C2 is used. For this 0.02 common background alert rate, we found an empirical detection probability of 0.55 for W2 vs 0.37 for C2. For p-values < 10^-4, the plots show that C2 and W2 have similar sensitivities, but C2 has a much higher background cluster rate.

In (a) and (b), the x-axis gives the estimated background cluster rates from the dotted curves of Figure 7, and the y-axis gives the estimated signal detection rates from the solid curves. Figure 8(a) was derived for 100 sets of 50 injected counts, while Figure 8(b) is the same for 100 injected counts. Thus, in Figure 8(b) for 100 added cases, the 0.02 background level corresponds to 0.55 and 0.37 for the W2 and C2 estimation methods, as discussed above. For the zip codes of the injected cases, these plots are analogous to ROC curves giving sensitivity as a function of false alarm rate. Thus, for the 50-inject signals, a 5% background alarm rate gives sensitivities of 38% and 28% for cluster detection using the W2 and C2 estimation methods, respectively. For the 100-inject signals, the analogous sensitivity estimates are approximately 70% for W2 and 50% for C2.