Environmental Injustice in Mexico City: A Spatial Quantile Approach

Abstract

The majority of studies on environmental justice show that groups with lower socio-economic status are more likely to face higher levels of air pollution. Most of these studies have assumed simple, linear associations between pollution and deprived groups. However, empirical evidence suggests that health impacts are greater at high-pollution concentrations. In this paper, we investigate the associations of extreme levels of particulate matter up to 10 micrometres in size (PM₁₀) and ozone with deprived conditions, children and elderly people at sub-municipal level in Mexico City, using Áreas Geoestadisticas Básicas (AGEBs) as the unit of analysis. We used spatial quantile regression to analyse the association for each quantile of the range of pollution values, while also addressing spatial autocorrelation issues. Across AGEBs, higher levels of PM₁₀ are significantly positively associated with deprived economic conditions and elderly people. These results demonstrate clear variations in the associations between PM₁₀ and vulnerable groups across the ranges of these pollutants. Ozone levels are positively associated with higher numbers of children. The findings reflect differences in the source and degradation of these pollutants and provide important evidence for decision-makers addressing air pollution inequalities and injustice in Mexico City and other cities.

Introduction

Environmental injustice refers to the unequal impact of environmental degradation on social groups depending on their social, economic, racial and ethnic background (Zimmerman 1993; Pulido 1996; Mohai et al. 2009; Raddatz and Mennis 2013; Laurian and Funderburg 2014).¹ Evidence has been accumulating on the unequal distribution of environmental risk across social groups, with people of low socio-economic status (SES) living in close proximity to hazardous facilities (Zimmerman 1993; Krieg 1995; Pastor et al. 2001; Saha and Mohai 2005; Mohai and Saha 2007; Schoolman and Ma 2012; Raddatz and Mennis 2013), and incinerators (Laurian and Funderburg 2014), having fewer green areas nearby (Johnson-Gaither 2011; Wolch et al. 2014), living in areas with a high risk of flooding (Grineski et al. 2015a), and being exposed to air pollution (Asch and Seneca 1978; Grineski et al. 2007; Downey and Hawkins 2008; Havard et al. 2009). However, other studies have not found environmental injustice (Hajat et al. 2013; Richardson et al. 2013; Padilla et al. 2014). The mixed nature of this evidence may be explained by differences in the type of hazard, geographical unit, methodology and local context (Briggs et al. 2008; Havard et al. 2009). This study focus on air pollution, where again, there are studies showing that people with lower socio-economic status are more exposed to air pollution (Grineski et al. 2007; Bell and Ebisu 2012; Carrier et al. 2014; Clark et al. 2014; Zou et al. 2014), while other studies debate whether such links exist (Branis and Linhartova 2012; Hajat et al. 2013; Richardson et al. 2013; Padilla et al. 2014). This literature emphasises the importance of targeting heterogeneous environments in policy-making. This paper aims to investigate spatial heterogeneity in the relationship between air pollution and social vulnerability in the urban setting of Mexico City, to provide further information that could facilitate the improved targeting of local policies to mitigate or adapt to pollution threats. Spatially targeted environmental programmes potentially have the advantage of focusing on small areas, where policy measures can have a higher impact, particularly for the most vulnerable groups, than is the case when public resources are dissipated across the city. Moreover, local programmes may increase confidence and capacity to incentivise community participation in policy initiatives (Smith 1999; Tunstall and Lupton 2003).

The analysis will focus on PM₁₀ and ozone, because of the serious impacts on human health at high levels of these pollutants (Maantay et al. 2009), which highlight the need to better understand where the social heterogeneity in exposure to air pollution may occur. PM₁₀ can cause heart disease, lung cancer, asthma, and acute lower respiratory infections, with more than 2 million people dying annually because of breathing tiny particles, present in indoor and outdoor air pollution (World Health Organization 2011b). Ozone is positively associated with daily mortality levels (World Health Organization 2006), and can cause reduction of lung capacity and serious lung damage (Levy et al. 2001). Moreover, Arceo et al. (2016) estimated that 1 μg/m³ increase of in 24-h PM₁₀ in our study area, Mexico City, results in an additional 0.24 deaths per 100,000 births. Thus, disproportional exposure of air pollution is strongly linked with health inequalities. Jerrett et al. (2001) defined this as ‘triple jeopardy’; people with deprived economic and social conditions are more likely to be exposed to high levels of contaminants and hence experience more negative impacts on their health. As a result, air pollution is a concern at the public health level due to its detrimental impact on human health and on the economy (Lagercrantz and Sundell 2000; Schwartz and Repetto 2000; Jerrett et al. 2004; Maantay et al. 2009; Nishimura et al. 2013; Parent et al. 2013; Beatty and Shimshack 2014). Hanna and Oliva (2015) estimated some of the economic costs of air pollution in Mexico City, showing that a 20% increase in air pollution can result in a reduction of 1.3 working hours in the following week. Filippini and Martínez-Cruz (2016) found that the individual’s willingness to pay for improved air quality in this city amounts to an average of US $262 (2008 US dollars) annually.

In this study, we apply a quantile regression approach to investigate the hypothesis that the association between vulnerable groups and air pollution grows stronger as pollution concentrations increase. That is, in highly polluted locations, we expect the evidence for environmental injustice to be stronger than in locations with lower levels of pollution. We, therefore, extend the work of Chakraborti et al. (2017) and Rissman et al. (2013), with the latter authors showing that higher airport-contributed PM_2.5 (fine particulate matter that is 2,5 microns in diameter and less) concentrations have different relationships with social minority indicators compared with the rest of the PM_2.5 distribution. Our approach contrasts with most empirical studies, which explore the association between socio-economic conditions and air pollution assuming a homogeneous air pollution pattern across the studied area, using mean levels of air pollution with standard regression (i.e., ordinary least squares regression, OLS), leaving aside its lower and upper values with respect to the mean. The use of quantile regression, examining different percentiles of the conditional air pollutant distribution, can better account for spatial heterogeneity in air pollution levels and identify changes in its relationship with deprived economic conditions, that may be missed by the application of conventional mean regressions. It can hence be more informative to policy-makers. Programmes that mitigate air pollution impose social consequences associated with compliance with new regulations, as well as health benefits; the greater insights from quantile regression can help to identify any distributional issues that may need addressing in spatially targeted policies. In this sense, it is of interest to analyse the sensitivity of the environmental justice hypothesis (i.e. differences in vulnerable groups’ exposures) at locations with extreme values of pollutants, as this may be the result of particular actions in these locations that do not necessarily occur elsewhere (e.g. the location of industrial facilities near age-vulnerable communities). In that case, policy-makers can consider spatially targeted emission-reduction policies (e.g. truck-rerouting, low-emission zones, industry re-allocation) in some locations to produce the strongest benefits to environmental justice (Nguyen and Marshall 2018).

This research takes the sub-municipality Áreas Geoestadisticas Básicas (AGEBs), the smallest administrative units in Mexico, as the spatial unit of analysis. This has the advantage that socio-economic characteristics are likely to be fairly homogeneous within these small geographical areas, which will enhance the reliability of results obtained (Bowen et al. 1995; Maantay 2002). Potential spatial autocorrelation will be accounted for to ensure robust hypothesis testing and the estimation of coefficients. Otherwise, the assumptions regarding the independence and identical distribution of the residuals would not meet. This may bias the estimators due to the inflated values of t statistics, and hence to reject the null hypothesis incorrectly, Type error I results (Anselin 2002; Dormann et al. 2007). That is, the estimators would appear significant when they are not.

Data and Methods

Area of Study

The choice of the case study of Mexico City is consistent with growing concerns about air pollution in urban areas in developing countries, where high population densities and low quality health services collide with high levels of harmful pollutant concentrations, impacting on residents’ health and well-being (Krzyzanowski et al. 2014; Marlier et al. 2016). In Mexico City, the annual average PM₁₀ concentration for 2014 and 2015 was 43.5 ug/m³, and the concentration of ozone is increasing, with an annual average of 27 and 29.5 parts per billion (ppb) in 2014 and 2015, respectively (Air Quality in Mexico City 2014). These values are higher than the World Health Organization (World Health Organization 2006) thresholds,² above which there are significantly increased risks to health. Moreover, there is also limited literature on environmental justice in developing countries (Pearce and Kingham 2008; Rooney et al. 2012). For the specific case of Mexico, there is only scant evidence around environmental injustice, which is focused on industrial contaminants in the north and border regions with US, and with emissions generally obtained by measures of proximity to industrial facilities (Blackman et al. 2003; Grineski and Collins 2008, 2010; Lara-Valencia et al. 2009; Grineski et al. 2015b; Chakraborti et al. 2017). Chakraborti et al. (2017) provide an exception, as they undertook a nation-wide analysis focusing on water disposal of toxic metals. They found a positive association between marginalisation (poorer communities) and pollution, with stronger evidence at locations with higher levels of water toxic pollutants.

Pollution Data

The pollution data, ozone and PM₁₀ for the year 2015, were obtained from the measuring stations operated by the Automatic Air Quality Monitoring Network of Mexico City (RAMA n.d.), which provide hourly records. We estimated the 24hr mean for PM₁₀ and ozone (from 10am to 6pm), each averaged into annual mean concentrations, following previous studies (e.g. Romieu et al. 2012). This analysis included data for all the stations that had at least 75% of information in the studied year. The numbers of measuring stations that met this criterion, and were therefore used to compute the 24hr values, were 23 and 31 available stations for PM₁₀ and ozone, respectively. The geographical coverage of the monitoring stations network contains some areas with sparse data (see Figs. 3, 4 on Appendix). The discussion further elaborates on this issue.

We applied a universal kriging algorithm to obtain interpolated values for each pollutant, at the AGEB level, from the monitoring stations data. This technique is considered one of the best interpolation methods because it deals better with erroneous local variability compared with other interpolation techniques such as inverse distance weight (IDW) (Jerrett et al. 2005a). We, therefore, complement previous work (Hanna and Oliva 2015; Arceo et al. 2016) which has used the IDW technique to carry out the interpolation, leaving aside the spatial variability which is common in pollutant datasets. We chose universal instead of ordinary kriging because this approach considers the global trend over the area of study and takes into account the spatial dependence (Burrough and McDonnell 1998). Moreover, universal kriging models have been used previously in the area of environmental justice along with epidemiological studies (Jerrett et al. 2001, 2005b; Finkelstein et al. 2003; Künzli et al. 2005; Su et al. 2011) to interpolate air pollution data.

Economic and Geographic Data

Economic information was obtained from the Population and Housing Census (Instituto Nacional de Estadística y Geografía - INEGI 2010) at the level of the AGEB, which includes the number of households with TV, car, computer, landline phone, mobile phone and internet; this information, households’ purchasing power, is used in this study to characterise the households’ SES. Demographic information from the same data source was obtained on the number of children and elderly people. Children were considered from 0 to 11 years old (US Food and Drug Administration 2014) and elderly people were from 65 years old onwards (World Health Organization 2011a). Those AGEBs for which this information was either not available or labelled in the dataset as confidential (n = 126) were excluded from the analysis, resulting in a total of 2287 AGEBs in the analysis. Interpolated pollution and economic–geographic datasets were merged, assigning a pollutant value to each of the economic–demographic variable at the AGEB level.

Statistical Analysis

A principal component analysis was used to generate a deprivation index as a proxy for households’ economic deprivation conditions. This approach follows the previous literature, where the method has been used to create socio-economic indices (Richardson et al. 2013; Rissman et al. 2013; Grineski et al. 2015a), particularly in a developing country context, due to a frequent lack of official data on income e.g. (Fiadzo et al. 2001; Fotso and Kuate-Defo 2005). A principal component analysis also controls for the high collinearity among the economic variables. This analysis identifies the components which explain a significant cumulative proportion of the variance of the data set; We extracted the components with Varimax rotation to simplify the expression and hence its interpretation.

A Global Moran Index (Anselin et al. 1996) was calculated to explore the spatial distribution of PM₁₀, ozone (the original values from the monitoring stations and the interpolated values), deprivation index and vulnerable-aged groups. This index takes values from − 1 to 1, where a large negative or positive value means that there is spatial autocorrelation, there are some clusters, where the values of the neighbouring AGEBs are dissimilar or similar, respectively. In contrast, when the value approaches zero, it means that there is random spatial pattern.

As an initial exploratory analysis to assess environmental injustice, We assigned the annual averages of PM₁₀ and ozone into five economic deprivation categories (quintiles) across Mexico City’s AGEBs and used one-way ANOVA and Tukey–Kramer test to evaluate whether there were any statistical differences in the mean pollution levels for the extreme quintiles (i.e., between AGEBs with households with the lower and higher levels of deprivation conditions). We then used regression analysis, as done in previous studies (Rissman et al. 2013; Carrier et al. 2014; Fecht et al. 2015) to assess the association of a given minority group with each pollutant, and determine its statistically significance, after controlling for the other groups. We first carried out standard ordinary least squares (OLS) regressions to quantify general associations of economic deprivation and vulnerable-aged groups with PM₁₀ and ozone concentrations. We also analysed the potential heteroscedasticity, in terms of different variance of the residuals across the distribution, of the regression with the Breusch and Pagan (1979) test. To better understand exposure to high pollutant levels, and to deal with the potential heteroscedasticity, we applied a quantile regression. This simple technique allowed us to assess the levels of association of the economic deprivation index and the proportion of children and elderly people with concentrations of PM₁₀ and ozone across the full range of concentration levels for each pollutant. This allowed us to examine how the relationship of pollution levels for vulnerable groups changes at different levels of the pollutants (for example at the highest and lowest pollution levels). The presence of residual spatial autocorrelation was examined using the Moran Index for the analysis of both PM₁₀ and ozone, which led to evidence of the potential biased estimators; therefore, a spatial regression was applied in the quantile analysis to obtain accurate coefficients.

Briefly, we describe below the quantile regression which estimates the conditional quantile functions in contrast with the conditional mean functions of ordinary least square (OLS). Quantile regression uses the full sample and allows to determine the effect of the determinants across the full distribution (quantiles) of the dependent variable. Unlike OLS, the quantile approach can deal with heteroscedasticity, outliers and unobserved heterogeneity (Koenker and Hallock 2001; Koenker 2005); in this analysis, this is convenient because it does not assume any distributional assumption (independent and identically distributed) of the residuals, allowing uneven distribution on the PM₁₀ and ozone. According to Koenker (2005) instead of minimizing the sum of the squared residuals as in OLS, quantile regression focuses on minimizing a weighted sum of the absolute deviations:

$$\mathop {{\text{min}}}\limits_{{\{ b_{j} \} _{{j = 0}}^{k} }} \mathop \sum \limits_{{i = 1}}^{n} |y_{i} - \mathop \sum \limits_{{j = 0}}^{k} \beta _{j} X_{{j,i}} |h_{{i}}$$

where $y$ = dependent variable, $X$ is the vector of the covariates and $\beta$ is the vector of the slopes. The weight is defined either as ~ $h_{i } = 2q$ when the residual for the ith observation is positive or as ~ $h_{i } = 2 - 2q$ if the residual is negative; and $q \in (0,1)$ denotes the quantile of the dependent variable to be estimated.

Spatial autoregressive models (SAR), lag model and spatial error, are commonly used to tackle the potential spatial autocorrelation in linear regressions (Anselin and Arribas-Bel 2013). In this paper, We applied the lag model approach due to the fact that the quantile analysis is not applicable to the spatial error model (Liao and Wang 2012) and because the dependent variable is highly clustered, with nearby AGEBs tending to have similar levels of pollution. Following McMillen (2012), the quantile regression with spatial lag model is defined as:

$$y = \lambda \left( q \right)Wy + X\beta \left( q \right) + u$$

where $W$ is the spatial weight matrix, which denotes the spatial relation between each value of $y$ and its neighbours; $\lambda \left( q \right)$ is the spatial lag parameter; $u$ is the error term. The $W$ (weighted matrix) was constructed to model the structure of the spatial lag component, as shown above, using the first contiguity method. This method was chosen because the size of AGEBs is highly heterogeneous (with mean 0.317 and standard deviation 0.324 square kilometres); it involves creating regions with links if AGEB i and AGEB j share one or more boundary points. Three isolated AGEBs were excluded from analysis because they did not have any neighbour.

There are different methods to handle the spatial lag component in a quantile regression model. Kim and Muller (2004) introduced the Two-Stage Quantile Regression (2SQR) which requires the estimation of two consecutive quantile regressions. In this paper, however, We used an Instrumental Variable Quantile Regression (IVQR) (Chernozhukov and Hansen 2006), where the same quantile is used just in one stage leading to more robust results (McMillen 2012). First, an instrumental variable is created for $Wy$ from the predictive values of an OLS regression of $Wy$ on a set of instruments $Z$ ($X {\text{and}} WX$). Then, a quantile regression is fitted (one regression for each $\lambda$ value) $y - \lambda Wy$ on $X$ and $\widehat{Wy}$, using the created instrumental variable for $Wy (\widehat{Wy}).$ The estimated value of $\lambda$leads to small coefficients (closest to zero) on $\widehat{Wy}$. Having the values of $\widehat{\lambda }$ a quantile regression $y - \widehat{Wy}$ on X is fitted to get the estimated values of $\beta$. It is expected that $\widehat{Wy}$ will be zero when the instruments are chosen properly (McMillen 2012). We analysed the estimated spatial lag variable to illustrate the level of spatial autocorrelation. Finally, an analysis of variance was applied to test whether the spatial coefficients across the different quantiles on the pollutants are statistically different from one another (Koenker 2006). The final dataset contained 2287 observations. We used R program version 3.2.3 for the analysis.

Results

Descriptive Analysis

Table 1 provides an overview of the descriptive statistics of the pollution, economic and demographic variables, of the households, in percentages (%) in Mexico City. More than half of the households lacked internet access (60%) or had no car access (53%). In addition, approximately a quarter of them did not have landline or mobile phone access (28% and 25%, respectively).

Table 1 Description and descriptive statistics of the pollution, and households’ economic and demographic variables in percentages (%) in Mexico City

Table 2 shows the results of principal component analysis (PCA) capturing the households’ economic deprivation conditions. It shows that the first component explains 89% of the cumulative in all collinear economic variables and its eigenvalue is greater than 1. This component comprises one cluster which describes the level of deprivation of car, computer, landline, mobile phone and internet for the households. All variables had high loading values which reflect the important contribution on this first component. We considered this as an economic deprivation index. Households in AGEBs with low values of this deprivation index had better purchasing power (in terms of the economic variables), while high values represent worse deprived conditions.

Table 2 Principal component analysis (PCA) results with component loadings for the economic variables for Mexico City AGEBs

Spatial Analysis

An analysis on the potential spatial autocorrelation within the original pollutant dataset, coming from the monitoring stations, showed a positive and significant Global Moran Index [Moran Index = 0.34 and 0.48 with a p-value < 0.001] for PM₁₀ and ozone, respectively. Similarly, the Moran Index results for the interpolated pollution data also showed that there is a significant positive spatial autocorrelation in the concentration levels of both pollutants (Moran Index = 0.99 with a p value < 0.001, for both pollutants). Thus, Fig. 1a and b show the spatial distribution of the PM₁₀ and ozone interpolated values, with a high level of clustering. The darker red shading shows the highest levels of concentration of PM₁₀ and ozone. Figure 1a illustrates that the higher levels of PM₁₀ were mainly found in the north of the city, whereas the south-faced lower levels. For ozone, it is the opposite: the south area presented higher levels of concentration, while in the north the levels were lower (Fig. 1b). These pollutants were thus found to be highly negatively correlated (r = − 0.77, Pearson correlation), with a p-value < 0.001. These different patterns reflect differences in the sources and chemical processes associated with particulate and ozone pollution. Particulate pollution is concentrated in cities due to its source in power generation, domestic heating and traffic. Ozone is not emitted directly into the city environment to any great extent but is formed through a series of photochemical reactions involving reactive organic compounds and nitrogen oxides emitted from combustion engines. However, the nitrous oxide (NO) produced from combustion engines also reacts with ozone itself to form nitrogen dioxide (NO₂), thus removing ozone from the city environment. In suburban or peri-urban areas (those areas which are on the edge of the city), there is less traffic, hence less available nitrous oxide to react with ozone, and ozone is therefore more persistent (Briggs et al. 2008; Fecht et al. 2015). Suburban or peri-urban areas are common to the south of Mexico City, and hence there is greater ozone in this area. The spatial distribution of the economic deprived index in Fig. 1c shows a relatively high level of clustering (Moran Index = 0.7 with a p-value < 0.001). In such figure, the red colour shows the areas with the more deprived conditions. The distribution of the deprivation index shows an economic gradient from the less-deprived AGEBs (green colour) to the most deprived AGEBs (red colour). The AGEBs with the less-deprived households form a big cluster located at centre-west, where the purchasing power is higher than in other locations in the city. Households with most deprived conditions form small clusters in the north, centre-east and south. Higher percentages of vulnerable groups (children and elderly people) are located in red-coloured AGEBs. The Moran Indices for children and elderly people were 0.45 and 0.37, respectively, indicating the presence of statistically significant clustering of AGEBs with similar proportions of vulnerable age groups (with a p-value < 0.001). Figure 1d and e, with the dark red and whiter colours show the areas with the higher and lowest proportion of children and elderly, respectively.

Fig. 1

a and b Spatial distribution of PM₁₀ and ozone across Mexico City in 2015. c–e Spatial units according to AGEBs. Spatial distribution of deprivation index, proportion of children and elderly people across Mexico City in 2010. Spatial units according to AGEBs

Table 3 reports the descriptive statistics of PM₁₀ and ozone concentrations according to the quintiles of the deprivation index. The category q1 denotes the AGEBs with less economic deprivation conditions in terms of purchasing power and q5 captures those with the most deprived conditions. The results in this table show that the most economic deprived AGEBs (q5) experienced higher PM₁₀ concentration levels (43.3 ug/m³) compared with the less-deprived AGEBs (q1) (40.3 ug/m³). In contrast, for ozone the most deprived AGEBs (q5) had lower concentrations (56.7 ppb) and the less-deprived AGEBs (q1) had higher concentrations of this pollutant (57.4 ppb).

Table 3 Distribution of PM₁₀ and ozone concentrations levels by deprived index in quintiles, in ug/m³ and ppb units, respectively

We carried out one-way ANOVA and Tukey–Kramer tests to compare the differences in the mean values of both pollutants across the different deprived categories (i.e., across all the quintiles and between each quintile, respectively). One-way ANOVA shows that the means for all the quintiles were statistically different (with a p value < 0.005) for both pollutants. The Tukey–Kramer test indicates that the differences for q1-q2, q1-q3, q1-q4 and q1-q5 were significant (with a p value < 0.005) for PM_10. Similarly, the differences for q2-q1, q2-q3, q2-q4 and q2-q5 were also significant (with a p value < 0.005) for ozone.

Regression Analysis

The regression analysis further assesses the relationship of pollution levels with the deprivation index, and the proportion of children and elderly. Table 4 illustrates the results of both OLS and the spatial quantile regressions, showing the relationship of PM₁₀ (table 4a) and ozone (Table 4b), the economic deprivation index, percentage of children and percentage of elderly people, at AGEB level. The heteroscedasticity identified in the OLS model (Breusch-Pagan with a p value < 0.001] and the positive spatial autocorrelation detected in the residuals of the quantile regression [Moran Index ≥ 0.8 with a p value < 0.001) justify the use of this spatial quantile regression. Both models in Table 4 show that the two pollutants were significantly related with the deprivation index and vulnerable-aged groups (p value < 0.001), but in different ways. In general, the analysis of the coefficients, for both models, shows that PM₁₀ (ozone) was positively (negatively) associated with economic deprivation conditions, after controlling for age groups. This provides evidence of environmental injustice for population with deprivation economic conditions residing in locations with higher PM₁₀ levels. In contrast, no such inequity was found for ozone; in fact, populations with lower deprivation conditions were associated with higher ozone exposure. Likewise, elderly people (children) were associated positively (negatively) with PM₁₀ concentrations, after controlling for SES conditions. In contrast, ozone levels were positively (negatively) associated with children (elderly).

Table 4 OLS and spatial quantile regression estimates assessing the relationship of (a) PM₁₀ and (b) ozone and the economic deprivation index, children and elderly within each AGEB

We found similar results between the mean of the OLS and the 50% quantile (of the spatial quantile regression) for each explanatory variable (except for elderly people where the difference was around 6 units for the case of ozone). However, the spatial quantile regression outcomes showed a clear variation in the relationship of PM₁₀ with the economic deprivation index and elderly people, except the 90% for elderly (see Table 4a and Fig. 2). Figures 2 shows the coefficients of elderly, children and deprivation index with the different quantiles of PM10 and ozone (each dot represent from the 10% to the 90% quantiles). Regarding the economic deprivation, the quantile estimators showed that within the 80% and 90% quantiles, PM₁₀ levels of pollution were more strongly positively associated with the economic deprivation index than the lower quantiles (the 10% and 20% quantiles). The variation was almost two times larger in the right tail (3.2) than the left one (around 1.7). Figure 2, left side, also shows a general increasing trend of estimators across the lower and upper levels of PM_10. With respect to elderly people, higher levels of PM₁₀ were more strongly positively associated with higher percentage of elderly people within the higher quantiles of the distribution with β = 43.9 (averaged for the 70% and 80% quantiles) than within the lower quantiles with β = 28.8 (averaged for the 10% and 20% quantiles). Therefore, these results showed that the average estimates, β = 40.2 and β = 38.4, provided by the OLS and the 50% of the spatial quantile regression, respectively, were lower with respect to the values of the higher quantile estimates but above the values of the lower quantile estimates of PM₁₀. Figure 2, left side, also illustrates a raised pattern from the 10% to 80% quantiles for elderly people and PM₁₀. Note that these results do not report an interesting heterogeneity for the different levels of PM₁₀ and its association with the percentage of children in each AGEB.

Fig. 2

Spatial quantile regressions with different quantiles of PM₁₀ (left) and ozone (right) as response variables

The variation between ozone higher and lower levels and deprived economic conditions and vulnerable age groups was not clear (see table 4b and Fig. 2, right side). For example, there was no significant variation in the negative association of ozone with the deprivation index or the percentage of elderly people, at higher and lower levels of concentration. Nevertheless, We can identify some patterns. At lower levels of ozone, β = − 38.8 and β = − 51 (20% and 30% quantiles) the elderly people had a stronger negative association with this pollutant than when ozone was at its upper levels, β = − 26 and β = − 25.7 (80% and 90%). Likewise, within the upper quantiles, of ozone pollution, the level of pollution, β = 26 (averaged at the 80% and the 90% quantiles), was more strongly positively associated with the percentage of children in the AGEB than in the lowest quantile, β = 9.9 (10% quantile).

The values of the spatial lag variable, ${\text{Wy}}$, were consistently small (all 0.01, except for the 90% quantile which was < 0.01) and significant, p value < 0.005 across all quantiles, suggesting that spatial autocorrelation is minimal for both pollutants over the full range of quantiles.

Importantly, the results of the spatial quantile model highlight the non-linearity of some associations across the all pollution levels, especially for PM₁₀ in relation to deprivation and percentage of elderly people; and for ozone in relation to the percentage of children. One-way ANOVA, which measures the precision of the different estimators across the quantiles of each pollutant, confirmed a high significant difference between such estimators (p value < 0.005).

Discussion

The analysis investigated spatial heterogeneity, comparing exposure to higher and lower levels of each pollutant, PM₁₀ and ozone, across vulnerable groups in Mexico City. Overall, our results show a positive association between deprived economic conditions and PM₁₀ and a negative association between lower socio-economic conditions and ozone. Even though the analyses focused on different levels of air pollution, which have been rarely studied, our findings are consistent with previous studies that focused on the mean of air pollution. The positive association of PM₁₀ with the deprivation index is congruent with the previous literature (Briggs et al. 2008; Fecht et al. 2015). Moreover, other studies analysing PM_2.5 (Gray et al. 2013; Hajat et al. 2013) found a positive association between better socio-economic conditions and lower exposure of this pollutant. With respect to ozone and its negative association with deprivation index, our findings are also similar to previous research (Briggs et al. 2008; Gray et al. 2013). Conversely, Grineski et al. (2007) found a positive relation of ozone with more deprived economic status. With respect to vulnerable age groups, the findings show mixed evidence as previous studies: higher concentrations of PM₁₀ were significantly associated with higher proportions of elderly people but with lower proportions of children. With respect to children our results are similar to Fecht et al. (2015) and Carrier et al. (2014), with the latter study using PM_2.5. Likewise, elderly people were to be found disproportionally exposed to other pollutants, SO₂ and NO₂ (Clark et al. 2014; Zou et al. 2014) which is consistent with our outcomes. When considering ozone, the outcomes illustrate that elderly people are not disproportionally exposed to this pollutant. Instead children were found to face disproportional exposure to ozone. Calderón-Garcidueñas and Torres-Jardón (2012) showed that children, living in the South of Mexico City, were highly exposed to ozone, which is congruent with our results.

The results highlight that the higher the PM₁₀ level is, the greater the level of disproportionate exposure of this pollutant to people in deprived economic conditions and with elderly people. Thus, the findings show that the association of the AGEBs with economic deprivation conditions was significantly heterogeneous on the different level quantiles of PM₁₀, especially for the lower and upper concentration levels. In general, our results verify the hypothesis of an increasing pattern of this association from the lower to the higher quantiles of PM_10. This is, higher levels of PM₁₀ were more strongly and positively associated with those AGEBs with deprived conditions than those with lower levels of this pollutant. This result is consistent with the findings of Rissman et al. (2013), who found a slightly decreased association between median income and concentrations of PM_2.5 pollution due to aircraft, from the 50 to 90% quantiles. In the case of elderly group, we also identified an increasing trend of PM₁₀ exposure, from lower to higher quantiles (excepting the 90% quantile). This would imply that the health of these groups (those with deprived economic conditions and elderly people) is at risk due to high levels of PM₁₀ concentration. In Mexico City, elevated levels of this pollutant are more than double WHO’s threshold levels, which were established to avoid health risks. Therefore, these specific groups should be targeted in pollution reduction policies at those locations.

The spatial distribution analyses partially explains the higher exposure of PM₁₀, where traffic and industry processes are their principal sources (Querol et al. 2008), on deprived conditions and elderly people. Clusters of elderly people were found in the municipalities of Cuauhtémoc, Miguel Hidalgo and Venustiano Carranza, where the high proportion of this age group is due to lower fecundity rates and better medical services (Negrete 2003) than in other areas of the city. These areas are also (particularly Cuauhtémoc and Miguel Hidalgo), where most of the public services and jobs in Mexico City are located (Instituto de Políticas para el Transporte y el Desarrollo ITDP México 2015), attracting much commuter traffic. From 2008 to 2012 the vehicular fleet increased by close to 11%, this figure was elaborated based on the information of ‘report of the quality of the air’ (Instituto Nacional de Ecología y de Cambio Climatico (INECC) 2013). Moreover, Cuauhtémoc has two of the main and busiest avenues: ‘Paseo de la Reforma’ and ‘Insurgentes’. Therefore, policies, such a congestion tax, in these two municipalities, that incentivise the use of low-emission public transport and less frequent vehicle usage would benefit the health of people living there by lowering the level of PM₁₀ emitted by vehicles. Such spatially targeted policy has been applied in cities like Stockholm, Gothenburg and London (Leape 2006; Börjesson and Kristoffersson 2018). Central London, after the introduction of congestion charge in February 2003, experienced a reduction of about 20 percent on automobile traffic (Litman 2005). This allowed to lower the pollution emitted by vehicles (Leape 2006). Note that these policies may need an improvement of public transportation in advance; as it was the case in London, where there was an expanded bus lane system and major renovations to the subway system (Litman 2005).

The association of deprived conditions and PM₁₀ can be explained using the arguments of Calderón-Garcidueñas and Torres-Jardón (2012) that less economically privileged people spend considerable time in the traffic or close to it, walk long distances to take the crowded transport or work on the streets. In that sense, an improvement in low-emission public transport, as mentioned above, would benefit the poorer communities as well. Moreover, note that the spatial analysis identified clusters of AGEBs with lower SES particularly in the north area, which includes the municipality of Gustavo A. Madero, which is recognised as one of the areas with the greatest concentration of people in poverty (CONEVAL n.d.). This northern area (mainly the municipality of Azcapotzalco and Gustavo A. Madero) is also characterised by having many industries and main roads (Air Quality in Mexico City 2014). This industrial character in the north is related to the availability of nearby facilities, new housing construction, and better quality infrastructure (Cruz and Garza 2014). Therefore, spatially targeted policies could be implemented, in this northern area, to reduce PM₁₀ pollution from the industries there. This could be done by either obligating and/or incentivising better housekeeping, material substitution, recycling or process innovations (Cairncross 1992; Willig 1994).

Our findings should be interpreted with some caution due to some methodological and data limitations. First, data on air pollution could be improved using different modelling methods such as Atmospheric Dispersion Modelling System (Havard et al. 2009), Land-use Regression Models (Ryan and LeMasters 2007) or other models such as Integrated Meteorological-Emission Models or Hybrid Models. This is because these models use more variables and information such as traffic volumes, land-use, meteorology, topography to accurately model air pollution. However, it is because of this extra information and special equipment and software (Jerrett et al. 2005a) that we could not use them in the study. We applied the universal kriging interpolation, which is based not only on the distance between the measured points but also on the overall spatial arrangement of the measured points to overcome this issue. One advantage of applying the kriging approach is the production of standard errors which quantify the degree of uncertainty of the spatial prediction, allowed us to identify the places with less reliable interpolation values (Mulholland et al. 1998). In that sense, we expect PM₁₀ pollution estimates in the south, where there are sparser data due to fewer monitoring stations, to have larger standard errors; meaning, that these errors may influence the results. As monitoring stations are not equally distributed across space, this problem is often acknowledged in the literature. For example, Künzli et al. (2005) obtained larger standard errors, as the result of universal kriging, on the periphery of Los Angeles metropolitan area with 23 monitoring stations. We followed Künzli et al. (2005) study by carrying out a sensitivity analysis to check the robustness of our regression outcomes, coming from the interpolated pollution values, especially for PM₁₀. This involves down-weighting estimates with larger errors, in weighted least-square models (the weights specified as the inverse of the standard errors) and comparing the results with the main models with the universal kriging estimates. Thus, we accounted for the standard errors, obtained from the kriging interpolation, in the regressions. The outcomes were robust and similar to what we found in the original regression model, especially for ozone (results available in Fig. 4 in Appendix). Figure 4 shows the coefficients of elderly, children and deprivation index with the different quantiles of PM10 and ozone (each dot represents a percentage quantile from 10 to 90%). There was a variation of six units for elderly people after controlling for the standard errors for PM₁₀, but the results for the other variables were quite similar to the original regressions. A second limitation of the approach is that we did not consider the mobility of people. It is difficult to measure the activity patterns of people, which is often ignored in the literature (Havard et al. 2009; Fecht et al. 2015); for example, where they spend more time, at home or in their jobs, and how far away they live from their jobs. This would require extensive data on behavioural patterns that were not available for our study site. Finally, while we recognised that alternative theories and approaches address the relationship between income and pollution (Martinez-Alier 1995; List and Gallet 1999; Yaduma et al. 2015; Stern 2017), we followed the existing literature on environmental justice (Rissman et al. 2013; Carrier et al. 2014; Fecht et al. 2015) by not making a specific attempt to explain the co-location of vulnerable groups and pollution. This would have required us to control for the problem of reverse-causality (i.e., income may affect pollution through greater production levels, or the amount of pollution may affect income as people of low SES live in cheaper, but often also more polluted, areas), and for important omitting variables such as political or regulatory efforts, strong enforcement institutions, research and development activities or infrastructure (on this point see Lin and Liscow 2012; Germani et al. 2014). An alternative, longitudinal approach may have allowed us to gain additional insights into the chronological causal relationships that contribute to environmental inequalities (Briggs et al. 2008; Havard et al. 2011; Rissman et al. 2013), but such an approach is dependent on a suitable time series of data. Here, we used a cross-sectional approach to analyse the evidence for environmental injustice across all AGEBs in Mexico City, more than two thousand, facing heterogeneous levels of air pollution (particularly for those at the edge of the distribution, lowest and highest values).

Aside from these caveats, this study provides some distinct advantages over much previous work. We used spatial quantile regression, which shows the heterogeneous spatial distribution of the link between air pollution and vulnerable social conditions, with stronger unequal exposures for SES and elderly people in locations with upper levels of concentrations of PM₁₀. Moreover, we used AGEBs, the smallest geographical units in Mexico City, and accounted for the spatial effect of clustering of the data set, and hence avoided producing biased estimators. These methodological aspects all contributed to enhancing the robustness of the results.

Appendix

See Figs.

Fig. 3

a PM₁₀ monitoring stations across urban AGEBs in Mexico City in 2015. b Ozone monitoring stations across urban AGEBs in Mexico City in 2015

3 and

Fig. 4

Spatial quantile regressions with different quantiles of PM₁₀ (left) and ozone (right) as response variables, controlling for the standard errors from the kriging interpolation

4.