Modelling geographical variations in fertility and population density of Italian and foreign populations at the local scale: a spatial Durbin approach for Italy (2002–2018)

Studies on fertility determinants have frequently pointed to the role that socio-economic, cultural and institutional factors play in shaping reproductive behaviours. Yet, little is known about these determinants at an ecological level, although it is widely recognised that demographic dynamics strongly interact with ecosystems. This research responds to the need to enhance the knowledge on variations in fertility across space with an analysis of the relationship between fertility and population density of Italians and foreigners in Italy at the municipal level for the period 2002–2018. Using global and local autocorrelation measures and a spatial Durbin model, we show that there is a negative association between the fertility and population density of the Italian population, while the density of foreigners is correlated with higher fertility. This second result poses new insights on the relationship between space and fertility. Moreover, we find that the features of neighbouring areas, measured by population density, contribute significantly to explaining spatial fertility variation, confirming the importance of the study of spatial diffusion in demographic processes.


Introduction
The geographical distribution of the human population is widely acknowledged to be a crucial variable in shaping the environment and social behaviours (Bocquier and Costa 2015;Boyle 2003;White 2016). Spatial variations in the distribution of individuals are the result of multiple interconnected factors engendering a complex system (Xu and Cohen 2019;Yang et al. 2015). Natural properties or resources partly explain the dissimilarities in population densities. Nevertheless, the socio-demographic characteristics of human groups also provide other equally relevant explanatory elements. A recurrent theme in the contemporary debates on sustainable development concerns the question of optimising human densities according to what a territory can bear in terms of its congestion risks and biocapacity (Costanza et al. 2014;DeFries et al. 2004;Rees and Wackernagel 1996). This issue strongly highlights the need to consider, on the one hand, the ecological footprint of the inhabitants (Ferguson 2002;Mancini et al. 2018;Wackernagel and Yount 2000), and, on the other hand, the demographic dynamics which can be markedly different depending on the country (Salvati et al. 2019). Therefore, density is not a univocal concept related only to the demographic dimension, but rather it features in important discussions among diverse, interconnecting disciplines. In the field of population studies, human density is identified as a fundamental measure for two main reasons. Firstly, it reflects the spatial models of the different population distributions that may pose the question of territorial equity regarding a more spatially balanced and therefore sustainable development (European Commission 1999). Secondly, it reveals trends of demographic phenomena within a territory through the events associated with them (births, deaths, migrations). While extensive literature has recently developed on the spatial variation of human populations, including research focusing on alternative measures (Benassi and Naccarato 2019;Cohen et al. 2013;Benassi 2018, 2020;Xu and Cohen, 2019), to our knowledge, studies devoted to understanding the relationship between population density and demographic phenomena are relatively scarce. This is particularly true with reference to fertility. In fact, despite a strand of literature dating back to the 1970s that points to the influence of population density on human fertility (Easterlin 1976;Firebaugh 1982;Janowitz 1971;Leet 1977), over recent decades few studies have investigated population density as a determinant (Lutz et al. 2005(Lutz et al. , 2006Lutz and Qiang 2002;Rotella et al. 2021;Testa 2003).
These last studies present two major aspects that can be discussed: the territorial level investigated and the results achieved. Concerning the first dimension, in all of these, the geographical scale of analysis is the macro one (i.e. country level). This scale is necessary when the aim of the study is an international comparison (i.e. Lutz et al. 2006;Testa 2003). Nevertheless, for an analysis conducted at national level it may represent a serious limitation in understanding the intimate relationship between population density and fertility. This is because both indicators are characterised by a very deep and important geographical variability that needs a control at a microscale territorial level (Burillo et al. 2020;Salvati et al. 2020). With regard to the second aspect concerning the results attained, all studies on this topic agree in finding a negative effect of population density on fertility and therefore, they confirm the importance of using this variable as a crucial determinant in modelling fertility processes and dynamics (Lutz et al. 2005(Lutz et al. , 2006Lutz and Qiang 2002;Testa 2003).
It should be noted that in recent years, research on fertility based on spatial approaches has also developed in Italy. Using geographically weighted regression models and spatial panel regressions, some studies have investigated the link between socio-economic factors and fertility, modelling their spatial dependence and spatial diffusion patterns at provincial level (Mucciardi and Bertuccelli 2013;Vitali and Billari 2017). Mucciardi and Bertucelli demonstrated the importance of local effects in modelling the spatial variation of total fertility rate, paying attention to the use of local regression models (i.e. geographically weighted regression model, GWR). In their pivotal paper on this topic, Vitali and Billari (2017) investigated the association between fertility and some cultural and economic indicators across provinces showing that it varies spatially and over time. Furthermore, in the framework of the diffusionist approach on the fertility decline in Europe, their work used local and global spatial regression models to prove that fertility in a given area depends not only on the economic, institutional and cultural characteristics of that area but also on those of neighbouring ones. In both studies, the territorial level of analysis was the provincial one, that is to say a territorial level intermediate between the local (municipality) and the regional. Another recent work using Morans' indices to measure spatial fertility dependence at different territorial levels  suggests that the spatial variation of fertility reflects its response at the local level to economic expansion and recession over the last decade. This study was conducted both at provincial and municipal (local) levels but without any 'causal' approach, i.e. without using any kind of regression model. This paper aims to fill these gaps, on one hand, and to complement and update the knowledge on this topic, on the other, by examining, at the local level, the relationship between the spatial variation of fertility and population density in Italy from 2002 to 2018, a particular period that includes pre-and post-economic crisis years. In doing this, we refer to the population density of the total population and separately, to that of Italians and foreigners. To the best of our knowledge, the originality of our work consists precisely of investigating, for the first time in Italy, this relationship adopting, as a territorial scale, the municipality one.
It is important to underline that we are not directly interested in detecting the main drivers of fertility (a topic also studied by many important papers in relation to Italy) but rather to investigate the spatial effects (direct and indirect) of population density on fertility at a local level controlling for the composition of population (Italians and foreigners). To pursue this purpose, we adopted a spatial approach using univariate and bivariate global and local measures of spatial autocorrelation and applied a spatial Durbin model that allowed us to investigate whether and to what extent the fertility rate of a certain area may also be associated with the characteristics of its neighbouring areas beyond its own socioeconomic or cultural features. The key conceptual framework underpinning our analysis is synthetised by the assumption that, 'Everything is related to everything else, but near things are more related than distant things' (Tobler 1970, p. 236), also known as the first law of geography.
Our findings show that the spatial density of Italian and foreign populations influences the fertility of the municipalities in Italy in a different way and they also reveal the role played by the neighbouring areas in shaping spatial fertility variation.
The paper is structured as follows: In the next section, theoretical reflections on the relationship between population density and demographic dynamics are proposed in a very brief overview; in Sect. 3, the data and methods are presented; finally, in Sect. 4, the discussion and conclusion are presented.

Understanding the relationship between population density and demographic dynamics
One of the most common aims of the current research on human population distribution is to understand and assess its environmental impact. It is commonly known that an increase in densities is accompanied by population growth, and this relationship can produce diverse implications that shape the features and configurations of territorial spaces (Polinesi et al. 1 3 2020). Over time, settlement patterns and the question of the spatial concentration of populations have stimulated an intense and unresolved debate on the 'weaknesses' of urban sprawl and the 'virtues' of dense cities (Johnson 2001;Nguyen 2010). Currently, population density triggers both a fear of an excessive concentration of human populations within some geographical areas and that of depopulation within others. Thus, in the global context of sustainable development, scholars and policy makers have highlighted the need to maintain or regenerate a sustainable socio-ecological system (Ostrom 2009;Zuindeau 2007;Bastianoni et al. 2012). Nevertheless, the relationship between human population distribution and ecosystems is not linear and univocal (de Sherbinin 2007;Hunter 2000); rather, it is mediated by demographic dynamics which are associated with population density. Behaviours related to fertility, mortality and territorial mobility are themselves endogenous drivers of density that, in turn, can imply significant changes in settlement patterns. Consequently, exploring the geographical and temporal distribution of individuals through density is a step toward better understanding the socio-ecological system in which demographic behaviours of intertwined human societies and ecosystems strongly interact (Xu and Cohen 2019). Certainly, population density is a measure that presents a high degree of spatial variability. Furthermore, the geography of spaces, in which densely populated and uninhabited areas alternate, often makes arduous the comprehension of the logic that governs the spatial distribution of the human species. From this point of view, the study of population density should focus on the relevance of the spatial variation of demographic patterns. For shaping trends and the intensity of these patterns, fertility is a crucial assessment of the demographic dynamics in an area. Persistent disparities in fertility rates along geographical spaces reveals a heterogeneity of reproductive behaviours that, investigated from both a diachronic and synchronic perspective, refer to a specific geographical context in which reproductive preferences can be influenced by multiple factors. Most of the literature on this topic has examined spatial variations in fertility on different geographical scales, highlighting the diffusion of its low level in Europe (Arpino and Tavares 2013;Burillo et al. 2020;Campisi et al. 2020). Some authors have focused on the determinants of a decline in birth rates, pointing to the role played across space by economic , socio-cultural and political factors (Campisi et al. 2020;Vitali and Billari 2017). Conversely, the research on the possible predictors of fertility at the ecological level remains inconsistent. To date, while an important effort has been devoted to the analysis of density-dependent mechanisms of population increase (Burillo et al. 2020;Cohen 2003;Lima and Berryman 2011;Sibly and Hone 2002), only a few works have examined the spatial relationship between demographic density and fertility. The findings of these studies demonstrate that density affects reproductive behaviours (Loftin and Ward 1983), reducing fertility preferences (Lutz et al. 2006) and also interacting with other covariates (Campisi et al. 2020). However, several questions regarding the need to employ this kind of analysis data at a finer geographic scale, as well as to perform a spatial model that considers the features of neighbouring areas, remain to be addressed. In this context of analysis, Italy represents without doubt an interesting and unique laboratory of analysis because of the remarkable heterogeneity of its territory and the strong territorial differences in the temporal and spatial distribution of the population and fertility dynamics at both regional and sub regional levels (Benassi et al. 2021;Reynaud et al. 2020;Salvati et al. 2020). Our study aims to contribute to the advance of the specific knowledge on this subject by answering to the following research questions:

RQ1.
Is fertility (measured here using general fertility rate, GFR hereafter) a phenomenon that is distributed randomly across space in Italy at the local level, or does it follow spatial patterns? RQ2. Are the spatial patterns of GFR spatially correlated with the density patterns? RQ3. What kind of relationships exist between GFR (y) and population density (x), controlling for spatial biases of y and x?
Additionally, it is important to note that Italy presents a significant geographical variation in fertility behaviour patterns that has been associated somewhat with the economic and cultural diversity observed at the regional level along the north-south axis (Kertzer et al. 2009;Zambon et al. 2020) but that, explored at a local level, reveals further 'fragmented and mixed micro-scale behaviours' ). In addition to these peculiarities, Italy is experiencing the lowest-low level of fertility in Europe (Billari and Kohler 2004;Kohler et al. 2002): the annual birth rate has fallen steadily since the 1960s and 1970s, despite a rebound recorded in the late 1990s and early 2000s. This downward trend is explained by a decline in the average number of children per woman (total fertility rate) and by a reduction in the percentage of the female population in their fertile years (Mencarini and Vignoli 2018). In Italy, women born in the 1960s who have completed their reproductive lives are now being replaced by a much smaller-sized cohort (Istat 2020). As a result, in Italy, like in all southern European societies characterised by an elderly population and a fast ageing process, the childbearing of immigrant women is recognised as a crucial factor underpinning the demographic functioning of the societies (Giannantoni and Strozza 2015), which can also strongly shape the spatial fertility patterns. In this framework, an important (or even crucial) role is played by the foreign population and its interplay with demographic dynamics at the micro and macro levels. There is strong evidence that the territorial distribution of the foreign population (more than 5 million in Italy currently) is quite different from the Italian population, which presents important spatial patterns Strozza et al. 2016). From this perspective (i.e. detecting differential effects and/or relationships between demographic density of foreign and autochthon populations and fertility levels), we modelled two kinds of population density: one related to foreign populations and another related to Italian ones in order to answer the fourth research question.

RQ4.
Is the relationship between GFR (y) and population density (x) the same for Italian and foreign populations?

Data and methods
In this section, we describe in detail the data used and methods applied in order to provide a comprehensive analysis of the level of spatial dependence of demographic density and fertility at a local level in terms of univariate and bivariate distribution, and to estimate the spatial variation of fertility as function of population density considering the spatial lagged effect of both the dependent variable and independent one. 1 3

Dependent and independent variables
As in the pioneering paper by Loftin and Ward (1983) on the effect of demographic density on fertility in the USA, the dependent variable is the s.c. general fertility rate (GFR). For each statistical unit (i), i.e. each Italian municipality, GFR was obtained as the ratio between the number of births and the mean female population aged (15-49 years) within a time interval. So that: GFR indicates the number of births in a municipality per 1000 women aged 15-49 years in that same area. In our case, the numerator of the ratio was the number of births that occurred in each Italian municipality between 2002-2018. The denominator of the ratio was obtained as the mean between women aged 15-49 years residing in each municipality in the two years of 2010-2011 (01.01). The latter can be considered as a sort of mean population for the whole period (2002)(2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010)(2011)(2012)(2013)(2014)(2015)(2016)(2017)(2018). As known demographic rate was computed using a mean population as the denominator. Here, the period of observation was from 2002 to 2018; therefore, 2010 and 2011 can be considered as the central (mean) years of the entire period. That is why we computed the mean population of 2010 and 2011 (01.01) as a proxy of the mean whole period population in each municipality. We were obviously aware that this indicator (GFR) has limitations. Indeed, it ignores who is at risk of having births and partially, in our case, because we were considering the female population aged between 15-49 years, the age structure of the population. Nevertheless, this indicator is the only one that can be computed at municipality level as other measures, although more accurate, are not available as this finer level ). The independent variable was the population density (D) that represented, for each municipality, the number of residents per square kilometre. It was obtained as: for which the numerator was the population residing in each municipality (i), and the denominator was the surface of the same municipality (i) expressed in square kilometres. As mentioned above, we computed three types of c densities: one related to the total population (Eq. 2), and the other two that considered whether the population was foreign or Italian. These two densities were simply obtained by counting, as the numerator, the foreigners (average population between 2010-2011) or Italians (average population between 2010-2011) separately. Italians and foreigners were identified using the variable of country of citizenship. The source of all demographic data (births and residential population broken down by gender, country of citizenship, and age) was the intercensal population reconstruction provided by the Italian National Institute of Statistics. The source of the geographical data (shape files) was the Italian National Institute of Statistics ('territorial bases'). Statistical units were the 7904 Italian municipalities. As 14 of them were neighbour less, i.e. islands, the spatial analysis was carried out on 7890 municipalities.

Methodologies
The methods applied belong to spatial demography tools and techniques (Matthews and Parker 2013). In the first part of the paper, we explored the spatial distribution of GFR by using global and local univariate versions of Moran's index I (Anselin 1995;Moran 1948). (1) We were interested in verifying whether the distribution of GFR was random or not and in identifying local spatial clusters. Next, we ran the same indices but in a bivariate way (Anselin et al. 2002) to detect the spatial correlation (global and local) between the two observed variables (i.e. between GFR and D). In more detail, a positive (and significant) global univariate Moran's I index indicates spatial dependence; conversely, a negative (and significant) global univariate Moran's I index indicates spatial dispersion (Moran 1948). Maps of the local version of univariate Moran's I (Anselin 1995) were used to identify any significant clusters of GFR and D that could be characterised as either homogenous spatial regimes (a combination of high-high or low-low attribute values among neighbouring areas as defined by the spatial weight matrices) or spatial outliers (derived from high-low or low-high patterns). Both of these indices were computed separately for the two variables, GFR and D. On the contrary, the bivariate version of the indices allowed consideration of the spatial distribution of both variables simultaneously. In particular, the global bivariate Moran's I measures the correlation between a variable and a spatial lag for another variable. It should be noted, however, that the index does not consider the inherent correlation between the two variables. In other words, the index measures the degree to which the value for a given variable at a location is correlated with its neighbours for a different variable (Anselin 2020). The local version of the index can reveal the spatial correlation (or even disparity) of the relationship between two variables, which in our case were: the GFR and D. For the attribution of spatial weights, we used a first-order 'Queen' based contiguity matrix so that two territorial units (municipalities in our case) were considered neighbouring if they shared a boundary or a vertex geographically. In all spatial autocorrelation analyses, the variables were expressed in a standardised form, such that their means were zero and their variance one. In addition, the spatial weights were row standardised. The hypothesis of the existence of a condition of univariate and bivariate spatial clustering was tested at a 5% level of statistical significance (p-value ≤ 0.05).
In the third part of the analysis, we modelled the relationship between y (GFR) and x (D) by using a spatial Durbin model (Anselin 1988). The spatial Durbin model (SDE, hereafter) has been proven to outperform more classic autoregressive spatial econometrics' models (like spatial error or spatial lag) because, as underlined in Yang et al. (2015), it is the only means of producing unbiased coefficient estimates (Elhorst 2010), it has no restriction imposed on the magnitude of spatial effects, and both global and local effects are produced (LeSage and Pace 2009; Yang et al. 2015). The relative scarcity of demographic studies that use this kind of model (Sabater and Graham 2019;Vitali and Billari 2017;Yang et al. 2015) is quite surprising and the fact that, to our knowledge, no application of this approach has been made for the study of fertility in Italy at the local (municipality) level. 1 As underlined in Sabater and Graham, who applied the SDE model for the study of fertility variation and international migration in Spain during the Great Recession, this kind of model is particularly useful in approaching demographic process as spatial (diffusion) process. In particular, the spatial Durbin model includes both spatially lagged dependent variables (like in the classic spatial autoregressive model, or SAR) and spatially lagged explanatory variables as well (LeSage Space 2009). This latter point is crucial because our task here was to estimate the spatial lagged effect of GFR on GFR but, especially, the spatial lagged effect of D on the dependent variable. Following Yang et al.

3
(2015), three components comprise a spatial Durbin model (LeSage and Pace 2009): a spatial lagged dependent variable; a set of explanatory variables of a spatial unit; and a set of spatial lagged explanatory variables, which can be expressed as Eq. 3: where y denotes an n × 1 vector of the dependent variable (i.e. GFR); W is the spatial weight matrix; W y represents the spatial lagged dependent variable (endogenous interaction relationships); and denotes the effect of W y , which is known as the spatial autoregressive coefficient. l n indicates an n × 1 vector of those associated with the intercept parameter α. X represents an nxk matrix of k explanatory variables (population density, in our case), which are related to the parameter . WX reflects the spatial lagged explanatory variables (exogenous interaction relationships), and denotes a kx1 vector of the effects of WX . The error term, , follows a normal distribution with a mean of 0 and a variance 2 I n , where I n is an nxn identity matrix. Equation 3 clearly demonstrates that the characteristics of a specific unit (the municipality, in this study), and its neighbours are simultaneously considered in the analysis (Yang et al. 2015). This paper used this approach to explore whether the GFR of a municipality was related to the features of its neighbours and, if so, to answer how they were related. Specifically, we investigated whether the GFR of a municipality was associated with the features of the surrounding municipalities after accounting for the characteristics, in terms of demographic density, of the specific municipality. As in the case of spatial autocorrelation analysis, the spatial weight matrix ( W ) was obtained as the 'Queen' contiguity matrix of the first order. One crucial aspect of SDE, and in general of all spatial autoregressive regression models, is the interpretation of the coefficients that cannot be interpreted as in an OLS model, but rather it is necessary to refer to direct and indirect (spatial spill overs) effects (Golgher and Voss 2016). The direct effect, 'represents the expected average change across all observations for the dependent variable in a particular region due to an increase of one unit for a specific explanatory variable in this region' (Golgher and Voss 2016, p. 185), while the indirect effect, 'represents the changes in the dependent variable of a particular region arising from a one-unit increase in an explanatory variable in another region' (Golgher and Voss 2016, p. 185). As usual in spatial regression analysis, we took an OLS model (aspatial regression model) as a benchmark and a tool of comparison of the results of the spatial model of regression (Benassi and Naccarato 2017;Burnham and Anderson 2002;Cutrini and Salvati, 2021;Yang et al. 2015) using the Akaike Information Criterion, AIC (Akaike 1974) as the parameter of comparison (i.e. the model with the smaller AIC should be preferred). From an interpretative point of view, the comparison between OLS and SDE was straightforward. The OLS model parameters were estimated under the explicit assumption that the observations were independent, which means that changes in values for one observation did not 'spill over' to affect values of another observation (Golgher and Voss 2016). Spatial regression models, alternatively, assume that observations are not independent and that they can exert a reciprocal influence, as the first law of geography by Tobler (1970) clearly states. Global and local univariate and bivariate spatial analysis was carried out with GeoDa (version 1.18 10.12.2020). Regression models (OLS and SDE) were estimated by using R Studio (Mendez 2020). Thematic maps were created using Qgis 'Odense', version 3.20.2.

Results
In the presentation of the results, we first describe the global and local univariate spatial autocorrelation patterns of selected variables and then we focus on the bivariate global and local spatial autocorrelation between both. Then, we show the results of the regression models.

Descriptive results
The spatial distribution of GFR was clearly not random because Moran's I was equal to 0.44 (RQ1). This means that the spatial distribution of GFR was characterised by a positive spatial autocorrelation: territorial units with the same values of GFR (similar level of fertility) tended to cluster and, therefore, to stay close in space. However, the spatial distribution of GFR was not, as in the past, simply linked to the north-south divide but rather to different types of heterogeneities ). This was evident in terms of local spatial autocorrelation ( Fig. 1. GFR Lisa maps). In fact, in the cluster we can clearly appreciate an urban-rural divide. Most of the large Italian cities that are the capitals of metropolitan areas were HH clusters: Milan, Turin, Bologna, Rome, Naples, Messina, Catania and Palermo, among others. A coastal-non-coastal divide emerged: coastal municipalities were more present in the HH cluster while non-coastal municipalities were more present in the LL cluster. Particularly evident were the s.c. 'inner areas' (groups of municipalities characterised by a low level of accessibility and located in inland areas) that were in most cases classified as LL clusters (i.e. comparatively low values of GFR were closed in space). Local clusters of positive spatial autocorrelation (LL and HH) were formed out of a total of 2447 municipalities (31% of the total territorial units) while negative spatial autocorrelation (LH and HL) was detected in 271 municipalities (3.4%). In the rest of the municipalities no local spatial association for GFR was detected; therefore, we can assume that there the spatial distribution is a random process. It is important to note that in absolute terms the HH cluster was the biggest group (1267 municipalities).
Population density presented a more marked geographic variability compared to GFR, and a higher level of global spatial autocorrelation (Moran's I now equal to 0.72). It is interesting to note that, contrary to GFR, most local clusters now belonged to the LL cluster (2108 municipalities), i.e. where comparatively low values of demographic density were spatially closer. The territorial distribution of this cluster of municipalities spanned across Italy, especially in the inner part of the s.c. 'Appennino' mountain chain and in the border areas (mountain areas across the 'Alpi'). A major portion of the LL clusters were concentrated in Sardinia (one of the two island regions of Italy), one of the less urbanised regions of Italy, except for the metropolitan municipality of Cagliari and its surrounding area, which are the red polygons on the southern part of the island of Sardinia. The HH clusters were instead quite concentrated and referred to the major urban municipalities like Milan and its surrounding area, Turin, Rome, Florence, Naples, Bari, Cagliari and Palermo. It is quite clear that, in some cases, HH clusters of GFR overlapped with HH clusters of D, but this occurred only for urban areas and large municipalities. In most cases, except for the northeast quadrant, an LL GFR corresponded to LL D. Obviously, as expected, important differential aspects both in the levels and in the variability characterised the demographic densities of Italians and foreigners (Table 1). From a spatial and statistical point of view, the complexity of space and the relevance of local scale analysis clearly emerged. In fact, bivariate local spatial correlations between GFR and D (total, Italian and foreign population) revealed important spatial patterns and territorial heterogeneities (RQ2). In the three cases, global bivariate Moran's I was positive, indicating the existence of a positive spatial correlation between GFR and D. The level was not particularly high (around 0.21), with the highest value recorded between GFR and the demographic density of the foreign population (0.23). The importance of the urban context clearly emerged with the HH clusters that were located precisely in the major urban municipalities and their surrounding areas, confirming the relevance of the spatial dimension (Fig. 2). This aspect was more evident regarding the foreign populations, which presented a sort of continuum of HH clusters from Milan to the north-eastern part of Italy (in the Veneto region). This is one of the most industrialised and urbanised areas of Italy .

Regression models
When modelling the GFR in function of D, important differences emerged between the OLS (i.e. not spatial models) and the SDE (spatial model), Table 2. In the first case, the relationship between D and GFR was positive, although low (RQ3). This means that where D increased, the level of GFR also increased. In the SDE, the effect of D on GFR was not statistically significant, while the lag D was positive and statistically significant. This  seems to indicate the relevance of the areas surrounding the municipality (in terms of demographic density level). In fact, the indirect effect of D was positive and statistically significant, underlying the existence of spatial spill overs. In terms of performance, the SDE was better than the OLS (because the AIC of the SDE was lower than that of the OLS). As a matter of fact, as mentioned before, it is important to make a distinction between the demographic density of Italian populations and foreign populations ( Table 3).
The results of Table 3 are extremely interesting (RQ4). The SDE model outperformed the OLS model, as the AIC for the SDE was lower than for the OLS (93,666 versus 95,868). In the OLS model, both demographic densities had a positive effect on GFR, but in the case of Italians it was very weak (0.01). Conversely, for the demographic density of foreigners, the impact was higher (0.5). Things changed significantly when the spatial dimension was introduced in the model. In the SDE model, in fact, the effect of population density of the Italian population on GFR was negative and comparatively quite weak. This means that high values of population density for the Italian population corresponded to low values of GFR. These were typically the urban areas when the density of the Italian population was comparatively higher. Obviously, this can depend on the different age structures and the propensity of the Italian population to live in an urban context as compared to the Italian population that live in semi-urban and rural contexts and, moreover, by the extremely difficult condition of the real estate market in urban areas. This is the intrinsic mechanism that has boosted suburbanisation in recent decades, and which explains the difficulty of living in Italian cities where the cost of living is comparatively high, and the quality of life is comparatively low. The opposite holds for foreign populations: the density of foreigners still had a positive effect on GFR, confirming the importance of this variable in detecting fertility levels and patterns. The direct and indirect effects, again, were different between the two groups of populations. In the case of demographic density of foreign populations, both direct and indirect effects were positive, while in the case of demographic density of Italian populations the direct effect was negative while the indirect effect was positive. This means that in both demographic densities, spatial spill overs (indirect effects) were positive while the effects were divergent in terms of direct impact. Differential aspects also arose between Italians and foreigners in terms of the effect of the lag covariate variable (lag D). In the case of Italians, the effect was positive and statistically significant, while in the case of foreigners the effect of the lag D on GFR was not significant. This implies a different spatial structure characterising the two populations, confirming the spatial matters in terms of spatial distribution (and attitudes) of the different populations that live in each territory (de Castro 2007). Finally, the Rho parameter ( ) was positive and statistically significant. This indicates the existence of a spatial diffusion pattern of GFR (in terms of contagion) because it is influenced (positively) by the variation of the same variable in the neighbouring spaces, which confirms the nature of the spatial vicious (or virtuous) circle of demography that characterises the Italian territory.

Discussion and conclusion
In this paper, we analysed the role that population density can play in shaping fertility patterns in Italy. Our findings show that the heterogeneity of the Italian territory and the different distributions of the population interact with fertility, which can be influenced in different ways. This study develops new knowledge by showing that the heterogeneity in fertility is also correlated to different population groups. Indeed, the analysis provides a deeper understanding of fertility in Italy by examining its spatial variation at the municipal level in relation to the population density of both Italian and foreign populations. We found that while there was a negative association between fertility and the population density of Italians, the density of foreigners was correlated with higher fertility.
The first result is in line with earlier ones (Lutz et al. 2006;Testa 2003) and is also confirmed by a recent study in which it was observed, for 116 of 174 countries, that increases in population density were associated with declines in fertility rates (Rotella et al. 2021). In explaining these findings, the authors suggest that, in more dense and safer environments, individuals must possess higher skills and knowledge to compete with their counterparts. This requires longer times to achieve higher education and suitable work which in turn determines a delay in family formation and the transition towards parenthood. Moreover, Rotella et al. (2021) argued that in high-density contexts the delay in family planning is also associated with the choice of couples to have fewer children to concentrate more time and resources per child who, in a more competitive environment, require greater investments. This assumption derives from the life history theory which, providing an interpretation of the strategies adopted by individuals and families to optimise their available resources, suggests that human behaviours may depend on ecological constraints.
The research conducted so far could offer some responses to why higher population densities of Italians lead to lower fertility rates. Nevertheless, in contrast to these results, our study also shows a positive relationship between population density of foreigners and fertility rates. Plausible explanations for this relationship could be identified in some cultural factors that would contribute to invert the effect of increasing density on fertility, influencing the life history strategies of foreigners independently of other environmental dimensions or constraints. Indeed, Rotella et al. (2021) demonstrated that, in some countries observed in their study, the links between density and fertility were attenuated by key moderator variables such as religiousness and social norms. In accordance with these findings, we can assume that where population densities of foreigners are higher, cultural values and religious traditions could be stronger and more decisive in orienting reproductive behaviours, especially during the immediate post-settlement period in the host territory (Carella et al. 2021;Mussino et al. 2015;Mussino and Strozza 2012). Furthermore, the strength of ties and cultural norms is even more relevant among foreigners who intentionally opt for concentrated settlement patterns to ensure mutual support and to be cohesive in addressing the challenges produced by the environment and the host society.
In performing the analyses, our study goes beyond what might be learnt from previous works by using a spatial Durbin model as an alternative approach to study the relationship between population density and fertility at the local level. This allows us to understand how the features of the neighbouring areas (assessed by population density) contribute to defining the spatial fertility variation.
The results of this article highlight the fact that it may be inappropriate to make generalisations about the quantum and tempo of fertility if it is not investigated using a spatial approach. As for all demographic phenomena, studies on fertility cannot have meaning outside a geographical context where the number of births or the birth rates are recorded or where reproductive preferences and behaviours manifest themselves. Conclusions about fertility differentials are very likely to depend upon the way in which fertility is measured and upon the groups that are being investigated. As a result, in assessing fertility, spatial structures should be considered as a leading variable to better explain reproductive preferences.
Additionally, our study leads to two important considerations. First, in evaluating the relationship between demographic dynamics and population density, it is also important to note that the latter shifts according to the level at which it is observed. The choice of a finer spatial scale may complicate the investigation for several reasons (e.g. a lack of data), but it implies a different perception of demographic phenomena. Indeed, a spatial analysis conducted at the local level can reveal forms of differentiation that otherwise would remain hidden within the overall values of indicators measured on national or regional scales.
Second, in the analysis of the interactions between a demographic phenomenon and ecological factors such as population density, it is equally important to use an adequate spatial approach that can include the greatest number of factors influencing these interactions.
Comparing our results to those of studies on the spatial patterns of fertility in Southern Europe (Burillo et al. 2020;Carioli et al. 2021;Sabater and Graham 2019), it is important to note similar evidence that relies on the important role played by geographical environments in shaping reproductive behaviours.
Recent research concerning Spain highlighted that fertility trends in this country, characterised by a sharp decrease and the delay of childbearing, are associated with important differences in the space (Carioli et al. 2021). In showing that the heterogeneity of fertility in Spain is also spatially driven, another work investigated the influence of the territorial context on local fertility using ecological variables measured at the municipal scale (Burillo et al. 2020). Indeed, the authors examined the spatio-temporal changes in reproductive behaviour over the periods of economic expansion and recession documenting diverse subnational patterns of fertility. They found a sharp increase in fertility in the municipalities with medium and high population densities and a persistent decline of its rates especially in depopulated rural districts. These results contradicted those of the previous and recent studies that investigated the relationship between ecological variables and fertility at the territorial macroscale level. Another study explored the impact of unemployment, immigration and emigration on fertility across Spanish provinces during the Great Recession (Sabater and Graham 2019). In their paper, using a spatial Durbin model, Sabater and Graham (2019) demonstrated that the spatial dependence of fertility decisions can derive not only from the features of a given area but also from those of neighbouring ones. Furthermore, employing the same spatial approach and confirming similar results, recent studies have also documented, in other European countries, that the influence of the variables related to neighbouring areas can influence reproductive behaviours. Vitali et al. (2015) provided evidence on the spread of childbearing within cohabitation in Norway while, in a current work, Brée and Doignon (2022) discuss the strong effect of birth control on the fertility in neighbouring districts of Paris during the last years of nineteenth century.
In this theoretical framework, the findings of our study match state of the art methods but, as mentioned above, they go beyond previous results, paying special attention to the relationship between fertility and population density in Italy at a municipal scale.
Future investigations are still necessary to understand more completely the key tenets of this topic. Generally speaking, the relationship between ecological variables and fertility should be identified as a spatial process where it originates partly as a result of the contextual factors and behaviours changing at various geographical scales. In a space studied at different levels of aggregation (global, national and local), the effects of demographic dynamics produce multiple pathways that should not be interpreted as homogeneous. This assumption strengthens the value of the spatial approach and requires further studies, conducted at a microscale level, to capture the variety of demographic behaviours and the relevance of the context in explaining them.
Looking forward, future research could better illustrate the relationships between population density and fertility over time using geographically weighted regressions to analyse distinctive local interactions and also investigate the determinants that affect the negative density-fertility relationship in the case of natives and the positive one in the case of foreigners. Moreover, in the framework of the diffusionist perspective on the fertility transition, a comparative analysis could be conducted on this relationship among southern European countries with the lowest fertility. Of course, the evidence based on the macro model can be enriched by a model based on micro (individual) data. In this perspective, availability of microdata on individuals that reside in specific contexts can play a fundamental role in better understanding the micro foundation of what comes out from the ecological model.