Geographic disparities and predictors of suicide mortality risk in Florida: spatial scan statistics and negative binomial modeling

Ethics review

This study was reviewed by the University of Tennessee Institutional Review Board (IRB) which determined that it was not human subjects research (IRB Number: UTK IRB-24-083540-XM) because the secondary data used in the study did not contain information on living individuals. Therefore, the board determined that IRB oversight was not required.

Study area

This study was conducted in Florida, a US state composed of 67 counties (Fig. 1). As of 2020, Florida’s population was estimated at 21.5 million, with the largest percentage residing in the southeast (32.1%), followed by southwest (20.8%), central (20.1%), northeast (11.0%), south (9.0%), and northwest (7.0%) regions (United States Census Bureau, 2020) (Fig. 1). Miami Dade, Broward, and Palm Beach, all located in the Southeast, are among the most populous counties in the state and collectively represent 65.3% of Florida’s population (Khan et al., 2022). Approximately 62.7% of the counties in Florida are considered mostly urban, 32.8% are mostly rural, and 4.5% are completely rural (UnitedStatesCensusBureau, 2010) (Fig. 1). This classification, based on the 2010 US Census Bureau County rurality levels, categorized counties as mostly urban if less than 50% of the population lived in rural areas, mostly rural if 50–99.9% of the population lived in rural area, and completely rural if 100% of the population lived in rural areas (UnitedStatesCensusBureau, 2010).

Data sources

The secondary data used in the study were obtained from the Florida Department of Health (Florida Department of Health, 2024). The data were aggregated to county geographical scale and split into three 3-year time periods: 2011–2013, 2014–2016 and 2017–2019. Death records of all suicide deaths in Florida were used in the study. Deaths resulting from suicide were identified using the International Classification of Diseases 10th revision (ICD-10) External Causes of Death Codes X60-X84, and Y87.0 (Saman et al., 2012). The county-level population estimates were obtained from the 2012, 2015, and 2018 American Community Survey and were used as denominators for 2011–2013, 2014–2016, and 2017–2019 time periods, respectively (UnitedStates Census Bureau, 2023).

The potential predictor variables used in the analyses and data sources for each are listed in Table 1. The main data sources were the Florida Department of Health (Florida Department of Health, 2024), County Health Rankings and Road Map (Peterson et al. 2024; University of Wisconsin Population Health Institute, 2024), and the American Community Survey Veteran Status S2101 (Lavrakas, 2013). The 2023 Geographic Boundary Files for cartographic displays were obtained from the United States Census Bureau TIGER Geodatabase (United States Census Bureau, 2024).

Computation of age-adjusted suicide mortality risks and cartographic displays

Suicide mortality risks were calculated and directly age-standardized to the 2010 US population in STATA version 8 using the “dstdize” command (Li et al., 2019; Mehmetoglu & Jakobsen, 2022). Using first-order queen contiguity spatial weights, spatial empirical Bayesian (SEB) smoothed risks were computed in GeoDa version 1.22.0.2 (Anselin, 2024). Risk smoothing was necessary to account for spatial autocorrelation and variance instability resulting from low population areas (Saman et al., 2012). Raw (or unsmoothed) age-adjusted suicide mortality risks and SEB risks were displayed as choropleth maps for each time period using Jenks’ classification scheme (Brewer & Pickle, 2002). All map visualizations were created using ArcGISPro version 3.2.0 (Environmental Systems Research Institute, 2024).

Investigation of high-risk clusters

To identify counties with significantly high suicide mortality risks, Tango’s Flexible Spatial Scan Statistics (FSSS) was used to detect both circular and irregegular shaped clusters using FleXScan version 3.1.2 (Takahashi, Yokoyama & Tango, 2010). The FSSS works by imposing a large number of overlapping, flexibly shaped windows of varying sizes over the study area to identify both circular and irregular shaped clusters (Tango & Takahashi, 2005). For this study, the specifications of Tango’s FSSS analysis were as follows: (i) Poisson probability model; (ii) restricted log likelihood ratio test (ensures that counties with low risk are not included in high-risk clusters); (iii) default alpha value of 0.2 (recommended by Tango & Takahashi (2005) and Tango & Takahashi (2012)); (iv) maximum scanning window size of 34 counties which is approximately half of the counties in the study area (adopted from Odoi (2020)); and (v) 999 Monte Carlo replications for statistical inference (Tango & Takahashi, 2005; Tango & Takahashi, 2012; Tango, 2021).

Investigation of predictors of county-level suicide mortality risk

Data for the most current available time period (2017–2019) were used to investigate predictors of county-level suicide mortality risks. Potential predictor variables used in the analysis were from 2018. Potential predictor variables considered in the investigation were selected based on a conceptual model (Fig. 2) (Posey, 2009; Kaplan et al., 2014; Blanchflower & Oswald, 2020; Kearns et al., 2020; Khader et al., 2020; Amiri, 2022; Graham & Ciciurkaite, 2023; Fabiano et al., 2024; Kareem et al., 2024; Blais et al., 2025). The conceptual model includes social determinants, behavioral risks factors, and mental health stressors that have been shown in past studies to be associated with suicide. Sociodemographic variables such as gender, race, age, education, veteran status, and marital status were investigated as potential predictors of exposure to adverse life conditions that impact risk of suicide. Previous studies have shown that veteran status, often characterized by effects of deployment stress, was associated with physical, psychological, and psychosocial outcomes which significantly elevates suicide risk particularly when combined with binge drinking (Posey, 2009; Blais et al., 2025). It has been reported that veterans who engaged in binge drinking were 33% more likely than non-veterans to report any suicide risk, and 72% more likely to plan a suicide without an attempt (Blais et al., 2025). This factor was also highlighted by Kaplan et al. (2014) who found that suicide decedents had significantly higher odds of acute alcohol use and intoxication prior to death. The sociodemographic factors included in the conceptual model have also been shown to be associated with unemployment, food insecurity, and single parenthood. Amiri (2022) reported that unemployment was strongly associated with suicide ideation, attempts, and mortality. Graham & Ciciurkaite (2023) showed that food insecurity increases suicide ideation among young adults, largely through elevated stress and social isolation. Similarly, Kareem et al. (2024) showed that single parenthood is linked to higher risk of depression and suicidality, emphasizing the emotional and physical toll of solo caregiving. Insufficient sleep, onset-insomnia, and sleep-related distress have been reported to be significant predictors of suicide ideation and mortality in adolescent and college student populations (Kearns et al., 2020; Khader et al., 2020). Additionally, a recent study showed that frequent mental distress has increased over time, particularly among economically disadvantaged populations (Blanchflower & Oswald, 2020), reinforcing its role as an important factor in suicide risk. Moreover, factors such as smoking, physical inactivity, obesity, and excessive drinking, often co-occur with frequent mental distress and may impact suicide risk (Fabiano et al., 2024). Thus, the conceptual model presented in Fig. 2 helped ensure that variable selection and model development were guided by evidence from existing scientific literature (Dohoo, Martin & Stryhn, 2012).

To reduce the likelihood of multicollinearity, pairwise correlations of all potential predictors were assessed using two-way Spearman rank correlation analyses. Only one of a pair of the highly correlated variables (r_s>—0.7—) was considered for further analysis. The choice of which of a pair of highly correlated potential predictors was retained for further analysis was based on literature review (biological plausibility) and statistical considerations. Potential predictors were further assessed for bivariate (unadjusted) associations with the outcome (county-level suicide mortality risk) using univariable negative binomial regression models with the number of suicide deaths as the outcome and the log of the total population per county as the offset. Negative binomial model was chosen because of presence of overdispersion in the data. Variables that were significantly associated with suicide mortality risk based on a liberal critical p ≤ 0.2 in the univariable analysis were included in subsequent multivariable modeling process (Dohoo, Martin & Stryhn, 2012). The assumptions of the negative binomial model include: linear relationship between the dependent and predictor variables, outcome is a count, no multicollinearity, presence of overdispersion (variance of the dependent variable is larger than its mean). This is a key feature that distinguishes a negative binomial model from a Poisson model. The form of the negative binomial model is: ${\displaystyle E\left(Y\right)=n\lambda \text{or}\frac{\mathrm{E}\left(\mathrm{Y}\right)}{n}=\lambda }$ where, E(Y) is the expected number of cases of suicide deaths, n is the measure of exposure, and λ is a function of the predictors.

The usual form of λ is derived from a linear equation on a natural log scale: ${\displaystyle 1n\left(\lambda \right)={\beta }_0+{\beta }_1{X}_1+{\beta }_2{X}_2+{\beta }_3{X}_3+\dots +{\beta }_k{X}_k.}$ This can also be expressed as: ${\displaystyle \lambda =e^{{\beta }_0+{\beta }_1{X}_1+{\beta }_2{X}_2+{\beta }_3{X}_3+\dots +{\beta }_k{X}_k}}$ where, β₀ is the intercept term representing the log count when all predictors are zero, β₁…β_k are the regression coefficients, and X₁…X_k are the predictors.

The mean and variance are represented as follows: ${\displaystyle \text{Mean}:E\left(Y\right)=\mu }$ ${\displaystyle \text{Variance}:Var\left(Y\right)=\mu +\alpha {\mu }^2}$ where, α is the dispersion parameter. If α = 0 it means the variance is equal to the mean and the model is then a Poisson model.

A manual backward elimination procedure was used to identify the most parsimonious multivariable negative binomial model using a critical p ≤ 0.05. Confounding was assessed by comparing the regression coefficients of variables in the model before and after the removal of a suspected confounder. A 20% or more change in coefficients of any variables in the model indicated confounding. Confounding variables were retained in the model regardless of their statistical significance. Biologically plausible two-way interaction terms of the variables in the final main effects model were also assessed and statistically significant ones retained in the final model. The coefficients from the final multivariable model were exponentiated and interpreted as risk ratios.

The goodness-of-fit of the final model was assessed using the chi-square goodness-of-fit test computed using both Pearson and the deviance residuals. Statistical analyses were performed using R statistical software version 4.3.2 (R Core Team, 2024), implemented in RStudio (R Studio Team, 2024). The geographical distributions of significant predictors were presented as choropleth maps created in ArcGISPro version 3.2.0 (Environmental Systems Research Institute, 2024).

Geographical distribution

There were slight increases in the overall state age-adjusted suicide mortality risks from 22.6 deaths per 100,000 persons during the time periods 2011–2013, to 23.2 during 2014–2016, and 24.3 during 2017–2019 time periods. Although the spatial patterns of the SEB smoothed suicide mortality risks mirrored those of the raw risks, the spatial patterns of the SEB risks were more apparent/identifiable (Figs. 3 and 4). It is evident that the highest suicide mortality risks were consistently observed in the northwest, northeast, southwest, and some parts of the central regions of Florida during all three time periods (Figs. 3 and 4). Based on the SEB suicide mortality risk maps, there was evidence that the number of counties with high risks were increasing over time (Fig. 4).

Suicide clusters

Consistent with the spatial patterns of both the raw (unsmoothed) and SEB smoothed suicide mortality risks, clusters of high suicide mortality risks were identified consistently in the northwest and central regions (Fig. 5 and Table 2). Significant high suicide mortality risk clusters identified during the time periods 2011–2013 and 2014–2016 were in the northwest, central and southwest regions, whereas the high-risk clusters during 2017–2019 were in the eastern side of the northeast, and some parts of central, southwest and south regions. Seven counties were persistently in high-risk suicide mortality clusters throughout the study period. These counties were Volusia and Brevard counties in the eastern region; Citrus, Hernando, Pasco, and Pinellas counties in the southwest region; and Marion County in the central region of the state (Figs. 1 and 5).

Predictors of geographic disparities in suicide mortality risks

A total of 16 potential predictors had significant univariable associations with suicide mortality risk based on a relaxed critical p = 0.2 (Table 3). Based on the final multivariable negative binomial model (Table 4), suicide mortality risks tended to be high in counties with high percentages of the population that: were aged 45–64 years (RR = 1.010; p = 0.017), were aged ≥ 65 years (RR = 1.008; p = 0.016), drank alcohol excessively (RR = 1.064; p < 0.001), had frequent mental distress (RR = 1.039; p < 0.001), and were veterans (RR = 1.022; p = 0.005). Although not statistically significant insufficient sleep (p = 0.264) was included in the model because it was an important confounder of age 65 years and above. The significant (p = 0.010) interaction between excessive drinking and insufficient sleep implies that the effect of the percentage of population involved in excessive drinking on county risk of suicide deaths depends on the percentage of population in the county that experience insufficient sleep (Table 4). The risk ratios seem relatively small because they represent the effect of 1% increase in each of the variables. Higher increases in the percentages of the significant predictors would result in much higher effects. For instance, a 10% increase in population reporting frequent mental distress would result in a 2.34 times higher risk of suicide deaths (e^(0.085∗10) = 2.34) at the county level when all other variables are kept constant. There was no evidence of either multicollinearity (variance inflation factor (VIF) were all below 1.9) or lack of fit (Pearson Chi-Square tests, p = 0.26; deviance chi-square test, p = 0.33).

Counties in the northwest and northeast regions of Florida tended to have high percentages of the population: aged 45–64 years, with frequent mental distress, and who were veterans (Fig. 6). Counties with high percentages of population who drank excessively were in the northeast, southwest, south, and the eastern side of central and southeast regions of the state while those with a high percentage of populations aged 65 years and above tended to be in the northwest, central, and southwest regions of the state.

Strengths and limitations

Use of SEB age-adjusted mortality risks ensure adjustment for both spatial autocorrelation and risk instability due to high variances in low population counties. The implementation of a FSSS with a restricted likelihood ratio enables the identification of both circular and non-circular clusters. Moreover, employing a restricted likelihood ratio test in place of the original log likelihood ratio minimizes false positives. The limitation of the study is that it was conducted at the county-level and hence patterns observed at this geographic scale may be different from those at lower geographic scales. Unfortunately, due to small numbers, the study could not be conducted at levels lower than the county spatial scale.