Metropolitan age-specific mortality trends at borough and neighborhood level: The case of Mexico City

Abstract

Karol Baca-López, Cristóbal Fresno Roles Conceptualization, Data curation, Formal analysis, Methodology, Supervision, Visualization, Writing – original draft Affiliation: Technology Development Department, National Institute of Genomic Medicine, Mexico City, Mexico Jesús Espinal-Enríquez Roles Investigation, Supervision, Writing – review & editing Affiliation: Computational Genomics Department, National Institute of Genomic Medicine, Mexico City, Mexico Miriam V. Flores-Merino Roles Methodology, Resources Affiliation: School of Chemistry, Autonomous University of the State of Mexico, Toluca, State of Mexico, Mexico Miguel A. Camacho-López Roles Methodology, Resources Affiliation: School of Medicine, Autonomous University of the State of Mexico, Toluca, State of Mexico, Mexico Enrique Hernández-Lemus Roles Conceptualization, Formal analysis, Methodology, Project administration, Supervision, Writing – original draft, Writing – review & editing * E-mail: [email protected] Affiliations Computational Genomics Department, National Institute of Genomic Medicine, Mexico City, Mexico, Centro de Ciencias de la Complejidad, Universidad Nacional Autonoma de Mexico, Mexico City, Mexico ORCID logo https://orcid.org/0000-0002-1872-1397 Introduction Identifying trends in Metropolitan Mortality (MM) is a challenging problem nowadays. The authors reported an increase in homicide mortality that surpassed positive outcomes from health care reforms oriented to promote life expectancy at the national level [34, 35]. [...]there is an urgent need to count on fine-grained data in order to better interpret spatio-temporal correlations in health outcomes. In this context, the spatial component is tackled at borough and neighborhood scale, whereas the time scale was used in a year and trimester levels. [...]using a multilevel modeling for longitudinal data provides a finer description of these phenomena that

About the Authors:

Karol Baca-López

Contributed equally to this work with: Karol Baca-López, Cristóbal Fresno

Roles Data curation, Software, Visualization, Writing – original draft

Affiliations School of Medicine, Autonomous University of the State of Mexico, Toluca, State of Mexico, Mexico, Computational Genomics Department, National Institute of Genomic Medicine, Mexico City, Mexico

Cristóbal Fresno

Contributed equally to this work with: Karol Baca-López, Cristóbal Fresno

Roles Conceptualization, Data curation, Formal analysis, Methodology, Supervision, Visualization, Writing – original draft

Affiliation: Technology Development Department, National Institute of Genomic Medicine, Mexico City, Mexico

Jesús Espinal-Enríquez

Roles Investigation, Supervision, Writing – review & editing

Affiliation: Computational Genomics Department, National Institute of Genomic Medicine, Mexico City, Mexico

Miriam V. Flores-Merino

Roles Methodology, Resources

Affiliation: School of Chemistry, Autonomous University of the State of Mexico, Toluca, State of Mexico, Mexico

Miguel A. Camacho-López

Roles Methodology, Resources

Affiliation: School of Medicine, Autonomous University of the State of Mexico, Toluca, State of Mexico, Mexico

Enrique Hernández-Lemus

Roles Conceptualization, Formal analysis, Methodology, Project administration, Supervision, Writing – original draft, Writing – review & editing

* E-mail: [email protected]

Affiliations Computational Genomics Department, National Institute of Genomic Medicine, Mexico City, Mexico, Centro de Ciencias de la Complejidad, Universidad Nacional Autonoma de Mexico, Mexico City, Mexico

ORCID logo https://orcid.org/0000-0002-1872-1397

Introduction

Identifying trends in Metropolitan Mortality (MM) is a challenging problem nowadays. Systematic approaches to discriminate the relevance of social, economic, demographic, educational, environmental or criminal factors in MM are matters of intense current research [1–3].

MM can be differentiated from rural mortality since the causes and risk factors are substantially different [3]. To develop a useful model that may explain the evolution in time and space of MM, is crucial to have a comprehensive information regarding the aforementioned factors. However, missing data is a common problem in developing and emerging countries, particularly for small-scale spatial and temporal level measurements. Diverse limitations and restricted access to public health data to collect information at the neighborhood and monthly level, lead to the modifiable areal unit problem (MAUP). It implies that a bias is generated affecting statistical hypothesis testing due to the combination of two or more spatial scales. To address this problem there have been several alternatives such as: correct the variance-covariance matrix using samples from individual-level data; focus on local spatial regression rather than global regression; design areal units to maximize a particular statistical result; developing statistics that change across scales in a predictable way; Bayesian hierarchical models as a general methodology for combining aggregated and individual-level data, or, using simulated data.

One of the most useful variables to explain MM is the age group, since the causes associated to risk factors, disease incidence and other several related variables are substantially different between such groups. Seminal works regarding age-specific MM can be found elsewhere [4–6]. Often, these age groups are separated into non-overlapping sets as follows:

1. Infant (x < 1 years old)

2. Pre-school (1 ≤ x < 4 years old)

3. School (4 ≤ x < 14 years old)

4. Productive (14 ≤ x < 64 years old)

5. Post-productive (x ≥ 64 years old)

In this work, we refer to the productive group as the economically productive population. Post-productive group is referred as economically dependent from the productive population. In addition, age-specific mortality can be influenced by social determinants at the individual level and at national or state level [7, 8]. Individual mortality is influenced by personal-level characteristics such as genetics, socioeconomic status and education [9].

However, there is growing evidence that collective or regional disadvantages can also be good predictors of individual mortality and population level mortality [1, 10]. Example of these regional disparities can be classified into social, demographic and environmental factors. In a recent study, Gavourova and Toth (2019) described how environmental factors influence changes in preventable mortality and how they impact differ from district to district in Slovakia [2]. In another context, individual and regional characteristics have been simultaneously analyzed for cardiovascular disease, to elucidate the effects of their interaction with air pollution, psychosocial stress, adverse childhood experiences and neighborhood deprivation index [11, 12].

Regarding the description of spatial mortality trends, two main approaches are used: i) All-cause or ii) Cause-specific mortality studies. In the first case, spatial and temporal variability have been measured, using different levels of granularity [13–17]. For cause-specific mortality several models have been developed for diseases such as cancer, diabetes, hypertension, chronic obstructive pulmonary disease, cardiovascular disease, hepatitis C and HIV/AIDS [6, 14, 17–26]. As an example of cause-specific mortality, Dwyer-Lindgren, et al. (2017), studied the variations in life expectancy, mortality rates and years of life lost from 152 causes of death at the county and neighborhood levels by age group and sex. They mainly conclude that county level estimations mask important local differences (between neighborhoods) [15, 27]. Studies with these characteristics have been performed in developed countries which generally have large, accessible and almost complete databases, thus, there is no MAUP.

In developing and emerging countries, it has been described an exacerbated problem in health disparities due to poverty, environmental threats, inadequate access to health care and educational inequalities that may lead to higher mortality rates [4]. In most cases, the conducted studies about mortality rates and other health outcomes rely on data that is not documented or incomplete. The last can be due to the lack of economic resources from public institutions or data restrictions for research purposes.

Mexico City, is one of the most populated cities worldwide. Its economy, employment rates and health services have improved over the last decades. Unfortunately, these resources are unevenly distributed, mainly affecting three risk groups: children, the elderly and the poor [28]. To worsen this situation, exposure to environmental risk factors derived from urbanization has increased, with its associated negative health effects gaining attention in recent years [29–32].

To set some context of mortality causes and trends in Mexico City, in a recent study, Aburto et al. (2018) analyzed lifespan and preventable mortality for Mexico City and the other 31 states of the country. They described changes in age groups. Particularly, they reported an increase in Diabetes and heart diseased-related mortality rates [33]. Furthermore, beyond diseased related mortality causes, an important number of deaths in Mexico City corresponds to homicides. Aburto et al. conducted different studies for the periods 2000-2010 and 2005-2015 regarding the impact of homicides in life expectancy and lifespan inequality. The authors reported an increase in homicide mortality that surpassed positive outcomes from health care reforms oriented to promote life expectancy at the national level [34, 35].

Regarding health outcomes, Gómez-Dantés, et al. (2016) report the leading risk factors for children and adults, supporting the fact that public health reforms and interventions should vary according to the specific risk factor for each age group. To give some examples of the latter, diarrhoeal diseases, undernutrition and poor sanitation were the leading risk factors for children, meanwhile for adults, chronic diseases and violence resulted the highest risk factors [36].

Despite the above mentioned previous background, at the borough and smaller scales, health, socioeconomic, educational or environmental disparities for Mexico City have not been formally quantified. In general, quality of life indexes of the urban and rural blocks are drastically different. However, these differences may be larger between urban neighborhoods than between rural and urban blocks.

Recently, new estimations suggested a high degree of social backwardness (term coined to refer to the lack of advancement of a group relative to the average) [37]. In Mexico City, health disparities such as illness and mortality rates can vary between close neighborhoods, just a few blocks apart, within the same borough. Taking this into account, assessment of health outcomes such as MM, using borough as a measurement unit, may result inadequate. A similar behavior arises from considering the temporal component: coarse-grained measurements may mask the variant behavior of MM. Once again, the MAUP is a matter to take into consideration at both spatial and temporal scales.

It is known that temporal variations in environmental factors (air pollutants for example) and socioeconomic variations (socio-economic level among boroughs and neighborhoods) can lead to biased results when not having measurements at the appropriate scales. Thus, there is an urgent need to count on fine-grained data in order to better interpret spatio-temporal correlations in health outcomes. To the best of our knowledge, there is no previous report in which time lapses for measurements were taken into smaller scales than in a year’s time. Nevertheless, both spatio-temporal data are available in Mexico City for MM. These data include, at the borough level, information for general mortality, age, gender and other descriptive variables from 2000 to 2016 in a year basis.

In this work, using available data for Mexico City, a two-level spatial and temporal description for MM was analyzed using geo-statistical interpolation. In this context, the spatial component is tackled at borough and neighborhood scale, whereas the time scale was used in a year and trimester levels. Thus, using a multilevel modeling for longitudinal data provides a finer description of these phenomena that may result in a more accurate spatio-temporal explanation of MM in Mexico City. This in turn may allow to better capture the multi-scale complexity of mortality patterns in large urban areas such as Mexico City aiming at improving integrative policy designs.

Materials and methods

Regarding the data collection in Mexico, the National Institute of Statistics, Geography (INEGI, Instituto Nacional de Estadística y Geografía) [38], each ten years develops a country-level census of population. In that national survey, economic, demographic and social data of all citizens in Mexico are collected and stored in a public database. Hence, INEGI provides the core of socioeconomic and demographic data, used as input for those databases that we used in this work. In addition, the mortality data was obtained from the Secretariat of Health in Mexico City (Secretaría de Salud de la Ciudad de México: SEDESA) [39]. Detailed description of both databases can be found in the following subsections.

Study area

The study area is the capital of Mexico, Mexico City which was known as the Federal District (Distrito Federal) until 2016 (Fig 1A). Mexico City belongs to the Valley of Mexico Metropolitan Area, the biggest metropolis in the central region of the country. The city is divided into 16 administrative boroughs (municipalities) and 2, 097 neighborhoods, according to the geo-spatial information, i.e., polygon shapefiles, obtained from the Geostatistic framework, December 2018 (Marco Geoestadístico, Diciembre 2018) [40].

[Figure omitted. See PDF.]

Fig 1. Mexico City study area.

A) The map of Mexico shows the location of Mexico City (in red), formerly known as Distrito Federal, one of the 32 states of Mexico which is located in the central area. B) Borough level description of the 16 cases (in numbers) corresponding to: 1. Álvaro Obregón, 2. Azcapotzalco, 3. Benito Juárez, 4. Coyoacán, 5. Cuajimalpa de Morelos, 6. Cuauhtémoc, 7. Gustavo A. Madero, 8. Iztacalco, 9. Iztapalapa, 10. La Magdalena Contreras, 11. Miguel Hidalgo, 12. Milpa Alta, 13. Tláhuac, 14. Tlalpan, 15. Venustiano Carranza and 16. Xochimilco. C) Neighborhood level description. Blue lines describe borough limits, whereas white lines at neighborhood areas (in gray). Red dots depict borough calculated centroids. In panels A and B, the scale bar and north arrow are also included. Notice that some boroughs have a dense neighborhood description in comparison.

https://doi.org/10.1371/journal.pone.0244384.g001

The neighborhood polygons were combined by borough identifiers to create the corresponding borough regions (Fig 1B), in order to have a multilevel hierarchical description of Mexico city. Borough centroids, were obtained using rgeos R package [41] (Fig 1C). In addition, neighborhood centroids were obtained by gCentroid or gPointOnSurface for convex or concave polygons respectively. All maps were created using R-software libraries: sf [42], rgeos [41], raster [43], geosphere [44], spacetime [45], sp [46], rgdal [47], ggplot2 [48], cowplot [49], gridExtra [50] and ggspatial [51].

Mortality database

A mortality database was created from available open data obtained at the Secretariat of Health in Mexico City (Secretaría de Salud de la Ciudad de México: SEDESA) [39]. The database comprehends a 16 borough, age-specific mortality rates follow up complete records, from the year 2000 up to 2016 and can be downloaded from [52].

The age-specific groups include the following five non-overlapping age descriptors:

1. Infant (x < 1 years old)

2. Pre-school (1 ≤ x < 4 years old)

3. School (4 ≤ x < 14 years old)

4. Productive (14 ≤ x < 64 years old)

5. Post-productive (x ≥ 64 years old)

In this work, we refer to the productive group as the economically productive population. Post-productive group is referred as economically dependent from the productive population. The mortality rates are expressed, in all cases, per 1, 000 inhabitants, but the Infant rate, which was calculated per 1, 000 liveborn according to the data obtained from INEGI and the Mexican National Population Council (Consejo Nacional de Población, CONAPO) [38, 52].

A five year time-step borough level evolution descriptive panel for every age-specific group is presented in Fig 2. Rows stand for age-specific mortality rate, and columns for a summarized time period (year) as reported by SEDESA. Naturally, with the available data at hand, the only alternative to visually explore the spatial data component, is to homogeneously color each borough with a single mortality rate for the reported year. Fortunately, there is discontinuity among boroughs, i.e., adjacent neighbor at using borough as area unit. Unfortunately, in this representation, neither time evolution nor spatial mortality age-specific rates are represented.

[Figure omitted. See PDF.]

Fig 2. Spatial age-specific mortality rates in Mexico City’s boroughs.

Each row corresponds to a particular age-specific mortality rate, i.e., Post-productive (x ≥ 64 years old), Productive (14 ≤ x < 64 years old), School (4 ≤ x < 14 years old), Pre-school (1 ≤ x < 4 years old) and Infant (x < 1 years old). Each column stands for the selected years 2000, 2005, 2010 and 2015 from the total yearly available period 2000 to 2016. Mexico City boroughs are treated as a the unit area and color coded according to the corresponding mortality rate, which make them comparable by row. Interestingly, notice the different mortality rate ranges (color bars) depending on the age-specific group of analysis. Polygon shapefiles files can be freely downloaded at INEGI’s website [40], whereas mortality data from SEDESA [52].

https://doi.org/10.1371/journal.pone.0244384.g002

To overcome this drawback, a visual description at borough level is presented in Fig 3. It can be seen that there are different age-specific time evolution patterns (linear trend, in blue), such as a decreasing mortality rate for Infant and Post-productive groups, in contrast to the increasing trend in the Productive case. Although the Cuauhtémoc borough (in green) is the smallest one in terms of its neighborhoods (38 in total), it seems to overcome the mortality rates for most age-specific cases in comparison to the rest of the boroughs. No apparent time-dependent correlation can be observed for School and Pre-school, where several borough vibrant mortality curves are presented.

[Figure omitted. See PDF.]

Fig 3. Temporal age-specific mortality rates in Mexico City’s boroughs.

Each panel corresponds to an age-specific mortality rate, i.e., Post-productive (x ≥ 64 years old), Productive (14 ≤ x < 64 years old), School (4 ≤ x < 14 years old), Pre-school (1 ≤ x < 4 years old) and Infant (x < 1 years old). All panels include complete age-specific mortality rates records by year from 2000 up to 2016 (in colour points) from SEDESA [52]. Color points stand for one of the 16 boroughs in Mexico City longitudinal mortality rate data. Borough data have been joined by their corresponding color dashed lines, whereas complete age-specific panel has been modeled by a linear regression (intercept and slope, blue line) with its respective standard deviation (grey area), to get a clear picture of the time evolution pattern. Interestingly, as the years pass, the mortality rate time evolution seems to diminish for Infant and Post-productive groups, whereas the Productive age-specific counterpart tends to increase.

https://doi.org/10.1371/journal.pone.0244384.g003

Spatio-temporal interpolation

In order to have a more robust input dataset to analyze trends in MM, the measured variables must have a temporal component as fine-grained as possible, since it is well known that several risk factors associated with mortality have a cyclic behavior, grounded on temperature, air pollution, seasonal pathogens or even individual social aspects [53–55]. Several studies on MM have taken into account yearly data to associate certain variables with the response outcome [10, 20–22]. Often, monthly data are not available for developing countries; in other cases, a great number of missing data results very common.

In general, although different area-level health outcomes might share a variety of explanatory variables such as socioeconomic, pollution, delinquency levels, health access, among others; by down scaling, spatio-temporal heterogeneity might arise [20, 21, 56–58]. To overcome the MAUP, a variety of geo-statistical procedures have been implemented to estimate mortality rates at different granularity levels. For instance, Population-Weighted Average, local and global Empirical Bayes and Poisson Kriging have been used to estimate disease-specific mortality rates from age-adjusted data. Accounting for spatial correlation patterns for low and high frequency rates, Poisson Kriging have shown better results [59, 60]. Kriging methods allow to estimate spatial risk considering heterogeneity among small areas from poorly reported databases [61].

In this context, area-to-point (AtP) Kriging provides instantaneous estimation of the spatial regression, which is valid for each time point, i.e., it is appropriate for a cross-sectional study, nonetheless, we have both longitudinal and spatial data points. In addition, whenever a time-wise progression estimate is needed, AtP falls short of considering spatio-temporal correlations unless an appropriate time-regression procedure is considered, i.e., one compliant with Gauss-Markov theorem. In this case, a sum-integrated joint method—minimizing error via least squares regression—will give rise to a valid measure; formally a best linear unbiased estimator. Doing this, it will however be formally equivalent to a joint (sumMetric) integrated spatio-temporal estimator.

To overcome data scarcity, spatio-temporal Kriging has been used to jointly interpolate spatial missing data at the county level as well as for temporal interpolation in mortality data [6, 13, 62]. Kriging estimates has been also used to describe changes in child mortality trends, evaluating between-countries and within-country sources of variation [5]. However, in current literature, a standard spatio-temporal methodology that properly addresses heterogeneity and data scarcity, has not been reported yet. The spatio-temporal description of the mortality database used in this work has been constructed on a yearly basis at the borough level.

Although, both spatial and temporal description levels can be accurate for open data summarized statistics, such description fails to represent the actual statistical metric unit, e.g., trimester measurements at the neighborhood level. In addition, it is known that in this context, data points have a spatio-temporal correlation which can be used to estimate values at unmeasured locations in space and time. Among the different geo-statistical interpolation methodological alternatives, the Kriging family provide unbiased estimates that minimize the variance-covariance spatio-temporal structure using an appropriate variogram modeling to tackle the MAUP [63]. Under this methodology, the choice of the variance-covariance spatio-temporal structure is not trivial and requires prior knowledge. As far as the authors know, we could not find in the literature reported mortality variance-covariance spatio-temporal structure. Thus, we tackle this issue by parsimony principle using computational force, in order to find the model structure with the least possible error among different mortality variance-covariance spatio-temporal structure available at gstat [64, 65] library: metric, separable, productSum, sumMetric and simpleSumMetric.

Briefly, a second data description level in both space and time, using neighborhood centroid locations and a trimester time basis, was carried out using Global Ordinary Kriging, for each age-specific mortality rate. All estimations were obtained using gstat R package implementation [64, 65]. Initial variogram parameter values were obtained from the complete mortality empirical (s)patio-(t)emporal rates variogram matrix γ(s, t) (Table 1 in S1 File):

* Nugget: The median value of the first three empirical variogram matrix row/column means, for the spatial or temporal initial guess respectively.

* Range: The spatial range is one third of the lagged maximum spatial value; for the temporal case, it corresponds to the maximum value.

* Sill: The median value of the last five empirical variogram matrix row/column means, for the spatial or temporal initial guess respectively.

* stAni: The spatio-temporal anisotropy was estimated using a linear model as implemented in gstat R package [64, 65].

* Joint spatio-temporal initial values, are based on the mean of the independent spatial and temporal values, respectively.

Using the initial variogram parameters, different spatial, temporal or joint covariance structures were tested to find the best parsimonious correlated data description, according to the available implementations in gstat (metric, separable, productSum, sumMetric and simpleSumMetric) [64, 65]. All possible single, double or triple variogram combinations (Exponential, Gaussian and Spherical) were tested according to the corresponding covariance structure. Hence, computational power was used in order to find the best variance-covariance spatio-temporal structure using the appropriate variogram in order to use the parsimony principle. All covariance models were fitted using a quasi-Newton box constrained method, where only the lower-bound was set to 0.001 for every parameter. The upper-bound was left to its default value (infinite), i. e., no box constraint was imposed for the maximum value. The selection criterion to choose a covariance structure model, was to minimize the weighted mean squared error. Complete age-specific mortality tested covariance fitted model results can be found in Table 2 in S1 File. Finally, the spatio-temporal mortality rate interpolation was obtained using all the available data points, under the best covariance model description, for the neighborhood centroids at a trimester level, in the period 2000 to 2016. Under this configuration the complete estimation root mean square error for the space-time Kriging data points was 8.16 × 10−10 (see Table 3 in S1 File for age-specific errors).

Time-evolution modeling

In the literature, there are many alternatives to assess the mortality time-evolution patterns using different types of models such as multilevel modeling for longitudinal data, longitudinal multilevel model, and longitudinal hierarchical linear model, among others. All of them target the different ways to model the variance-covariance structure, taking advantage of the inter/intra data structure and/or modifiable area unit.

Raudenbush & Bryk 1986 [66] took advantage of observed unit (individuals) when considering inter and intra school effects. Hence, their proposal can bee seen as a hierarchical or multilevel approach, where first, the within-group model is estimated by a separate regression equation for each school (referenced as Eq 1 in the original work). Then, the between-group model uses the regression coefficients as dependent output and try to model the within-school structural relationships (Eq 2 in [66]). Finally, the Raudenbush & Bryk can be coupled into a single equation by substitution of Eq 2 into Eq 1. The resulting equation allows to model the error term with many degrees of freedoms associated with the studied variables of interest.

However, the model estimation itself is not explained in detail. On the other hand, Kwok et al. 2008 [67] modeled variable time data points, but, they were focused on time correlation, hence, they introduce first-order auto-regression structure (AR(1)) to model the variance-covariance structure using SPSS (MIXED) and SAS (PROC MIXED) procedures. The two previous procedures are in essence, Linear Mixed-Effects Models. Recently, Anaya & Al-Delaimy 2017 [68] and Green et al. 2019 [69] did implement multilevel modeling for longitudinal data using linear mixed-effects models with R software [70] using the lme4 package.

In order to assess the mortality contribution of the age-specific group, the spatial unit area (borough or neighborhood) and polynomial time evolution, we used a linear mixed-effects model to account for data constraints and lack of error independence using the definitions of Eqs (1)–(4) [71]. The model was specified using Infostat software version 2018, which is an R [70] front-end as follows [72]: (1) (2) (3) (4)where, yijk is the mortality rate for the i-th age-specific mortality group (αi), at the j-th borough or neighborhood level (βj), for the k-th time (tk); μ is the global mean; γ, δ, ζ are the corresponding third order time polynomial coefficients; in addition, the double and triple complete fixed effects model interactions and the error term εijk were also specified. Indeed, the error term εijk in Eq (1) is modeled in Eq (2) by a two-level model, to account for the lack of independent errors; where, denotes the unobserved time effect and νijk is the idiosyncratic error term. The correlated errors were tackled using a first order auto-regressive model for time dependence as described by Eq (3), where is the corresponding coefficient and the individual effect. Finally, the heteroscedasticity in Eq (4) was modeled as a multiplicative effect of the residual variance σ2 times the variance error function g(.) using a varIdent definition of the different age-specific mortality groups αi [72].

The model was fitted using the R language with the nlme package under restricted estimation of the maximum likelihood [70, 73]. When possible, back-step model selection strategy was applied to remove the least significant fixed-effect term, one at the time, until no difference was found using a maximum likelihood test between competitor models. Type III sum of squares was used to assess an Analysis of Variance (ANOVA) table for marginal hypothesis tests for the fixed effects. Posterior Fisher’s Least Significance Difference (LSD) tests were applied over statistical significant terms, using a multiple comparison Bonferroni p-value correction. When possible, bilateral test was used and Type I error was set to 0.05.

Results

Spatio-temporal variogram estimation

Regarding spatio-temporal age-specific mortality estimations, the initial guesses obtained from the sample spatio-temporal variograms are shown in Table 1 in S1 File (see Material and methods section). Depending on the age-specific mortality group, the initial guesses are different for the nugget, range, sill and spatio-temporal anisotropy (stAni).

In this context, the nugget is the model intercept attributable to the smallest error measurements or spatial sources of variation. Interestingly, these sources of variations are negligible for Pre-school and School, in contrast to the other age-specific groups, with a wide range of nugget values (0.01 − 21.31). In addition, the correlation extends between measurements, also known as range, in all cases is exactly the same for all age-specific groups and last about 12 years for as far as 13.27 km.

The variogram values obtained at the range, a.k.a. the sill, is as close to the nugget for the School and Pre-school age-specific groups, and departs from it at most to double its value for Infant group. The anisotropy remains the same for Post-productive, School and Infant groups, but differ for Productive and Pre-school counterparts. Final covariance model weighted mean square error for all the tested variogram permutations can be found in Table 2 in S1 File.

It is worth to mention that the lowest error for the different covariance structure methods included the metric case for School and Pre-school; sumMetric for Productive and Infant and simpleSumMetric for Post-productive was the best one. Within these covariance models, there was no apparent pattern in the winner variogram model tested permutation (temporal + spatial + joint). The Gaussian + Gaussian (Gau + Gau) was the winner’s choice of temporal and spatial in the Post-productive age-specific group.

Moreover, Infant mortality also followed this pattern, with an additional Gau component for the joint variogram. The Exponential + Gaussian (Exp + Gau) was the winning choice for Productive mortality and Gau joint variogram. Indeed, this was a data-driven approach that required to explore the complete permutation grid, in order to reach a parsimonious spatio-temporal correlation model. A visual comparison of each winner covariance model and sample variogram can be found in File 1 in S1 File. Finally, under these configurations the complete root mean square error estimation for the space-time Kriging data points was 8.16 × 10−10 (see Table 3 in S1 File for age-specific errors).

Two-level mortality rate spatial description in Mexico City

Productive mortality rate.

The spatio-temporal mortality description in Mexico City starts at the raw data presented in Fig 2. Let us consider the whole picture, using, for example, productive mortality rate as presented in Fig 4. Using the borough spatio-temporal granularity top panel in Fig 4, there is a clear spatial pattern, no matter the selected year (2000, 2005, 2010 or 2015). It seems that there is a global maximum mortality rate value (hotspot) at Cuauhtémoc borough (number six in Fig 1C). There, the mortality rate values radially decrease as long as we depart from this location to the outer boundaries of Mexico City. Indeed, the decrease is not homogeneous, i. e., it is dictated by a spatial anisotropy where the north and northeast direction have a less marked decline in comparison to the south and south-west direction. Moreover, at this temporal description, there is also an increment in the global maximum rates as we move from 2000 to 2015, in agreement with the temporal patterns presented in Fig 3 for the Global case.

[Figure omitted. See PDF.]

Fig 4. Three levels of granularity for productive mortality rate in Mexico City.

First one, productive mortality rate original data description at borough level (panel A). Second one, productive mortality rate kriged data description at the neighborhood level is presented in panel B. In panel C, it is presented a zoom-in at a second description level for Cuauhtémoc borough. Interestingly, this is the borough with the highest productive mortality rate no matter the selected year (2000, 2005, 2010 or 2015), according to panel A and B (central borough in red). However, the mortality rate is not homogeneous at the neighborhood level, as depicted by the kriged values presented in panel C for the different years. Polygon shapefiles files can be freely downloaded at INEGI’s website [40].

https://doi.org/10.1371/journal.pone.0244384.g004

Moving towards a deeper data exploration, the spatio-temporal interpolation obtained by Ordinary Kriging provides a productive mortality rate smooth surface (see neighborhood level at Fig 4). For a fair comparison, the same time scale was used (years), whereas the spatial description considered the neighborhood centroids grid. With this zoom-in, a more realistic geographical-continuous mortality transitions can be observed, unlike the discrete phenomena at borough level, for adjacent boroughs, presented at borough level in Fig 4. The latter reflects a far less abrupt change in mortality rates from one borough to another, thus giving a continuity between neighborhoods.

Furthermore, the productive mortality rates increase throughout neighborhoods and boroughs consistently from 2000 to 2015 in both spatial-scales. Also, notice that the mortality rate scale bar has also decreased from a maximum of 6 (at the borough level) into a 4.2 when we move towards into the fine grain spatial description for this age-specific group.

If we zoom-in even further, we can see the neighborhoods at a single borough, e.g. Cuauhtémoc borough Fig 4, where we can distinguish neighborhood level trends. Here, Cuauhtémoc borough was selected, since it has the highest mortality rate. Although apparently imperceptible, there are distinguishable differences in mortality values among Cuauhtémoc neighborhoods. Similar differences are obtained for the rest of the boroughs as seen at borough level description in Fig 4.

School mortality rate.

Analogously to Fig 4, in Fig 5 we may observe the mortality rate through 15 years in the school group in three different levels of granularity. As it can be appreciated in the figure, mortality patterns in the school age are not homogeneously distributed in the 16 boroughs of Mexico City. Instead, each borough has its own pattern.

[Figure omitted. See PDF.]

Fig 5. Three levels of granularity for school mortality rate in Mexico City.

First, school mortality rate original data description is presented at borough level (panel A). Second, school mortality rate kriged data description at the neighborhood level is presented at panel B. In panel C, it is presented a zoom-in at the neighborhood level for the Cuauhtémoc borough. Interestingly, the school mortality rate for this borough increments from the year 2000 to 2005, while it reduces from 2005 to 2015 according to panels A and B, whereas the city’s southeast area shows an increment from 2005 to 2015. However, the mortality rate is neither homogeneous nor constant at the neighborhood level, as depicted by the kriged values presented in panel C. Polygon shapefiles files can be freely downloaded at INEGI’s website [40].

https://doi.org/10.1371/journal.pone.0244384.g005

Unlike the previous case regarding productive age, in the school group the spatial trend is not clear. However, a visible and measurable decrease is observed from 2010 to 2015 in practically the whole city.

It is worth noticing that the scales in both figures are different. In Fig 4 the upper value for mortality rate is close to 6, meanwhile for school age, the top value mortality rate is around 0.4.

The comparison between those different groups in Cuauhtémoc borough and its neighborhoods is also remarkable. Meanwhile, for school group, 2015 is the year with the lowest mortality rate, that year was the highest for the productive age. Additionally, the differences in mortality rate at the neighborhood level is more visible in the productive group than the schoool one.

Similar two-level mortality rate spatial description in Mexico City can be found for the rest of the age-specific groups (Infant, Pre-school and Post-productive) in Figs 1-3 in S2 File respectively.

Time-evolution mortality rate modeling in Mexico City

So far, we have filled the spatio-temporal gaps for the missing data points for both scales within the different mortality age-specific groups. Now, that we have overcome the summary time period for the mortality rate constraint, which is usually in a year base, we can study temporal/seasonal trends within the temporal range (12 years). In addition, we can explore the between mortality rate variability at borough level or even within its borough, a.k.a. neighborhood level data description. The later is a must, in order to address the MAUP, as we are changing the spatial scale, from several kilometers long (borough level), to a very different area unit scale at neighborhood (couple of blocks), where the prior estimated spatial correlation range is as far as 13.27 km. Finally, the mortality rate itself can be decomposed into both fixed-effects and random variance-covariance structure contribution using Eqs (1)–(4), according to the data level description as follows.

Borough level mortality rate contribution.

Borough-level results are presented in Table 1. The ANOVA results showed that the only non-significant effect at the borough level, is the triple interaction that includes the time to the third power (p = 0.22). Hence, the mortality rate time-evolution pattern in Fig 3 can be parsimoniously captured by our methodological proposal. In addition, the auto-correlation parameter had an impact not as high as one would expect at borough level (ϕ = 0.12). On the other hand, the variance function did address the different age-specific groups where Infant was the one with the highest value (11.59) followed by Post-productive (10.34), Pre-school (0.85), Productive (0.82) and School (0.43).

[Figure omitted. See PDF.]

Table 1. Analysis of the variance at two spatial data description level.

https://doi.org/10.1371/journal.pone.0244384.t001

To further explore the mortality behavior, posterior Fisher’s LSD tests results were obtained over the age-specific, αi, Mexico City’s borough, βj, and double borough times age-specific interaction terms αi × βj (Fig 6A–6C, respectively). The first remarkable result is the mortality contribution to the main effects αi and βj. Mexico City’s model results evidence, for the time-period studied, suggests that the age-specific term (αi) is the main responsible for the mortality rate when compared to the spatial borough contribution (βj). Indeed, the mortality rate means are definitely empowered by age-specific groups (upper bounded at 50 [x1000]) rather than spatial borough locations, which is upper bounded at 17[x1000].

[Figure omitted. See PDF.]

Fig 6. Posterior mortality rate test results at borough level for Mexico City.

Fisher’s Least Significant Difference (LSD) tests were performed over the estimated age-specific mean mortality rate, according to model description in Eqs ((1)–(4)). In panels, Mexico City Fisher’s LSD test results for: A) Age-specific mortality, B) Borough level, and, C) Borough × Age-specific interaction. The LSD group mean test results are presented as bar-plots with their corresponding standard deviation bar and mean group letters (A, B, C …), where bars that share at least a single letter, are not statistically different after Bonferroni multiple test correction (p > 0.05). Interestingly, mortality rates is mainly composed by age-specific contribution when compared to Borough impact. In addition, Borough × Age-specific interaction retains the mortality age-specific rate pattern, but, it is modulated by the borough contribution.

https://doi.org/10.1371/journal.pone.0244384.g006

Secondly, the age-specific mortality rate groups do not overlap between each other, due to the different Fisher’s LSD letters (A-E) in Fig 6A. Moreover, the Post-productive group (A) outruns any other group, but also doubles its following competitor, the Infant group (B). Interestingly, this trend also remains when we compare two consecutive groups, i.e., B vs. C, C vs D and so on. In this age-specific context, the school is the group with letter E, i.e, is the one with the least mortality rate. Interestingly, previous results could be, to some extend, the explanation that models the Mexico’s actual population pyramid shape data [74]—pencil like figure, i.e., long bar with sharp-pointed end.

Here, the Post-productive people are the most underrepresented group with a sharp-pointed shape at the tip of the population pyramid, due to the high mortality rate evidence of our borough level model results. Following this rationale, the second highest mortality rate group is B (Infant), which could be probably the reason why the typical pyramid shape is broken into a uniform/bar shape shared by the Infant, Pre-school, School and Productive age-specific mortality groups, influenced by their low mortality rate.

Thirdly, taking into consideration Fisher’s LSD results for the 16 borough’s mortality rate contribution βj (Fig 6B), it can be seen three important aspects: i) Among group letters from A to I, Cuauhtémoc (A) and Tlalpan (I) are the boroughs with the highest and lowest model estimated mean values respectively; ii) The LSD mortality rate group mean letters cluster up to five boroughs per cluster, i.e., borough that share a single letter, are not statistically different after Bonferroni multiple test correction (p > 0.05), thus, belong to the same cluster; iii) The borough Fisher’s LSD cluster structure is correlated to spatial proximity.

Fourthly, moving towards αi × βj group results (Fig 6C), i.e., age-specific times borough model’s interaction LSD group means, it can be seen that the original and predominant age-specific mortality rate contribution in panel A is modulated by the borough contribution of panel B. Thus, if we average LSD group mean results in panel C, by age-specific groups, we should return into panel’s A results. At this level, we are dealing at “between” borough mortality rate level description, where the number of Fisher’s LSD group letters has increased proportional to the number of possible αi × βj levels.

Neighborhood level mortality rate contribution.

In order to model “within” borough mortality rate data description, we need to change our attention into another spatial scale representation, a.k.a. neighborhood level. If we pick the borough with the highest mortality rate (Cuauhtémoc), we found that the model presented in Eqs (1)–(4) is not well-suited for this level description data. Hence, a model selection process was considered in order to obtain the best parsimonious data description.

The current model results are presented in Table 1, where it can be seen that some cells are empty due to terms discarded from the analysis. Interestingly, the third and second order time triple interaction terms have been excluded from the analysis, in addition to the borough times time to the third power. Hence, in this context the modeling complexity has been reduced, at the expense of a higher autocorrelation (ϕ = 0.88) and different variance function parameters and ranking, i. e., Productive (3.80), Post-productive (2.43), Pre-school (2.22), School (1.35) and Infant (0.04).

The Cuauhtémoc borough results can be found in Fig 7. Its 38 neighborhoods were numbered according to the high mortality rate downwards (Fig 7A). This result is complemented by Fig 7B, where the Fisher’s LSD means were used to describe the mortality rate landscape. Interestingly, at the neighborhood level, there is also a radial mortality rate decay starting at the central neighborhood with number one and the letter A, which correspond to Tabacalera neighborhood. Conversely, the two neighborhoods with the lowest mortality rate (letter D) are situated at the opposite borough extremes—south-west (Hipódromo de la Condesa) and northeast (Valle Gómez) borders. Moreover, the age-specific mortality patterns have changed from Post-productive and Infant mortality at the borough description in Fig 6A into Post-productive and Pre-school as seen in Fig 7C. Moreover, the Fisher’s LSD test results in this borough, but at the neighborhood level, showed some regions with up to four possible over-lapping letters (Fig 7D).

[Figure omitted. See PDF.]

Fig 7. Posterior neighborhood mortality test results zoom-in at Cuauhtémoc borough in Mexico City.

A) Cuauhtémoc borough is divided into its 34 neighborhoods. Numbers are ordered from the highest (1) to the lowest (34) model estimated neighborhood fixed effect mortality mean according to Eqs (1)–(4) description: 1. Tabacalera, 2. Centro, 3. Juárez, 4. Doctores, 5. Buenavista, 6. Guerrero, 7. San Rafael, 8. Roma Norte, 9. Obrera, 10. Cuauhtémoc, 11. Santa María la Ribera, 12. Tránsito, 13. Centro Urbano Benito Juárez, 14. Esperanza, 15. Unidad Hab. Nonoalco Tlatelolco, 16. Morelos, 17. Buenos Aires, 18. Atlampa, 19. Vista Alegre, 20. Algarín, 21. Paulino Navarro, 22. Roma Sur, 23. Ex-hipódromo de Peralvillo, 24. San Simón Tolnáhuac, 25. Santa María Insurgentes, 26. Hipódromo, 27. Maza, 28. Ampl. Asturias, 29. Felipe Pescador, 30. Condesa, 31. Asturias, 32. Peralvillo, 33. Valle Gómez and 34. Hipódromo de la Condesa. B) Neighborhoods are filled according to the Fisher’s Least Significant Difference (LSD) group mean obtained at this level of representation. C) Age-specific Fisher’s LSD results within Cuauhtémoc borough. D) Fisher’s LSD neighborhood contribution. In all cases, capital letters stand for Fisher’s LSD groups, where bars that share at least a single letter, are not statistically different after Bonferroni multiple test correction (p > 0.05). Results are presented as ordered mean ± standard error estimation according to model description of Eqs (1)–(4). Notice that at this data level decomposition, the LSD group means are at the same mortality rate order, i.e., the age-specific group mean values are comparable to the different neighborhood contribution. In addition, age-specific mortality values do not overlap whereas, most of the neighborhoods within this specific borough, share a common LSD group.

https://doi.org/10.1371/journal.pone.0244384.g007

Taking about the different mortality rate model description, we have to keep in mind that Fig 6 panel B “between borough” mortality rate data order of magnitude, is going to be fine grained (“within borough”) modeled using the same framework depicted in Eqs (1)–(4). The first result to be discussed is that, unlike borough level, neighborhood mortality rate data does not cope with third order neighborhood interaction and the only triple significant interaction has a linear time tendency for age-specific times neighborhood mortality, as shown in Table 1.

The second aspect, is the new insight of the mortality rate at the neighborhood level. Here, the first mortality rate level description (at borough) leaves Cuauhtemoc’s borough near 15 [x1000]. Now, Fig 7 panels C and D, decomposed into a finer grain considering the same age-specific groups, but, now Cuauhtemoc’s neighbors are included, respectively. At this data level description, mortality age-specific groups within Cuauhtemoc are closer to each other (same magnitude order, units) unlike borough level (one magnitude order, tenths). In addition, between the neighborhood mortality rate is almost shared by all the 34 neighborhoods (above 7 and below 8), is we consider the shared letters of the LSD results in Fig 7 panel D. Finally, the rest of the mortality rate contribution to sum up to 15 (Cuauhtemoc’s borough mortality description), it is distributed upon the different neighborhood model terms of Eqs (1)–(4).

Discussion

Understanding metropolitan mortality is not only relevant in terms of a mere descriptive statistical approach, but, to take into account variables that could determine, with the highest accuracy, this crucial health outcome. This knowledge, in due time, could lead to the establishment of appropriate public policies to improve citizen life quality in metropolitan areas.

In the aforementioned terms, by means of the development of a modeling approach, based on a systematic interpolation of missing data, one may observe spatio-temporal dynamics of mortality in the urban areas with higher precision. Kriging family methods have proven to be useful to achieve this goal.

In this work we have demonstrated how the improvement of the granularity level at both spatial and temporal definition, could explain some of the socio-demographic variables underlying the changes in mortality rates between boroughs. This task was achieved by statistically interpolating those data points into a fine-grained level. In what follows, we will comment on the findings using spatio-temporal Kriging-based methods.

It can be easily noticed that Cuauhtémoc borough is the one with the highest mortality rate in all age groups, but for pre-school –which ranks in third place– as depicted in Fig 3. This borough has several particularities that should be commented, in order to unveil some hints towards plausible explanations, for the unique behavior observed there along time and space.

As observed in Fig 6C, Cuauhtémoc is the borough with the highest mortality mean value in all age groups, but in Pre-school, which shares places with Cuajimalpa and Iztapalapa boroughs. Except in that case, Cuauhtémoc has an outstanding mortality behavior over time. The case is particularly dramatic in the productive age (14 to 64 years old).

In addition, in the Productive time-evolution pattern in Fig 3 can be observed the different temporal behavior of Cuauhtémoc (solid green dashed line), compared to the rest of the boroughs. Cuauhtémoc is the economic center of Mexico City. There, the executive and legislative powers are placed, as well as the most important commerce hotspots. Indeed, Cuauhtémoc is the most densely populated borough of Mexico City.

In economic terms, Cuauhtémoc concentrates 4.6% of the gross domestic product of the entire country [75]. Around 5 million people pass through this borough every day, despite its population oscillates only around 500, 000 inhabitants. This is the place in Mexico City with the highest number of public transport stations. Cuauhtémoc also concentrates the largest markets of informal commerce of the city (Tepito Market, in the Morelos neighborhood). The high density and the flux of money and services may help to explain, to some extent, the different behavior of the mortality rate in the productive age in this borough.

Another point to take into account when mortality rate in Cuauhtémoc is observed, lies on high levels of insecurity for the aforementioned reasons regarding population density and economic concentration. Additionally, the Drug War launched at the end of 2006 by former President Felipe Calderón, affected mortality rates in a large part of the country [76], being the capital of the country also upset, in particular, the city downtown, Cuauhtémoc.

By taking into account the spatio-temporal Kriging, it has been possible to modeled mortality rates at the neighborhood level. With this spatio-temporal kriging model, the four places with the highest mortality rates were Tabacalera, Centro, Juárez and Doctores neighborhoods, part of Cuauhtémoc borough.

According to reports of the Executive Secretariat of the National System of Public Security System (Secretariado Ejecutivo del Sistema Nacional de Seguridad Pública, SESNSP), Centro is one of the most insecure boroughs of Mexico City and the homicide rate is the highest in the city [77]. By integrating transit and other accident-related issues, with employment determinants and crime-associated mortality, it may be possible to present an explanation for which the productive age mortality presents a consistent increase in Cuauhtémoc, compared to the rest of the boroughs.

Possible intrinsic biases due to possible errors in the registry must be taken into account. For example, the record of a death is registered once a certificate of death has been provided. In some cases, it could last days, depending on diverse factors. However, in the case of Mexico City, that under-registry is extremely low since social and health services are guaranteed in practically the whole area of the city.

The model presented here also has caveats, since it shows an interpolation of the coarse-grained data at the borough level of description. For instance, according to the model, Tabacalera is the neighborhood with the highest mortality rate, however, by looking at the data, the homicide rate as an example is not as large as Centro or Morelos [77].

At this stage is not possible to disambiguate whether these inconsistencies are due to the interpolating strategy or indeed reflect different causes of death, such as the ones related to environmental factors like air pollution and other contaminants. Despite this caveat, the model shows a radial decrease of mortality in Mexico City, starting from Cuauhtémoc downtown, similar to the general behavior observed at the borough level (Fig 2).

It is relevant to notice that some of the caveats and limitations just discussed are actually instances of a well known issue known as the modifiable areal unit problem (MAUP). MAUP states that a bias is generated affecting statistical hypothesis testing due to the combination of two or more spatial scales on a given geostatistical analysis. The reason is that data aggregation may become dependent on the choice of modifiable areal unit (MAU) used as a primary source in the analysis. MAUP will in this case induces statistical biases that may lead to a form of ecological fallacy. In the particular case we presented here, since Kriging analysis is a form of regression, there are some ways to solve the MAUP or at least alleviate some of its consequences. There are a number of methodological choices to do so. We decided to follow a mixed approach by incorporating several of these as follows:

1. One of the alternatives is to correct the variance-covariance matrix using samples from individual-level data

Indeed, to this first end, we have tackled the MAUP in a hierarchical manner:

* Spatio-temporal kriging itself was modeled first, selecting the best combination of spatio-temporal variance-covariance structure (metric, separable, productSum, sumMetric and simpleSumMetric) with the appropriate single, double or triple variogram combinations (Exponential, Gaussian and Spherical), to minimize the weighted mean squared error (see Tables 2 and 3 in S1 File). Once the complete variance-covariance has been fixed, the neighborhood spatio-temporal grid was used to obtain the kriged values.

* Over the original data, at borough level, a linear mixed-effect model using the definitions of Eqs (1)–(4). In these definitions, as in Kwok et al. 2008 [67], we did include a first-order autoregression (AR(1)) structure. In addition, data heteroscedasticity was modeled as var(εijk) = σ2 g2(αi), i.e., a multiplicative effect of the residual variance σ2 times the variance error function g(.) using a varIdent definition for the different age-specific mortality groups αi, as described in Material an Methods section.

* Over the kriged data, at the neighborhood level, also the model of Eqs (1)–(4) was fitted to cope with the data, i.e. considering a different MAU to meet our needs.

2. A second way to correct for the MAUP is to focus on local spatial regression rather than global regression

3. Here, we used a global regression for the kriging process, but, considering the best variance-covariance structure possible for our data, letting the model to adjust its contribution instead of using a fixed local spatio-temporal regression. The RMSE for each mortality age-specific group is presented in Table 3 in S1 File.

4. A third strategy to correct the MAUP bias is by resorting to design areal units to maximize a particular statistical result. In this manuscript, we first:

* Maximized the model likelihood while minimizing the weighted mean squared error when selecting the best spatio-temporal variance-covariance structure.

* Maximized the linear mixed-effect model likelihood for both borough and/or neighborhood data level description.

Related approximations to solve the MAUP have been discussed in the specialized literature for decades. For instance, Raudenbush & Bryk 1986 [66] took advantage of observed unit (individuals) when considering between and within group (in their case, schools) effects. Hence, their proposal can bee seen as a hierarchical or multilevel approach, where first, the within-group model is estimated by a separate regression equation for each school (referenced as Eq 1 in the original work). Then, the between-group model uses the regression coefficients as dependent output and try to model the within-school structural relationships (Eq 2 in [66]). Finally, the Raudenbush & Bryk can be coupled into a single equation by substitution of Eq 2 into Eq 1. The resulting equation allows to model the error term with many degrees of freedom associated with the studied variables of interest. However, the model estimation itself is not explained in detail.

Kwok et al. 2008 [67] on the other hand, resort to modeling variable time data points (not necessary at fixed intervals), unlike repeated-measurements Analysis of the Variance (ANOVA) or ANOVA polynomial trend analysis. They propose to model the variance structure by means of linear growth models using multilevel models. They start with a simple random intercept model and simple linear growth model, to show how the variance-covariance structure can tackle within-class correlation. Then, they move to models with time-invariant covariate and pseudo-R2 statistic to evaluate the model’s effectiveness. Finally, they present how to model covariance structure for the within-individual random errors were, the proximal autocorrelation in longitudinal data, was addressed using a first-order autoregression (AR(1)) structure.

Recent models by Anaya & Al-Delaimy 2017 [68] explicitly mention to use a “mixed effect linear model with random effect estimation and repeated measurements” in their work. Additionally, the hierarchical linear model in [78] and the multilevel modeling of longitudinal data in [79] were fitted using SPSS (MIXED) and SAS (PROC MIXED) procedures respectively. On the other hand, Green et al. 2019 [69] used linear mixed-effect models with R software [70] by means of lme4 package, to fit multilevel models with random intercepts for schools and random intercepts and random slopes for time at the individual level since students are nested within waves.

These and other methods approach the MAUP from different angles. It is clear, though, that no definite nor one-size-fits-all solution to the problem has been found. Different problems and datasets will probably call for tailored analytical perspectives determined by their specific constraints. At the same time, further theoretical developments to the MAUP may shed some light on what kind of features a generalistic approach may have. The spatio-temporal representativeness problem (one of which instances have been discussed here) may need to be re-visited whenever further advances on how to effectively tackle the MAUP are developed.

Conclusions

The granularity enhancement in mortality and health outcomes data will further improve its usage in a variety of public policies, such as urban development, security, creation of health centres, hospitals, public transport and even water re-usage. In these terms, efforts in gathering data at the lowest possible source (time and space), are highly desirable for future works.

Finally, in order to build predictive models that allow to tailor public policy design, the spatio-temporal dynamics of MM should incorporate socio-demographic, environmental, economical factors and relevant covariates, as detailed as possible. Future directions include, but are not limited to, utilizing environmental data, such as air pollution, quality of water, disposition of water supplies among others. Some of the above mentioned variables are actually available as open data resources. The ultimate understanding of MM by using the data presented here, as well as environmental and other risk factors, will help in the searching for the improvement of life quality of the metropolitan areas in Mexico City.

Supporting information

[Figure omitted. See PDF.]

S1 File. Variogram modeling for age-specific mortality rates as well as root mean square errors for kriged data.

https://doi.org/10.1371/journal.pone.0244384.s001

(PDF)

S2 File. Three levels of granularity for infant, pre-school and post-productive mortality rates in Mexico City.

https://doi.org/10.1371/journal.pone.0244384.s002

(PDF)

Acknowledgments

Authors want to thank to the Subdirección de Bioinformática at the Instituto Nacional de Medicina Genómica. KBL thanks to the Consejo Mexiquense de Ciencia y Tecnología. This work is part of her PhD Thesis.

Citation: Baca-López K, Fresno C, Espinal-Enríquez J, Flores-Merino MV, Camacho-López MA, Hernández-Lemus E (2021) Metropolitan age-specific mortality trends at borough and neighborhood level: The case of Mexico City. PLoS ONE 16(1): e0244384. https://doi.org/10.1371/journal.pone.0244384

References

1. Gaskin DJ, Roberts ET, Chan KS, McCleary R, Buttorff C, Delarmente BA. No man is an island: the impact of neighborhood disadvantage on mortality. Int J Environ Res Public Health. 2019;16(7):1265.

2. Gavurova B, Toth P. Preventable Mortality in Regions of Slovakia—Quantification of Regional Disparities and Investigation of the Impact of Environmental Factors. Int J Environ Res Public Health. 2019;16(8). pmid:30999607

3. Rodríguez-Sanz M, Gotsens M, dell?Olmo MM, Borrell C. Trends in mortality inequalities in an urban area: the influence of immigration. Int J Equity Health. 2019;18(1):37. pmid:30808362

4. Ayele DG, Zewotir TT. Childhood mortality spatial distribution in Ethiopia. J Appl Stat. 2016;43(15):2813–2828.

5. Burke M, Heft-Neal S, Bendavid E. Sources of variation in under-5 mortality across sub-Saharan Africa: a spatial analysis. Lancet Public Health. 2016;4(12):e936–e945. pmid:27793587

6. Tlou B, Sartorius B, Tanser F. Space-time patterns in maternal and mother mortality in a rural South African population with high HIV prevalence (2000–2014): results from a population-based cohort. BMC Public Health. 2017;17(1):543. pmid:28578674

7. Montez JK, Zajacova A, Hayward MD. Explaining inequalities in women’s mortality between US States. SSM Popul Health. 2016;2:561–571. pmid:27722192

8. Shin J, Cho K, Choi Y, Lee S, Park EC, Jang SI. Combined effect of individual and neighborhood socioeconomic status on mortality in patients with newly diagnosed dyslipidemia: A nationwide Korean cohort study from 2002 to 2013. Nutr Metab Cardiovasc Dis. 2016;26(3):207–215. pmid:26895648

9. Braveman P, Gottlieb L. The social determinants of health: it’s time to consider the causes of the causes Public health reports. 2014; 129, p. 19–31

10. Kindig DA, Cheng ER. Even as mortality fell in most US counties, female mortality nonetheless rose in 42.8 percent of counties from 1992 to 2006. Health Aff. 2013;32(3):451–458.

11. Hazlehurst M, Nurius P, Hajat A. Individual and Neighborhood Stressors, Air Pollution and Cardiovascular Disease. Int J Environ Res Public Health. 2018;15(3):472. pmid:29518012

12. Casas I, Delmelle E, Delmelle EC. Potential versus revealed access to care during a dengue fever outbreak. Journal of Transport & Health. 2017;4:18–29.

13. Vaughan AS, Kramer MR, Waller LA, Schieb LJ, Greer S, Casper M. Comparing methods of measuring geographic patterns in temporal trends: an application to county-level heart disease mortality in the United States, 1973 to 2010. Ann Epidemiol. 2015;25(5):329–335. pmid:25776848

14. Bethea TN, Palmer JR, Rosenberg L, Cozier YC. Neighborhood socioeconomic status in relation to all-cause, cancer, and cardiovascular mortality in the Black Women’s Health Study. Ethn Dis. 2016;26(2):157. pmid:27103765

15. Dwyer-Lindgren L, Stubbs RW, Bertozzi-Villa A, Morozoff C, Callender C, Finegold SB, et al. Variation in life expectancy and mortality by cause among neighbourhoods in King County, WA, USA, 1990–2014: a census tract-level analysis for the Global Burden of Disease Study 2015. Lancet Public Health. 2017;2(9):e400–e410. pmid:29253411

16. Manda SO, Abdelatif N. Smoothed temporal atlases of age-gender all-cause mortality in South Africa. Int J Environ Res Public Health. 2017;14(9):1072. pmid:28914783

17. Vaughan AS, Schieb L, Quick H, Kramer MR, Casper M. Before the here and now: What we can learn from variation in spatiotemporal patterns of changing heart disease mortality by age group, time period, and birth cohort. Soc Sci Med. 2018;217:97–105. pmid:30300762

18. Parise CA, Caggiano V. Regional variation in disparities in breast cancer specific mortality due to race/ethnicity, socioeconomic status, and urbanization. J Racial Ethn Health Disparities. 2017;4(4):706–717. pmid:27604380

19. Dwyer-Lindgren L, Bertozzi-Villa A, Stubbs RW, Morozoff C, Kutz MJ, Huynh C, et al. US county-level trends in mortality rates for major causes of death, 1980-2014. JAMA. 2016;316(22):2385–2401. pmid:27959996

20. Cho KH, Lee SG, Nam CM, Lee EJ, Jang SY, Lee SH, et al. Disparities in socioeconomic status and neighborhood characteristics affect all-cause mortality in patients with newly diagnosed hypertension in Korea: a nationwide cohort study, 2002–2013. Int J Equity Health. 2016;15(1):3. pmid:26743664

21. Cho KH, Nam CM, Lee EJ, Choi Y, Yoo KB, Lee SH, et al. Effects of individual and neighborhood socioeconomic status on the risk of all-cause mortality in chronic obstructive pulmonary disease: a nationwide population-based cohort study, 2002–2013. Respir Med. 2016;114:9–17. pmid:27109806

22. Li C, Hu S, Yu C. All-Cause and Cancer Mortality Trends in Macheng, China (1984–2013): An Age-Period-Cohort Analysis. Int J Environ Res Public Health. 2018;15(10):2068.

23. Ford MM, Desai PS, Maduro G, Laraque F. Neighborhood inequalities in hepatitis c mortality: spatial and temporal patterns and associated factors. J Urban Health. 2017;94(5):746–755. pmid:28623451

24. Sifuna P, Otieno L, Andagalu B, Oyieko J, Ogutu B, Singoei V, et al. A Spatiotemporal Analysis of HIV-Associated Mortality in Rural Western Kenya 2011–2015. J Acquir Immune Defic Syndr. 2018;78(5):483. pmid:29683994

25. Roth GA, Abate D, Abate KH, Abay SM, Abbafati C, Abbasi N, et al. Global, regional, and national age-sex-specific mortality for 282 causes of death in 195 countries and territories, 1980–2017: a systematic analysis for the Global Burden of Disease Study 2017. The Lancet. 2018;392(10159):1736–1788.

26. Rawla P, Sunkara T, Gaduputi V. Epidemiology of Pancreatic Cancer: Global Trends, Etiology and Risk Factors. World J Oncol. 2019;10(1):10. pmid:30834048

27. Dwyer-Lindgren L, Bertozzi-Villa A, Stubbs RW, Morozoff C, Mackenbach JP, van Lenthe FJ, et al. Inequalities in life expectancy among US counties, 1980 to 2014: temporal trends and key drivers. JAMA Intern Med. 2017;177(7):1003–1011. pmid:28492829

28. Ramos GCD, et al. Real Estate Industry as an Urban Growth Machine: A Review of the Political Economy and Political Ecology of Urban Space Production in Mexico City. Sustainability. 2019;11(7):1–24.

29. Cohen AJ, Brauer M, Burnett R, Anderson HR, Frostad J, Estep K, et al. Estimates and 25-year trends of the global burden of disease attributable to ambient air pollution: an analysis of data from the Global Burden of Diseases Study 2015. The Lancet. 2017;389(10082):1907–1918. pmid:28408086

30. Lelieveld J, Evans JS, Fnais M, Giannadaki D, Pozzer A. The contribution of outdoor air pollution sources to premature mortality on a global scale. Nature. 2015;525(7569):367–371. pmid:26381985

31. Bel G, Holst M. Evaluation of the impact of bus rapid transit on air pollution in Mexico City. Transport Policy. 2018;63:209–220.

32. Davis LW. Saturday driving restrictions fail to improve air quality in Mexico City. Scientific reports. 2017;7:41652. pmid:28151487

33. Aburto JM, Riffe T, Canudas-Romo V. Trends in avoidable mortality over the life course in Mexico, 1990–2015: a cross-sectional demographic analysis. BMJ Open. 2018;8(7):e022350. pmid:30068622

34. Aburto JM, Beltrán-Sánchez H, García-Guerrero VM, Canudas-Romo V. Homicides in Mexico reversed life expectancy gains for men and slowed them for women, 2000–10. Health Affairs. 2016;35(1):88–95. pmid:26733705

35. Aburto JM, Beltrán-Sánchez H. Upsurge of homicides and its impact on life expectancy and life span inequality in Mexico, 2005–2015. American journal of public health. 2019;109(3):483–489. pmid:30676788

36. Gómez-Dantés H, Fullman N, Lamadrid-Figueroa H, Cahuana-Hurtado L, Darney B, Avila-Burgos L, et al. Dissonant health transition in the states of Mexico, 1990–2013: a systematic analysis for the Global Burden of Disease Study 2013. The Lancet. 2016;388(10058):2386–2402. pmid:27720260

37. Benita F. Social backwardness in Mexico City metropolitan area. Soc Indic Res. 2016;126(1):141–160.

38. Instituto Nacional de Estadística y Geografía. Instituto Nacional de Estadística y Geografía; 2019. https://www.inegi.org.mx/.

39. Secretaría de Salud de la Ciudad de México. Secretaría de Salud de la Ciudad de México; 2019. http://data.salud.cdmx.gob.mx/portal/.

40. Instituto Nacional de Estadística y Geografía. Instituto Nacional de Estadística y Geografía: Marco Geoestadístico, diciembre 2018; 2019. http://en.www.inegi.org.mx/app/biblioteca/ficha.html?upc=889463674658.

41. Bivand R, Rundel C. rgeos: Interface to Geometry Engine—Open Source (’GEOS’); 2018. Available from: https://CRAN.R-project.org/package=rgeos.

42. Pebesma E. Simple Features for R: Standardized Support for Spatial Vector Data. The R Journal. 2018;10(1):439–446.

43. Hijmans RJ. raster: Geographic Data Analysis and Modeling; 2019. Available from: https://CRAN.R-project.org/package=raster.

44. Hijmans RJ. geosphere: Spherical Trigonometry; 2017. Available from: https://CRAN.R-project.org/package=geosphere.

45. Pebesma E. spacetime: Spatio-Temporal Data in R. Journal of Statistical Software. 2012;51(7):1–30.

46. Pebesma EJ, Bivand RS. Classes and methods for spatial data in R. R News. 2005;5(2):9–13.

47. Bivand R, Keitt T, Rowlingson B. rgdal: Bindings for the’Geospatial’ Data Abstraction Library; 2019. Available from: https://CRAN.R-project.org/package=rgdal.

48. Wickham H. ggplot2: Elegant Graphics for Data Analysis. Springer-Verlag New York; 2016. Available from: http://ggplot2.org.

49. Wilke CO. cowplot: Streamlined Plot Theme and Plot Annotations for’ggplot2’; 2019. Available from: https://CRAN.R-project.org/package=cowplot.

50. Auguie B. gridExtra: Miscellaneous Functions for “Grid” Graphics; 2017. Available from: https://CRAN.R-project.org/package=gridExtra.

51. Dunnington D. ggspatial: Spatial Data Framework for ggplot2; 2018. Available from: https://CRAN.R-project.org/package=ggspatial.

52. Secretaría de Salud de la Ciudad de México. Secretaría de Salud de la Ciudad de México: Mortalidad; 2019. http://data.salud.cdmx.gob.mx/portal/index.php/informacion-en-salud/103-informacion-salud/354-mortalidad.

53. Gouveia N, Junger WL, Romieu I, Cifuentes LA, de Leon AP, Vera J, et al. Effects of air pollution on infant and children respiratory mortality in four large Latin-American cities. Environ Pollut. 2018;232:385–391. pmid:28966023

54. Salto-Quintana JN, Rivera-Alfaro G, Sánchez-Ramos EL, Gómez-Gómez A, Noyola DE. Post-pandemic influenza-associated mortality in Mexico. Pathog Glob Health. 2019; p. 1–8.

55. Wang N, Mengersen K, Tong S, Kimlin MG, Zhou M, Wang L, et al. Lung Cancer Mortality in China: Spatial and Temporal Trends Among Subpopulations. Available at SSRN 3315845. 2019.

56. Wändell P, Carlsson AC, Gasevic D, Sundquist J, Sundquist K. Neighbourhood socio-economic status and all-cause mortality in adults with atrial fibrillation: a cohort study of patients treated in primary care in Sweden. Int J Cardiol. 2016;202:776–781. pmid:26474468

57. Yang W, Olson DR, Shaman J. Forecasting influenza outbreaks in boroughs and neighborhoods of New York City. PLOS Comput Biol. 2016;12(11):e1005201. pmid:27855155

58. Ryvicker M, Sridharan S. Neighborhood environment and disparities in health care access among urban Medicare beneficiaries with diabetes: A retrospective cohort study. NQUIRY J Health Car. 2018;55:0046958018771414. pmid:29717616

59. Goovaerts P. Geostatistical analysis of disease data: estimation of cancer mortality risk from empirical frequencies using Poisson kriging. Int J Health Geogr. 2005;4(1):31. pmid:16354294

60. Lin WC, Lin YP, Wang YC, Chang TK, Chiang LC. Assessing and mapping spatial associations among oral cancer mortality rates, concentrations of heavy metals in soil, and land use types based on multiple scale data. Int J Environ Res Public Health. 2014;11(2):2148–2168. pmid:24566045

61. Tadesse S, Enqueselassie F, Gebreyesus SH. Estimating the spatial risk of tuberculosis distribution in Gurage zone, southern Ethiopia: a geostatistical kriging approach. BMC Public Health. 2018;18(1):783. pmid:29940913

62. Chien LC, Yu HL, Schootman M. Efficient mapping and geographic disparities in breast cancer mortality at the county-level by race and age in the US. Spatial Spatio-temporal Epidemiol. 2013;5:27–37.

63. Oliver M, Webster R. A tutorial guide to geostatistics: Computing and modelling variograms and kriging. Catena. 2014;113:56–69.

64. Pebesma EJ. Multivariable geostatistics in S: the gstat package. Computers & Geosciences. 2004;30:683–691.

65. Gräler B, Pebesma E, Heuvelink G. Spatio-Temporal Interpolation using gstat. The R Journal. 2016;8:204–218.

66. Raudenbush S, Bryk AS. A hierarchical model for studying school effects. Sociology of education. 1986; p. 1–17.

67. Kwok OM, Underhill AT, Berry JW, Luo W, Elliott TR, Yoon M. Analyzing longitudinal data with multilevel models: an example with individuals living with lower extremity intra-articular fractures. Rehabilitation psychology. 2008;53(3):370. pmid:19649151

68. Anaya G, Al-Delaimy WK. Effect of the US-Mexico border region in cardiovascular mortality: ecological time trend analysis of Mexican border and non-border municipalities from 1998 to 2012. BMC public health. 2017;17(1):400. pmid:28477629

69. Green EP, Cho H, Gallis J, Puffer ES. The impact of school support on depression among adolescent orphans: a cluster-randomized trial in Kenya. Journal of Child Psychology and Psychiatry. 2019;60(1):54–62. pmid:30055002

70. R Core Team. R: A Language and Environment for Statistical Computing; 2018. Available from: https://www.R-project.org/.

71. Pinheiro J, Bates D. Mixed-effects models in S and S-PLUS. Springer Science & Business Media; 2006.

72. Di Rienzo J, Casanoves F, Balzarini M, Gonzalez L, Tablada M, Robledo C. InfoStat; 2011.

73. Pinheiro J, Bates D, DebRoy S, Sarkar D, R Core Team. nlme: Linear and Nonlinear Mixed Effects Models; 2018. Available from: https://CRAN.R-project.org/package=nlme.

74. Instituto Nacional de Estadística y Geografía. Instituto Nacional de Estadística y Geografía: Población 2015; 2020. https://www.inegi.org.mx/temas/estructura/.

75. Gobierno de la Ciudad de México. Programa Delegacional de Desarrollo en Cuauhtémoc 2016-2018; 2019. http://www.data.seduvi.cdmx.gob.mx/portal/index.php/programas-de-desarrollo/programas-delegacionales.

76. Espinal-Enríquez J, Larralde H. Analysis of Mexico’s narco-war network (2007–2011). PLoS One. 2015;10(5):e0126503. pmid:25993657

77. Gobierno de la Ciudad de México. Secretariado Ejecutivo del Sistema Nacional de Seguridad Pública; 2019. https://www.gob.mx/sesnsp.

78. Nesterko Y, Turrión CM, Friedrich M, Glaesmer H. Trajectories of health-related quality of life in immigrants and non-immigrants in Germany: a population-based longitudinal study. International journal of public health. 2019;64(1):49–58. pmid:29947823

79. Dolev T, Zilcha-Mano S. The role of the therapeutic relationship in the association between interpersonal behaviors and outcome: Comparison of two competing models. Psychotherapy Research. 2019;29(5):553–564. pmid:29310548