Study design and population

In this ecological study, we utilized a spatiotemporal hierarchical Bayesian model to analyze mortality data related to CAD, focusing specifically on women aged 50 to 79 across 399 municipalities in Paraná, southern Brazil. We specifically focused on this age group due to the significant physiological changes associated with menopause, in which the prevalence of CAD increases after the 5th decade of life [18]. It is well known that estrogen provides a protective effect against cardiovascular diseases, and including pre-menopausal women could introduce a confounding factor into our analysis [7]. However, this protective effect ceases during the menopausal transition. Regarding the location, Paraná is the fifth most populous state in the country, with 11,597,484 inhabitants in 2021 [19]. Paraná is one of the three states that constitute the southern region of Brazil and is home to over 11.46 million residents across an expanse of 199,315 km², as reported by the Brazilian Institute of Geography and Statistics (IBGE) (Figure 1) [19]. This places its population in proximity to that of notable countries like Belgium, Greece, and Bolivia [20]. In terms of land area, Paraná shares similarities with Senegal and surpasses Greece and South Korea [20].

Paraná is divided into 22 health regions based on geographical proximity and population coverage, with the primary goal of ensuring equitable access to healthcare services. They play a pivotal role in coordinating and integrating healthcare initiatives, as well as in resource allocation and the implementation of health policies. Furthermore, health regions serve to understand the health situation and living conditions of populations in delimited territories [21]. Moreover, one health region is administratively independent of the other and generally the population of a municipality within a health region is not referred to or cared for in the neighboring health region, even if they are neighboring municipalities. Utilizing health regions as the analytical units in our research provides a comprehensive insight into the spatial distribution of health outcomes and allows for the assessment of healthcare interventions within specific geographic contexts [22].

Data sources and variables

We gathered data on the number of CAD-related deaths among women in the state of Paraná between 2010 and 2019. This data was obtained from the Mortality Information System of the Brazilian Ministry of Health. Our analysis focused on the following International Statistical Classification of Diseases and Related Health Problems (ICD) codes: Angina pectoris (I20), Acute myocardial infarction (I21), Recurrent myocardial infarction (I22), Certain current complications following acute myocardial infarction (I23), Other acute ischemic heart diseases (I24), and Chronic ischemic heart disease (23, I25).

Sociodemographic and healthcare services factors were chosen considering their potential to influence CAD mortality. Therefore, the following variables were collected: primary healthcare service coverage, IPARDES municipal performance index, exercise stress testing rate, echocardiography rates, scintigraphy rates, heart catheterization rates, cardiologist availability rates, accessibility to the cardiology center, mortality, and population of individuals aged above 15 years in each municipality (Table 1).

Data treatment

Primary healthcare service coverage is the percentage of the population that, potentially, would be covered by a specific primary care service, considering the reference of 3,000 individuals per Family Health Team (ESF). This reference is based on the average number of individuals to be served by each ESF, following the recommendation of the Ministry of Health [27].

The IPDM (IPARDES Municipal Performance Index) is an indicator that assesses the performance of the 399 municipalities in the State of Paraná, considering three key dimensions: income, employment, and agricultural production; as well as health and education. Its compilation relies on diverse administrative statistics made available by public entities. This number varies from 0 to 1, where values closer to 1 indicate better living conditions for the population considering the 3 dimensions [28].

Rates for exercise stress testing, echocardiography, scintigraphy, heart catheterization, and cardiologist availability were determined by calculating the number of respective procedures or practitioners between 2010 and 2019, dividing this figure by the population aged 15 and above, and then multiplying the result by 1,000. We used a divider for the population over 15 years of age to avoid pediatric cardiologists being included in the analysis.

To assess the accessibility to the cardiology center, we employed the Two-Step Floating Catchment Area (2SFCA) method as described by Luo et al. [29]. This analysis was performed in two steps. In the first step, we begin by defining a service area around each hemodynamic center. In this case, we used a 60-minute travel radius by car.

The 60-minute travel radius was chosen based on European Society of Cardiology for the management of acute myocardial infarction with ST-segment elevation (STEMI), which recommend a total reperfusion time of less than 120 minutes for STEMI patients. Studies indicate that longer transport times lead to longer door-to-balloon (D2B) times. By limiting transport to 60 minutes, it’s possible to stay within this time window, taking into account the initial care and preparation for percutaneous coronary intervention (PCI) at the referral hospital. An American study on patients with myocardial infarction with ST-segment elevation shows that, for ground transfers, the average door-to-device (D2D) time was 93 minutes for travel times of up to 30 minutes, 117 minutes for travel times of 31 to 45 minutes, and 121 minutes for travel times over 45 minutes. This highlights that longer transfer times result in longer D2D times. Thus, a 60-minute travel radius is justified to ensure that, after initial care and transport, there is sufficient time for the admission and preparation of the patient for PCI, keeping the total time within the recommended limit of 120 minutes [30, 31]. Within this service area, we summed up the number of people residing in that space. In the second step, we calculated the overall accessibility to hemodynamic centers for the entire population. To do this, we sum up the center’s capacity of 100 beds by demand within the area. This accessibility index provides a measure of the accessibility of the population to hemodynamics within the specified travel time radius. The higher the value of the accessibility index, the better the accessibility [32, 33].

Mortality was created using the number of CAD-related deaths, as described in data sources, and multiplied by 100,000.

Exploratory analysis

The initial exploratory measure in this analysis involved raw CAD mortality rates. Subsequently, we constructed standardized mortality rates (SMR) to assess CAD mortality risk for each municipality in Paraná on an annual basis. To calculate these rates, we divided the observed mortality (the actual number of deaths) by the expected mortality (predicted number of deaths obtained through indirect standardization, adjusting the mortality rate for our study population aged 50 to 79 years). This computation was performed using the “dplyr” package in RStudio software, version 4.1.0. SMR values below one indicates a low risk of death from CAD, while values above one indicates a high risk [34].

We conducted a collinearity assessment using the Variance Inflation Factor (VIF), an index that measures the degree to which collinearity affects the variance of a regression coefficient. Detecting multicollinearity is crucial as it helps identify highly correlated relationships among independent variables in which could lead to misinterpretation. This excessive correlation can make interpreting results challenging, lead to instability in regression coefficients, and diminish the predictive capacity of the model [35]. Furthermore, it can complicate the selection of relevant variables when values exceeding five are indicative of potential bias [36]. Therefore, all independent variables utilized in this study demonstrated VIF values close to one, signifying a minimal level of correlation between these independent variables and the outcome. Importantly, while we have minimized the risk of collinearity, it does not necessarily mean that our study is absent of bias.

In the next step, we employed Poisson’s Log-Linear Model, a regression analysis suitable for modeling variables related to count data or rates [37]. We utilized a 95% confidence interval and a significance level of 5% to assess the spatiotemporal correlation between CAD mortality and other variables. Moreover, the dependent variable was number of deaths, and the independent variables were primary healthcare service coverage, IPARDES municipal performance index, exercise stress testing rate, echocardiography rates, scintigraphy rates, heart catheterization rates, cardiologist availability rates, accessibility to the cardiology center. We tested several combinations of independent variables and chose the best set based on the performance according to the lowest values of DIC (Deviance Information Criterion) and WAIC (Watanabe Akaike Information Criterion), and highest values of loglikelihood.

The presence of spatiotemporal autocorrelation was determined by residuals of models, extracting then, the data were extracted year by year to calculate the Moran’s Index (Moran’s I). The values of Moran’s Index vary between –1 (negative autocorrelation) and +1 (positive autocorrelation), with a value of zero corresponding to spatial independence [38, 39].

Hierarchical Bayesian Spatiotemporal Modeling

The first-order autoregressive process was used to quantify the evolution of the spatial pattern in mortality risk over time, designating the random effects vector for all time units. The error vector was modeled as spatially autocorrelated [40].

We applied hierarchical Bayesian spatiotemporal modeling using the CARBayesST package within RStudio 4.1.0 software to achieve posterior spatial distribution smoothing and inference [41]. This method combines prior results with posterior inferences based on Bayes’ theorem.

The model was fit using Markov Chain Monte Carlo (MCMC) simulations, which are iterative algorithms for sampling from a probability distribution. Subsequent risk estimates were obtained from 220,000 simulations, with 20,000 discarded during the burn-in period. We used Gelman-Rubin diagnostic to verify the convergence of the MCMC estimated values as an additional check, and considered good quality when a value is less than 1.1 [42]. Moreover, Geweke diagnostics was also used in order to verify the convergence based on values between –2 and 2 [41]. After verifying both convergence diagnostics, we can proceed with the inferences.

CAD mortality relative risk per year were estimated after dividing the mean of expected from MCMC value by the actual value. Relative risk was reported using exponential of the mean while adopting credible interval of 95% [43].

The same analysis was conducted to assess the independent variables’ ability to globally explain CAD mortality, utilizing the estimated MCMC analysis to determine the mean and standard deviation of each variable [43]. Values below one indicate negative autocorrelation (a protective factor), 1 corresponds to spatial independence, and values above 1 indicate positive correlation. The average trends in disease risk over time and space were calculated based on the average risk for each unit of analysis, for each year, and graphically illustrated using QGIS software version 3.32 [44].

Ethical aspects

Given the data is secondary and publicly available, an exemption was granted by the Committee on Ethics in Research Involving Human Beings of the State University of Maringá, in accordance with Resolution 510/2016 of the National Health Council [45].

Descriptive analysis

Between 2010 and 2019, there were 14,603 deaths due to coronary artery disease (CAD) among women aged 50 to 79 years. This corresponds to an average of 1,460 ± 61.30 deaths per year during this period. The highest number of deaths was recorded in 2014, with 1,531 deaths, while the lowest occurred in 2019, with 1,338 deaths. Despite a noticeable decline in the number of deaths over time, specifically a reduction of 159 deaths (10.62%) since 2010, the overall variation in death rate was less than 2% (see Figure 2).

The average number of deaths from CAD per 100,000 inhabitants (raw rates) was 111.27 ± 61.30. These rates decreased from 132.59 deaths per 100,000 inhabitants in 2010 to 87.36 deaths in 2019. A spatiotemporal analysis revealed a heterogeneous distribution of these rates during the study period, with most regions exhibiting a raw rate lower than 114 deaths per 100,000 inhabitants. Only specific areas within the state, notably in the central-northeast, northwest, and southern regions, recorded values exceeding 200 (Figure 3).

CAD mortality standardized ratio (SMR) and relative risk

An analysis of spatial autocorrelation for CAD mortality revealed a Moran’s Index of 0.63, with p < 0.05, indicating a significant spatial dependency associated with this disease. Consequently, we assessed the standardized mortality ratio and relative risk within each municipality to identify covariates related to CAD. The spatiotemporal map depicting the average SMR values reveals heterogeneity in CAD mortality rates across the state of Paraná. Generally, municipalities with an SMR greater than 1 indicate that the observed mortality exceeded the expected rates. Notably, SMR values exceeding 1 are prevalent in the majority of the Northwest, Central-West, and South-Central regions of the state (see Figure 4).

The comprehensive temporal analysis indicates a notable decrease of approximately 22.6% in the relative risk of CAD mortality among women between 2010 (with a RR of around 1.06) and 2019 (RR of around 0.82). This declining trend became particularly prominent after 2014 (see to Figure 5).

However, in Figure 6, we present the Relative Risk (RR) and the Probability of Exceedance (PEP) every five years for greater clarity. When examining the RR for each municipality individually, it becomes evident that the majority of RR values are below 1, but with some nuances between 1 and 1.5, mainly in certain municipalities in the Central-West, Central-South, Central-East, and Southern regions of the Central-North region. PEP is a metric used to assess the occurrence of a specific event, meaning that the higher the PEP value, the greater the probability associated with the event. For example, a PEP of 0.2 (or 20%) indicates a 20% probability that the event in question will occur. Consequently, the risk associated with the event, with a value of 1.2, indicates that if the event occurs, its impact is 1.2 times greater than a reference situation with a risk of 1.0. This implies that the event has a significantly greater impact than the reference situation. The PEP shows us that there was a dispersion of excessive values in the state, and at the same time, there was also a decrease in municipalities presenting lower values over time. However, the risk still existed in 2019, primarily concentrated in the same regions mentioned above.

Covariates related to CAD mortality

Out of eight covariates tested in our model, only five presented the best performance. IPARDES municipal performance index had a RR of 0.46, representing a 34% reduction in CAD mortality. Similarly, the accessibility to the cardiology center presented a RR of 0.92, corresponding to 8% decrease, while the exercise stress testing rate had a RR of 0.96, signifying a 4% reduction. These results indicate their protective effects (see Table 2).

To comprehend the factors associated with a protective effect on CAD mortality, we have plotted the average values for each variable (Figure 7). From this, it becomes evident that primary healthcare service coverage and the IPARDES municipal performance index are relatively high across the state. Conversely, the availability rate of cardiologists and the rate of exercise stress testing showed lower values. Regarding accessibility to cardiology centers, there are clusters of regions with high values in the Northwest, West, Southwest, and South-Central areas. However, there are also regions with lower values, such as Central-West, North pioneer, East-Central, and the southern part of Southeast.