Urban environment and mortality differentials in Spain

1 Institute for Economy, Geography and Demography, Center for Humanities and Social Sciences (CCHS‐CSIC), Madrid, Spain Department of Production and Dissemination of Demographic and Social Statistics, ESRI Spain, Madrid, Spain Center for Human Geography, Faculty of Geography and History, Complutense University of Madrid (UCM), Madrid, Spain Dissemination of Demographic and Social Statistics, Institute of Statistics and Cartography of Andalusia, Seville, Spain

unfavourable diets in China's urban population raises concerns about future health disadvantages. However, vaccination rates, hygiene, and access to health care are significantly better in China's urban centres than in most rural areas (Gong et al., 2012).
Traditionally, rural and urban subpopulations are compared by using aggregated data sources and dichotomous measures to distinguish between the rural and urban area types. This, however, fails to capture existing heterogeneity and interactions between environmental area features and the variable of interest, which are often masked by differences in the population composition . Considering such hidden clustering effects and accounting for positive health aspects of urban life like the closer proximity to hospitals and specially qualified doctors, it may appear as if urban residents are more likely to have a health advantage over their rural counterparts. However, urban dwellers or those who spend most of their day in a highly urbanised area are exposed to specific environmental risk factors like high levels of air pollution and the lack of access to green space (Guan, Zheng, Chung, & Zhong, 2016;Pearce, Shortt, Rind, & Mitchell, 2016). Numerous studies also link unfavourable health effects with increasing spatial proximity to areas with high exposure to negative environmental features, like natural gas wells or high relative amounts of nanomaterial (Colvin, 2003;Fergusson, 1990;Rabinowitz et al., 2015).
Thus, health and well-being of individuals nested in neighbourhoods or other forms of small areas are influenced by shared exposure to certain environmental stress factors. Such effects may occur independently of their individual characteristics and may be mediated through the degree of "urbanicity," referring to the degree or extent of urban features within residential areas Xu & Wang, 2015). Even if the majority of health and survival disparities can be traced to behavioural, socio-economic, or biological differences, accounting for exposure to environmental hazards complements the analysis of individual risk factors (Waller & Gotway, 2004).
To our knowledge, only a handful of studies have addressed disparities between rural and urban populations in large regions of high-income countries, possibly due to much lower growth rates when compared with cities in China and other vastly urbanising countries.
Even if these changes are modest, attractive job markets in or near urban centres and the growing demand for service work have led to a continuously changing population composition between rural and urban settings with regard to age distribution, education, and other wealth parameters. Consequently, population movements and area developments lead to environmental changes, which in turn affect the population's health (Green, Subramanian, Vickers, & Dorling, 2015;Morris, Manley, & Sabel, 2018;Riva, Curtis, & Norman, 2011).
Although comprehensive health insurance coverage in most highincome countries prevents large-scale health and survival disparities, the examination of indirect environmental risk factors, some of which may be specifically urban or rural, can help to explain and ideally prevent the evolving health and mortality gap between social groups. Particularly in places like Southern Europe, it is necessary to analyse phenomena like urban heat islands where the asphalt and other artificial surfaces store and accumulate summer heat during the day and thereby create a substantially warmer environment during the night for individuals residing and working in urban centres (Oke, 2011).
In the context of this work, we aim to contribute to the debate on urban-rural differences particularly in the field of small area analysis by estimating the effects of environmental impacts and urban characteristics on individual-level survival over time. First, we introduce the conceptual framework focusing on the measurement of "urbanicity." Second, we explain in more detail the construction of our index to measure the latent concept "urbanicity" and introduce the data infrastructure. Third, we fit a Cox proportional hazards (PH) model with mixed effects to estimate individual-level mortality risks and the effects of exposure to urban environment over time. We account for possible stratum homogeneity by including random area effects.
Finally, we compare our results to alternative models and discuss the main findings and limitations.

| Measuring "urbanicity" and small area environment features
The increasing availability of spatially referenced data, greatly improved computing power, and the development of more advanced statistical approaches helps generate new strategies for comparing rural and urban subpopulations. Classic studies often apply dichotomous indicators to distinguish between urban and rural areas. They commonly rely on a core set of characteristics aggregated at differently specified areas, most commonly based on administrative boundaries .
A dichotomous classification according to administrative boundaries appears to be straightforward but fails to capture part of thebetween-area variation and relevant urban characteristics related to infrastructure, geographic position, and distribution of active space (Cyril, Oldroyd, & Renzaho, 2013). In fact, cities are often surrounded by heavily populated areas that might not be part of the same administrative unit but are still exposed to similar conditions.
To address calls for a more nuanced approach to the subject matter (McDade & Adair, 2001), Vlahov and Galea (2002) proposed a refined conceptualisation for measuring urban space. The crucial step of their approach was to disentangle the two related and often synonymously used concepts urbanisation and "urbanicity." While they defined urbanisation as a process of growth of cities in terms of area and population over time, they related the term "urbanicity" to a current state of an area that can differ by degrees of certain urban characteristics like the proportion of built in environment. In other words, their concept captures the "nature of urban environments" (Dahly & Adair, 2007). "Urbanicity" is a concept strongly dependent on the regional context and changes over time (Leon, 2008). As such, it is difficult to define. The lack of a universal definition, however, provides the opportunity to propose new measures that can then be compared over time.

| Indicator of urban environment
The degree of "urbanicity," or in other words the degree of how urban an area is, refers to a multidimensional and latent concept. Accounting for the complexity of such a concept, we choose a mix between a theoretical and data-driven approach to construct a multicomponent index that allows us to distinguish between different degrees of urban environments (cf. Dahly & Adair, 2007;Jones-Smith & Popkin, 2010;Rey, Jougla, Fouillet, & Hmon, 2009;Vlahov et al., 2005). By examining graphical tests and correlation coefficients between all accessible environmental variables related to urban settings, we identify four main contributors to the latent concept "urbanicity," depicted in Figure 1.
First, we calculate the population density per census tract as a measure of relative crowdedness and standardised observations based on the overall deciles to make them comparable with other scale components. Second, we calculate and add the average coverage with medical service based on estimated service area polygons that represent the distance that can be covered between a health facility and any point on the map within 30-min driving-to-facility time. 2 Third, we use satellite imagery data from the CORINE land-cover database to estimate the artificial surface area per census tract for the year 2006. 3 Fourth, we obtain and add road density by estimating the total length of line objects (transportation networks) within area units (kilometer of road per square kilometer of land surface).
The weights with which single components enter the index are estimated through maximum likelihood factor analysis incorporating standardised single-component values (cf. Brown, 2014;Jreskog, 1967). Factor weights are represented in Table A1. The resulting index variable is further normalised and centred around zero. A Crohnbachs alpha score of 91% suggests a sufficient internal consistency of all indicator components. A graphical quality test of the index is depicted in Figure 2, which displays a choropleth map of Seville, the biggest city, population-wise, in Andalusia (INE, 2017). Scores for the FIGURE 1 Components of the "urbanicity" indicator FIGURE 2 Municipality of Seville-Scores for the urban environment indicator by census tract multicomponent indicator for urban environment are represented by census tract, and darker shades are associated with a higher degree of "urbanicity." Areas identified as urban are rather small and located in the city centre, which indicates a good graphical fit for our measure.

| Environmental and population-based area features
The rather physical "urbanicity" indicator does not allow us to capture area-specific heterogeneity regarding environmental hazards and potentially harmful population or social features unrelated to "urbanicity." As we attempt to distinguish between different kinds of urban and rural areas, it is important to account for further heterogeneity between areas similar to what we can observe in reality. To capture between-area differences regarding the latent socio-economic situation and the exposure to area-specific environmental hazards,

| Study population and individual-level information
There are two main sources for individual-level information. The BDLPA is a mortality follow-up, which is semiannually updated and corrected. By using an individual identifier, we are able to link all subjects in the study to their answers from the population and housing  (Pla & Cabrerizo, 2004).
In general, individuals between ages 35 and 80 have resided in their house or apartment for relatively long periods and are, as they grow older, more likely to own their dwelling. The probability of moving is rather low for those age groups, and individuals are more likely to be exposed to the same environment for the time of our study.
Limiting the age range was further motivated by the distribution of the event of interest as more than 90% of all deaths in Andalusia occur after age 35, but before age 92, the highest age individuals reach at the end of the follow up period in 2014 (cf. Ocaña Riola & Mayoral-Cortés, 2010;Viciana-Fernández, Ruiz-Ramos, & Pujolar, 2008). A consequence of selecting these age groups is that the sample size decreases from 723,234 to 351,769 individuals. We tested all models with different age ranges but did not observe greatly different patterns in population composition or general model outcomes. To assure that neither the observed population nor the area-specific "urbanicity"

| Statistical approach
The incorporation of area effects into an analysis of individual-level mortality differences requires statistical testing for potential impacts of cluster-specific effects and, in case of geographical data, the spatial distribution. The graphical representation of the multicomponent "urbanicity" index and statistical tests regarding the environmental variables suggest that observations are more likely to be similar if they are geographically closer to each other. To assess if observations are spatially autocorrelated, intrinsic stationarity is assumed before we calculate a row-standardised matrix of spatial weights that is based on the list of contiguous neighbours. At least one point of the boundary of a spatial polygon that represents a census tract has to be within snap distance of at least one point of a neighbouring polygons' border to meet our contiguity condition (cf. Anselin, 2013;Bivand, Pebesma, & Gomez-Rubio, 2013). We then calculate the product-moment correlation coefficient (Moran's I) as statistical test for spatial autocorrelation (Moran, 1948;Sokal & Wartenberg, 1983). Although spatial autocorrelation of mortality indicators would justify the analysis of area differences in the first place, we also estimate Moran's I for other central variables including the index for "urbanicity" and aforementioned area-specific environmental features.

| Statistical model
Following the descriptive tests for spatial autocorrelation, we estimate mortality disparities by degree of "urbanicity" and environmental impact with an extended version of the Cox PH model. The original model is the most commonly used approach to model censored time to event data, particularly when the main interest is to obtain relative effects of covariates (Mills, 2011). Such effects are assumed to be proportional over time and enter the model multiplicatively as expressed in the following equation (cf. Kleinbaum & Klein, 2010).
where h 0 is the baseline hazard and exp (β i X i ) is the nonnegative function of covariates. Hazard ratios are obtained through the maximisation of the partial log likelihood with respect to β i X i (Allison, 2014;Therneau & Grambsch, 2000). Because only the right-hand side of the formula is maximised, the Cox PH model does not require you to specify the underlying baseline distribution. Due to our assumption that individuals are nested in small areas where they are exposed to similar environmental hazards and the same degree of "urbanicity," we choose to apply an extended version of the Cox PH model, which allows us to account for such homogeneity within clusters. Thus, a stratum-specific frailty term is added to the original model (cf. Austin, 2017). The resulting Cox PH model with mixed effects can be considered as shared frailty model with a normally distributed stratum specific frailty term Z j as follows (Therneau, 2015).
where Z j is the design matrix for random effects, which captures homogeneity within clusters. The model can be interpreted as multilevel survival model with shared frailties. The added random effect term can be understood as relative effect of given covariate patterns on the baseline hazard, which varies across census tracts (Pankratz, De Andrade, & Therneau, 2005). Given the set-up of our analysis, it is necessary to account for left truncation (Cain et al., 2011). This adjustment affects survival estimates for everybody in the sample because their time under risk of dying before the start date of the study remains unobserved. In other words, we select individuals based on their survival upon the start year of the examination. To account for left-truncation and assure we measure agespecific mortality differences, we choose to use person years as the time scale in our models. Cohort effects are accounted for by including birth cohort effects as covariate (Canchola et al., 2003).

| RESULTS
In order to analyse area differences, spatial variation must be present in the variable of interest. We determine to what extent these differences are spatially associated by estimating Moran's I for the variables in our analysis, which can be interpreted as the correlation between a variable and its spatial neighbours (Anselin, Syabri, & Smirnov, 2002).
These estimates and associated p values are presented in Table 1. All variables of interest exhibit significant spatial autocorrelation, justifying our analysis of area differences in Andalusia.
As the analysis aims to highlight the impact of shared environmen- After including individual-level differences (Model 2), the effect of urban environment appears more pronounced than in Model 1. In Model 2, every unit increase in the "urbanicity" scale of a census tract increases the estimated hazard of dying by three percentage points.
Although changes in the individual-level impact factors are negligible between different models, the effect of the degree of "urbanicity" on survival varies with the incorporation of additional area-level characteristics. In Model 3, we incorporate the effects of perceived cleanliness, noise, air, and water pollution in an attempt to control for different kinds of heterogeneity between urban areas with the same degree of "urbanicity." The estimates suggest that including such environmental area features reduces the effect of the degree of "urbanicity" to 1.7 percentage points for every unit increase. Both cleanliness and pollution are found to have a highly significant but small effect on survival. Perceived cleanliness of the area is estimated      Table A2, which we also provide as a summary of the likelihood ratio tests between all models and their counterparts without the additional random effects. The test

| DISCUSSION
In this work, we aim to extend classic approaches for analysing mortality disparities between rural and urban subpopulations. To capture the rural-urban gradient and estimate "urbanicity" effects on survival, we disentangle population and physical features from environmental impact factors in residential areas. We use satellite imagery data and census information to identify four universal predictors of "urbanicity," the multidimensional latent concept that describes the "urban nature" of an area (Dahly & Adair, 2007). Our index represents an improvement over classic binary measures that are based solely on population density. The use of census tracts as a clustering unit increases the precision with regard to area size and reduces the risk of misclassifying large areas as urban if, for example, only a part of the overall area exhibits urban features. Therefore, our approach offers an advantage over comparable measures, such as the "rurality index" proposed by Ocaña Riola and Sánchez-Cantalejo (2005) in which data were aggregated at the municipality level.
We incorporate our index in mixed effects Cox PH models to estimate long-term survival according to different degrees of "urbanicity" The most recent financial and debt crisis of 2008 hit Andalusia particularly hard and led to significant job loss, a continuously increasing at-risk-of-poverty rate, and a high number of evictions (Eurostat, 2017;Romanos, 2014). Thus, the assumption of residential continuity may not hold among economically disadvantaged groups. Moreover, there is no information measuring the average exposure to the estimated environmental and social features of the residential area. The average amount of time someone spends in his or her residential area and is therefore exposed to its environment probably differs by age, employment status, and other unmeasured area features such as access to "third places" (Mehta & Bosson, 2010;Oldenburg & Brissett, 1982).
In spite of these data constraints, to our knowledge, this analysis is the first study to combine detailed small area (census tract) information, individual-level variables, and survival follow up in Southern Europe. Our results differ from a previous analysis on a larger provincial scale for Spain (Regidor et al., 2015). This other study suggested a negative association between per-capita income and average survival times, although we find that environmental features and, to a greater extent, population composition affect survival probability.