# Modelling the age pattern of fertility: an individual-level approach

Daniel Ciganda, Nicolas Todd

<figure><img src="image_1.png" /></figure>

## To cite this version:

Daniel Ciganda, Nicolas Todd. Modelling the age pattern of fertility: an individual-level approach. Royal Society Open Science, 2024, 11, 10.1098/rsos.240366. http://hal-04866069

HAL Id: hal-04866069

https://hal.science/hal-04866069v1

Submitted on 6 Jan 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>
<figure><img src="image_3.png" /></figure>

Research

<figure><img src="image_4.png" /></figure>

<figure><img src="image_5.png" /></figure>

Check for
updates

**Cite this article:** Ciganda D, Todd N. 2024
Modelling the age pattern of fertility: an
individual-level approach. *R. Soc. Open Sci.* **11**: 240366.

https://doi.org/10.1098/rsos.240366

Received: 4 March 2024

Accepted: 17 September 2024

**Subject Category:**

Science, society and policy

**Subject Areas:**

computer modelling and simulation, statistics,
behaviour

**Keywords:**

agent-based modelling, fertility, reproductive
process, demography, approximate Bayesian
computation

**Author for correspondence:**

Daniel Ciganda

e-mail: ciganda@demogr.mpg.de

Electronic supplementary material is available
online at https://doi.org/10.6084/
m9 figshare.c.7477930.

# Modelling the age pattern of fertility: an individual-level approach

Daniel Ciganda<fnref n="1" />,<fnref n="2" /> and Nicolas Todd<fnref n="3" />

DC, 0000-0002-6181-9947

Fitting statistical models to aggregate data is still the dominant approach in many demographic and biodemographic applications. Although these macro-level models have proven useful for a variety of tasks, they often have no demographic interpretation. Individual-level modelling, on the other hand, offers a deeper understanding of the mechanisms underlying observed patterns. Their parameters represent quantities in the real world, instead of pure mathematical abstractions. However, estimating these parameters using real-world data has remained a challenge. The approach we introduce in this article attempts to overcome this limitation. Using a likelihood-free inference technique, we show that it is possible to estimate the parameters of a simple but demographically interpretable individual-level model of the reproductive process by exclusively relying on the information contained in a set of age-specific fertility rates. By estimating individual-level models from widely available aggregate data, this approach can contribute to a better understanding of reproductive behaviour and its driving mechanisms, bridging the gap between individual-level and population-level processes. We illustrate our approach using data from three natural fertility populations.

## 1. Introduction

Modelling the age distribution of fertility rates is an essential step in a number of demographic applications. When only a tight fit to an observed schedule is required, as in the generation of single-age rates from grouped data, nonparametric models, typically based on splines, tend to produce the best results [1]. No particular price is then placed on whether each model parameter may be interpreted in any meaningful way. Other

<footnote n="1">¹Max Planck Institute for Demographic Research, Rostock, Germany</footnote>

<footnote n="2">²Instituto de Estadistica, UDELAR, Montevideo, Uruguay</footnote>

<footnote n="3">³UMR7206 ‘Eco Anthropologie’, Musée de l’Homme, CNRS, Paris, France</footnote>
applications require models that can fit the data well while also providing a well-defined, ideally small set of parameters that represent quantities that can be interpreted in demographic terms [2–6].

The interpretability of parameters is particularly relevant in a forecasting context, where it may prove especially useful to associate a change in the value of the parameters with an underlying behavioural process, such as fertility postponement or the diffusion of contraceptive methods.

A difficulty associated with most existing parametric models is that the interpretation of their parameters can be elusive [7]. Even in the best cases, the relationship between mechanisms and parameters is ambiguous and indirect. In essence, macro-level models always operate one level above the processes that drive fertility.

A natural solution to this problem is to model at the individual level. This approach may be decomposed into three steps: (i) develop a micro model of the reproductive process, (ii) use this model to generate synthetic data and compute a set of simulated age-specific fertility rates (ASFRs) and (iii) estimate the parameters of the model by minimizing some measure of distance between simulated and observed rates.

The framework required for step (i) has long been developed by some of the pioneers of mathematical demography [8]. The ideas that emerged from this initial stage were later expanded to more realistic settings using simulation methods. The objective of this second wave of models was the incorporation of basic characteristics of the reproductive process such as age-dependent fecundability, or to represent change over time in fertility rates [9,10]. The use of simulation methods allowed researchers to complete step (ii), by generating synthetic ASFRs from individual-level models, drawing a direct connection between individual-level behaviour and observed, aggregate, quantities.

These contributions, however, remained largely theoretical. One key issue was that parameter values were borrowed from previous studies or calibrated through trial and error until the simulated data would appear to fit a given target distribution. In other words, estimation procedures remained too fuzzy, if not arbitrary, to decide whether individual-level models could be of use to analyse actual data.

Interest in these type of models dwindled after the 1980s. Ironically, the first statistical methods that would enable inference on complex simulation models started to emerge at around this time [11]. Likelihood-free inference methods, such as approximate Bayesian computation (ABC), have been extensively researched since then, providing a robust statistical framework to estimate the parameters of computational models [12]. In short, likelihood-free methods have reconciled statistical and mechanistic modelling by allowing the rigorous statistical estimation of sophisticated mechanistic models.

What we attempt to show in the remainder of this article is that these advances in statistical computing, that have proved so useful in fields such as ecology and population genetics, also provide the missing piece to the programme outlined above, finally completing step (iii). Specifically, these advances make it possible to infer individual-level quantities even when only aggregate fertility rates are available. More generally, the approach we advocate here provides a clear roadmap to fully integrate behavioural mechanisms in the modelling and forecasting of fertility trends.

To illustrate our approach, we use high-quality data from three natural fertility populations: Hutterite communities from 1860 to 1914, the French-Canadian population from 1700 to 1750 and a subsample of the Louis Henry Enquête in eighteenth-century France. The three datasets provide full maternity histories as well as information on the dates of other relevant events in a woman's life course like marriage, death or death of spouse. However, we use the available individual-level data only for the validation of our results (see §4.3) and not for the estimation of our model because our objective is to show that individual-level quantities can be recovered when only aggregate data are available.

All the results presented in this article are fully reproducible using the code and data available at https://github.com/dciganda/individual_fx.

## 2. Model

The simplest form of the reproductive process is one in which the outcomes of the process (exact timing and number of births) are not a function of individual preferences. Louis Henry coined the term 'natural fertility' to refer to this type of process. A model of the reproductive process in a natural fertility setting requires three basic inputs: the moment when the process starts, typically marriage; the probability of a conception, known as fecundability [13]; and the length of the period in
<figure><img src="image_6.png" /><figcaption>**Figure 1. Fecundability model.** Although it depends on only two parameters, the model used to represent the evolution of the risk of conception across the reproductive life course can capture a wide range of patterns. The dashed lines represent the two Bernstein polynomials. The solid red lines represent different fecundability profiles associated with different combinations of $\phi_1$ and $\phi_2$ values.</figcaption></figure>

which a woman is not able to conceive following a childbirth, colloquially identified to post-partum amenorrhea.

We follow this approach and model the reproductive experience of a cohort of women to age 50 using a discrete-time simulation model, in which time advances at fixed increments of 1 month. The age (in months) at marriage for a woman *i* is simulated from a lognormal distribution, parametrized via its mean $\mu_m$ and standard deviation $\sigma_m$.

Every month, married women who are neither pregnant nor amenorrheic are exposed to a probability of conceiving (her fecundability), which depends on her age. Whether susceptible woman *i* conceives in a given month when aged $x^*$ is the outcome of a Bernoulli trial with probability $\phi(x^*)$.

It is well established that fecundability progressively increases from puberty, reaches a peak during the early adult years and decreases progressively afterwards until menopause [14–16]. In order to capture this pattern, we model fecundability as a polynomial of degree 3. We define $\phi(x)$ as the fecundability at age $x$ in the reproductive window (assumed 10–50 years), and $x_s = (x - 10)/(50 - 10)$, the age scaled to belong to [0, 1]. Instead of using the standard basis $(1, x_s, x_s^2, x_s^3)$, we use the Bernstein polynomials $B_1: x_s \mapsto 3x_s(1 - x_s)^2$; $B_2: x_s \mapsto 3x_s^2(1 - x_s)$; $B_3: x_s \mapsto (1 - x_s)^3$; $B_4: x_s \mapsto x_s^3$, decomposing $\phi(x)$ as $\sum_{j=1}^4 \phi_j B_j(x_s)$. This enables us to naturally restrict our search to a subspace of dimension 2: since at age 10 ($x_s = 0$): $\phi(10) = \phi_3 B_3(0) = \phi_3$ and at age 50 ($x_s = 1$): $\phi(50) = \phi_4 B_4(1) = \phi_4$, the fact that fecundability is known to be 0 at ages 10 and 50 immediately sets $\phi_3$ and $\phi_4$ to 0. What remains is therefore:

$$ \phi(x) = \phi_1 \cdot B_1(x_s) + \phi_2 \cdot B_2(x_s). \qquad (2.1) $$

Different combinations of $\phi_1$ and $\phi_2$ generate different age profiles for the risk of conception, as seen in figure 1. Given the age profile we expect, $\phi_1$ will have a stronger influence on the level at which fecundability peaks and $\phi_2$ will largely control the pace at which the risk of conception decays, until permanent sterility is reached. While the specific values of $\phi_1$ and $\phi_2$ are not informative on their own, the resulting age profile of fecundability produced by a given combination reflects both the behavioural (frequency of sexual intercourse) and biological mechanisms that influence the risk of conception.

A conception brings about a state of non-susceptibility (no risk of conception) equal to the duration of the pregnancy, 9 months, plus the duration of post-partum amenorrhea, $\delta$, that must be estimated together with the other parameters of the model.

The model outputs individual reproductive trajectories, from which a set of simulated ASFRs $s_f$ can be computed. We then leverage the distance between simulated and observed rates, $\delta_f$, to estimate the parameters of our model, namely ($\mu_m, \sigma_m, \phi_1, \phi_2, \delta$). This procedure is explained in more detail in §3.5, after we describe the data used in our analysis in the next section.
# 3. Data

## 3.1. Hutterites

The Hutterites are an Anabaptist community, originating in the fourteenth century in the Tyrolean Alps. After a long history of migration in Europe and Asia, they relocated to North America at the end of the nineteenth century where they still live today in self-sustained, largely isolated colonies.

Like other religious communities, the Hutterites oppose the use of birth control methods [17,18]. What made them stand out for demographic analysis was how strictly they adhered to these beliefs, at least until the second half of the twentieth century. The average marital fertility of the Hutterite cohorts born until the early 1900s remained around 10 children per woman, providing researchers with an exceptional opportunity to study the reproductive process under natural fertility conditions [19].

Another reason the Hutterites became the gold standard in natural fertility research was their practice of keeping detailed family records. These records, personally checked for consistency by colony preachers, were made available for various scientific studies in the 1950s and 1960s [19,20]. In the course of one of these studies, the original data were extended through follow-up interviews, which resulted in complete maternity histories for 562 families [21].

## 3.2. Historic Quebec

BALSAC is a longitudinal population database that contains information on individuals and families who have lived in Quebec from the seventeenth century to the present. The information used in this article, that concerns cohorts born before 1750, was gathered using family reconstitution methods from parish registers [22]. Given the high quality of these records, the information on seventeenth-century Quebec populations has become an important reference in the literature on natural fertility [15,23].

## 3.3. Historic France

Another important dataset for natural fertility research has been the Enquête Louis Henry. This dataset was created by Louis Henry and his co-workers in an effort to draw a nationally representative sample of the eighteenth-century French population. It contains data on life events for a sample of 378 parishes for the period 1670-1829 from parish and civil registers—the ‘anonymous sample’. Family reconstitutions were carried out for a subsample of 40 rural parishes—the ‘nominative sample’. For a detailed description of the Enquête, see [24]. We use the ‘nominative sample’.

A well-known problem in historical demography is that a small fraction of births may be missing from the parish and civil registers used. Specific procedures were developed by Henry to systematically assess and correct for these unregistered births [25]. Given our objectives here, we make no attempt to compute fertility rates corrected for these unregistered births. However, doing so would not pose any additional difficulties.

## 3.4. Sample sizes

To keep things simple, we model a process without death or union dissolution. Therefore, we restrict our data to ‘intact’ marriages, i.e. marriages that were neither lost to follow-up nor interrupted by the death of one of the spouses before the woman reached age 50. Although our model simulates the reproductive experience of a single birth cohort of women, we include information from multiple cohorts in the computation of our observed ASFRs in order to obtain reasonable sample sizes. To keep our hypotheses realistic, we analyse a range of birth cohorts that are as homogeneous as possible with respect to age at first marriage and total fertility rate (table 1).

Figure 2 shows the observed ASFRs for our three natural fertility populations. As expected, the Hutterites show consistently higher fertility from age 20, with an average number of children per woman in these cohorts of 10.8. The French-Canadian cohorts have the second largest fertility in the three populations, with an average of 7.7 children per woman, followed by the French cohorts with an average of 5.7 children per woman. The distribution of ASFRs for the French cohorts also peaks at older ages, which is not surprising given their higher mean age at marriage, 25.8 years, compared with 20.6 years for Hutterites and 22.1 years for French-Canadian cohorts.
<figure><img src="image_7.png" /><figcaption>**Figure 2.** Age-specific fertility rates for three natural fertility populations. Dashed lines show the fit of a cubic smoothing spline.</figcaption></figure>

<table><caption>**Table 1.** Sample information for the three natural fertility datasets.</caption><thead><tr><th></th><th>no. marriages</th><th>no. births</th><th>cohorts</th></tr></thead><tbody><tr><td>Hutterites</td><td>161</td><td>1726</td><td>1860–1914</td></tr><tr><td>Historic Quebec</td><td>14 303</td><td>110 772</td><td>1700–1750</td></tr><tr><td>Historic France</td><td>3235</td><td>18 623</td><td>1680–1760</td></tr></tbody></table>

## 3.5. Estimation

As previously mentioned, our primary goal is to demonstrate that the lack of individual-level data is no impediment to understanding the behaviours and mechanisms underlying fertility processes.

When only aggregate data are available, the probability of individual observations cannot be calculated, and model parameters cannot be estimated using a maximum likelihood approach. Instead, we rely on ABC, a popular likelihood-free approach for the estimation of computational models.

The objective of ABC is to get a sample drawn from an approximation of the posterior distribution, without ever explicitly calculating the likelihood function. Instead, ABC algorithms simulate data from the model for a set of parameter values drawn from a prior distribution, compute summary statistics from each of these simulated datasets and keep only those parameter values that yield summary statistics that are close to those computed on the observed data, according to some predefined distance function. We naturally use the vector of ASFRs as our summary statistics and the Euclidean distance as our distance function.

The basic ABC rejection algorithm we use is straightforward to implement, but it requires a large number of simulations to obtain a sufficiently large sample from the posterior distribution [26,27]. One must accept keeping parameter values yielding simulated statistics only reasonably close to the observed statistics. To correct for this approximation and thus improve the efficiency of the method, one natural approach is to model the parameters as a function of the summary statistics in order to correct the parameter values that were accepted in the first step. This procedure, known in the ABC literature as 'regression adjustment', leads to a better approximation of the posterior distribution [28]. For this second step, we use a random forest (RF), as described in the pseudocode below.
**Algorithm: ABC + Regression Adjustment**

**Input:** Prior distribution $\pi$ for the $J$ parameters, $\theta = (\theta_1, ..., \theta_J)$, model $\mathcal{M}(\theta)$ generating simulations, observed age-specific fertility rates $o_f$, distance threshold $\epsilon$.

**Output:** Sample from the posterior distribution $\{\tilde{\theta}_1, \tilde{\theta}_2, ..., \tilde{\theta}_N\}$.

begin

Draw samples from $\pi$

For each sampled $\theta$, simulate data $y$ from $\mathcal{M}(\theta)$

Compute simulated fertility rates $s_f$ from $y$

Compute Euclidean distance between $s_f$ and $o_f$

Accept $\theta$ if the distance is less than $\epsilon$

With accepted parameter values $\{\theta_1^*, \theta_2^*, ..., \theta_N^*\}$ and corresponding simulated rates $\{s_f1, s_f2, ..., s_fN\}$:

begin

For each dimension $j\theta$ of $\theta$, train a Random Forest $m_j$ that maps $s_f$ to $j\theta$

Adjust accepted parameters values:

$$j\tilde{\theta}_i = m_j(o_f) + (j\theta_i^* - m_j(s_fi))$$

end

Return $\{\tilde{\theta}_1, \tilde{\theta}_2, ..., \tilde{\theta}_N\}$ as a sample from the adjusted posterior distribution.

end

## 3.6. Sample size

As shown in table 1, the number of individual trajectories used to compute ASFRs varies widely between the three populations. Smaller sample sizes naturally produce noisier distributions of observed rates (as seen in figure 2) and therefore provide less information on the data generating process: having fewer data must result in a larger degree of posterior uncertainty about the parameters. To take this source of uncertainty into account, we match the number of marriages in our simulations to the number of marriages observed in each population, with the only exception of the French-Canadian sample. For sample sizes above 5000, the sample-to-sample variation in estimated fertility rates becomes extremely small: all vectors of rates simulated with 5000 marriages or more are similarly smooth. Since the computational cost of simulating the actual number of marriages in the French-Canadian sample (14 303) is prohibitive, in this case we simulate cohorts of 5000 marriages.

## 3.7. Prior choice

We use uniform priors for all our parameters and define a single prior distribution for the three populations analysed. In order to define lower and upper limits for these distributions, we reviewed the literature in a search for available estimates for our quantities of interest, as explained below. In all cases, this information was used to define a wide, plausible range of values for our search.

A compilation of estimates for peak effective fecundability (the maximum of function $\phi$) in natural fertility populations can be found in [29] and [30]. Taken together, these estimates range from about 0.15 to 0.30. We set the bounds of $\phi_1$ and $\phi_2$ to comfortably include these numbers, leading to peak fecundability values in the 0.12–0.35 range.

The duration of post-partum amenorrhea cannot be computed directly from maternity history information alone. However, a rough estimate of post-partum amenorrhea duration after the first birth can be obtained by subtracting the mean duration between marriage and first births from the mean duration between first and second births. Given that the duration of non-susceptibility periods increases with age and parity, this procedure offers a good lower bound estimate for mean amenorrhea duration across the reproductive life course. Both [21] and [31] estimate this lower bound at 6 months in the case of the Hutterites. Leridon [29] provides estimates of 7.2, 8.6 and 10.2 months for three eighteenth-century
<figure><img src="image_8.png" /><figcaption>**Figure 3.** Results of model cross-validation. The location of the red dots displayed in the figure is defined by the mean of the estimated posterior distribution in the y-axis, and the value of the parameter that generated the pseudo-observed data in the x-axis. (a) $\phi_1$, (b) $\phi_2$, (c) $\delta$, (d) $\mu_m$, (e) $\sigma_m$.</figcaption></figure>

French samples, which are not far from the estimate of 7.5 months we obtain by applying the method described above to the French sample (6.6 months for the French-Canadian sample). Taking these figures into account, we set the lower limit for the $\delta$ prior at 6 months and the upper limit at 18 months, which appears to be an upper bound for populations practicing extended breastfeeding [32].

Finally, we define the prior range for $\mu_m$ at ages 19–27 years, and for $\sigma_m$ at 2–7 years, both comfortably including the empirical means and standard deviations of the age at marriage distributions in our samples.

# 4. Results

## 4.1. Model validation

We use cross-validation to assess whether the parameters in our model can be correctly identified from the information contained in a set of ASFRs. For this, we randomly select the ith simulation as our validation set. The associated simulated rates are taken as pseudo-observed data. With all other simulations, the model parameters are then estimated as previously described, yielding a mean posterior $j\tilde{\theta}_i$ for parameter j. Repeating this procedure for 100 randomly chosen simulations, we compute the prediction error for parameter j as

$$ERR_j = \frac{\sum_i (j\tilde{\theta}_i - j\theta_i)^2}{\mathrm{Var}(j\theta_i)}.$$

The results of the cross-validation exercise are shown in figure 3. The three sets of results correspond to the different sample sizes used in the simulations. The results obtained with the largest sample size show that the estimation procedure consistently returns parameter values that are very close to the true
<figure><img src="image_9.png" /><figcaption>**Figure 4.** Model fit, 95% credible intervals. (a) Hutterites, (b) eighteenth-century Quebec, (c) seventeenth–eighteenth-century France.</figcaption></figure>

parameters that generated the data. As expected, simulation results obtained with a smaller number of marriages are associated with large prediction errors.

Within a given sample size, the degree of uncertainty with which a parameter is estimated depends on to what extent the available data allow us to assign each region in the parameter space a distinct probability of having generated the data. When the effect of a change in a given parameter has an effect on the outcome that is similar to the effect associated with another subset of parameters, different regions of the parameter space will have similar probabilities of having generated the data, leading to wider posteriors.

Given this, it is not unexpected that the parameter with the smallest associated prediction error is $\mu_m$ and the parameter with the largest associated error is $\delta$. Indeed, a change in the value of $\mu_m$ produces a shift of the entire distribution of rates that is very distinct from the effect of a value change in any other parameter. A change in the value of $\delta$, on the other hand, produces an effect on the rates that overlaps to a certain extent with the effect of a change in the value of $\phi_1$ (see the electronic supplementary material for some animations illustrating the effect of each parameter on the schedule of ASFRs).

## 4.2. Model fit

Figure 4 depicts the model’s posterior predictive distribution against observed ASFRs for the three populations analysed. The model accurately captures the general characteristics of the data, such as the shape, location of peaks and rate of decline of the three distributions. In fact, the fit appears very satisfactory at most ages, except for the rates at the very end of the reproductive age window (ages 47–50), which the model tends to slightly underestimate. Credible intervals reflect the higher degree of uncertainty with which parameters are estimated in the Hutterites’ case.
<figure><img src="image_10.png" /><figcaption>**Figure 5. Estimated fecundability and post-partum amenorrhea durations.** The distributions of the probability of conception by age (a) were obtained by sampling a set of values for $\phi_1$ and $\phi_2$ from the posterior distribution (see appendix A) and applying equation (2.1). The lines represent the mean values of the resulting distributions and the shaded areas the quintiles that contain 95% of the distribution. (b) The posterior distributions for the length of post-partum amenorrhea in our three populations.</figcaption></figure>

### 4.3. Individual-level mechanisms

Figure 5 shows the estimated age pattern of fecundability (figure 5a) and the estimated duration of post-partum amenorrhea for our three populations (figure 5b), while figure 6 displays the estimated versus observed age at marriage distributions. Taken together, they indicate that the higher ASFRs of the Hutterites are the result of a combination of higher fecundability, earlier marriage and shorter non-susceptibility periods. They also show that the model correctly explains timing differences in fertility in these populations not as a consequence of differences in the age pattern of fecundability, which is very similar in the three cases, but as a result of differences in marriage patterns.

Although the model captures well the highest density region of the age at marriage distribution, when the right tail of the distribution is more pronounced, as in the case of our French-Canadian sample, it tends to slightly underestimate its variance. The difficulties in correctly representing the right tail of these distributions are most likely related to the fact that we do not explicitly model any parity or marriage duration effects. By ignoring these effects, we are effectively computing the average fecundability and average length of post-partum amenorrhea over parity and marriage duration, making it harder for the estimation process to identify the exact proportion of marriages at each age. Nevertheless, the misalignment is relatively small considering the simplicity of the model and the data used for its estimation.

The peak fecundability that is implied by our estimates for $\phi_1$ and $\phi_2$ for the French and French-Canadian populations clearly falls on the lower end of the available estimates range (see §3.7). This is due to the fact that our model does not account for differences in couples' conception probabilities. In the real world, couples are heterogeneous with respect to their fecundability. Given that the expected waiting time to conception for each couple is the inverse of its fecundability, the mean waiting time in the population as a whole is simply the harmonic mean of fecundability [29], on which couples with a very low fecundability exert a large influence. Thus, given an observed mean waiting time, not accounting for heterogeneity can lead to a fecundability estimate that is significantly lower than the mean fecundability effectively found in the population.

Figure 7 compares fecundability as estimated with our original model with the estimates obtained with a model that accounts for heterogeneity (see below). As expected, the estimate for mean fecundability increases as the variance of the distribution of fecundability increases. Therefore, depending on the objectives of the analysis and characteristics of the data analysed, assumptions about the distribution of fecundability might be necessary.

To obtain the estimates under different variance assumptions, we model heterogeneity in fecundability assuming that each couple $i$ has fecundability $\phi_i(x) = z_i \cdot \phi(x)$, where $z$ is Gamma distributed in the population: $z \sim \Gamma(s, s)$. Setting the shape and rate equal ensures that $E(z) = 1$. Different
<figure><img src="image_11.png" /><figcaption>**Figure 6.** Estimated versus observed age at marriage distributions. Shaded red areas represent densities for lognormal distributions with mean and standard deviation equal to the mean posterior values for $\mu_m$ and $\sigma_m$ (see Appendix A, figures 8–9). Purple bars represent observed ages at marriage in the three samples analysed. (a) Hutterites, (b) eighteenth-century Quebec, (c) seventeenth–eighteenth-century France.</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>**Figure 7.** Fecundability estimates. Model comparison for French cohorts. Triangles represent the estimated probability of conception at each age in the original model (as displayed in figure 5a). The rest of the curves show estimated fecundability by age for a model that accounts for differences in the probability of conception between women. Each line corresponds to an assumption about the variance of the distribution of fecundability at the age fecundability is at its peak. The shaded areas represent 95% prediction intervals.</figcaption></figure>

heterogeneity scenarios can be obtained by fixing $\sigma^2 = s^{-1}$ at different, credible, values. Estimating the model for each scenario then gives posterior medians for $\phi_1$ and $\phi_2$. The variance values given in the legend of figure 7 are those for peak fecundability associated with these posterior median values, as variance at peak fecundability is the measure of heterogeneity that is most often found in the literature.

# 5. Discussion

The solutions at the aggregate level have predominated in dealing with the need of representing the age distribution of fertility rates. These solutions traded parameter interpretability in exchange for statistical formalization. Thanks to the approach we introduced in this article, this compromise is no longer necessary.

The individual-level model we introduce here is undeniably a simplistic, if not crude, description of the reproductive process. Nonetheless, it provides important insights into the behavioural and
biological mechanisms underlying differences in ASFRs between populations, with minimal data requirements.

In particular, our results consistently show that it is possible to obtain a very good approximation of the timing of marriage for a cohort of women by exclusively relying on the information contained in a schedule of ASFRs. With respect to fecundability and the length of post-partum amenorrhea, the quality of the approximation will depend on the degree to which the assumptions implied in the structure of the model are met in the population analysed. In the case of fecundability, we showed that the estimation obtained depends on the variance of fecundability in the population analysed. In other words, as originally implemented, our model recovers the age pattern of fecundability implied by the data, under the assumption of homogeneous fecundability. When heterogeneity with respect to fecundability in the population of interest is high, or the objective of the analysis is to obtain precise fecundability estimates, our approach can be easily extended to incorporate assumptions about the variance of its distribution.

More broadly, our approach can be seen as a tool to transform a simpler representation of some data into a much richer representation of the same information. In our particular application, this new representation allowed us to understand differences in aggregate fertility rates in terms of differences in individual-level mechanisms across populations.

This approach is particularly useful in contexts where a process of interest is only observed through simpler, aggregate data. Such a scenario arises in biomography, for example, when the demographic dynamics of wild animal populations are measured through census counts, instead of longitudinal monitoring data. For some of these populations, our model can be readily applied by simply adjusting the interpretation of age at marriage to age at the onset of reproduction.

The most obvious scenario, however, is contemporary human fertility, given that most vital statistics systems only collect count data on births. In fact, one of the objectives of our article is to spur the exploration of more complex models with the aim of understanding the microfoundations of present and future fertility patterns. These models will have to take into account the effect of individual preferences on the reproductive process, expressed both through the notion of a desired family size, as well as through the notion of short-term reproductive intentions that control the spacing of births across the reproductive lifespan.

The estimation of such models will very likely need additional information regarding, for example, the planned/unplanned nature of the observed births. In this regard, an important benefit of our approach is that such information can be incorporated as an additional target in the estimation procedure but can also be seamlessly integrated via prior distributions.

In light of the growing importance of data-driven, machine learning approaches in the study of natural and social systems, our emphasis on mechanisms may seem outdated. Ever larger amounts of data and computational resources combined with new algorithms seem to be shifting the balance in favour of pure-prediction versus explanatory approaches. However, as a number of recent studies suggest, the key to solving some of the most challenging problems in various disciplines may lie in a combination of mechanistic and machine learning approaches [33]. We hope the ideas presented here help bring together the mechanistic and statistical tools needed to further the study of demographic dynamics.

**Ethics.** This work did not require ethical approval from a human subject or animal welfare committee.

**Data accessibility.** Data are available at GitHub [34].

Electronic supplementary material is available online [35].

**Declaration of AI use.** We have not used AI-assisted technologies in creating this article.

**Authors' contributions.** D.C.: conceptualization, data curation, formal analysis, funding acquisition, investigation, methodology, project administration, resources, software, supervision, validation, visualization, writing—original draft, writing—review and editing; N.T.: conceptualization, data curation, formal analysis, investigation, methodology, software, supervision, validation, writing—original draft, writing—review and editing.

Both authors gave final approval for publication and agreed to be held accountable for the work performed therein.

**Conflict of interest declaration.** We declare we have no competing interests.

**Funding.** No funding has been received for this article.

**Acknowledgements.** We would like to acknowledge the valuable assistance provided by Diego Araujo and Ignacio Campón in the development of this project.

## Appendix A

Figures 8 and 9.
<figure><img src="image_13.png" /><figcaption>Figure 8. Adjusted posterior distributions.</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>Figure 9. Adjusted posterior distributions. Pink dashed lines show kernel density estimators for each sample.</figcaption></figure>
1. Schmertmann CP. Calibrated spline estimation of detailed fertility schedules from abridged data. *Rev. bras. estud. popul.* **31**, 291–307. (doi:10.1590/S0102-30982014000200004)

2. Coale AJ, Trussell TJ. 1974 Model fertility schedules: variations in the age structure of childbearing in human populations. *Popul. Index* **40**, 185–258.

3. Coale AJ, Trussell TJ. 1978 Technical note: finding the two parameters that specify a model schedule of marital fertility. *Popul. Index* **44**, 203–212.

4. Kostaki A, Paraskevi P. 2007 Modeling fertility in modern populations. *Demogr. Res.* **16**, 141–194. (doi:10.4054/DemRes.2007.16.6)

5. Chandola T, Coleman DA, Hiorns RW. 1999 Recent European fertility patterns: fitting curves to 'distorted' distributions. *Popul. Stud.* **53**, 317–329. (doi:10.1080/00324720308089)

6. Schmertmann CP. 2003 A system of model fertility schedules with graphically intuitive parameters. *Demogr. Res.* **9**, 81–110. (doi:10.4054/DemRes.2003.9.5)

7. Hoem JM, Madsen D, Nielsen JL, Ohlsen EM, Hansen HO, Rennermalm B. 1981 Experiments in modelling recent Danish fertility curves. *Demography* **18**, 231–244. (doi:10.2307/2061095)

8. Sheps MC, Menken JA, Radick AP. 1973 *Mathematical models of conception and birth*. Chicago, IL: University of Chicago Press.

9. Barrett JC. 1971 Use of a fertility simulation model to refine measurement techniques. *Demography* **8**, 481–490. (doi:10.2307/2060684)

10. Le Bras H. 1993 Simulation of change to validate demographic analysis. In *Old and new methods in historical demography* (eds D Reher, R Schofield), pp. 259–280. Oxford, UK: Clarendon Press.

11. Rubin DB. 1984 Bayesianly justifiable and relevant frequency calculations for the applied statistician. *Ann. Statist.* **12**, 1151–1172. (doi:10.1214/ aos/1176346785)

12. Beaumont MA. 2019 Approximate Bayesian computation. *Annu. Rev. Stat. Appl.* **6**, 379–403. (doi:10.1146/annurev-statistics-030718-105212)

13. Gini C. 1924 Premières recherches sur la fécondabilité de la femme. In *Proc. Int. Mathematical Congress, Toronto, Canada, 11–16 August 1924*, vol. 2. Toronto, Canada: University of Toronto Press.

14. Bendel JP, Hua C. 1978 An estimate of the natural fecundability ratio curve. *Soc. Biol.* **25**, 210–227. (doi:10.1080/19485565.1978.9988340)

15. Larsen U, Yan S. 2000 The age pattern of fecundability: an analysis of French Canadian and Hutterite birth histories. *Soc. Biol.* **47**, 34–50. (doi:10.1080/19485565.2000.9989008)

16. Weinstein M, Wood JW, Stoto MA, Greenfield DD. 1990 Components of age-specific fecundability. *Popul. Stud.* **44**, 447–467. (doi:10.1080/0032472031000144846)

17. Lee SC, Brattrud A. 1967 Marriage under a monastic mode of life: a preliminary report on the Hutterite family in South Dakota. *J. Marriage Fam.* **29**, 512. (doi:10.2307/349594)

18. Ingoldsby BB, Stanton ME. 1988 The Hutterites and fertility control. *J. Comp. Fam. Stud.* **19**, 137–142. (doi:10.3138/jcfs.19.1.137)

19. Eaton JW, Mayer AJ. 1953 The social biology of very high fertility among the Hutterites; the demography of a unique population. *Hum. Biol.* **25**, 206–264.

20. Mange AP. 1964 Growrh and inbreeding of a human isolate. *Hum. Biol.* **36**, 104–133.

21. Sheps MC. 1965 An analysis of reproductive patterns in an American isolate. *Popul. Stud.* **19**, 65–80. (doi:10.1080/00324728.1965.10406005)

22. Vézina H, Bournival JS. 2020 An overview of the BALSAC population database: past developments, current state and future prospects. *Hist. Life Course Stud.* **9**, 114–129. (doi:10.51964/hlcs9299)

23. Eijkemans MJC, van Poppel F, Habbema DF, Smith KR, Leridon H, te Velde ER. 2014 Too old to have children? Lessons from natural fertility populations. *Hum. Reprod.* **29**, 1304–1312. (doi:10.1093/humrep/deu056)

24. Séguy I. 2001 *La population de la France de 1670 à 1829: l'enquête Louis Henry et ses données*. Paris, France: INED.

25. Henry L. 1968 The verification of data in historical demography. *Popul. Stud.* **22**, 61–81. (doi:10.1080/00324728.1968.10405526)

26. Tavaré S, Balding DJ, Griffiths RC, Donnelly P. 1997 Inferring coalescence times from DNA sequence data. *Genetics* **145**, 505–518. (doi:10.1093/ genetics/145.2.505)

27. Pritchard JK, Seielstad MT, Perez-Lezaun A, Feldman MW. 1999 Population growth of human Y chromosomes: a study of Y chromosome microsatellites. *Mol. Biol. Evol.* **16**, 1791–1798. (doi:10.1093/oxfordjournals.molbev.a026091)

28. Blum MG. 2018 Regression approaches for ABC. In *Handbook of approximate Bayesian computation*, pp. 71–85. London, UK: Chapman and Hall/CRC.

29. Leridon H. 1977 *Human fertility: the basic components*. Chicago, IL: University of Chicago Press.

30. Baird DD, Strassmann BI. 2000 Women's fecundability and factors affecting it. In *Women & health* (ed. MB Goldman), pp. 126–137. New York, NY: Academic Press. (doi:10.1016/B978-012288145-9/50014-0)

31. Bongaarts J. 1980 Does malnutrition affect fecundity? A summary of evidence. *Science* **208**, 564–569. (doi:10.1126/science.7367878)

32. Todd N, Lerch M. 2021 Socioeconomic development predicts a weaker contraceptive effect of breastfeeding. *Proc. Natl Acad. Sci. USA* **118**, e2025348118. (doi:10.1073/pnas.2025348118)

33. Procopio A, Cesarelli G, Donisi L, Merola A, Amato F, Cosentino C. 2023 Combined mechanistic modeling and machine-learning approaches in systems biology: a systematic literature review. *Comput. Methods Programs Biomed.* **240**, 107681. (doi:10.1016/j.cmpb.2023.107681)

34. Ciganda D, Todd N. 2024 Individual_fx. GitHub. See https://github.com/dciganda/individual_fx.
35. Ciganda D, Todd N. 2024 Supplementary material from: Modelling the age pattern of fertility: an individual-level approach. Figshare. (doi:10.6084/m9.figshare.c.7477930)