Spatial Relationship Quantification between Environmental, Socioeconomic and Health Data at Different Geographic Levels

Spatial health inequalities have often been analyzed in terms of socioeconomic and environmental factors. The present study aimed to evaluate spatial relationships between spatial data collected at different spatial scales. The approach was illustrated using health outcomes (mortality attributable to cancer) initially aggregated to the county level, district socioeconomic covariates, and exposure data modeled on a regular grid. Geographically weighted regression (GWR) was used to quantify spatial relationships. The strongest associations were found when low deprivation was associated with lower lip, oral cavity and pharynx cancer mortality and when low environmental pollution was associated with low pleural cancer mortality. However, applying this approach to other areas or to other causes of death or with other indicators requires continuous exploratory analysis to assess the role of the modifiable areal unit problem (MAUP) and downscaling the health data on the study of the relationship, which will allow decision-makers to develop interventions where they are most needed.

approach yielded smaller prediction errors and more precise and accurate probability intervals and that it allowed for better discrimination between counties with high and low mortality risks.
Poisson kriging, in this context, presents a spatial methodology that allows for filtering the noise caused by the small number problem and enables the estimation of mortality risk and the associated uncertainty at different spatial scales. This approach has been implemented to modeling cancer risk by a number of authors: Oliver et al. [11] studied cases of cancer in children under fifteen years of age, and Goovaerts and collaborators considered lung cancer [12,13], breast cancer [14,15], prostate cancer [16], cervical cancer [17], and pancreatic cancer [18], and all found it to be relevant for this particular problem.
Selection of scale is perhaps the most important factor in creating and analyzing a relationship between environmental exposure and health outcomes [19]. This issue is similar to the modifiable area unit problem (MAUP), a term introduced by Openshaw [20,21]. The MAUP can cause differences in the analytical results of the same input data compiled under different zoning systems [22,23].
The present study aims to evaluate spatial relationships at three levels of aggregation: the IRIS level, an intermediate scale (the grid level), and the county level between health outcomes (mortality attributable to cancer) initially aggregated to the county level, district socioeconomic covariates, and exposure data modeled on a regular grid. The approach is illustrated using age-adjusted lip, oral cavity and pharynx, and pleural cancer mortality rates over the period 2000-2009 for the Picardy region. The deprivation index and trace metal exposure indicators are used as putative risk factors.

Study Area
The region of Picardy covers an area of roughly 19,500 km 2 and is located between the North Artois, the Ile-de-France in the south, the Bay of the Somme to the west and the East Champagne. It covers the departments of Somme, Oise and Aisne. The urbanization rate in this region is far below the national average (60.4% versus 74% for the entire country). The agricultural sector provides more than 4% of French agricultural production. It also has significant industrial activity through which fine chemicals and specialty chemicals account for nearly 15% of jobs and the vehicle industry comprises 40% of industrial employment (26.5% of assets employed in industry versus 19.5% at the national level). Three administrative or statistical spatial units, of different sizes, were considered: IRIS (the smallest administrative units in Picardy, 2,129 units) with irregular sizes and shapes, 64 km 2 grid cells (308 units) that are all squares of same size, and counties (112 units) with irregular sizes and irregular shapes. Figure 1 shows the counties of the study area.

Exposure Indicators
The environmental indicators (inhalation and ingestion) used were those described in Caudeville et al. for building GIS-based modeling platforms for quantifying human exposure to chemical substances [24]. Trace elements (nickel-Ni, cadmium-Cd, and lead-Pb) were modeled within the Picardy region [25] for various exposure pathways: atmospheric contaminant inhalation and ingestion of soil, vegetation, meat, eggs, milk, fish and drinking water.

Deprivation Index (SE)
The deprivation index used was developed by Rey [26] and was built at the French census block (IRIS) using the following socioeconomic data: the median household income, the percentage of high school graduates in the population aged 15 years and older, the percentage of blue-collar workers in the active population, and the unemployment rate. The deprivation index was also constructed for the county. For each county, the deprivation index was calculated as the population-weighted average score for all of the IRISes in the county.

Health Data
The health data came from the Regional Health Observatory of Picardy [27], where the age-adjusted mortality rates are calculated for each county from 2000 to 2009. Table 1 shows the cumulative, maximum and minimum number of mortality and age-adjusted rates per 100,000 person-years by county from 2000 to 2009.  Figure 2 shows the spatial distribution age-adjusted lip, oral cavity, pharynx, and pleura cancer mortality per 100,000 person-years. It should be noted that (1) the population is not evenly distributed throughout the study area and (2) the age-adjusted rate calculated from the less-populated counties tend to be less reliable. This implies that the interpretation of the map must be carried out with caution. The scatter plot at the bottom of Figure 1 illustrates this effect, well-known as the -small number problem‖. Table 2 presents the different spatial scales of measurement and the approaches used to homogenize spatial coverage.  To correct for the instability attributable to the small number problem, a number of algorithms have been developed that aim at estimating risk. The geostatistical approach, in this context, presents an interesting alternative; it conducts the noise filtering and allows for risk estimation along with the associated uncertainty at different scales. This section provides a brief overview of the geostatistical methodology for estimating risk values. See Goovaerts [17] for more details about this approach.
The cancer mortality count d(v  ) within a county v  is interpreted as the realization of a random variable D( v  ) that is Poisson distributed with a parameter (expected number of counts) that is the product of the population size n( v  ) by the local risk R( v  ). R( v  ) might be thought of as a noise-filtered mortality rate for area v  , which we also refer to as the mortality risk. It is estimated by using a variant of kriging with nonsystematic errors known as Poisson kriging [28].
The mortality risk and the associated kriging variance for an area v  are estimated as: Kriging variance is computed as follows: x represents either an area v  (ATA kriging). The kriging weights () i  and the Lagrange parameter μ ( v  ) are computed by solving the Poisson kriging system of equations: where ij  = 1 if i = j and 0 otherwise. The -error variance‖ term, m * /n( i v ) leads to smaller weights for rates measured over smaller population sizes. The covariance ( , ) i s is approximated as the population-weighted average of the point-support covariance () R C h computed between any two discrete locations between the areas i v and j v .

Spatial Autocorrelation
Global Moran's I was calculated for all of the explanatory variables as well as for the dependent variables within three spatial structures to determine the role of spatial representation using global spatial autocorrelation. The Global Moran's I spatial autocorrelation statistic measures the self-similarity of a spatial variable's value as a function of adjacency [29], using a first-order Queen's case spatial weight matrix and 999 permutations.

Exploring the Relationships between Health, Environment and Socioeconomic Factors
Analyses of correlations between health data and putative factors are traditionally performed using a global or -aspatial‖ regression model, under the implicit assumption that the impact of variables is constant over the entire study area. This assumption is likely unrealistic for large areas, which can display large geographic variations. Fotheringham and colleagues developed Geographically Weighted Regression (GWR) to explore spatial non-stationarity and map statistics to visualize the spatial patterns of the relationships between dependent and independent variables [30][31][32].

Aspatial Regression
The explanatory power of SE and exposure indicators was first investigated using the following multiple linear regression model: where is the kriging risks estimate for observation i, 0 is the intercept, is the regression coefficient (slope) of each factor , and is the error term. To account for the reliability of the kriged risks in the regression, each observation receives a weight that is the reciprocal of the kriging variance [33].

Geographically Weighted Regression.
In geographically weighted regression, the regression is conducted within local windows centered around each observation. The regression coefficients are thus location-dependent: Within each window, observations are weighted according to their proximities to the center of the window. A variety of distance decay functions are available. In this paper, we used the XX function, which is characterized by a bandwidth that corresponds to the distance beyond which the weight rapidly approaches zero.
The bandwidth is estimated by minimizing the AICc value: where n is the number of observations in the dataset,  is the estimate of the standard deviation of the residuals, and tr(S) is the trace of the hat matrix. For more information on the theory and practical application of GWR, the reader is referred to Fotheringham et al. [34]. All maps are substantially smoother than the original rate map because the noise caused by small population sizes has been filtered. These maps allow a better visualization of areas of higher risks: the lip, oral cavity and pharynx cancer mortality rates vary between 2.81 and 37.40 per 100,000 person-years. After the application of Poisson kriging, the minimum rate increased from 2.81 to 8.87 deaths/100,000 person-years, and the maximum rate of 37.40 decreased to 25.14 deaths per 100,000 person-years. We can note, for instance, that the high rates recorded in sparsely populated counties such as Sains-Richaumont (37.40 deaths/100,000 person-years), north of the Aisne department, are strongly smoothed (24.46 deaths/100,000 person-years). The highest rate recorded in a densely populated county (i.e., Abbeville North county-26.60 deaths/100,000 person-years) remained nearly the same after smoothing (24.90 deaths/100,000 person-years). Zero pleural cancer mortality rates recorded in sparsely populated counties were also smoothed, leading to minimum values of 1.00 deaths/100,000 person-years.

Poisson Kriging for Health Indicator
The maps of the kriging variance indicate the higher reliability of risks estimated in densely populated areas such as Amiens, Beauvais, Saint Quentin, and Abbeville. The variance of the risk estimates decreased as the geographic unit area increased: from the IRIS level to the grid level and then to the county level ( Table 3).
The risk estimates are characterized by positive spatial autocorrelation within the three spatial scales (p ≤ 0.05) but display low levels of statistically significant spatial autocorrelation at the IRIS level in comparison with the grid and county levels ( Table 3). In this case, the counties are internally homogeneous in terms of mortality according to the risks estimated by kriging.

Spatial Aggregation for the Explanatory Variables
The mean values of the explanatory variables under each of the three spatial structures are similar ( Table 4). The variables F1 (exposure inhalation indicator) and F2 (exposure ingestion indicator) display greatly reduced variability under different spatial structures. The variable SE has the highest levels of variability. In contrast, the SE index is affected by the use of different spatial structures. The lowest average was −4.5 at the county level in comparison with −7.32 at the IRIS level, and the variance decreased with increasing aggregation, unlike with F1 and F2, for which the lowest and strongest averages were somewhat similar for the three spatial scales.
All of the variables were characterized by positive spatial autocorrelation within the three spatial scales at levels of p ≤ 0.05. For F1 and F2, the Moran's I values of the three spatial structures were similar (Table 4); however, the F1 and F2 were not affected by the use of different spatial structures. This is explained by the fact that the exposure indicator presented a homogeneous distribution within each county (the original resolution of exposure indicator data was a grid (15 × 10 km); see Caudville et al. [24,25]. Figures 5 and 6 show maps of the SE index and the F1 aggregated at (a) county level, (b) grid level and (c) IRIS level.   Table 5 shows the correlation coefficient between each covariate and mortality rates before and after noise-filtering using Poisson kriging. Filtering the noise because of the small number problem clearly enhanced the explanatory power of the covariates: the proportion of variance explained (adjusted R 2 ) increased by nearly one order of magnitude: lip, oral cavity and pharynx cancer mortality, 0.11 to 0.26, and pleural cancer mortality, 0.11 to 0.25. The uncertainty attached to the risk estimates can be accounted for by weighting each estimate according to the inverse of its kriging variance, leading to correlation coefficients and adjusted R 2 values that were slightly larger for pleural cancer mortality and a slightly lower mortality rate for lip, oral cavity and pharynx cancers.

Aspatial Regression
Application of the linear model for lip, oral cavity and pharynx cancer mortality data explained a moderate proportion of the total variance (adjusted R 2 = 0.22 and 0.19) at the county level and grid level, respectively. This proportion was lower when the analysis was conducted at the IRIS level (adjusted R 2 = 0.11). It is noteworthy that the correlation coefficients for the SE factor were always significant for the different aggregation levels and were higher at the county level than at the IRIS level, which was the expected result because aggregation is known to increase the strength of correlation [20,35]. Linear association between SE index and cancer mortality has been demonstrated in other works [36][37][38][39].  The model explains a moderate proportion of the total variance when the dependent variable y is pleural cancer mortality for different levels of aggregation. The adjusted R 2 ranges between 0.24 and 0.28, with significant correlation coefficients of up to 0.5 for F1. The results showed the consistency of the association between trace metal inhalation exposure and pleural mortality across the different aggregation levels. This is explained by the fact that the pleural mortality presented a homogeneous distribution within each county and the exposure indicator was not affected by the use of different spatial structures. Pleural mesothelioma is one of four types of mesothelioma, but it accounts for approximately 75 percent of all diagnoses of asbestos-related cancers. The disease starts in the pleura, the soft protective tissue surrounding the lungs, which can be attributed directly to its cause: repeated or heavy occupational exposure to airborne asbestos fibers. However, Peterson suggests that a significant number of cases of this cancer are apparently not asbestos-related and that even in the absence of asbestos, other agents can induce malignant mesothelioma [40]. Some types of nanoparticles, especially those containing nickel, could also be a risk of pleural diseases [41].

Geographically Weighted Regression (GWR)
In the aspatial analysis, we implicitly assumed that the impact of covariates was constant across the study area. This assumption is likely unrealistic for large areas, which can display substantial geographic variation in socioeconomic and environmental conditions. Therefore, the global statistics reported by this traditional regression model could potentially hide a number of interesting local relationships. The question is then to examine whether there are any meaningful spatial variations in these relationships.
County-level data were used to optimize the bandwidth of the GWR distance decay function. Figure 7 shows the relationship between AICc and bandwidth size for the two types of cancer. The following bandwidths were found to be optimal: lip, oral cavity and pharynx cancer mortality (47 km), and pleural cancer mortality (54 km). The local regression was based on the following number of closest observations, which represented 15%-20% of the data: 22 for the county level, 39 for the grid level, and 426 for the IRIS level. A comparison of the AICc values suggests that all of the GWR models outperformed the global model (Table 6).
These results strongly suggest that the relationships between cancer (lip, oral cavity and pharynx-pleural) mortality and the environmental and deprivation indexes are not stationary but instead vary over the study area. The application of GWR clearly enhances the explanatory power of the covariates for the three spatial levels: the proportion of variance explained (adjusted R 2 ) is almost doubled ( Table 6).  The three spatial scales share the same southeast-northwest trend in the explanatory power of the local regression models for lip, oral cavity and pharynx cancer mortality: the lower mortality values in the south are better explained by the SE index than is the higher risk recorded in the northwest. As a recall, the largest R 2 observed in the south (Figure 8 a) corresponds to areas with low SE index variation. Similar to the global model (Table 5), GWR led to local correlation coefficients and R 2 values that were higher at the county level than at the IRIS level.
The lower pleural cancer mortality values are better explained in areas of low F1 variation (Figure 9a). The maps of the local correlation coefficients and the R 2 values also present the same spatial structure over the different spatial coverages of aggregation: lower in the west оf tһе Aisne department and higher in the north of the Oise department. Very similar results were obtained for the different spatial scales.

Discussion
Our results substantiate the work on noise filtering described in the introduction section from Oliver et al. [11] and Goovaerts et al. [12,18]. Indeed, we found the following: (1) the mean risk values estimated under each of the three spatial structures were similar; (2) the IRIS risk estimates were non-negative; (3) and their sums were equal to the original county risk. Although the disaggregation of cancer data on a small scale is somewhat arbitrary, in particular when it does not take into account secondary information to guide this downscaling, the approach should facilitate the analysis of the relationships between health data and the putative covariates (i.e., environmental, socioeconomic, behavioral or demographic factors) that are typically estimated for different spatial scales. These covariates can potentially be subsequently used as secondary information in the kriging, leading to more detailed risk maps at finer scales [42].
The other issue was the so-called modifiable areal unit problem (MAUP), for which different geographic scales can lead to inconsistent results for relationships analysis. For example, the mortality rate reported at (1) the county level requires an aggregated deprivation index at the same resolution, and this aggregation obscures the intra-county variation and thus the relationship and (2) the IRIS level, at which the disaggregation leads to a large variance in estimated risk. Exploratory methods, such as the univariate Moran's can serve as indications of the potential effect of the MAUP in the study how relationships based on the homogeneity and heterogeneity of spatial data are affected by the study level and may affect the ability of the study to detect a relationship [43].
In this study, very similar results were obtained for the different spatial scales between:  pleural cancer mortality and the exposure indicator F1.  lip, oral cavity, pharynx cancer mortality and the SE index.
Whereas other studies on the relationship between heath and deprivation showed that the use of spatial representations other than the census tract produced different analytical results [35,44], the significant association between lip, oral cavity, pharynx cancer mortality and the SE index estimated using the county structure were stronger than they were under the IRIS structure.
This difference between the significant correlation coefficients is the result of aggregation because the aggregation level is known to increase the strength of correlation [20]. Future work could provide tools to exhaust all possible aggregations and generate empirical frequency distributions of the statistical estimates that could be used to evaluate the sensitivity of results to aggregation effects.
Based on the results obtained, we can confirm that the presence of this significant statistical association was likely not induced by the use of a particular geography. At the three spatial scales, the strongest correlation coefficients were found where low deprivation was associated with low lip, oral cavity and pharynx cancer mortality and where low environmental pollution was associated with low pleural cancer mortality.

Conclusions
This paper presents an approach for evaluating spatial relationships between health outcomes (mortality attributable to cancer) initially aggregated at the county level, district socioeconomic covariates, and exposure data modeled on a regular grid. The approach was illustrated using age-adjusted lip, oral cavity and pharynx, and pleural cancer mortality rates measured over the period 2000-2009 for the Picardy region. The deprivation index and trace metal exposure indicators were used as putative risk factors. For the different spatial scales, the strongest associations were found where low deprivation was associated with low lip, oral cavity and pharynx cancer mortality and where low environmental pollution was associated with low pleural cancer mortality. However, applying this approach to other areas, for other causes of death, or with other indicators always requires exploratory analysis to assess the role of the MAUP and downscaling health data in the study of the relationships that will allow decision-makers to develop interventions where they are the most needed.