## Abstract

Most work on the human fertility transition has focused on declines in mean fertility. However, understanding changes in the *variance* of reproductive outcomes can be equally important for evolutionary questions about the heritability of fertility, individual determinants of fertility and changing patterns of reproductive skew. Here, we document how variance in completed fertility among women (45–49 years) differs across 200 surveys in 72 low- to middle-income countries where fertility transitions are currently in progress at various stages. Nearly all (91%) of samples exhibit variance consistent with a Poisson process of fertility, which places systematic, and often severe, theoretical upper bounds on the proportion of variance that can be attributed to individual differences. In contrast to the pattern of total variance, these upper bounds increase from high- to mid-fertility samples, then decline again as samples move from mid to low fertility. Notably, the lowest fertility samples often deviate from a Poisson process. This suggests that as populations move to low fertility their reproduction shifts from a rate-based process to a focus on an ideal number of children. We discuss the implications of these findings for predicting completed fertility from individual-level variables.

## 1. Introduction

The recent transition from high to low fertility in many human populations is one of the best documented demographic regularities of the modern era, a finding supported by volumes of data across diverse countries and historical periods [1]. In contrast to widely available data on mean fertility, much less is known about how other aspects of fertility distributions—most notably their variance—change as means decline [2–4]. Although the vast majority of demographic research has neglected changing variances in the study of fertility transitions, an evolutionary perspective demands attention to them [5–7]. Heritable variation is the raw material that natural selection works on, and efforts to estimate the heritability of fertility at different stages of the transition rely on the magnitude of phenotypic variance [2,3,8–13]. The determinants of variance in fertility and reproductive skew is an active area in evolutionary theory in its own right [14–16]. Moreover, demographers from both traditional and evolutionary backgrounds frequently aim to identify individual determinants of fertility [17]. Identifying the main sources of variance in fertility can shed light on observed correlations between individual fertility and differences in personality [18], physiology [19,20], ecology [21,22] and family history [8,11].

To identify key sources of variance in fertility, it is useful to examine the process giving rise to it. Completed fertility in humans arises through a counting process of mostly singleton births punctuated by variable wait times. One of the simplest classes of such a counting process, a Poisson process, illustrates how substantial variation in lifetime fertility can arise without any underlying differences between individuals [23]. Imagine 100 individuals all of whom have the same window of opportunity for reproduction (i.e. 20 years) and the same constant rate of giving birth (i.e. one every 4 years). Although the rate is constant, each individual's reproductive history will differ in actual timing of conception and the number of losses before a successful pregnancy [23]. This exceedingly simple process leads to a distribution of completed fertility across the 100 identical individuals showing substantial variability—roughly 10 people having seven children, three having two and only 18 people having the expected value of five children. Many complexities one might add to such a Poisson process, such as allowing interbirth intervals to increase over time or assuming some offspring will die with a certain probability, will not attenuate the stochastic variability introduced by the counting process [24]. In short, starting with a sample of identical participants, a straightforward counting process can create substantial variance in completed fertility that has nothing to do with individual variability in rates or windows of opportunity and everything to do with the stochastic timing of events [25–29].

Individual-level variables, such as family background, education and wealth, can only account for non-stochastic variation. For that reason, stochastic variation potentially introduced by a counting process can attenuate the kinds of associations between completed fertility and underlying individual differences often investigated by demographers. This can also limit the ability to detect such associations [25–28]. Thus, if we are given a sample's mean completed fertility without its variance, we have no way of judging how much of the variance in completed fertility is potentially attributable to true variation between individuals or is equally consistent with a Poisson-like process where everyone reproduces at a similar underlying rate. In the latter case, we should expect vanishingly small correlations between individual-level variables and completed fertility since there are minimal true individual differences. In the former case, even if there were some true variation in rates, the amount of variance potentially due to individual differences might be small (e.g. 5% of variation) or large (e.g. 60% of variation) and this will affect the expected size of correlations and the power required to detect such correlations. Substantial stochastic variation in completed fertility might account for the difficulties demographers have often faced in accounting for large proportions of variation in fertility based on between-individual variables (e.g. education, wealth, occupation, ethnicity, long-term physiological status and other family attributes) [30].

A few studies have used rich historical data to track changing variance in fertility across time in select European populations and have shown both increasing and declining variance over time [2–4]. In this paper, we analyse data from a much wider range of worldwide populations—200 Demographic and Health Surveys (DHSs) in 72 countries at varying levels of fertility (mean live offspring born from 1.8 to 8.7)—to map how the variance in completed fertility changes with declining mean fertility. First, we track changes in variance with changing mean fertility within countries that have data at multiple time points. Second, we examine whether samples have variance consistent with a Poisson process with varying rates. Third, for those samples with variance consistent with a Poisson process with varying rates, we estimate the theoretical upper bound on the proportion of variance due to individual heterogeneity. Fourth, we compare these theoretical upper bounds with the amount of variation empirically accounted for by three of the demographic variables most frequently associated with fertility differences—education, wealth and rural residence. Finally, we examine how these theoretical upper bounds depend on a sample's completed fertility.

## 2. Models and methods

### (a) Sample

We use data from 200 DHSs across 72 low- to middle-income countries that have collected systematic, comparable data on women's reproductive histories [31]. We focus on women for two practical reasons. First, the DHSs have primarily aimed their data collection efforts at women and children, and relatively few have collected comparable data for men. Second, even for DHS datasets that have information about male marital fertility, varying opportunities for extra-marital fertility make it impossible to estimate male fertility in the same way we assess female fertility. We consider DHSs collected between 1990 and 2014 that include information about fertility and household wealth [32].

### (b) Measures of fertility

We examine estimates of completed fertility for women aged 45–49 years as total live born children. In the most comprehensive study to date, 90% of women ended reproduction by 45 years and nearly 100% by 50 years [33]. Thus, focusing on 45–49 year olds will capture most live births.

### (c) Demographic variables

Education, rural residence and wealth have been shown to be some of the most important demographic predictors of fertility-related outcomes [34–36]. Thus, we include these three variables in regressions predicting completed fertility—education as a four-category dummy variable (no education, primary education, secondary education and higher education), rural residence as a dichotomous variable, and wealth as a 14-category dummy variable indicating different levels of household wealth *per capita* [37] based on 2011 international units, purchasing power parity [38].

### (d) Counting models of fertility

The simplest model for offspring counts is a Poisson process that assumes all women have the same rate of offspring production and the same window of opportunity for reproduction. According to this admittedly simple model, any variation in realized fertility is due to stochastic timing of events, and the variance of the distribution should be very close to the mean. Figure 1*a* illustrates a sample of completed fertility—Bangladesh 2004 DHS—with a close fit to the distribution expected from a simple Poisson process.

In a real population, the variance in completed fertility can be higher than expected if every woman followed a Poisson process with the same rate (e.g. over-dispersion). Figure 1*b* shows how the distribution of completed fertility in the Honduras 2006 DHS is wider than the best fitting distribution from a simple Poisson process. One potential reason for such over-dispersion is heterogeneity among women in their underlying rates due to inheritance from parents, physiological status, education, wealth, occupation, ethnicity or some other individual-level variable. A Poisson regression can be used to identify how women vary based on observed variables, such as wealth and education, but does not permit additional unexplained heterogeneity. A generalization of the Poisson model, the Gamma-Poisson (or negative binomial), permits such variation in underlying rates (following a gamma distribution). In addition to parameters for the mean lifetime rate of reproduction, the Gamma-Poisson includes a parameter for this rate variance across women [23,39–41]. Figure 1*c* is the same empirical distribution from Honduras 2006 DHS, but now compared with the best fitting distribution from a Gamma-Poisson process which reflects the wider shape of the empirical distribution.

It is also possible that a sample has lower variance than expected from a simple Poisson process, which is called under-dispersion. Figure 1*d* shows the distribution of completed fertility for the Ukraine 2007 DHS, compared to the best fitting distribution from a Poisson process. In this case, the theoretical distribution is too wide and has far too low expectations for parities two and three. In such cases, at least some of the sample is likely following a process of reproduction that is inconsistent with a Poisson process.

For those samples that are consistent with a Poisson process, a negative binomial modelling framework allows us to estimate what part of the variation in a population is due to stochastic Poisson processes, and what parts are due to differences in underlying individual rates. As long as the variance is larger than the mean then this estimate of the proportion of variance potentially due to individual heterogeneity can closely be approximated by 1 − Mean/Variance. For example, in a population with an average of 5.6 offspring, and a fertility variance of 6.8, then about 17% (1–5.6/6.8) of the observed variance can be attributed to stable between-individual differences in rate. This approach to characterizing variance in fertility is different from other common approaches, in that it relies on linear rather than quadratic scaling between the mean and variance. For interested readers, we present results in the electronic supplementary material for other measures of variance of potential interest, including Crow's Index of Selection due to Fertility (Variance/Mean^{2}) [2,6] and the coefficient of variation (s.d./Mean) [10,11,42], which assume quadratic scaling between mean and variance.

### (e) Analyses

We conduct analyses to answer five questions. For all tests run across the 200 samples, we correct inferences for multiple tests using the Holm–Bonferroni method with an initial alpha of 0.05.

#### (i) How does variance change with changing means?

We plot the overall sample means and variances for fertility across the 200 samples, with lines of least square fit for countries with multiple surveys over time. We also compare these values to a theoretical lower bound on variance from a simple Poisson process.

#### (ii) To what degree are samples under- and over-dispersed?

We estimate the degree of under- and over-dispersion in each of the 200 samples using a quasi-Poisson regression (glm in R, and dispersiontest in AER package) [43,44]. For each sample, we run: (i) an intercept-only null model to assess raw dispersion, and (ii) a model including the three demographic variables—education, rural residence and wealth—to assess dispersion after accounting for variation due to these key individual-level variables. A raw dispersion parameter from the intercept-only model is consistent with a mixture of Poisson processes for any value greater than or equal to unity. For subsequent analyses, we focus on those samples whose variance is consistent with a Poisson process (i.e. a raw dispersion parameter ≥ 1).

#### (iii) What is the theoretical upper bound on proportion of variance explained by individual-level variables?

For each of the samples that had a dispersion parameter greater than unity, we fit a negative binomial regression in R (glm.nb) to the sample's distribution of completed fertility [45]. We assess whether variance in rates is different from zero using a likelihood ratio test. We estimate the proportion of observed variation in fertility due to individual rate heterogeneity as the excess variance due to rate heterogeneity (excess variance = *μ*^{2}/*θ*) divided by the total variance of the fitted negative binomial distribution (total variance = *μ* + (*μ*^{2}/*θ*)). This provides a theoretical upper bound on the proportion of total variance that might be accounted for by individual-level variables.

#### (iv) To what degree is the theoretical upper bound consistent with empirical analyses?

For each of the samples that had a dispersion parameter greater than unity, we also estimate the proportion of variance in completed fertility accounted for by three demographic variables—education, wealth and rural residence. To do this, we fit two negative binomial regressions in R—one with no predictors (i.e. the null model) and one with education, wealth and rural residence predicting each individual's rate (i.e. the full model). We calculate the proportion of variance explained as 1 − SSE_{full}/SSE_{null}, where SSE is the sum of square errors from the full and null models, respectively. We then plot this proportional reduction in model residuals versus the theoretical upper bound for the samples.

#### (v) How does the theoretical upper bound depend on completed fertility?

For these analyses, we categorize samples consistent with a Poisson process into classes based on their mean completed fertility (completed fertility less than 3.5 = 3, 3.5 to 4.5 = 4, 4.5 to 5.5 = 5, 5.5 to 6.5 = 6, 6.5 to 7.5 = 7, greater than 7.5 = 8). We summarize: (i) variance expected from a single-rate Poisson process (equal to the mean rate), and (ii) excess variance due to individual heterogeneity in rates (excess variance = *μ*^{2}/*θ*). We also summarize how the theoretical upper bound on proportion of variance explained by individual variables (defined in §2e(iv)) depends on the mean completed fertility across these categories.

## 3. Results

### (a) Sample descriptives

We analysed 200 surveys collected from 72 countries between 1990 and 2014. These included 234 261 women (45–49 years) who had non-missing data on total liveborn children (fertility). Survey samples ranged from 115 (Dominican Republic 1999) to 9806 (India 2005) and mean fertility of 8.7 (Afghanistan 2010) to 1.8 (Ukraine 2007). The average fertility across all surveys was 5.8. Figure 2*a* illustrates the relationship between mean fertility in each survey and mean household wealth *per capita*. The outlying points in the upper right corner represent surveys from Jordan between 1990 and 2012, suggesting a late start to the fertility transition relative to the country's current household wealth *per capita*.

### (b) How does variance in completed fertility change with changing mean?

As expected from a Poisson process, nearly all countries showed a downward trend in variance with declining mean fertility (figure 2*b*). Despite declining variance, nearly all samples maintain a variance above the theoretical minimum variance of a Poisson process (figure 2*b*). The few samples falling below the variance expected from a Poisson process are at very low (Armenia, Kyrgyzstan) and high (Burkina Fasa, Rwanda) ends of fertility.

### (c) To what degree are samples under- and over-dispersed?

Of the 200 samples, 182 (91%) had dispersion greater than unity in the intercept-only quasi-Poisson regression. *In subsequent sections using the negative binomial model, we focus on these 91% of samples*. The 18 samples with dispersion less than unity were concentrated primarily in low fertility populations (45% of 22 samples with fertility less than 4), with some in high fertility populations (19% of 32 samples with fertility greater than 7). Only two cases of under-dispersion occurred in intermediate fertility populations (1% of 146 samples).

Dispersion tests indicate that most of the 200 samples (147 or 74%) are significantly over-dispersed. A small percentage (8 or 4%) were significantly under-dispersed, with all but one (Burkina Faso 2010) coming from relatively low fertility countries in Europe and Central Asia (Ukraine 2007, Armenia 2000, 2005 & 2010, Moldova 2005, Albania 2008 and Kyrgyzstan 2012). Adding three demographic predictors—education, wealth and rural residence—to the quasi-Poisson model regression makes three additional low fertility samples (1.5%) significantly under-dispersed (Vietnam 2002 and 2005, Dominican Republic 2013).

### (d) What is the theoretical upper bound on proportion of variance explained by individual-level variables?

Of the 182 samples with dispersion greater than unity, the negative binomial regression results indicated that 91% (165) had variances in completed fertility that were significantly greater than would be expected if women were following the same underlying rate. The excess variation attributable to individual differences in lifetime fertility rates ranged from zero to 68% of the total observed variation (mean = 34%).

### (e) To what degree is the theoretical upper bound consistent with empirical analyses?

The negative binomial regression estimated that the proportional reduction in variance in completed fertility due to the three demographic variables ranged from zero to 40% (mean = 12%, s.d. = 8%, *n* = 182). Figure 3 plots the relationship between the theoretical upper bound on variance explained by individual-level variables and the actual variance explained by three key demographic variables. Samples with observed values below the theoretical upper bound (and thus consistent with the upper bound) are represented by points below the diagonal. Very few cases (5%) go above the diagonal defining the theoretical upper bound, and these were on average only 7% above that bound (range = 1–7%). On the basis of chance alone we would expect some samples to fall above the line, but this could also reflect that in these populations, a Poisson process was not an accurate description of fertility patterns.

### (f) How does the theoretical upper bound depend on completed fertility?

Figure 4*a* depicts how the total variance in completed fertility changes with changing mean fertility along with its two theoretical components. Total variance increases substantially with increased mean completed fertility. However, examining the two components of variance, we see that a large part of this increase is likely due to increased stochastic variation. Figure 4*b* illustrates how these trends in total variance and its components affect the upper bound on the proportional reduction in variance due to individual-level variables. Notably, we see that the largest upper bounds occur among middle fertility samples rather than the highest fertility samples.

## 4. Discussion

In this paper, we outline how variance in women's completed fertility changes as completed fertility declines across 200 samples from 72 low- to middle-income countries. Most of the 200 samples (91%), and nearly all of the samples in middle and high fertility settings were consistent with a model of fertility based on a Poisson process. If something roughly similar to a Poisson process underlies reproduction in these samples, then this puts systematic upper bounds on the amount of variation that can be attributed to individual differences in variables such as personality, physiology, ecology and family background. We further show that this theoretical upper bound is consistent with empirical findings about the proportional reduction in variance when predicting completed fertility with education, wealth and rural residence—three demographic variables most commonly found to be associated with fertility (figure 3).

The upper bounds identified in this way are also consistent with recent analyses in other samples from mid and high fertility settings that focus on predicting completed fertility from another individual-difference variable—the completed fertility of genetic kin. Among high fertility Hutterite females (completed fertility = 7.1), heritability is estimated as 0.20 which is just below the estimated proportion of between-individual variation of 0.26 [46]. Among Finnish populations at different stages of the fertility transition, heritability estimates fall below the theoretical upper bound by much more (pre-1880 completed fertility = 5.1, *h*^{2} = 0.11 of 0.55 available; post-1880 completed fertility = 2.5, *h*^{2} = 0.17 of 0.60 available) [3].

This approach permits us to partition variance in completed fertility into two components—that potentially due to individual differences and that due to stochastic events. With this partition, we further show that although total variance uniformly decreases with declining mean completed fertility, the proportion of variance that is potentially attributable to individual differences is on average highest in intermediate fertility settings (figure 4*b*) [47]. This may reflect increasing individual heterogeneity in these countries in terms of education and wealth. While this makes it potentially easier to detect the effect of individual-level variables on observed fertility in such populations, it also means comparisons of the effects of individual-level variables on fertility across low, middle and high fertility regimes should be interpreted in light of these upper bounds [8].

While the theoretical upper bound changes systematically with changing mean fertility, figure 3 also illustrates that individual samples vary considerably in their upper bound—ranging from close to zero to 68% of total variance. The formula 1 − mean/Variance gives a rough guideline to the maximum possible correlations expected with individual-level variables in cases where the variance is greater than the mean. The upper bound identified here is similar to the partitioning of variance into between- and within-group components using multilevel models that are then used to identify the sources of variation at these two levels. In the same way, the upper bound indicates how much of the variance is potentially due to individual differences in a given sample, and gives a guideline for the maximum associations one might expect to find with specific variables—genetic, cultural and ecological.

Despite quite severe upper bounds on the proportion of within-population variance attributable to individual differences, nearly all samples had variances in completed fertility that were significantly greater than would be expected if women were following a single underlying rate. We show here (as expected from prior demographic work) that an often substantial proportion of this excess variation can be attributed to individual differences in household wealth, education and rural residence. As shown in figure 2*a*, these within-population results are also consistent with associations of mean household wealth and completed fertility across: (i) populations from different countries and (ii) populations from the same country over time.

Our analyses have focused mainly on low- and middle-income samples with middle to high mean completed fertility. While a Poisson process model is consistent with completed fertility in most of these samples, the few samples that showed under-dispersion point us to periods in the transition where a Poisson model needs to be substantially modified or even replaced. For example, there were a few very high fertility samples that showed under-dispersion (in Rwanda and Burkina Faso, mean completed fertility approx. 7–8). Further examination of the distribution of completed fertility in these high fertility samples confirms established studies of fecundity that show completed fertility should hit a wall at the highest levels of fertility due to natural limits on female fertility [23]. This may be resolved in part by fitting a counting process that takes some time (e.g. 12 months) to restart after a birth event, which would put a natural upper limit on the number of children produced in a normal reproductive lifespan.

The more common cases of under-dispersion in low fertility settings may need a complete change in models. These cases of under-dispersion are striking for two reasons. First, they are more severely under-dispersed than the high fertility cases. Second, in many of the low fertility samples with under-dispersion, we also see substantial amounts of variance accounted for by individual-level variables. For example, in five of the 18 samples with dispersion parameters of less than unity, we see that the three demographic variables, education, wealth and rural residence, account for more than 10% of the variance in completed fertility (Moldova, Albania, Kyrgizstan, Vietnam 2002 and 2005). These findings suggest that at the late stages of the fertility transition, there may be a shift from a rate-driven Poisson process to one based on stopping at a limited set of family sizes (e.g. 0, 1 or 2 children). This would be consistent with demographic theories proposing a fundamental shift in decision-making rather than a simple change in how many children women are choosing to have [47–49]. This would also give rise to a new kind of variation between populations. In addition to variation in mean rates, we would also begin seeing variation in the proportion of individuals who are stopping at specific parities. Newly developed multiple-hurdle models likely provide one way forward in modelling this potentially important shift in decision-making [50]. Models in this vein could permit, for example, estimates of how much samples over-prefer one, two or three children relative to a Poisson model, and could determine if under-dispersion is due to such targeted preferences in the low range of fertility.

More generally, a Poisson process is a reasonable first-order approximation to completed fertility in human populations from low- and middle-income countries, and thus serves as useful starting point. However, future work should assess key assumptions of this base Poisson model and make necessary refinements—such as declining rates with age, twinning, excess proportions of nulliparous individuals, varying windows of opportunity for reproduction, different patterns of interbirth intervals, as well as mixtures with other distributions in low fertility settings. This will hopefully help determine what aspects of the fertility process are most important for observed changes in fertility distributions.

## 5. Conclusion

An evolutionary approach to reproduction naturally draws attention to the variance in reproduction. Here, we show how expanding the focus from means to variances (and distributions more generally) can: (i) identify the upper bounds to variation in completed fertility that might be accounted for by individual-level variables and (ii) track how this upper bound changes over the course of the fertility transition. These upper bounds can provide important context for interpreting the effect of individual-level variables across diverse samples. Moreover, in places where simple process models break down—such as the inability of a Poisson process to account for the distributions of completed fertility in very high and low fertility settings—an approach that explicitly specifies process models leads us to ask new questions and to propose testable refinements that can potentially shed light on shifting patterns of reproductive decision-making.

## Ethics

The study was approved by ASU's Office of Research and Integrity IRB (Protocol no. 1302008836).

## Competing interests

We declare we have no competing interests.

## Funding

We received no funding for this study.

## Acknowledgements

We thank Mary Towner, Gert Stulp, Joe Hackman, Heidi Colleran, Paul Hooper and members of both the NESCENT Evolutionary Demography of Fertility Working Group and the National Science Foundation-funded research seminar at the School for Advanced Research for helpful comments and conversations.

## Footnotes

One contribution of 14 to a theme issue ‘Understanding variation in human fertility: what can we learn from evolutionary demography?’

- Accepted December 29, 2015.

- © 2016 The Author(s)

Published by the Royal Society. All rights reserved.