Based on individual-level administrative data, a prospective research design is used in which level of educational attainment and region of residence are measured in 2001, and subsequent fertility hazards for women aged 15–49 are observed in the 2002–2005 follow-up period. ¹ During the follow-up period, women at risk of a birth are observed until they (i) experience a birth; or (ii) are censored due to emigration, death, a 50th birthday or the end of the observation period (i.e. 31 December 2005).

For Belgium, individual microdata from the 2001 census, which cover the complete residential population (Deboosere and Willaert, 2004), are linked to information regarding household composition from the National Register for 2002–2005. This probabilistic linkage based on cross-source household and individual identifiers in combination with registered relations within households has been validated and used in previous research (Wood and Neels, 2017, 2019; Wood et al., 2017, 2020). For Denmark, microdata from Statistics Denmark’s population register for 2001 are linked to Statistics Denmark’s databases for fertility, mortality, educational attainment, student enrolment and migration. For Finland, a nationally representative 11% random sample of the population present in Finland in 1987–2007 who were aged 15 and older in 2007 is drawn from administrative registers, and is linked to different registers by Statistics Finland (license number TK-53-339-13). For France, we use the 2016 database of the Permanent Demographic Sample (EDP), which consists of a large-scale socio-demographic panel (all people born from 1 to 4 October ² (Couet, 2006)) that includes official birth, marriage and death registration along with census information. Although the EDP is currently the most reliable source of data for studying the education-fertility nexus by subnational region from a period perspective, previous research has provided evidence of a potential underestimation of fertility, although mostly among men (Robert-Bobée, 2006). This underestimation has been attributed to limited information on emigration and a lack of information on events occurring outside of the French metropolitan area (Couet, 2006). For the Netherlands, the features of the Belgian dataset are replicated using the System of Social statistical Datasets (SSD), which is the most important source for official social statistics in the country (Bakker et al., 2014). As the SSD consists of a sample in which older registrations of higher education are overrepresented, weights are used to correct for this selectivity, an approach that has been routinely adopted in previous research (e.g. van Gaalen, 2016). For Norway, residence, age and childbirth data are extracted from annual updates of the Central Population Register (Statistics Norway, 2021) and linked to information from the National education database (Statistics Norway, 2017). For Sweden, we link the registered population at the end of 2001 to population registers for the 2002–2005 period.

To capture educational attainment in a comparable way for all countries, we use the International Classification of Education 1997 (ISCED97), which is a standard international classification developed by UNESCO to further the comparability of educational groups across countries (OECD, 1999) (cf. supplementary material, available online at https://doi.org/10.1553/p-4cgn-nz7f, for more info). In order to maximise the comparability of our findings to those of previous research, we followed the common practice of distinguishing between low educational attainment (primary and lower secondary: ISCED codes 0–2), medium educational attainment (upper secondary: ISCED codes 3–4) and high educational attainment (tertiary: ISCED codes 5–6) (e.g. Jalovaara et al., 2019; Wood et al., 2014). In the four northern European countries – Denmark, Finland, Norway and Sweden – the information on educational attainment is based on educational registers. In the Belgian 2001 census as well as the French EDP data, educational degrees are self-reported. In the Netherlands, educational information is based on registers combined with a large set of pooled Labour Force Surveys (LFS), and is made representative via a weighting procedure for the entire population. We use a time-constant indicator for level of education (measured at the start of the observation window), as our dataset does not include time-varying information on level of education.

As the fertility behaviour of individuals enrolled in education has been shown to differ substantially from that of those who have completed their education (Kravdal, 2007; Lappegård and Ronsen, 2005), and this study aims to address gradients by educational attainment, individuals enrolled ³ in education at the start of the observation window are excluded. This approach would be problematic when investigating the causal effects of education on fertility due to reverse causality between fertility and education, as enrolled women who become mothers may drop out of school early ⁴ (Baizan and Martín-García, 2006). However, this report does not focus on the causal relationship between education and fertility, but rather on their association. In addition, excluding individuals who are still enrolled in education also limits problems regarding the usage of a time-constant indicator of level of education. The problem of the measured level of education in 2001 becoming outdated during the observation period would be greater if enrolled subgroups were included in the analyses using their highest level attained at the start of the observation. We also dropped all women for whom no educational attainment information is available. The resulting sample sizes by educational attainment, region and parity are provided in the supplementary material (cf. Tables S1–S2).

For the measurement of the region of residence, we use NUTS level 2 (Nomenclature of Territorial Units for Statistics) categories. These categories generally yield more regional units for selected countries with higher population numbers, and fewer regional units for countries with lower population numbers. The only exception is the use of NUTS level 1 regional units for France, as the French sample size does not allow us to use the NUTS level 2. In all countries, the applied NUTS categories capture strong differences in the degree of urbanisation (measured by population density (inhabitants per square kilometre)), as illustrated by Table 1. For Belgium, NUTS 2 regional units are used that distinguish the extremely densely populated capital, Brussels, from five Dutch-speaking northern provinces with higher population densities (Antwerp, Flemish Brabant, East Flanders and Limburg), and five French- or German-speaking southern provinces with lower population densities (Hainaut, Walloon Brabant, Liège, Namur and Luxembourg). It should be noted that the Brussels capital region is an outlier in our data, as Belgium includes some of the most densely populated areas in Europe, but also because the Brussels region includes only the capital city, in contrast to some other regions that include large and/or capital cities in addition to surrounding areas. For Denmark, the NUTS 2 categories distinguish five regions: the relatively densely populated capital region of Hovedstaden and four other regions with lower population densities (Sjælland, Southern Denmark, Midtjylland and Nordjylland). For France, we adopt NUTS 1 categories according to the classification until 2016 that distinguish the densely populated capital region Région parisienne from the relatively densely populated area Nord and the increasingly rural regions of Est, Méditerranée, Centre-Est, Ouest, Bassin parisien and Sud-Ouest. For Sweden, NUTS 2 level categories distinguish the Stockholm region from seven increasingly rural regions: South Sweden, West Sweden, East Middle Sweden, Småland and islands, North Middle Sweden, Middle Norrland and Upper Norrland. For the Netherlands, NUTS 2 units distinguish 12 regions referring to the 12 provinces that vary considerably in population density, ranging from the densely populated Zuid Holland, Noord Holland and Utrecht to the progressively rural provinces of Limburg, Noord Brabant, Gelderland, Overijssel, Groningen, Flevoland, Zeeland, Friesland and Drenthe. For Finland, NUTS 2 categories distinguish four large regions: the relatively densely populated capital region Helsinki-Uusimaa and the three rural regions of South Finland, West Finland and East and North Finland. Finally, for Norway, the division into NUTS 2 regions is expanded with an extra category for Oslo, creating eight regions. The strongly urbanised Oslo region is distinguished from seven more rural regions: the Akershus region surrounding the capital, Agder og Rogaland, Sør-Østlandet, Vestlandet, Trøndelag, Hedmark og Oppland and Nord-Norge in the north.

Finally, the register data also include an indicator for migration background, more specifically a variable indicating whether a person was born in their country of residence.

To assess the relationship between female educational attainment in 2001 and subsequent birth hazards in the 2002–2005 follow-up period, discrete-time hazard models with a logit link function are used. We adopt late-entry models as women enter the observation window in 2002 at different exposures, and thus contribute to different partitions of the hazard function.

As a result of user regulations related to the register data for the six countries studied (e.g. stipulating that the aggregate tables can be shared only if the cell counts are higher than a given threshold), the aggregation of person-period files to tables that are collapsed over all indicators used in the models allows us to estimate models using these tables. These tables include occurrence (i.e. events) as well as initial risk sets (with cell counts as weights), which allow us to attain estimates identical to those of individual-level discrete-time hazard models based on individual data. This approach is similar to the usage of event-occurrence tables in previous research (Kulu et al., 2017). Models are estimated separately for first, second and third births.

Model 1 includes (I) exposure, (II) region and (III) educational attainment. For first births, exposure distinguishes age categories 15–21 (reference category), 22–28, 29–35 and 35–45; whereas for second and third births, the indicator distinguishes duration categories 0–3 (reference category), 4–9 and 9–30 years after the previous birth. Hence, this variable reproduces relatively low first birth hazards among the young and old age groups and particularly high second and third birth hazards shortly after the previous birth. Region distinguishes the 56 regions under consideration, and is included to control for differences in fertility levels between all regions within and across countries. As birth hazard functions are likely to depend on the context considered, the effect of exposure is allowed to vary by region. Education captures differential birth hazards between women with medium or high education and women with low education (reference category). In line with the widely supported finding that the shape of first and higher order birth functions varies by level of education (e.g. Kreyenfeld, 2002; Wood et al., 2014), and as a result of significant improvements to null models without such interactions models (i.e. log likelihood ratio tests), model 1 interacts the effects of exposure with level of education.

In order to assess regional variation in the female education-fertility nexus, model 1 also includes an interaction between level of education and region, allowing the educational gradient in first, second and third births to vary across 56 subnational regions in seven countries. ⁵

The results of these models will be illustrated using average marginal effects of medium and high levels of education, with corresponding significance tests, by exposure. The advantage of using average marginal effects – representing the additive difference in discrete-time hazards between educational groups – is that these measures provide information on educational differentiation while also reflecting the general level of the birth hazard in the part of the hazard function considered. As such, educational differences for less relevant parts of the hazard function with very low birth hazards (e.g. second birth hazards after more than 10 years after the first birth) are less likely to be (mis)interpreted as important in terms of their contribution to educational differentiation in fertility. In addition, as the use of odds ratios can be considered complementary to the use of average marginal effects by providing an indicator of differentiation regardless of the general level of the birth hazards, results in terms of odds ratios are provided in the online supplementary material (cf. Tables S3–S4).

Following previous research (Neels, 2006; Wood et al., 2014), we use estimated birth hazards by level of education, region and exposure to calculate synthetic parity progression ratios (SPPRs) using equation 1. In general, parity progression ratios (PPR) can be interpreted as the proportion of women with a certain number of children who go on to have another child. The SPPR for birth order i and time t reflects the proportion of women in the period fertility table having a birth of order i within 30 years after entry into the risk set (i.e. age 15 for first births and year of the previous birth for second and higher order births). SPPRs are based on a period, or “synthetic”, life table. Consequently, SPPRs indicate the share of women who would progress to the next parity if they were exposed to the age- or exposure-specific birth hazards estimated for the 2002–2005 follow-up period (Hinde, 1998). Hence, in contrast to cohort parity progression ratios (CPPR), they do not represent the proportion of women who progress to the next birth within an actual birth cohort. However, it is well-documented that SPPRs approximate cohort fertility indictors considerably more strongly than period total fertility rates (e.g. Neels, 2006). SPPR i t = 1 − ∏ x = 0 30 ( 1 − q ( i , x , t ) ) (1) This report focuses on two educational contrasts in estimated SPPR i : (I) between medium and low educated groups, and (II) between highly educated and low educated groups. As the magnitude of the absolute difference in SPPR 1 – 3 between different educational groups is, by definition, influenced by the overall levels of parity progression in a given region, the supplementary material also includes differential SPPR for medium and highly educated groups as a percent of the SPPR among the low educated group (cf. Tables S5–S6).

In a final step of the analyses, SPPR 1 – 3 differentials between educational groups are correlated with regional-level population density in order to not only explore regional variation as such, but also to provide statistical tests for the association between urbanisation and the educational gradient in fertility. This is done using the estimates of the models as discussed in section Hazard models (i.e. model 1), but also using alternative model specifications (models 2–5) to address the robustness of this association. Whereas the first model includes all the aforementioned indicators and interactions, model 2 additionally includes country fixed effects by level of education, implying that only within-country cross-regional variation is used in the association with the degree of urbanisation. Subsequently, model 3 additionally excludes the Brussels capital region, which is a strong outlier in terms of urbanisation. Model 4 adopts the same model specifications as model 2, yet excludes all foreign-born women from the analyses to address whether the association between educational gradients and the degree of urbanisation changes when taking foreign-born groups out of the equation. Finally, model 5 combines the aforementioned specifications by including country fixed effects by level of education, excluding the Brussels capital region and excluding all foreign-born women.

In line with well-accepted findings regarding differential first birth schedules, Figure 1 shows clear indications of postponement and recuperation, with first birth hazards for medium and highly educated women being lower in the younger age categories, and higher in the higher age categories. However, large variation in these patterns can be observed between countries around the turn of the century. Northern European countries stand out, with more positive associations between medium or high levels of education and first birth hazards in the 15–21 and 22–28 age categories, and stronger positive associations between high levels of education and first birth hazards in the 29–35 age category. Results indicate that the educational gradients in first births vary significantly between regions within countries. Model comparison tests indicate that allowing varying educational gradients by subnational region in model 1 significantly improves model fit compared to a restricted model allowing varying gradients by country only ( Δ -2LL: 1014.95; Δ df. 98; p<0.000).

Figure 2 illustrates how the aforementioned variation in first birth schedules yields considerable regional variation in the educational gradient in SPPR 1 , and allows us to examine regional and urban-rural differences in the educational gradient in the transition to a first birth. While our results generally show a positive relationship between level of education and SPPR 1 , in most countries, more densely populated areas exhibit stronger negative or weaker positive differential SPPR 1 for medium and highly educated women, as illustrated by strongly negative or weakly positive gradients in capital regions. Regional-level correlations between educational gradients in first births, as estimated by model 1, and population density (Table 2) are −0.40 for both medium and highly educated women. When controlling for between-country differences in model 2, these correlations weaken, with corresponding values of −0.25 (p<0.01) and −0.27 (p<0.001) for medium and highly educated women, respectively. However, when excluding the extremely densely populated Brussels region (model 3), which is an outlier (cf. Table 1), only the negative correlation between population density and mid-low educational differentials remains significant. The negative correlation between educational gradients and population density disappears when excluding foreign-born women (models 4 and 5), which illustrates the impact of population composition in terms of country of birth on urban-rural variation in the educational gradient in fertility.

Furthermore, a visual inspection of the results in Figure 2 at least suggests that our conclusion that there was significant regional variation during the early 2000s, with stronger negative and weaker positive educational gradients in more densely populated regions, might not hold for all countries independently. Results for Belgium and the Netherlands indicate considerable variation in both the direction and the magnitude of the educational gradient in first births, yet the only clear patterning by population density can be observed in the relatively strong negative gradients of the Brussels and Zuid-Holland capital regions. Findings for Denmark, Norway and Sweden show limited variation in the positive educational gradient in SPPR 1 , yet clearly weaker positive gradients in more densely populated regions. The variation in the positive educational gradient in Finland does not seem to go hand-in-hand with population density, and France is characterised by limited variation in a neutral educational gradient in SPPR 1 .

Figure 3 indicates that in all regions considered, second birth hazards (predominantly in the 1–3 years after the birth of the first child) varied significantly by women’s level of educational attainment in the early 2000s, with mostly positive educational gradients. This result is consistent with the routinely documented finding that higher educated women space their second birth closer to the first (Kreyenfeld, 2002). These effects are followed by a mixture of positive and negative associations between medium or high education and second birth hazards for the 4–9 exposure category, and less relevant educational differentials during later exposure periods with extremely low second birth hazards.

Figure 4 demonstrates considerable regional variation in the resulting SPPR 2 levels within countries. Model comparison tests indicate that allowing the educational gradient in second births to vary by subnational regions in model 1 in addition to between-country variation significantly contributes to the model fit ( Δ -2LL: 952.54; Δ df. 98; p<0.000).

In Belgium, Denmark, Finland and Sweden, the capital regions exhibit the weakest positive or even strongest negative differential SPPR 2 for medium educated women, highly educated women or both. Regional-level correlations between educational gradients in SPPR 2 and population density (Table 2) are −0.31 and −0.46 for medium and highly educated women when controlling for between-country differences (model 2). Similar to the findings for SPPR 1 , this correlation weakens for medium educated differential SPPR 2 and disappears for highly educated women after excluding the Brussels outlier (model 3).

However, we find a positive correlation between population density and the educational gradient in SPPR 2 when excluding foreign-born women from the analysis, in addition to excluding the Brussels outlier (model 5). This implies that educational gradients in the transition to a second birth – which are generally positive – are actually more positive (and less likely to be negative) when only native-born groups are considered. This finding suggests that the differential composition in terms of foreign-born populations in less urbanised and more urbanised regions was again responsible for the less positive gradients in more urbanised regions shortly after the turn of the century.

With respect to cross-country differences in regional variation and urban-rural patterning, it should be noted that the degree of variation in the educational gradient in second births seems to vary considerably by country. In contrast to the strong variation in the direction and the magnitude of educational gradients in France, and the relatively large variation in the positive educational gradient in Belgium and Finland, Denmark, the Netherlands, Norway and Sweden exhibit a more consistent positive educational gradient in the progression to a second birth. Thus, largely in line with the findings for first births, northern European countries exhibit less regional variation in educational gradients compared to the selected western European countries.

Figure 5 shows that medium educational attainment is associated with lower third birth hazards compared to those for low educational attainment, with the exception of positive associations in some Norwegian, French and Dutch regions. Regarding the relationship between high educational attainment and third birth hazards, positive effects compared to those for lower educational attainment (in the 1–3 years following the second birth) are found not only in Norwegian and Dutch regions, but also in many Danish regions, a considerable number of Belgian regions and a few French regions. Differential third birth hazards by educational attainment exhibited considerable subnational geographical variation in the early 2000s. Model comparison tests indicate that allowing the educational gradient to vary by subnational regions in model 1 significantly contributes to the model fit ( Δ -2LL: 564.61; Δ df. 98; p<0.000).

The resulting SPPR 3 levels illustrated in Figure 6 indicate a consistently negative gradient of varying magnitude in Sweden, Denmark and Finland, and greater variation in the shape and magnitude of educational differentials in Belgium, France, the Netherlands and Norway. Such large variation in educational gradients in third births compared to those for lower birth orders and the occurrence of (reverse) U-shaped educational gradients in third births have also been identified in previous literature (Neels and De Wachter, 2010; Wood et al., 2014).

With respect to urban-rural patterns in the educational gradient in SPPR 3 , a visual inspection of Figure 6 shows signs of more negative gradients in more densely populated regions in Belgium, Denmark and Norway. A more formal test of the urban-rural patterning of educational gradients in SPPR 3 is provided through regional-level correlations between educational gradients in third births and population density (Table 2). We find relatively strong negative correlations between regional-level population density and the educational gradient in SPPR 3 , which strengthen when controlling for between-country differences (model 2) and after excluding the Brussels outlier (model 3). However, when also excluding foreign-born groups (model 5), such negative correlations are considerably weaker and no longer reach statistical significance. Consequently, these findings indicate that, in addition to the generally negative gradients, the disproportionate presence of foreign-born groups in more densely populated regions drives the stronger negative gradients in SPPR 3 in these regions. When controlling for this compositional difference by excluding foreign-born women, educational gradients in third births are not more positive in more urbanised regions. This implies that foreign-born women play an important compositional role in the study of educational gradients in fertility in urban contexts, given that foreign-born groups have relatively low educational attainment and higher fertility (upon arrival).

Finally, similar to our findings for first and second births, it is noteworthy that the degree of within-country regional variation in the educational gradient in third births seems to depend on the country and educational contrast considered. For instance, particularly large regional variation in the direction and magnitude of associations in some country-specific educational differences (e.g. high versus low education in Belgium) contrasts with very limited variation in others (e.g. medium versus low education in Sweden).