1. Introduction

In 2015, the World Health Organization (WHO) recognized ambient air pollutant exposure and its effects on human health as a public health priority requiring further study (^{WHA68, 2015}). The rapid expansion of urban areas worldwide has increased the number of individuals exposed to pollutant concentrations above the upper limits recommended by the WHO (^{Riojas-Rodríguez et al., 2016}). For instance, the Mexico City Metropolitan Area (MCMA) is the third-largest metropolitan area of the member states of the Organization for Economic Cooperation and Development (OECD) and is considered one of three megacities in Latin America (^{OECD, 2016}). Historically, this area has registered high rates of air pollution (^{INE-SEMARNAT, 2011}), posing a constant risk to the health of its nearly 21 million inhabitants (^{INEGI, 2015}).

Epidemiological studies have shown that chronic exposure to ambient PM_2.5 and NO₂ is associated with health outcomes including increased mortality (^{Beelen et al., 2014}), cardiovascular disease (^{Dockery, 2001}; ^{Hoek et al., 2013}), lung cancer (^{Pope et al., 2002}; ^{Chen et al., 2008}; ^{Hamra et al., 2014}), cardiopulmonary conditions (^{Krewski et al., 2009}) and diabetes (^{Li et al., 2014}). These studies have used a variety of methods to estimate air pollutant exposure and have developed deterministic and/or probabilistic models as the cornerstone of their analyses. In Mexico, the majority of epidemiological studies on the health effects of air pollutants have estimated pollutant exposure using proximity methods, such as those based on data from the nearest air quality monitor (^{Rojas-Martinez et al., 2007}; ^{Barraza-Villarreal et al., 2008}; ^{Escamilla-Nuñez et al., 2008}; ^{Hernández-Cadena et al., 2009}), averages by city or municipality (^{Téllez-Rojo et al., 2000}; ^{Carbajal-Arroyo et al., 2011}), or Inverse Distance Weighting (IDW) interpolation (^{Riojas-Rodríguez et al., 2014}; ^{Trejo-González et al., 2019}). One disadvantage of these approaches is their failure to capture the spatial variability of exposures caused by local pollution sources, urban topography and local meteorological factors (^{Jerrett et al., 2005}), which leads to considerable uncertainty regarding accurate exposure estimation. To the best of our knowledge, the available literature offers few epidemiological studies which have employed more sophisticated, high-resolution methods (e.g., based on satellite data) for estimating exposure to ambient air pollutants in the Mexican population (^{Rosa et al., 2017}; ^{Gutiérrez-Ávila et al., 2018}; ^{Téllez-Rojo et al., 2020}). Furthermore, previous studies have noted that human exposure to air pollutants can differ between individuals due to variable factors such as time spent outdoors and mobility patterns, thereby introducing bias into individual-level exposure estimates in epidemiological analyses (^{Setton et al., 2011}; ^{Nyhan et al., 2019}; ^{Yu et al., 2020}).

To offset these limitations, our study has geocoded the home and work addresses of the study population and fitted generalized additive models (GAMs), which predict PM_2.5 and NO₂ concentrations at both locations. The principal advantage of the GAMs is their use of semiparametric methods to model non-linear functions through penalized splines (^{Yanosky et al., 2008}).

This study represents a refinement of exposure models previously developed and applied in the MCMA (^{Just et al., 2015}; ^{Rivera-González et al., 2015}; ^{Son et al., 2018}). In order to craft a feasible model adaptable to the limited information available in developing countries, we used the highest-quality and largest amount of secondary source data available on air pollutant predictors in Mexico. This study therefore provides a novel methodological tool to be used in the Mexican Teachers’ Cohort (MTC) study, as well as for subsequent epidemiological studies which seek to assess the effects of PM_2.5 and NO₂ exposures on health. Our aim was to determine the level of monthly exposure to outdoor PM_2.5 and NO₂ of approximately 16 500 MTC participants living and working in the MCMA between 2004 and 2019.

2. Material and methods

We developed spatiotemporal prediction models for PM_2.5 and NO₂ using the pollutant concentrations at monitoring locations as the dependent variable, and meteorological and geographic variables, as well as other pollutants, as predictors. For statistical modeling, geostatistical modeling, geoprocessing and spatial analysis, we used R software version 3.5.3 (^{R Core Team, 2019}).

2.1 Study population and geocoding

This study is nested within a larger study known as the MTC, an ongoing prospective observational study established in 2006-2008 that includes 115 314 female public-school teachers from 12 Mexican states, with a median age at enrollment of 44 years. A detailed description of the full cohort has been published elsewhere (^{Lajous et al., 2015}). In brief, recruitment was divided into two phases: the first (2006) took place in two states (Jalisco and Veracruz, n = 27 979), and the second (2008) incorporated 10 additional states (Baja California, Chiapas, Mexico City [CdMx], Durango, Guanajuato, Hidalgo, Mexico State [EdoMex], Nuevo León, Sonora and Yucatán, n = 87 335). Participants completed a reference questionnaire with sociodemographic characteristics (including home address [HA] and work address [WA]), reproductive history, diet, lifestyle and health status. Follow-up takes place every 3-4 years to obtain information on any disease diagnoses and update participant risk factor profiles. The first follow-up cycle (2011-2013) obtained a response rate of 83%, and the second cycle (2014-2020, ongoing) of 61%.

For the present study, we selected a subsample of the MTC which excluded participants with the following characteristics: (1) from states other than CdMx or EdoMex (n = 92 571), (2) without a valid HA or WA on file (n = 6071), or (3) with an HA or WA unable to be geocoded inside the area of influence of the MCMA environmental monitoring network (n = 265, methods described below). The analysis included a total study population of n = 16 407 MTC participants.

The original MCMA is composed of 16 municipalities of CdMx, 59 of EdoMex and one of Hidalgo. We did not include the entire MCMA due to the low/under-representation of pollutant measurements at locations far from monitoring stations, which do not achieve full coverage of the area (Fig. 1). Based on this determination, we drew the boundaries of the analysis area using the extreme coordinate position of a 5 km radius around each of the most peripheral monitoring stations. Within these boundaries, we defined a grid composed of 1 × 1 km cells (in gray, Fig. 1).

Fig. 1 Area of analysis and geocoded addresses. (a) Location of work addresses at baseline (2008) and first follow-up cycle (2011-2013). (b) Location of home addresses and number of teachers per zip code. 

We geocoded the HA and WA of each participant within the study population. For HA, we geocoded zip codes (ZC) as reported at baseline. Based on the Mexican Postal Service (SEPOMEX) criteria for valid ZCs, we selected only those containing five digits and verified that the first two corresponded to the codes of home municipalities also reported by teachers. Finally, we used the geographical layer of SEPOMEX’s ZCs (^{Gobierno de México, n.d.}) to assign the coordinates of the ZC’s centroid to each HA. Since information on HA changes was unavailable for before and after the baseline, we assumed that HA were constant throughout the follow-up period. To geocode WAs at baseline and follow-ups, we used the workplace codes provided by the National Teaching Career Program. The geographic coordinates (longitude and latitude) were matched from the National School Information System (^{SEP, n.d.}).

3. Results

3.1 Study population and geocoding

The national-level MTC included 115 314 total participants, of which 22 743 lived either in the CdMx or EdoMex. From this subsample, we were able to geocode HA and WA within the MCMA of 73% (n = 16 672). In terms of HAs, the median of teachers per ZC was 23, with a range of one to 175 (Fig. 1). The areas of these ZCs ranged from 0.005 to 17.003 km². Therefore, the final study population included 16 407 MTC participants who lived and worked in the area of analysis (Fig. 2), covering (totally or partially) 72 municipalities in the CdMx or EdoMex.

Fig. 2 Flow chart of geocoding process and identification of study population. (MTC: Mexican Teachers’ Cohort; MCMA: Mexico City Metropolitan Area; SEPOMEX: Mexican Postal Service; HA: home address; WA: work address; WC: workplace code). 

3.3 Generalized additive model for PM _2.5 and predicted values

The PM_2.5 model predictors were altitude, relative humidity, wind speed, season and year, in addition to the pollutants PM₁₀ and NO₂. Only PM₁₀, relative humidity and wind speed were included as continuous variables with smoothing functions, using penalized splines of up to seven degrees of freedom (Fig. 3).

Fig. 3 Smooth plots from the final models showing non-linear effects of selected predictors and associated predicted PM_2.5 (uptrend) and NO₂ (downtrend) concentrations. CI is set at two standard errors above and below the mean. 

The maps in Figure 4 illustrate the average spatial distribution of the monthly PM_2.5 levels predicted by the model per year on the surface of the area of analysis (the grid). In general, predicted levels were higher further north within the MCMA and trended downwards until 2014. The 2019 map considers only measurements until July; therefore, the rainy season, which yielded the lowest PM_2.5 levels (see Table I), was underrepresented. The percentage of deviance explained by the model was high (DE = 87.3%) for the 2004-2019 period (Table I). The RMSE in the sample indicated high predictive accuracy with a low level of error (RMSE = 0.100), and showed no overfitting of data, as determined by its proximity to the cross-validation error (CV-RMSE = 0.102).

Fig. 4 Raster of the annual average of the monthly PM_2.5 and NO₂ predicted concentrations on the surface of the area of analysis (grid made up of 1 × 1 km cells) (2004-2019). 

3.5. Exposure assignment adjusted for mobility

All teachers reported a HA different from their WA, which strengthens the case for considering mobility when estimating air pollutant exposure. For the time period and population analyzed, the ranges of differences in estimated exposure at the homes and workplaces of participants were 0-6.01 µg m^-3 for PM_2.5 and 0-14.76 ppb for NO₂. Only 0.11% (n = 1799) of teachers changed WAs between the first and second follow-up cycles. On average, they traveled 6.285 km (SD = 6.216) between their homes and their workplaces (min: 0.0098 km, Q1: 1.846 km, Q2: 4.282 km, Q3: 8.746 km, max: 56.155 km).

We obtained 187 weighted monthly averages of PM_2.5 and 187 of NO₂ assigned to each teacher as a proxy for long-term exposure. Figure 5 illustrates the annual distribution and principal statistics for the weighted monthly averages (HWA _weig ) assigned to the 16 407 MTC participants during the 2004-2019 period. For both pollutants, we determined the expected exposure variability attributable to the spatial distribution of pollutants and the location of HAs and WAs within the study area. On average, for the population and period analyzed, weighted exposures were 24.38 µg m^-3 (SD = 6.78) to PM_2.5 and 28.21 ppb (SD = 8.00) to NO₂.

Fig. 5 (b, d) Violin plot and descriptive statistics for assignment of monthly PM_2.5 (B and D) and NO₂. (a, c) weighted exposures () by year. First, second and third quartiles estimated for n = 16 407 teachers. 

Average exposures to both pollutants were lower at HAs than at HWA _weig (PM_2.5: HA 24.37 µg m^-3, HWA _weig 24.39 µg m^-3; NO₂: HA 28.15 ppb, HWA _weig 28.21 ppb). Average differences were -0.016 (range -2.00-1.91) for PM_2.5 and -0.063 (range -4.92-4.60) for NO₂. In both cases, the differences between HA and HWA _weig did not equal zero and were statistically significant (p < 0.0001).

4. Discussion

The models developed in our study provide a useful method for predicting monthly exposure to outdoor PM_2.5 and NO₂ at any location within large areas throughout the MCMA, with high spatial (1 × 1 km) and temporal (16 years: 2004-2019) resolution. These models performed well and delivered high predictive accuracy, with cross-validation errors slightly higher than those of the models (CV-RMSE = 0.102 for PM_2.5 and CV-RMSE = 4.497 for NO₂). A large difference between these values would have suggested that GAMs overfit the data, as performance would have been lower in data not included in the adjustment during cross-validation. Because of their high resolution, these models satisfy the permanent need of epidemiological studies to predict and assign exposures to ambient air pollutants in numerous MCMA locations with a high degree of precision. Furthermore, they cover a long time period (16 years) on a small temporal scale (monthly), allowing the assessment of effects of medium- and long-term exposures.

To our knowledge, only three published studies have proposed models for determining exposure to ambient air pollutants in Mexico (^{Just et al., 2015}; ^{Rivera-González et al., 2015}; ^{Son et al., 2018}), of which all focus on Mexico City and/or the MCMA. ^{Rivera-González et al. (2015)} provide a comparison of four basic methods based on proximity and interpolation for all criteria pollutants. The principal limitation of several of these approaches is that they do not capture the spatial variability of exposure (^{Jerrett et al., 2005}). Our study overcomes this limitation by using high spatial resolution, which captures much greater variability. Son et al. (2018) proposed a mixed-effects regression model based on land use for all criteria pollutants. Although this study follows the traditional assumption of a linear relationship between predictors and the dependent variable, it is innovative in that it explores the statistical model-fitting techniques of least absolute shrinkage and selection operator (LASSO) method, in addition to analyzing various spatiotemporal scales (hourly, daily, monthly, semi-annual and annual). Instead of this model, our study proposes GAMs, which maximize the quality of dependent variable prediction by accounting for its non-linear relationship with the predictors, and incorporate smoothed (non-parametric) terms using penalized splines. However, we explored a monthly exposure scale for only two criteria pollutants. Finally, ^{Just et al. (2015)} developed a more sophisticated mixed-effects regression model, including satellite-based aerosol optical depth measurements as primary PM_2.5 predictor, over a period of 10 years. While innovative, this method is limited in that it cannot be applied where little or no PM_2.5 monitoring data exist for validation, and also because no satellite measurements can be taken on cloudy days (^{Paciorek and Liu, 2009}; ^{Bravo et al., 2012}). We overcame this first limitation, an obstacle faced by many developing countries, by using information on predictors from secondary sources which are easily accessed and interpreted, and are permanently monitored. Furthermore, we covered a longer time span (16 years) and also modeled NO₂ exposure. Despite important distinctions from previous studies, we obtained largely similar results; predicted spatial patterns for PM_2.5 were comparable and, on average, predicted concentrations for both pollutants (PM_2.5 and NO₂) were only slightly lower.

The results of our analyses demonstrate the relevance of using meteorological variables and other pollutants as primary PM_2.5 and NO₂ predictors. This could indicate the feasibility of generating highly predictive models using only secondary monitoring data from local networks. Our study therefore confirms that, despite limited availability of information, it is possible to design feasible methods for predicting outdoor exposure to air pollutants, including for other pollutants (such as PM₁₀) when no data is available on PM_2.5, as is the case in multiple cities within developing countries. In terms of spatiotemporal prediction of PM₁₀, we expect to observe predictive accuracy as high as that of PM_2.5, since in urban areas both fractions of PM have similar potential sources, though in different proportions (^{Karagulian et al., 2015}). Furthermore, geographical predictors are a key consideration given the nature of the contaminants of interest and the geographical characteristics of the study area. In the two models we developed, the geographical predictors showed smaller coefficients than the other predictors. In addition to altitude and the vehicle motorization index, we generated and tested other geographical variables. The linear meters of different types of avenues in buffers of 100-, 300- and 500-m diameters around the monitoring stations did not prove to be significant and therefore were not included in the final models. We tested the UV-B radiation variable and found a significant relationship with NO₂; however, it was excluded from the final model due to insufficient monitoring sites in the study area (only 0-6 stations during the 2004-2019 period) which did not allow a semivariographic analysis, causing an increase of 0.583 units of CV-RMSE in the NO₂ model (CV-RMSE including UV-B = 3.914).

A recent study which analyzed and compared the home-work mobility patterns of approximately 400 000 individuals, affirmed that ignoring mobility causes erroneous classification when estimating health effects (^{Nyhan et al., 2019}). The principal strength of our study was in considering home-work mobility to estimate PM_2.5 and NO₂ exposures. Given the findings of ^{Nyhan et al. (2019)}, we expect that this consideration will help reduce measurement errors that could lead to null association measurements in subsequent epidemiological studies. Another strength of our study was the method used to ascertain the spatial distribution of predictor variables at unmeasured sites. Although no existing method can offer an exact estimate of an unknown reality, GCS considers the dispersion of the phenomenon analyzed and allows the obtention of one among the possible spatial of a random function (surface). The simulated variables thus present the same variability and correlation characteristics as the observed data (measured according to their mean, variance and semivariogram) (^{Goovaerts, 1997}; ^{Deutsch and Journel, 1998}). Other potential methods considered included IDW interpolation or kriging; however, despite their extensive application in earth and environmental sciences, these provide only a smooth image of reality, whereas GCS assumptions include seasonality and a normal data distribution (^{Goovaerts, 1997}; ^{Deutsch and Journel, 1998}).

We are aware that inherent precision and measurement errors exist in predictive models such as the one we developed, which depend not only on mobility but also on the residential history of participants and precise geocoding of addresses. Prior studies have suggested that the geocoding technique can influence health effect estimates when using fine-scale exposure models (^{Jacquemin et al., 2013}). Therefore, we recognize three sources of error. First, we assumed that the HAs of cohort participants were permanent throughout follow-ups, supported by a previous study by ^{Pimienta-Lastra and Toscana-Aparicio (2019)} that registered low inter-municipal migration within the MCMA (4.8%) from 2010 to 2015 (Pimienta-Lastra and Toscana-Aparicio, 2019), thereby suggesting that our study population was similarly unlikely to change addresses. Second, we assessed HAs by geocoding the centroids of their corresponding ZCs as opposed to using exact HA coordinates, and ZC areas varied widely (0.005-17.003 km²) among participants. Previous studies have documented that estimating air pollution exposure by certain zones such as zip codes implies greater uncertainty, as compared to estimation by exact home addresses (^{Kinnee et al., 2020}); however, since home addresses were not available, we used zip code centroids as a proxy of participant residence. Finally, we assumed that the long-term outdoor exposure levels determined were also a reasonable representation of interior and microenvironment exposure.

The models developed were subject to various limitations. First, determination of the study area and population was restricted by the quality of secondary available data (the geographical layer of ZCs) to geocode the HAs of the MTC participants. Although the cohort was extended across 12 Mexican states, we were able to partially include only two of them in this study. Second, given the suboptimal quality of available data, we did not include the characterization of pollution sources attributable to land use and vehicular traffic emissions. Nonetheless, as a proxy for the latter, we used a vehicle motorization index created with official data on registered vehicles at the municipal level. This data was used under the assumption that it would stay constant within each municipality, and inter-municipal mobility was not considered. The vehicle motorization index variable proved significant in the NO₂ model, but not the PM_2.5 model, which is consistent with our theoretical framework given that NO₂ is primarily generated by fuel combustion in vehicles. The result of including aggregated variables to the models may have been a misclassification of predictions’ spatial variability, as a result of the dilution of all possible values within each aggregated area. Finally, another important limitation to consider was the lack of monitoring data, which limits the type of statistical and spatial analyses possible. Air pollution models require datasets with sufficient data for training models and validation, therefore less available data limits the statistical techniques than can be applied to achieve prime model performance or cross-validation.

5. Conclusions

In conclusion, it is crucial for environmental epidemiology to apply methods for the estimation of ambient air pollutant exposure that can overcome limited data availability, which is a particularly relevant factor in developing countries. Our results suggest that the model presented herein offers a promising alternative to be used in epidemiological studies for predicting PM_2.5 and NO₂ exposure with high spatiotemporal resolution within the MCMA in the 2004-2019 period. Beyond the MCMA context, these models can be replicated in other cities using similar information sources. Our results highlight the importance of developing pollutant exposure models with high predictive accuracy based on available secondary data from official sources. Further research should focus on expanding these models to other Mexican cities.