Nitrogen dioxide levels estimated from land use regression models several years apart and association with mortality in a large cohort study

The setting

Rome is the largest Italian city with a population of 2.8 million inhabitants on a surface of 1,290 km². It is a radial city, and the urban development in Rome took place gradually from the centre to the suburbs, with significant urbanization in the 1930s, after the second world war, and in the 1990s [21]. The city is divided into 5,500 census blocks with an average population of 470 inhabitants.

Air pollution measurements in 1995–96

Nitrogen dioxide (NO₂) was measured at 70 schools in a cross sectional study on children conducted in 1995 (Italian Study on Respiratory Diseases and Environment), the Italian part of the International ISAAC study [22]. The objective of the measurement campaign was to estimate average NO₂ exposure at schools in Rome and the 70 schools were a random sample of all primary and junior high schools of the city (stratified by the 20 city districts). Three Palmes diffusion tubes, manufactured by the Emilia–Romagna Environmental Protection Agency, measured outdoor air pollution in each location simultaneously over three one-week periods in June 1995, November 1995, and March 1996 [23]. The tubes were placed outside the school, two meters above ground, and the mean of the measurements for the three tubes was considered. The school mean NO₂ concentration over the three periods was considered as an estimate of the annual mean level. Only sites that had complete data coverage over the three periods were included, thus 67 sites remained for the analysis [23]. The locations were 44 urban background sites, and 23 traffic sites.

Air pollution measurements in 2007

In 2007, we measured NO₂ concentrations using Ogawa passive samplers at 78 sites during three one-week periods in February, May, and October. As in the previous survey, the samplers were positioned within two days and collected exactly after one week. We selected the same locations used in 1995/96, and we added 11 sites to increase the variability of the measurements (3 urban background, 7 traffic sites, and 1 regional background site). In fact, we realized that schools selected in 1995/96 were not enough to represent the high variability of NO₂ in Rome, and a previous LUR developed with 1995/96 measurements tended to underestimate high values of NO₂[23]. In 27 locations, we used the same Palmes tubes and Ogawa samplers to compare the different measures and calculate a correction factor between the two methodologies. An Ogawa badge is a passive sampler that collects nitrogen oxides through diffusion on a pre-coated filter [19,24]. The filter is removed in the laboratory and analyzed with a spectrophotometer for nitrites using a Saltzman reaction. Five lab blanks from each measurement campaign were kept at our department. All the analyses of the Ogawa samplers were performed at the Institute for Risk Assessment Science of Utrecht University [19]. Field blanks and field duplicates were taken to document detection limits and precision of each measurement campaign. The Palmes tubes were analyzed by the Emilia-Romagna Environmental Protection Agency.

We assigned to each sampling location a single NO₂ level, the mean of the three measurements. The map of Rome with the sites of measurement is presented in Figure 1.

Land use and traffic variables

Several land-use and traffic variables were available to characterize the locations used for air pollution measurements. Traffic flows of 2005 (only for roads with >4,000 vehicles per day) were provided by the Municipality of Rome and were used for both models. For each location, we used the variables already chosen to develop a previous LUR model for Rome [23]: the proximity to green urban areas (a binary variable, which indicated if the census block was adjacent to a green urban area, selected using CORINE land cover data), the altitude, the size of the census block, the number of residents in the census block, and the inverse of the population density for the census block (m² per person). To increase our predictive ability, we added as linear terms the distance from the city centre and the geographical coordinates (X and Y). The city center has the highest traffic and population density levels and it is expected that air pollution tends to be higher in the center and lower in the suburbs so that the distance variable should capture this aspect. On the other hand, the orography of the area is peculiar (the Mediterranean Sea at South-West and the hills towards East and North-East ) and the prevalent wind direction tend to have a West- East direction. Therefore, the X and Y coordinates may capture the influence of orography and wind. In addition, we considered three GIS variables to characterize the traffic that surrounds each location: meters of high traffic road (more than 10,000 vehicles/day) within a buffer of 150 meters, traffic density (number of cars multiplied by the length of roads within 150 metres divided by the buffer area), distance from the closest high traffic road. For sake of simplicity, only one buffer was chosen (150 metres) considering the decay gradients of NO₂ with the distance to high traffic roads [15].

We collected and stored all geographical variables using ArcGis 9.1 (ESRI, Redlands, California, USA). We used the Word Geodetic System of 1984 with the Universal Transverse Mercator 33 N (WGS84 UTM33N) as the coordinate system and map projection.

The Rome Longitudinal Study (RoLS)

The Rome Longitudinal Study is a cohort based on the Rome Municipal Register’s data, described in detail elsewhere [25]. Briefly, it is a closed prospective cohort of residents of Rome on the 21^st October 2001 with a mortality follow up to the end of 2006. All 2001 residential addresses were geocoded using the Italian road network (Tele Atlas, Italy). Several individual socio-demographic variables at baseline were available for each subject of the cohort (including the residential history, education, occupation, place of birth, marital status). In addition, we used a small-area (census block, average population: 470 inhabitants) composite index of socioeconomic position (SEP), built with 2001 Census data [26], to better characterize social deprivation.

We selected all adults aged 45–80 years in 2001 who had not changed their residence in the five years preceding enrolment and during the follow-up. We linked subjects to the Regional Mortality Registry, and collected the date of death and the underlying cause of death (coded to the International Classification of Diseases (ICD) revision 9). We evaluated their vital status until December 2006 and analysed the association between NO₂ exposure (estimated with the two LUR models) and natural mortality (ICD-9 < 800), taking into account individual characteristics (age, gender, marital status, place of birth, education, occupation) and area-based socioeconomic position, as a surrogate measure of lifestyles and access to health services.

Statistical analysis

We used data from the 27 sites with both Palmes tubes and Ogawa samplers in 2007, to compare the measures of ambient NO₂ performed with different methods. Therefore, we applied a correction factor derived from a linear regression to all 1995/96 Palmes tube measurements in order to align them with the measurements collected with the Ogawa passive samplers.

We evaluated the correlation between NO₂ measurements in 1995/96 and 2007 with the Pearson correlation coefficient, and their agreement using the Bland-Altman plot. To build the period-specific LUR models, the association between the logarithmic transformation of NO₂ concentrations and each land-use or traffic variable was assessed by univariate and multiple linear regressions. We used the logarithmic transformation because the log-normal NO₂ distribution. The final model was constructed using a manual backward stepwise elimination procedure (p > 0.20): we eliminated from the model the predictors one-by-one on the basis of the highest p-value. When one of the coordinates was selected, we kept into the model the other coordinate independently on its p-value. Variables with a regression slope of opposite sign to what was expected were not included; in other words, we excluded all the variables known to be positively associated to air pollution, which seemed to have a negative effect instead.

The validity of the regression models was evaluated by leave-one-out cross-validation, i.e. by systematically subtracting each of the data points from the model one by one, and then comparing the predicted value for each specific point with the measured level at the location without that value in the model [27]. We plotted measured and predicted levels from the cross-validation in 1995/96 and 2007, and we calculated the root mean square errors (RMSE) and R² of regression analysis between measured and estimated concentrations. We checked each model for influential observations, performing the Cook’s distance, for multicollinearity by estimating the variance inflation factor (VIF), and we calculated the Moran’s Index of residuals, to check their spatial autocorrelation.

We assessed the ability of LUR models to predict concentrations both prospectively (using the 1995/96 model to predict measured concentrations in 2007), and retrospectively (using the 2007 model to predict 1995/96 measured concentrations).

We estimated the NO₂ concentrations for all residence addresses of the participants of the Rome Longitudinal Study using the two LUR models. We used the Pearson correlation coefficient to measure the correlation between the results of the LUR models from 1995/96 with those from 2007 at all the residential addresses of subjects enrolled in the cohort. We used a test for trend across ordered groups for comparing levels of NO₂ among age and socioeconomic categories, and a t-test for comparing the exposure levels between genders.

Finally, we applied Cox proportional hazard regression models to analyse the association between NO₂ exposure (from both LUR models) and natural mortality, while adjusting for age, gender, marital status, place of birth, education, occupation and area-based socioeconomic position. We categorized the predicted NO₂ levels into quintiles and into four pre-specified categories (<=35 μg/m³, 35–45 μg/m³, 45–50 μg/m³, >50 μg/m³). We calculated hazard ratios of natural mortality for 10 μg/m³ NO₂ and for interquartile range.

Land use regression models

The univariate and multivariate association between land-use and traffic variables and the log of NO₂ concentration at the 67 sites in 1995/96 and at the 78 sites in 2007 is shown in Table 2. Multiple regression models indicated that the following variables predicted measured NO₂ in 1995/96: altitude, geographical coordinates, distance from the city centre, inverse population density of the census block, traffic density in the 150 m buffer around the site, and distance from the closest road with more than 10,000 vehicles daily. The best-fitting multiple regression model resulted in a determination coefficient (R²) of 0.737, and a root mean square error (RMSE) of 0.076 (i.e. 1.09 μg/m³).

The following variables were associated with NO₂ concentrations in 2007: altitude, the X-coordinate, distance from the city centre, size of the census block, and meters of high traffic road within 150 m of the site. The best-fitting multiple regression model resulted in a determination coefficient (R²) of 0.704, and a RMSE of 0.132 (i.e. 1.14 μg/m³).

For both models there was no evidence of multicollinearity among the chosen variables or of influential observations: the mean variance inflation factor was 1.55 (the maximum was 2.08 for Y-coordinate) in 1995/96 and 1.39 (the maximum was 1.92 for Y-coordinate) in 2007, the Cook’s distance was less than 0.3 in both cases. There was no evidence of spatial autocorrelation of residuals: the p value of Moran’s Index was 0.300 in 1995/96 and 0.218 in 2007.

Cross-validation of the two models

For the first period, the determination coefficient (R²), the adjusted-R² and the root mean square error of regression analysis between measured and estimated concentrations were 0.67, 0.66, and 3.24 respectively.

In 2007, the determination coefficient (R²), the adjusted-R² and the root mean square error of regression analysis between measured and estimated concentrations were 0.61, 0.61, and 5.38 respectively.

Comparison of the 1995/96 and 2007 observed and predicted levels

Figure 2 shows the scatter plots of the comparison of the measurements of NO₂ concentrations taken at the same locations in the two study periods (A), the ability of 1995/96 model to predict 2007 measurements (B), the ability of 2007 LUR model to predict 1995/96 measurements (C), and the comparison of predicted values of the two LUR models at the addresses of the study population (D). The Pearson correlation coefficient between the two measurements of ambient NO₂ was 0.79, it was higher for urban background sites (0.83) than for traffic sites (0.68).

The predictions of 1995/96 LUR model and measurements in 2007 had a correlation of 0.73, while the 2007 LUR predicted levels had a correlation with 1995/96 measurements of 0.83. When we applied the estimates of NO₂ concentration derived from the two LUR models to the residential addresses of the 684,204 subjects enrolled into the study, the correlation between the estimated measures in the two periods was 0.96.

NO₂ exposure for cohort members

Figure 3 shows the maps of levels of NO₂ from the two land use regression models in Rome. The socio-demographic characteristics of the study population and the levels of exposure according to the 1995/96 and 2007 models are presented in Table 3. According to the 1995/96 model, the 684,204 cohort members were exposed to a mean level of 45.7 μg/m³ (SD 5.9), the level was slightly higher in women than in men (p < 0.001), and was higher in the elderly than in the 45–64 year old population (p-trend < 0.001), and in subjects living in high compared to low socioeconomic position census blocks (47.2 vs. 42.9 μg/m³, p-trend < 0.001). According to the model developed with 2007 data, similarly to the 1995/96 model, the mean exposure was 43.9 μg/m³, and levels of exposure increased with increasing age class (p-trend < 0.001) and increasing socioeconomic position (p-trend < 0.001).

Association with natural mortality

During the study period, a total of 45,006 deaths from natural causes occurred in the cohort. Table 4 shows the results of the Cox regression analyses adjusted for several confounders. When we considered the quintiles of the two distributions as the relevant exposures, we obtained very similar results of increased mortality from natural causes with increasing quintiles of NO₂ at residence (HR = 1.10, 95%CI:1.06-1.13 in the 5^th quintile compared to the 1^st for 1995/96 estimates, and HR = 1.11, 95%CI:1.07-1.14 in the 5^th quintile compared to the 1^st for 2007 estimates). There was a similar association with natural mortality even when we used a 4-category variables with a 12 % increased risk of dying (95%CI: 8–18 %) for those with a concentration of NO₂ at residence of more than 50 μg/m³, compared to those exposed to less than 35 μg/m³. The increased mortality risk for an increase of 10 μg/m³ of NO₂ concentration in 1995/96 was 6 % (95%CI:4–8 %), while using NO₂ estimates for 2001 it was 4 % (95%CI:3–5 %). This result was related to the wider variability of NO₂ concentrations in 2007; in fact the hazard ratios of mortality for an increase of NO₂ corresponding to the inter-quartile range was 1.04 (95%CI: 1.03-1.06) for both periods.