Religion and Death in the United States: A Meta-Regression Comparative Assessment of Between-County Mortality Heterogeneity in the United States

Religion can have a favorable impact on individual-level health. The influence of religion on population health, however, remains less clear. This study investigated the association between religion and mortality at the population-level. Using county data, a meta-regression was performed to examine between-county mortality heterogeneity. The percent heterogeneity associated with religion variables were compared to demographics (i.e., place, race, language, age, and gender) and health factors (i.e., individual behaviors, clinical care, social and economic, and physical environment) as predictors of mortality. Religion was measured in terms of adherence (i.e., prevalence attending/belonging to a congregation), congregation density, and the diversity of adherents and congregation by denominations. Results showed counties with lower mortality were associated with higher proportions of religion adherents and a greater diversity of adherents and congregations. Counties with higher mortality were associated with higher religion congregation density. Religion, as a parsimonious multivariate model with all demographic and health factor predictors, had less added value when controlled for individual variables or constructs. The direction of association between religion and mortality was consistent, even when controlling for demographics and health factors, and thus merits further consideration as a population health determinant, as it may play a critical role in understanding other population health outcomes.


Introduction
Religion, however defined, is both ubiquitous and variant across populations; a simple road trip is often sufficient to reveal how ubiquitous and variant religion is in the United States. A road trip would also reveal considerable demographic variation such as population densities, race, and poverty/wealth; it would only be natural to wonder if these variances are related. There is considerable evidence of an individual-level relationship between religion and health, mostly with respect to mental and physical health and from a clinical perspective. However, the relationship at a population level has not been studied as extensively as it has at the individual level. The extant individual-level evidence and the ubiquity/variance of religion across populations warrants further study of religion as a determinant of population health [1]. and health data (e.g., U.S. Religion Census and the University of Wisconsin Population Health Institute's County Health Rankings and Roadmaps).
Conceptual models are necessary but not sufficient in establishing the nature of any relationship between population health and religion; empirical data, the what, is also required. The what will answer questions such as: (1) is religion ubiquitous and variant across populations; (2) how does population health vary with religion, and then how well does religion variation explain population health variance; and (3) does any relationship between religion and population health outcomes remain when controlling for more traditional population health determinants? These are complex questions for research as they lead to more questions such as how best to measure population health and religion.
This article focuses on population health measured as mortality. Religion is measured in terms of adherence (prevalence attending/belonging to a congregation), congregation density, and the diversity of adherents and congregation by denominations. The nature of any relationship between mortality and religion will be compared to well-established relationships between mortality and demographic and health factors at the county level.
Meta-regression with county data (in lieu of study data) is used to measure mortality heterogeneity across U.S. counties and the percent heterogeneity explained by religion, demographic, and health factor variables. Religion data were obtained from the 2010 U.S. Religion Census; mortality, demographic, and health factor data were obtained from the PHI County Health Rankings & Roadmaps. Demographic constructs include measures of place, race, and gender. Health factor constructs include individual behaviors, clinical care, social and economic factors, and physical environment. Heterogeneity explained is examined as single and multivariate covariates of mortality. The research questions include: 1. What is the direction and strength of the association, if at all, between religion variables and mortality, uncontrolled and controlled for county demographics and health factors? 2. Do religion variables explain less, more, or comparable levels of population mortality heterogeneity as do demographic and health factor variables? 3. Does the inclusion of religion variables have added value when studying population mortality controlling for demographics and health factors?

Data
Data obtained from the 2010 U.S. Religion Census included the number of adherents, the number of congregations, and county population size [7]. The census data were reported by state, by county, and by individual congregation. Summary statistics were also reported as an overall total (any denomination) and by six major denominations: Evangelical, Mainline Protestant, Black Protestant, Catholic, Orthodox, and Other [8]. Number of adherents were reported by individual congregations within each county and defined as including those with an affiliation to the congregation (children, members, and attendees who are not members); congregation was defined as a group of people who meet regularly (typically weekly or monthly) at a pre-announced time and location and may be churches, mosques, temples, or other meeting places [9].
Mortality, demographic, and health factor data were obtained from the County Health Rankings and Roadmaps website, which is a program of the University of Wisconsin Population Health Institute, funded by the Robert Wood Johnson Foundation [10]. County data are compiled from multiple sources over time and used to rank the health of individual counties within each state (a detailed description of individual data used for this research is provided as an Appendix A; additional information is available at their website). Demographic and health factor data include:

Variables for Descriptive Analyses
The proportion of adherents for each county was estimated as the number of adherents divided by the county population size, as reported in the 2010 U.S. Census. When the proportion was greater than 1 (29 counties, particularly for congregations near a county line reporting adherents from an adjacent county), the county's proportion of adherents was capped at 1. Congregation density was estimated as the number of congregations per 100,000 population. Diversity of religion was estimated both for adherents and for congregations using formulas reported elsewhere as Religious Diversity Index and Denominational Pluralism, respectively [11,12]. Both formulas are variants based on the Herfindahl-Hirschman Index, and values range from 0 to 1, with 0 reflecting no diversity and 1 an equal distribution of adherents and denominations across the six denomination major groups.
Mortality was estimated for each county as the average numerator (number of deaths over two time periods) divided by the average denominator (population for two time periods) per 100,000 population. All demographic and health factor variables were analyzed descriptively using their reported metrics (averaging when reported over multiple years-see Appendix A).

Variables for Meta-Regression Analyses
The mortality numerators and denominators were used to estimate double arcsin effect sizes and standard errors [13,14]. County population size and congregation density were converted to natural log values. Arcsin (angular) transformation was performed on all covariates measured as proportions to "stretch out" proportions close to 0 and 1 [15]. Adherent diversity, congregation diversity, and all remaining variables were analyzed in the meta-regressions in their original metric without any further transformation.

Analyses
Descriptive statistical analyses of data were conducted using Excel (Microsoft Office 365); correlations and collinearity were examined using SPSS v27 (IBM, Armonk, NY, USA). Meta-regression analyses were conducted using Stata v17. The meta-regression analyses were analyzed as random effects models using the Knapp-Hartung variance estimator [16].
The percent of heterogeneity explained by the addition of covariates was estimated as an adjusted R 2 (expressed as a percent) using Stata's "meta" suite of commands. The "added value" of religion data was calculated as the difference between the R 2 of parsimonious meta-regression models with religion variables and the R 2 of a model without religion variables.
After assessing correlations and collinearity, individual variables were explained as single covariate random effects meta-regression models; multivariate models were then assessed for religion, place (population size and rurality), race and language, age and gender, all demographics, individual behaviors, clinical care, social and economic, physical environment, and all health factors. These full models were then reduced to parsimonious multivariate models. The assessment of added value for religion for demographic and health factors were carried out using parsimonious models. Table 1 reports the descriptive statistics for the county variables in this study as number of counties with data, means, medians, and minimum and maximum values. Also included in this table are the single covariate meta-regression model results reported as percent heterogeneity explained (meta-regression adjusted R 2 values) and the direction of association and statistical significance (p < 0.05) of the covariate reflected by color: yellownot statistically significant; red-statistically significant in counties with higher mortality; and green-statistically significant in counties with lower mortality). Counties with higher mortality were statistically significantly (p < 0.05) associated with higher values of: congregation density; rurality; African American, American Indian, or Alaskan Native; age 65 years or older; female; smoking; obesity; physical inactivity; sexually transmitted infections; uninsured; preventable hospital stays; unemployment; children in poverty; income inequality; inadequate social support; single parent households; injury deaths; air pollution; drinking water violations; driving alone to work; and long commute when driving alone. Counties with lower mortality were statistically significantly associated (p < 0.05) with higher values of: adherent diversity; congregation diversity; population size; Asian, Hispanic, Native Hawaiian, or other Pacific Islander; not-proficient in English; age less than 18 years; food environment index (healthier foods); access to exercise opportunities; excessive drinking; primary care physicians; dentists; diabetes monitoring; mammography screening; high school graduation; some college; and severe housing problems. Variables that were not statistically significant (p ≥ 0.05) included: adherents (proportion); White; alcohol-impaired driving deaths; and violent crime.

Descriptive Statistics
The single variable explaining the highest between-county mortality heterogeneity was injury deaths (52.5%). Four variables explained between 40 to 50%: children in poverty (49.3%); physical inactivity (49.1%); some college (41.0%); and smoking (40.4%). Three variables explained between 30 to 40%: teen births (39.5%); congregation density (33.2%); and obesity (30.0%). All remaining variables individually explained less than 30%; adherent and congregation diversity explained 11.0% and 9.8%, respectively. Figure 1 portrays scatter plots for mortality and the four religion variables, all in the original (not transformed) metrics. Three of the variables (adherents, adherent diversity, congregation diversity) are negatively correlated with mortality, which reflects counties with lower mortality tend to have higher proportions of adherents (r = −0.02) and are more religiously diverse both in terms of adherents (r = −0.33) and congregations (r = −0.31). Counties with higher mortality tend to have greater congregation density (r = 0.39) or more congregations per 100 K population. Figure 1 portrays scatter plots for mortality and the four religion variables, all in the original (not transformed) metrics. Three of the variables (adherents, adherent diversity, congregation diversity) are negatively correlated with mortality, which reflects counties with lower mortality tend to have higher proportions of adherents (r = −0.02) and are more religiously diverse both in terms of adherents (r = −0.33) and congregations (r = −0.31). Counties with higher mortality tend to have greater congregation density (r = 0.39) or more congregations per 100K population.  Table 2 reflects substantial correlations (r ≥ 0.60) for the "top ten" (from Table 1) variables explaining heterogeneity as single covariate meta-regression models, and other substantial variable correlations (top half of table), as well as other variables (bottom half). Positive correlation coefficients reflect variables associated with counties having higher mortality rates, and negative with counties have lower mortality rates. Injury deaths, which explains 52.5% of between-county mortality, was not substantially correlated with any other variable. Children in poverty had substantial correlations with four other variables (teen births, college, healthy food, and single parent). Congregations (density) had  Table 2 reflects substantial correlations (r ≥ 0.60) for the "top ten" (from Table 1) variables explaining heterogeneity as single covariate meta-regression models, and other substantial variable correlations (top half of table), as well as other variables (bottom half). Positive correlation coefficients reflect variables associated with counties having higher mortality rates, and negative with counties have lower mortality rates. Injury deaths, which explains 52.5% of between-county mortality, was not substantially correlated with any other variable. Children in poverty had substantial correlations with four other variables (teen births, college, healthy food, and single parent). Congregations (density) had substantial correlations with three variables (population size, rural, and Asian). Physical inactivity was the third highest in explaining heterogeneity and was substantially correlated with obesity. Hispanic was substantially correlated with non-English proficiency. STI was substantially correlated with Black, White, and single parent. Adherent diversity was substantially correlated with congregation diversity.

Meta-Regression Results
The variable White (non-Hispanic) had a variance inflation factor (VIF) greater than 10 (as did African American and Hispanic), but when White was removed, all VIF values were less than 10. White was therefore not included in any multivariate model. Table 3 reflects the progression of individual variables into full and parsimonious multivariate random effects meta-regression models (in this table, the exact value for individual covariates is reported, and color coding, as in Table 1, is again used for all covariates). Variables that were reported by less than 3000 counties were also not included in the multivariate analyses in preference to having multivariate models with at least 3000 counties; these variables were: smoking (n = 2711), excessive drinking (n = 2225), inadequate social support (n = 2471), and violent crime (n = 2917).  The transition of specific single covariate variables to parsimonious multivariate models of note (direction of association change) included adherents, congregation diversity, Pacific Islander, age less than 18, sexually transmitted infections, and high school graduation. Specifically, adherents as a single covariate became statistically significant and associated with lower morality in the parsimonious model, while congregation density became statistically not significant. Pacific Islander (Native Hawaiians or other Pacific Islander) changed from being statistically significantly associated with lower mortality to higher mortality; age less than 18 and high school graduation had a similar change. Sexually transmitted infections changed from a statistically significant association with higher mortality to lower mortality. Variables that were statistically significant as single covariates and became statistically significant in the multivariate models included: adherents, and alcohol-impaired driving deaths. The construct than explained that the highest percentage of mortality heterogeneity was social and economic (74.9%), followed by individual behaviors (60.9%), religion (43.2%), race and language (42.6%), age and gender (24.4%), place (19.9%), and the least by physical environment (10.5%). Figure 2 reflects the consistency of the religion-mortality direction of association when controlling for separate demographic and health factor constructs, as well as when controlling for all demographic and health factor constructs simultaneously (using parsimonious models). A higher proportion of adherents and higher values of adherent diversity are consistently associated with lower county mortality (negative regression coefficients). Higher congregation density is consistently associated with higher county mortality (positive regression coefficients).  Figure 3 portrays the added value of religion when assessing between-county mortality heterogeneity when controlling for individual demographic and health factor constructs, as well as when controlling for all demographics; all health factors; all demographics and health factors (using parsimonious models). The religion construct alone (Table 3) explains 43.2% of between-county mortality. Depending on how much is explained by separate demographic and health factor constructs, the added value of religion ranges from 7.3% (when added to Social and Economic) to 45.7% (when added to Physical Environment). The added value of religion decreases when added to combined models: 12.9% (when added to all Demographics), 1.0% (when added to all Health Factors), to only 0.1% (when added to all Demographic and Health Factors).   Figure 3 portrays the added value of religion when assessing between-county mortality heterogeneity when controlling for individual demographic and health factor constructs, as well as when controlling for all demographics; all health factors; all demographics and health factors (using parsimonious models). The religion construct alone (Table 3) explains 43.2% of between-county mortality. Depending on how much is explained by separate demographic and health factor constructs, the added value of religion ranges from 7.3% (when added to Social and Economic) to 45.7% (when added to Physical Environment). The added value of religion decreases when added to combined models: 12.9% (when added to all Demographics), 1.0% (when added to all Health Factors), to only 0.1% (when added to all Demographic and Health Factors).  Figure 3 portrays the added value of religion when assessing between-county mortality heterogeneity when controlling for individual demographic and health factor constructs, as well as when controlling for all demographics; all health factors; all demographics and health factors (using parsimonious models). The religion construct alone (Table 3) explains 43.2% of between-county mortality. Depending on how much is explained by separate demographic and health factor constructs, the added value of religion ranges from 7.3% (when added to Social and Economic) to 45.7% (when added to Physical Environment). The added value of religion decreases when added to combined models: 12.9% (when added to all Demographics), 1.0% (when added to all Health Factors), to only 0.1% (when added to all Demographic and Health Factors).

Discussion
Our research asked questions regarding religion as a population health determinant of mortality variation between U.S. counties. We assessed data from up to 3141 U.S. counties, examining: 1. The direction and strength of any association between religion and mortality, uncontrolled and when controlled for known demographic and health factor mortality determinants. 2. Whether religion explained more, less, or was comparable to demographics and health factors. 3. If the data suggest religion has added value to population health research of mortality.
The best answer to these questions has to be "it's complicated". As single covariates or as a single construct, religion is consistently associated with mortality, regardless of what these variables are controlled for. Counties with lower mortality rates are associated with a higher proportion of adherents, greater adherent diversity, and greater congregation diversity. This is consistent with the religion-health evidence at the individual level where religion is a source of social support. Counties with higher mortality rates are associated with higher congregation densities. The only logical explanation for this it that congregations are fulfilling their mission of caring for those who need support. Counties with higher mortality also experience higher rates of poverty, obesity, smoking, teen births, etc.; when life in the present life is dismal at best, it is perhaps natural to hope for something better in the next, and religion may offer that, and thus a greater demand is responded to by greater supply.
However, religion is a label, and much like labels of race, it does not provide an accurate reflection of the complexity of individuals that fall under the label. Therefore, it is not surprising that, upon controlling for health factors that more directly connected with mortality, the strength of association between religion and mortality diminishes. In fact, in this study, once religion was controlled for by demographics and health factors, the "added value" of religion to population health studies of mortality decreases to where religion has an added value of less than one-tenth of a percent. But depending on what data is available to a research project, and if religion data is available, there may still be added value.
There are two important limitations to our research. First, this is an ecological study where counties are the unit of analysis. Readers may be tempted to generalize these findings to individuals, but that is not possible, and we would be remiss in our responsibility if we failed to make note of this. But our interest was in counties, and as an "individual entity", these data are generalizable to counties or other similar geographical units.
A second limitation is that the adherent and congregation density religion data reflected "any denomination" and thus may not apply to individual denominations. This is supported by the two diversity indices (adherents and congregation) that were both associated with lower mortality. Additional research examining specific denominations is warranted.
In the progression towards parsimonious meta-regression models, several variables dropped out due to statistical non-significance; even when in a single covariate model, they were statistically significant. Highly correlated variables are essentially proxies for each other, and when one drops out in a multivariate model, it does not lose its importance. Smoking and excessive drinking were not included in the movement towards a parsimonious multivariate model because these variables had less than 3000 counties reporting: 2711 and 2225, respectively. As a sensitivity analysis leaving these variables in (and thus reducing the multivariate model number of counties), this parsimonious model (2134 counties) explained 68.5% of between-county mortality, 7.63% points more than the individual behaviors parsimonious model reported in Table 2. In this analysis, smoking was statistically significantly associated with counties with higher mortality, while excessive drinking was significantly associated with counties with lower mortality. Excessive drinking was also significant at the single covariate level, and this association is more likely due to underreporting.
Last, mortality is a single measure of population health. These findings therefore do not generalize to other population health measures. Additional research of religion and other outcomes is also warranted.

Conclusions
Religion is a complex variable. The measure "adherents" does not assess the depth of one's faith, nor does it assess the frequency of attending religious services. More importantly, "religion" or "being a religious person" is a label that is not homogeneous. One "religious person" will most likely be substantially different from another "religious person" in terms of factors that are more directly linked to mortality and other health outcomes. The same is true for counties and any other population group. This study was not designed to produce definitive truths but rather to simply ask if religion was worthy of a population health study.
In that regard, this study is a success. The consistency of association between mortality and each of the religion variables, even when controlling for demographic and health factors, is intriguing. The hope of this research was to identify areas worthy of additional study. Why is the relationship consistent? Does the same relationship exist for other population health outcomes? Are the relationships similar between religious denominations? Can population health improvement strategies be more effective when engaging and collaborating with the faith communities?
The study of religion as a determinant of health, either at the individual level or at a population level, has been less extensive in public health than in the clinical sciences; this avoidance recalls the old maxim that cautions any discussion of politics and religion at the dinner table [17]. This avoidance of religion (and politics) in public health circles is ironic given that public health and religion have actually been connected for a long time; this avoidance may be changing. Former Dean of the Emory School of Public Health, James Curran, shares a face of history not well-known regarding John Snow and the removal of the infamous pump handle. Little reference has been given to public health students about Snow's collaboration with a local parish curate, Reverend Henry Whitehead, whose intimate involvement in the life of the people which led to "understanding how they walked to get water, offering critical insights to how disease spread" [18]. The current Dean of the Boston University School of Public Health, Sandro Galea, reminds us that religion and public health overlap through social justice and the safeguarding of vulnerable populations, and that given the influence of religion and its direct effect on billions of people, public health leaders today recognize how partnerships with religious organizations would help create healthier populations, and public health would be remiss to ignore this opportunity [19].  Data Availability Statement: The data of this study can be provided upon reasonable request to the corresponding author.

Conflicts of Interest:
The authors declare no conflict of interest.