Detecting Spatial Patterns of Disease in Large Collections of Electronic Medical Records Using Neighbor-Based Bootstrapping (NB2)

TL;DR

This paper presents NB2, a neighbor-based bootstrapping method for quantifying and analyzing geospatial disease variation using large-scale electronic medical records, aligning well with existing spatial autocorrelation techniques.

Contribution

The paper introduces NB2, a novel, flexible method for detecting spatial patterns in disease incidence from EMR data, capable of identifying variations at multiple scales.

Findings

NB2 agrees with Moran's I and kriging in detecting spatial autocorrelation.

NB2 can identify both large and small area geospatial variations.

Applied to 100 million records, it effectively ranks diseases by spatial variation.

Abstract

We introduce a method called neighbor-based bootstrapping (NB2) that can be used to quantify the geospatial variation of a variable. We applied this method to an analysis of the incidence rates of disease from electronic medical record data (ICD-9 codes) for approximately 100 million individuals in the US over a period of 8 years. We considered the incidence rate of disease in each county and its geospatially contiguous neighbors and rank ordered diseases in terms of their degree of geospatial variation as quantified by the NB2 method. We show that this method yields results in good agreement with established methods for detecting spatial autocorrelation (Moran's I method and kriging). Moreover, the NB2 method can be tuned to identify both large area and small area geospatial variations. This method also applies more generally in any parameter space that can be partitioned to consist of regions and their neighbors.

Equations

$$Y_{crude}^{age,gender}=\frac{\text{mbox}casesofICD-9}{\text{mbox}totalcases}\times \frac{100000\text{mbox}persons}{8\text{mbox}years}.$$

$$Y_{adjusted}={\Sigma }_{age,gender}{Y}_{crude}^{age,gender}\times \frac{\text{mbox}grouppopulation}{\text{mbox}totalpopulation}$$

$$I=\frac{n}{{\Sigma }_i{\Sigma }_j{w}_{ij}}\frac{{\Sigma }_i{\Sigma }_j{w}_{ij}(y_i-\bar{y})(y_j-\bar{y})}{{\Sigma }_i(y_i-\bar{y}{)}^2}$$

Full text

Detecting Spatial Patterns of Disease in Large Collections of Electronic Medical Records Using Neighbor-Based Bootstrapping (NB2)

Maria T Patterson1, 2*, Robert L Grossman1, 3, 4, 5**

**1 Center for Data Intensive Science, University of Chicago, Chicago, IL

2 Department of Astronomy, University of Washington, Seattle, WA***

3 Computation Institute, University of Chicago, Chicago, IL

4 Section of Computational Biomedicine and Biomedical Data Science, Department of Medicine, University of Chicago, Chicago, IL

5 Institute for Genomics & Systems Biology, University of Chicago, Chicago, IL

**

***Email: maria.t.patterson@gmail.com

**Email: robert.grossman@uchicago.edu

***MP currently works at the University of Washington

**

keywords: geospatial variation of disease incidence, geospatial correlation, electronic medical records

Abstract

We introduce a method called neighbor-based bootstrapping (NB2) that can be used to quantify the geospatial variation of a variable. We applied this method to an analysis of the incidence rates of disease from electronic medical record data (ICD-9 codes) for approximately 100 million individuals in the US over a period of 8 years. We considered the incidence rate of disease in each county and its geospatially contiguous neighbors and rank ordered diseases in terms of their degree of geospatial variation as quantified by the NB2 method.

We show that this method yields results in good agreement with established methods for detecting spatial autocorrelation (Moran’s I method and kriging). Moreover, the NB2 method can be tuned to identify both large area and small area geospatial variations. This method also applies more generally in any parameter space that can be partitioned to consist of regions and their neighbors.

INTRODUCTION

As the number of sources and the volume of electronic medical records (EMR) and electronic health records (EHR) increases, there is a growing ability to aggregate these data and extract information about population and public health [1, 2]. Over the past several years, new sources of digital health data, such as web searchers, social media, mobile phones, and personal health sensors have increased the number of sources and data volumes even more [3, 4, 5]. Much of these data can be geocoded with location information so that techniques from spatial epidemiology can be used to explore geospatial variation in disease, health outcomes, and population health using disease mapping, disease cluster analysis, and related techniques [6, 7, 8, 9]. With this much geocoded digital health data, there is a need for simple tools and algorithms that can be used by researchers across disciplines for identifying the presence of spatial autocorrelation in disease incidence data, especially in large datasets [10].

An initial starting point for evaluating the presence of patterns in disease or other geocoded data is determining whether the data are spatially autocorrelated, that is, whether the disease rates or values of interest are similar in nearby areas and fall off with distance, which could indicate the presence of core areas of disease risk [11, 12].

We introduce a Monte Carlo based algorithm that we call Neighbor-based Bootstrapping (NB2) that can be used to quantify geospatial autocorrelation. We apply this algorithm to approximately 100 million geocoded EMR and rank order 548 diseases as determined by ICD-9 codes from those with the strongest geospatial autocorrelation to those with the weakest geospatial autocorrelation. We compare this method’s results to Moran’s I statistic [13] and to kriging [14, page 44], two other techniques that have been used to quantify geospatial autocorrelation. The spatial size scale of disease patterns may range widely from small, localized affected regions to larger affected areas, depending on the nature of the underlying factors. We have developed two versions of the NB2 ranking, one favoring patterns of tight clusters and the other favoring broader less peaked patterns.

Applying geospatial analysis and visualization techniques to geocoded health data has long been understood to be important for identifying risk factors from the physical environment and for providing insights into the transmission of infectious and vector-borne diseases [15, 16, 17, 18, 19, 20, 21]. For example, spatial analysis of health data can be used to identify and manage risk associated with proximity to potentially harmful environmental exposures, such as chemical toxins or air pollutants [22, 23, 18]. More generally these techniques are also important for understanding a broader range of risk factors, including risk factors from the demographic, economic, social, cultural, regulatory, or legal environments [24, 25, 26, 27, 28, 29, 30, 31, 32].

MATERIALS AND METHODS

Data

The dataset consists of electronic medical records data from the Truven Health MarketScan Commercial Claims and Encounters Database, which includes approved inpatient and outpatient insurance claim information for a total of approximately 100 million unique and de-identified individuals across the United States for the time period from 2003 to 2010. The records include 1.3 billion diagnostic ICD-9 subdivided codes (12.89 unique codes per person), geotagged by county FIPS code. Here we restrict to using the approximately 800 non-subdivided ICD-9 codes from 001-799, which excludes injuries, poisonings, and accidents. We refer to ICD-9 codes by the 3 digit integer group. (For example, “005: Other bacterial poisoning” includes “005.0 Staphylococcal food poisoning”.)

We also restrict our analysis to the 3,109 counties in the continental US and to the ICD-9 codes that have data for two-thirds or more of the counties. This leaves 548 ICD-9 codes.

For each of these 548 codes, we adjust for age and gender by using standard populations [33] as follows. We determine crude incidence rates for the standard 19 groups of age populations for each gender by taking raw counts for each group and dividing by the population at risk, which in this case we take to be the total number of records for each county for each age/gender group converted to 100,000 person-year units:

[TABLE]

The age and gender adjusted rate is calculated by multiplying the crude rate for each group by the appropriate weight using the Census 2010 standard population and summing the products [34]:

[TABLE]

Neighbor-based bootstrapping (NB2) method

The NB2 method uses resampling to evaluate in this example whether or not the incidence rate of a disease can be accurately estimated from the incidence rate of the disease in counties that are neighbors. The first step in this method is to define regions and neighbors of regions. Here we define regions as counties and neighboring counties as counties that are geospatially contiguous to the county’s polygon border, including vertices (Queen style), though it is important to note that there are many options to consider when defining neighbor relationships (contiguity, distance, spatial weights) that have varying effects on results [35, 36, for example]. In this paper, we focus on geospatially defined neighbors, but an advantage of this method is that it is applicable without change to neighbors in any space of features, not just neighbors in 2 or 3 dimensional physical space.

We compute a bootstrapped estimate as follows. Fix a county Y. For each ICD-9 code, we sample with replacement a set of neighboring counties and a set of random counties and compare the normalized disease incidence (from Eq. 2).

More explicitly, fix a county $Y$ and assume that it has $n^{Y}$ neighbors. We estimate the log incidence rate $Z_{neighbor}$ for county $Y$ as the average log incidence rate of a list of $n^{Y}$ randomly chosen (with replacement) neighboring counties. We also estimate for each county the log incidence rate $Z_{random}$ for county $Y$ as the average log incidence rate of $n^{Y}$ randomly chosen (with replacement) counties from the full set of all counties. These counties may or may not be neighbors.

We compare the two estimates (neighbors vs random) to the known log incidence for each of the drawn counties in two separate ways. See Algorithm 1 and Algorithm 2.

In the first implementation, we take the difference from actual of the estimates of $Z_{neighbor}$ vs $Z_{random}$. We then use a paired Student’s t-test to evaluate whether the neighbor based predictions are a significant improvement over the random prediction. For ICD-9 codes with significant underlying spatial patterns, we expect that the $Z_{neighbor}$ estimates will be significantly closer to actual than the $Z_{random}$ estimates.

We repeat this process to obtain 1000 estimates of the neighbor-based vs random differences, and for each of these compute the paired t-test. We then take the median t-test value from these 1000 estimates. This gives us one t-test statistic value per ICD-9 code, describing how closely related incidence rates of that ICD-9 are in neighboring counties as compared with a random selection of counties.

In the second implementation, we compare the neighbors vs random estimates by counting, for each pair of bootstraps, the number of samples where the neighbor estimate is closer to actual than the random estimate. We then repeat this 1000 times, take the median number, and using this to calculate the log odds that the neighbor estimate is more accurate than the random estimate.

RESULTS

Performance

We first evaluated the impact of varying the number of times $M$ that we resampled. Running the entire procedure and resampling $M=1000$ times for all 548 diseases takes just over 28 hours on a virtual machine with 8 Xeon cores running at 2.00 GHz with 16 GB of RAM. This is about 25 minutes per disease using a single core.

For 100 bootstraps, the run time for 548 diseases on 8 cores takes about 220 minutes, or a little over 3 minutes per disease when using a single core. Comparing the NB2 statistic values for 1000 vs 100 simulations, the difference on average is 0.2% and at maximum 1.6%. For 10 bootstraps, the total run time is about 30 minutes, or about 30 seconds per disease using a single core. The mean difference between NB2 statistic values for 1000 and 10 simulations is 0.2%, and the maximum difference is 3.8%. The results are summarized in the table below.

[TABLE]

In the analysis that follows, we are primarily focused on the rank ordering of the ICD-9 codes according to these two implementations of the NB2 method. There is no significant difference in the rank orderings between 1000 and 100 or 10 repeated bootstraps.

Comparison with Moran’s I statistic

We compare the neighbor based bootstrapping results to the global Moran’s I statistic for detecting spatial autocorrelation, which is based on the sum over weights between units multiplied by the mean-adjusted outcome of interest divided by the squared mean difference of each point. Moran’s I is defined as [13]:

[TABLE]

where n is the total number of spatial polygons (counties), $y_{i}$ is the value of interest of the $i$th polygon, $\bar{y}$ is the global mean, and $w_{ij}$ is the spatial weight of the link between polygon $i$ and $j$.

Moran’s I ranges from $-$1 (perfect dispersion, as in a black and white checkerboard pattern) to 1 (black squares on one side, white on the other). A random distribution would have I close to 0. We compare the values of Moran’s I for the set of log incidence rates across counties for each ICD-9 code to both the implementation of the NB2 method using the paired Student’s t-test evaluation and the implementation using the log odds evaluation. If the geospatial variation that the NB2 method detects is similar to Moran’s I, then the NB2 statistic values should increase as Moran’s I goes to 1. In Figure 1, we show the NB2 t-test statistic estimate (left) and the log odds estimate (right) plotted against Moran’s I statistic for all ICD-9 codes tested. Generally, the NB2 statistic values increase as Moran’s I statistic increases, though there is noticeable scatter.

We rank ordered the ICD-9 codes using the two NB2 method implementations and the Moran’s I statistic to produce three ordered lists of ICD-9 codes from the strongest spatial correlation (largest Moran’s I statistic, largest NB2 t-test, largest NB2 log odds test) to the weakest. Tables Comparison with Moran’s I statistic - Comparison with Moran’s I statistic contain the top 25 ICD-9 codes for both NB2 implementations and Moran’s I. Here we will compare the properties of the spatial distributions for ICD-9 codes ranked highly by the two NB2 procedures with those ranked highly by Moran’s I statistic.

[FIGURE:]