https://orcid.org/0000-0002-1484-5296
Abstract
This article introduces bayesian spatial smoothing models for disease mapping—a specific application of small area estimation where the full universe of data is known—to a wider audience of public health professionals using firearm suicide as a motivating example. Besag, York, and Mollié (BYM) Poisson spatial and space–time smoothing models were fitted to firearm suicide counts for the years 2014-2018. County raw death rates in 2018 ranged from 0 to 24.81 deaths per 10 000 people. However, the highest mortality rate was highly unstable, based on only 2 deaths in a population of approximately 800, and 80.5% of contiguous US counties experienced fewer than 10 firearm suicide deaths and were thus suppressed. Spatially smoothed county firearm suicide mortality estimates ranged from 0.06 to 4.05 deaths per 10 000 people and could be reported for all counties. The space–time smoothing model produced similar estimates with narrower credible intervals as it allowed counties to gain precision from adjacent neighbors and their own counts in adjacent years. bayesian spatial smoothing methods are a useful tool for evaluating spatial health disparities in small geographies where small numbers can result in highly variable rate estimates, and new estimation techniques in R software have made fitting these models more accessible to researchers.
Introduction
Bayesian spatial smoothing models are useful tools for investigating rare outcomes or small geographies. These methods allow researchers to more precisely evaluate health outcomes to identify areas of higher or lower risk. While the value of this approach for small area estimation using complex survey data,1,2 or with incomplete data collection,3 has been well-recognized, it is also useful for disease mapping, a specific application of small area estimation where the full universe of data is known. Employing spatial smoothing models for disease mapping is valuable to (1) more robustly evaluate geographical health disparities when small numbers are an issue, and (2) present reliable rates for all geographies without suppression.
Firearm suicide deaths are a motivating example. Nationally, firearm suicide deaths are not rare. Suicide is the tenth leading cause of death in the United States and has risen to become the fourth leading cause among those in midlife aged 35-54 years.4 Firearms are the most common means of suicide death, accounting for half of these deaths. Firearms are by far the most lethal means of suicide5; although only about 4% of suicide attempts involve a firearm, approximately 90% of these attempts result in death.4 However, small numbers present challenges when examining county-level rates, particularly in rural counties, which typically experience higher rates of firearm suicide compared with urban areas.6,7 County-level analyses are crucial for firearm injury research to understand heterogeneity below the state level in firearm policy effects, take advantage of local natural experiments, and isolate risk and protective factors varying across space.
Traditional approaches to dealing with small numbers have been to combine or increase the level of geography under investigation,6,-9 aggregate over multiple years,10,11 or simply focus on metropolitan areas.12,-14 While these strategies circumvent the issue of small numbers, they introduce other sources of bias, including but not limited to measurement error or pure specification bias for data pooled across populations, potential ecological bias when large areas are used to make inference at a smaller scale, selection bias and a loss of generalizability when conditioning on areas included in analysis, and an inability to understand local variability or investigate rural areas.15
Mathematically, small numbers result in highly variable rates. Rates calculated from small numbers can be more extreme than rates calculated with larger denominators because the addition or subtraction of a single event has a much greater impact on the rate estimate. The Centers for Disease Control and Prevention (CDC) considers rates calculated from fewer than 20 occurrences to be unreliable and prohibits the public reporting of rates calculated with fewer than 10 occurrences due to the potential identifiability risk.16,17 In 2019, 2535 counties had fewer than 10 firearm suicide deaths, meaning a firearm suicide mortality rate could not be reported for 80.7% of US counties.18 Notably, to overcome these issues the CDC has integrated spatial smoothing methodologies into their Web-based Injury Statistics Query and Reporting System (WISQARS) fatal injury data dashboard making spatially smoothed estimates easily obtainable.
This article introduces the utility of spatial smoothing methods for mapping rare outcomes or small geographies using county-level firearm suicide as an illustrative example. These methods are becoming more common in research and public health practice. Users or consumers of the resulting modeled estimates must be aware of how the data were produced to avoid misinterpretation or misuse. This article presents the theory behind spatial smoothing for disease mapping, specification of the model, when they are appropriate to implement, and how to interpret the resulting estimates.
Methods
Data
Mortality data for 2014-2018 were obtained from the CDC National Vital Statistics System restricted Detailed Mortality files.19 Firearm suicides were classified by International Classification of Diseases, Tenth Revision, codes, and deaths were included if the primary or any of the additional cause of death fields indicated a firearm suicide (Appendix S1, Table S1, available at https://doi.org/10.1093/aje/kwae005). The count of firearm suicide deaths in each year was summed using county of occurrence—where the death took place. We included deaths in the contiguous United States (excluding Alaska, Hawaii, and US territories). Population data were obtained from the CDC Race-Bridged population files; the distribution of population counts across counties in 2018 can be found in Figure S1. Raw mortality rates were calculated for each year by dividing the count of firearm suicides by the county population to understand the spatial distribution of real-world values, suppressing rates for areas with fewer than 10 deaths.
Statistical approach
We first present an example of a spatial smoothing model for county-level firearm suicide mortality rates in a single year, and then extend this method to illustrate space–time smoothing of rates over several years. The following models were fitted on the full, complete mortality data for each county, without suppression. Accompanying code used to present these examples is available in Appendix S2.
Bayesian vs frequentist statistics
The spatial smoothing method outlined here uses a bayesian approach. In bayesian statistics, as opposed to frequentist statistics, the observed data—here, firearm suicide counts—are considered known or fixed, and the parameters being estimated—such as latent county firearm suicide rates—are random unknowns. The model uses the observed data inputs, their spatial structure, a specified likelihood, and model prior specifications to estimate posterior probability distributions of the mortality rates for each area, from which one can extract both the most likely estimate as well as its 95% credible interval as a measure of uncertainty.
Spatial structure and the neighbor matrix
Spatial smoothing methods are attractive because geography acts as a proxy for unmeasured covariates which make adjacent areas more similar to each other than to those farther away; they do not require the inclusion of carefully selected covariates to produce valid estimates. This follows Tobler’s First Law of Geography: “Everything is related to everything else, but near things are more related than distant things.”20(p236) Individual counties’ firearm suicide mortality estimates gain precision by leveraging outcome data from their neighbors. This analysis followed a common approach wherein county neighbors were defined by their adjacency using a binary indicator of 1 if the counties share a common border and 0 otherwise. The neighbor matrix is the array of 0 and 1 values classifying every county pair in the analysis as either neighbors or not neighbors.
The neighbor matrix informs the spatial random effects—we used intrinsic conditional autoregressive (ICAR) random effects as they are computationally efficient, depending only on local spatial autocorrelation from immediate neighbors.21 These spatial random effects allow the raw county estimate for each area to be weighted towards the mean of the rates in adjacent counties, thus smoothing the estimate towards rates in neighboring areas. The degree of smoothing is dependent upon the variability of the estimate in that county based on the number of expected outcomes, so counties with small numbers are more affected by the influence of their neighbors. ICAR random effects can be thought of as a spatial version of the temporal first order random walk.
Likelihood, prior distributions, and hyperparameters
Bayesian inference requires specification of two main components: (1) the likelihood, which is a probability distribution from which the data are assumed to be drawn, and (2) prior probability distributions for all model parameters. Prior distributions can be informative—based on previous data or expert opinion—or weakly/noninformative, which forces the model to depend almost entirely on the data observed. Weakly/noninformative priors are therefore most appropriate with large datasets. Prior distributions rely on modeler-specified hyperparameters to incorporate this prior knowledge (or lack thereof). Hyperparameters should be carefully chosen as they can have a potentially large effect on model results, especially when using small datasets.
The prior distributions are combined with a specified likelihood for the observed data to calculate a posterior distribution for the random parameters. It is the resulting posterior distributions that are used for inference via point estimates and measures of uncertainty. Point estimates are commonly calculated as the mean or median of the posterior distributions for each parameter of interest. Uncertainty is typically summarized with posterior credible intervals, which, for a specified probability (eg, 95%), give a range of values for which the posterior probability of a parameter being in this range is equal to the specified probability. Different prior distributions lead to different shapes of posterior distributions, which can affect both point estimates and credible interval widths.
Priors are one of the most challenging elements of the model to understand and implement—readers should consult a statistician if they have information from prior studies or trials to inform the model. Most tools for statistical modeling include recommended priors in their documentation, which can be used as a default, but sensitivity analyses should be done using different priors to examine their impact on results.
Model specification
The combination of spatial ICAR and IID normal random effects is known as the Besag, York, and Mollié (BYM) model22 and is argued to have the most robust spatial structure of commonly used spatial models.23 We used the BYM2 reparameterization of the BYM model,24 which specifies the total variance of the spatial and nonspatial random effects along with the proportion of this variance that is spatial. This aids computational efficiency and model interpretability by ensuring the random effects are on the same scale.
This analysis used an uninformative prior on the intercept, along with penalized complexity (PC) priors for the BYM2 random effects.25 PC priors are commonly recommended to aid interpretation of hyperparameters.26 For modeling, PC prior hyperparameter values were chosen to be relatively uninformative and favor a simpler model fit, which means the posterior probability distribution used for model inference will be informed mostly by the observed data.
After exponentiating, the model output provides an estimate of the firearm suicide count for each county in the dataset with corresponding posterior distributions, from which a rate can easily be estimated using the offset.
Measuring uncertainty
Smoothed county-specific firearm suicide mortality rates were estimated as the median of the posterior distributions for each county. Uncertainty was measured by calculating 95% credible intervals, which are similar to but have a different interpretation compared with frequentist confidence intervals. The 95% credible intervals were calculated by identifying the 2.5th and 97.5th percentiles of the smoothing model’s posterior distribution for each county. The posterior probability that the latent mortality rate is in this interval is equal to 95%; thus, the credible interval gives a plausible range of values for the county-specific firearm suicide mortality rate. The width of the 95% credible interval is a useful metric for understanding the magnitude of the uncertainty in the modeled estimates to compare across counties.
Extending to a space–time framework
While it is possible to calculate annual smoothed mortality rates separately for each year of data, a more efficient method for estimating a time-series is to model all years by including an additional random effect for year. This allows borrowing strength over time to improve yearly estimates. Here we model time as a random walk of order 1,27 which allows each year to be influenced only by the data from the previous year.
This model builds on the spatial smoothing equation above with the addition of two temporal terms: an unstructured temporal error term |${\mathrm{\omega}}_t$| and the |${\mathrm{\phi}}_t$| temporal first-order random walk term for each time t. |${\mathrm{\omega}}_t$| represents the deviation at time t from the overall mean. |${\mathrm{\phi}}_t$| is conditional on its immediate preceding neighbor in time and the amount of smoothing is governed by the smoothing parameter |${\mathrm{\sigma}}_{\mathrm{t}}^2.$| Similar to the spatial random effects, the magnitude of temporal smoothing depends on the variance of values over time. It is important to note that this model only has separate spatial and temporal effects and does not include a spatiotemporal interaction, which would allow spatial effects to vary across time and/or temporal effects to vary across space. Models with spatiotemporal interactions can be computationally prohibitive to fit and are beyond the scope of this paper, but further information can be found from other sources.27
Computation
Estimation of the posterior distribution can be done via various statistical methods. While certain models allow the prior distribution to be chosen such that the form of the posterior distribution can be calculated directly, this is often not the case in practice due to the complexity of models used for small area estimation. Thus, modelers must use other methods to approximate and sample from the posterior distribution. One of the most common general methodologies is Markov chain Monte Carlo (MCMC).28 MCMC methods sample a chain of autocorrelated realizations from the posterior distribution in a specific way allowing for full exploration of the posterior. Traditional MCMC approaches can be quite computationally intensive for complex spatial models, which is often a barrier for many practitioners. Another approach, and the one used in this paper, is integrated nested Laplace approximation (INLA),29 which is typically much faster and less computationally intensive because it approximates the necessary integrals for calculating posterior distributions in strategic ways. The INLA package26,30 in R (R Foundation for Statistical Computing, Vienna, Austria) implements this method for a large class of models—latent Gaussian models—which includes the vast majority of spatial smoothing models most practitioners might want to fit. Furthermore, it has excellent documentation for new users, making it accessible to a wide audience.
All analyses were performed in R, version 3.6.2,31 using the INLA,26,30 dplyr,32 sf,33 sp,34,35 and spdep35,36 packages. This analysis is considered nonhuman subjects research and did not require IRB review. For a deeper look into the statistical underpinnings and mathematical structure behind these bayesian spatial smoothing methods, please see prior work.26,37,-40
Results
All 3084 US contiguous counties were included in the analysis; 114 108 firearm suicide deaths were recorded during the 5-year study period of 2014–2018. In 2018, county-level raw (unsmoothed) death rates ranged from zero to 24.81 deaths per 10 000 people (Figure 1). However, this highest mortality rate was based on only 2 deaths in a population of approximately 800; 2484 counties experienced fewer than 10 firearm suicide deaths in 2018—80.5% of counties—and thus were required to be suppressed. An additional 324 counties experienced fewer than 20 deaths, meaning that based on CDC guidelines it would not be possible to report a reliable firearm suicide mortality rate for 91.1% of contiguous US counties in 2018.
Raw county-level firearm suicide mortality rates with counties experiencing fewer than 10 deaths suppressed, United States, 2018.
After the spatial and space–time smoothing models were employed, it was possible to report a reliable firearm suicide mortality rate for every county for the years 2014-2018. The 2018 spatially smoothed county-level firearm suicide rates ranged from 0.06 to 4.05 deaths per 10 000 people (Figure 2). Both the highest and lowest raw rates were smoothed towards the center of the data. Smoothing is more apparent for extreme high rates as they are often highly variable, but even low rates are smoothed towards the middle, particularly counties with zero deaths. These smoothed firearm suicide mortality rates can be interpreted as the estimated “true” risk, or the underlying firearm suicide rate for a hypothetical infinite superpopulation living in each county, and thus even if a county experienced zero deaths, the smoothed estimated death risk can never be zero since we assume it is impossible to have zero risk.
Spatially smoothed county-level firearm suicide mortality rates, United States, 2018.
Smoothed firearm suicide mortality rate estimates may be more precise than the raw rates but are still calculated with uncertainty. The widths of the 95% credible intervals, a measure of this uncertainty, ranged from 0.06 up to 9.46 deaths per 10 000 people, with a median of 1.09 (Figure 3). As expected, the results from the space–time smoothing model (Figure S2) were more precise than spatial smoothing alone because county estimates gained information both from adjacent spatial neighbors and their own mortality counts in the prior year. The median width of the credible intervals for the 2018 space–time smoothed mortality rate estimates was 0.77 and ranged from 0.3 to 3.9 deaths per 10 000 (Figure S3). Rate estimates between the spatial and space–time models differed slightly for 2018 due to the influence of additional data from adjacent years, with a median difference of 0.3 and ranging from −1.06 to 1.5 deaths per 10 000. Estimated space–time smoothed firearm suicide rates from 2014 through 2018 are shown for a random sample of 12 counties (Figure 4).
Width of the 95% credible intervals of spatially smoothed firearm suicide mortality rate estimates, United States, 2018.
Annual estimated space–time smoothed firearm suicide mortality rates with 95% credible intervals, in a sample of 12 US counties.
Total running time for the R script used to calculate the raw and smoothed firearm suicide rate estimates for every contiguous US county in 2018 was 3 minutes and 12 seconds on a 2018 MacBook Pro (Apple Inc., Cupertino, CA) with 16GB of SDRAM. Run time may vary based on computer infrastructure and resolution of geographic information system data used for the neighbor matrix, but with the R-INLA package, computation time is no longer a major barrier for this type of complex bayesian analysis.
Discussion
This firearm suicide example illustrates the two primary benefits of spatial smoothing methods for disease mapping. First, it allowed us to more reliably compare smoothed firearm suicide mortality rates across different geographic contexts. Second, we were able to estimate and display precise mortality rates for areas experiencing 20 or fewer deaths—the majority of US counties. Since this model-based smoothing approach pools and weights data from surrounding counties, the rate estimates do not pose the same concern for identifiability in reporting.
Smoothed firearm suicide mortality rates illustrate the geographic variability in risk of firearm suicide across the United States (Figure 2), as well as their general trends over time (Figure 4). The smoothed estimates clearly show the lowest rates of firearm suicide in the Mid-Atlantic Coast, parts of the Midwest near Wisconsin and Iowa, and most of mid and Southern California, while the highest rates tend to occur in the Mountain West and Southwest. Almost all counties show a gradually increasing rate of firearm suicide mortality over time. This pattern matches the steady increase in all suicide rates, which have continued to climb over the past 20 years.41,42
While spatial and space–time smoothed models have many favorable qualities, they require informed decision making and have several potential drawbacks. The decision to smooth rates is an example of the well-known statistical trade-off between variance and bias; bias is introduced by smoothing over some real local heterogeneity to increase precision and make comparisons across geographies more reliable. In space–time models, rates over time are also smoothed by borrowing information from adjacent years and thus real, informative spikes or dips across time can be flattened out as well. A misuse of spatially smoothed estimates would be for detection of statistically significant hotspots either over space or time. Some of the apparent clustering of firearm suicide mortality rates by geography is in fact due to the spatial smoothing approach. By design, these methods produce estimates in adjacent areas that are more similar to each other, and thus should never be used as inputs for cluster or hotspot detection analyses.
The goals of spatial smoothing for disease mapping are related to variance correction and not prediction. In this application the full universe of data is known, as is the case when using vital statistics or comprehensive registry data. The raw mortality rates are accurate, as they reflect real counts experienced in each geography, but they are not precise or reliable because the potentially high variance implicit in small number calculations creates unstable rates. While these methods allow us to compare rates of firearm suicide deaths across the entire country over time, individual areas should supplement their understanding of firearm suicide risk with local knowledge whenever possible.
Unlike other nonspatial small area estimation methods, the use of covariates in these spatial smoothing models is not necessary as the influence of the neighbors accounts for unmeasured covariates. However, their inclusion may aid in increasing precision, particularly when data are sparse or when there are strong known outcome-covariate relationships that do not follow a spatial pattern and can be measured well ecologically. This example of firearm suicide deaths includes sparse data, sometimes clustered in particularly rural areas such as the Mountain West, where highly variable-rate counties are unable to gain much precision from their equally variable neighbors, leading to higher uncertainty reflected in wider credible intervals. This level of uncertainty may be unacceptable depending on the intended application for the smoothed estimates. Including additional years of data in a space–time model can help increase precision by incorporating more information for these sparse estimates. Additionally, given the heterogeneous firearm-related legal landscape across states and the known association between permissive firearm laws and higher mortality rates, future studies may consider including indicators of state policies or simply state fixed effects. However, the spatial smoothing approach is still preferred over nonspatial approaches as adjacency to states with permissive firearm laws and even the occurrence of nearby gun shows has been found to increase risk of firearm-related deaths across state borders, even within states with more a more restrictive legal landscape.43,44
How to classify neighbors is also a crucial consideration as adjacency can be defined in several ways, which may be more or less appropriate for each specific research question. This analysis used a straightforward method of classifying neighbors by shared administrative boundaries. The resulting sparse, symmetrical, zero/one neighbor matrix allows for faster computation17 but may not reflect the true nature of between-county interaction as the simple structure allows only directly adjacent counties to share information. Neighbors can also be assigned using a distance threshold from the center of each area—either Euclidean or road network travel—or even by manually and painstakingly assigning neighbors using multiple criteria such as distance, contiguity, landscape features influencing contact, known social/economic networks, or other measures relevant to a particular analysis. Without local knowledge of the complex interdependencies among counties, the continuity approach is the simplest and most used. However, it may not work well for all counties, particularly for adjacent counties with very dissimilar contexts. It is also incapable of incorporating noncontiguous areas such as Alaska, Hawaii, and Puerto Rico, which is particularly unfortunate for firearm injury studies as these communities have unique legal and cultural landscapes and understudied vulnerable populations.
Spatial disease mapping is more accessible to researchers than ever. Traditional bayesian spatial smoothing techniques required time and computationally intensive MCMC iterations to create the posterior distributions used for inference. However, INLA has been found to perform well in fractions of the time and computational effort. Spatial smoothing using the R-INLA package can be a simple, straightforward, and computationally feasible method for most researchers familiar with R. Public health scientists may consider adding these methodologies to their toolkit, particularly when interested in rare outcomes or areas smaller than states.
In conclusion, with improvements made to data storage and computation, the ability to assess local mortality trends on a wide scale has increased and the demand for estimates below the state level has grown. Smaller areas have less susceptibility to ecological bias than larger aggregations because there is less within-area heterogeneity,15 but the reliability of these local estimates is a concern when based on small numbers. A bayesian approach leveraging the inherent spatial structure of outcomes can better identify true underlying risks in small geographies to identify genuine geographical disparities for more robust research.
Supplementary material
Supplementary material is available at the American Journal of Epidemiology online.
Funding
This work was supported by funds from the State of Washington to the Firearm Injury and Policy Research Program.
Conflict of interest
The authors declare no conflicts of interest.
Data availability
The mortality data used in this analysis is restricted-use and cannot be shared publicly. It was obtained from the National Center for Health Statistics, National Vital Statistics System. The R code for this analysis can be found in the manuscript supplement, and is also posted on GitHub: https://github.com/Epi-Emma/Firearm_Suicide_Disease_Mapping
