Study population

In the nineteenth and twentieth centuries, Switzerland underwent large-scale social and economic changes in connection with industrialisation, urbanisation and the expansion of the service sector. After 1850 the city of Basel experienced rapid economic and population growth (see Appendix, Figure S1A) and developed from a typical trading town to a major chemical and pharmaceutical centre. As early as the end of the nineteenth century, Switzerland's gross domestic product put it among the wealthiest nations in Europe (Maddison, Reference Maddison2001). Between 1850 and 1950, the real wages of skilled workers in Basel increased by a factor of two (although temporarily interrupted by World War I, Figure S1B; Studer, Reference Studer2008; Studer & Schuppli, Reference Studer and Schuppli2008). Infant mortality rate (Figure S1C), which was generally higher in Swiss cities than in rural parts until the end of the nineteenth century, decreased at a steady rate for Basel, as did the stillbirth rate (Figure S1D) (Floris & Staub, Reference Floris and Staub2019). Deaths owing to infectious diseases had been decreasing since the second half of the nineteenth century, also owing to large improvements in sanitation: access to clean water and flush sewage systems were provided in Basel from the 1860s and constantly expanded since then (Floris & Staub, Reference Floris and Staub2019). Consequently, life expectancy increased (Floris, Höpflinger, Stohr, Studer, & Staub, Reference Floris, Höpflinger, Stohr, Studer and Staub2019a; Floris, Staub, Stohr, & Woitek, Reference Floris, Staub, Stohr, Woitek, Juan, Hürlimann, Lorenzetti and Schiedt2019b). The average height of 19-year-old conscripts in Basel (covering 90–100% of male birth cohorts with Swiss citizenship) increased from 168.9 cm for young men born in 1894 to 173.8 cm for those born in 1938 (Figure S1E; Koepke et al., Reference Koepke, Floris, Pfister, Rühli and Staub2018; Schoch, Staub, & Pfister, Reference Schoch, Staub and Pfister2012). Thus conventional measures of living standards, monetary as well as biological, describe a steady improvement in overall well-being during this study period.

The demographic transition was in progress in Basel from the last quarter of the nineteenth century, when first deaths and later also births per 1000 inhabitants started to decrease (Figure S1G). Mechanical and chemical contraception and birth control were increasingly used in Switzerland and Basel from 1900, with higher usage being suspected among women of more advantageous socioeconomic positions (Labhardt, Reference Labhardt1930). An examination of the fertility of all Basel marriages concluded in 1921/1922 until the end of their reproductive phase in the 1940s showed that a total of 82% of those marriages had 0–2 children (25% of marriages remained childless; Labhardt, Reference Labhardt1943). Basel thus had a comparatively low number of children per mother (2.00). Nevertheless, parities of up to 19 are present in the birth records (see Table 1). The average marriage age of the women was 25 years, and the first child was born on average in the third year after marriage (in 78% of the marriages within the first 6 years). Furthermore, the percentage of children breastfed increased from 90% in 1911 to a constant 95% after 1920 (Jann, Reference Jann1934).

Birth record data

The present data (see Table 1) come from the birth records of the university maternity hospital (‘Frauenspital’) of the canton Basel-Stadt. Detailed data on each childbirth have been routinely recorded since 1896. The register books are carefully maintained and incomplete records are very rare, although only a third of the records were archived owing to space constraints. Access to the protected individual data was allowed by the Staatsarchiv Basel-Stadt upon signed contractual agreement. After linking the sources, the data have been fully anonymised. The maternity hospital was founded in 1868 and by 1897 around 30% of all births in Basel took place there. This increased to 50% by 1912 and reached approximately 75% by 1930 (Figure S1F). Compared with other cities at that time the proportion of missed home births and births in other hospitals is therefore relatively low. Furthermore, birth records in Basel cover births from women of both high and low socioeconomic position, as well as complicated and problem-free childbirths. Nevertheless, each individual woman's birth history is likely to be incomplete owing to the missing record books.

The birth records contain information about birth outcomes as well as maternal condition (Table 1). Specifically, in addition to maternal age and marital status, height was routinely measured, age at menarche was recorded and the physician provided a crude assessment of the mother's nutritional status and body shape; both of the latter were coded as three-point scales reflecting increasing health and embodied capital (see ‘Data issues’).

Based on a continuous and unique entry number for each year, nearly 100% of the birth records could be linked to the birth register (inventory Sanität X8). From this second register, additional variables regarding the socioeconomic background of the father (if known) were available. Our measure of socioeconomic position (SEP) was derived from the occupations of the father or the mother using the classification in Schüren (Reference Schüren1989). The results of this classification are six categories, which we aggregated to three groups: low (1, 2), medium (3, 4) and high SEP (5, 6). We use the SEP of the father, if available, and otherwise the SEP of the mother. However, records where mothers indicated ‘housewife’ and the SEP of the father was not available were treated as missing SEP data. Comparing the shares of father's occupation with the census shares in 1920 for male occupation groups shows that our sample closely reflects the SEP distribution in this population (Supplementary Table S1).

Data analysis

We used a Bayesian multivariate, multilevel model implemented in the ‘brms’ package (Bürkner, Reference Bürkner2017) in R. A multivariate model allowed us to simultaneously analyse multiple related outcomes (e.g. birth weight, birth length, etc.) in a single model while accounting for any residual correlations among them not captured by our primary predictors; for instance, all else being equal, heavier babies are also expected to be longer, have larger heads, etc. Accounting for such general correlations among outcome variables thus made our analysis more accurate and conservative, because parameter estimates were pooled across outcomes and extreme estimates shrunk towards the mean (McElreath, Reference McElreath2020), compared with treating all outcomes as independent in separate models. Fitting these models in brms also facilitated appropriate modelling of nonlinear effects using splines for variables with probably fluctuating effects (year) and monotonic effects (Bürkner & Charpentier, Reference Bürkner and Charpentier2020) for ordinal variables (SEP), as well as imputing missing data and propagating the associated uncertainty into final parameter estimates (McElreath, Reference McElreath2020; Zhou & Reiter, Reference Zhou and Reiter2010); specifically, we imputed 10 complete datasets and ran our models on each of them and then combined the posterior samples (Zhou & Reiter, Reference Zhou and Reiter2010). Based on our experience transcribing the data and inspecting missingness patterns, we have no reason to believe that missingness was non-random, and Figure S2 shows that imputed values have high plausibility as evidenced by density plots overlapping with the observed data (van Buuren & Groothuis-Oudshoorn, Reference van Buuren and Groothuis-Oudshoorn2011). Imputation makes no additional assumptions on patterns of missingness compared with complete-case analysis, but has the huge advantage of using all of the valuable information rather than discarding incomplete cases, thus generally leading to less biased estimates (McElreath, Reference McElreath2020). Further, Bayesian estimation allowed us to use regularising priors to impose conservatism on parameter estimates and reduce our risk of inferential errors (McElreath, Reference McElreath2020).

Our approach also accounted for unobserved variation at various levels that may affect both parity and maternal/offspring condition and is likely to mask tradeoffs. As noted in the introduction, unobserved variation in energetic condition or access to resources commonly masks life-history tradeoffs (van Noordwijk & de Jong, Reference van Noordwijk and de Jong1986), including in humans (Sear, Reference Sear2020; Sear et al., Reference Sear, Mace and McGregor2003). Technically this problem is referred to as endogeneity bias, which arises from the residuals of a model being correlated with a predictor (Steele, Vignoles, & Jenkins, Reference Steele, Vignoles and Jenkins2007); here, condition correlates with parity owing to unobserved variation in factors affecting both. This bias can be resolved by jointly estimating the processes that determine parity and maternal/offspring condition and letting the residuals from their respective linear predictors be correlated in a multivariate model (Steele et al., Reference Steele, Vignoles and Jenkins2007). If these correlations are different from 0 this indicates that endogeneity was indeed present, and the sign of the correlation may indicate underlying variation in energetic condition (if parity and outcomes are positively correlated) or provide further evidence of tradeoffs not accounted for by the parity fixed effect (for negative correlations). The same logic applies to random intercepts for mothers, which were also allowed to be correlated across all outcomes in our model. Specifically, the mother random effect captures consistent individual differences in birth outcomes among mothers, for women appearing more than once in the data, independent of our fixed effects. Such individual differences are expected to result from unobserved factors such as overall phenotypic quality (often referred to as ‘frailty’ by demographers) that would otherwise mask underlying tradeoffs within mothers’ reproductive careers (Doblhammer & Oeppen, Reference Doblhammer and Oeppen2003). More specifically, akin to the correlated residuals mentioned above, a positive correlation of mother ID intercepts across outcomes, e.g. reflecting higher than average offspring condition and (age-specific) parity, would provide evidence of individual variation in frailty, while a negative correlation would provide further evidence for tradeoffs that is not accounted for by the parity fixed effect. Place of residence (Basel, n = 10,515; others, n = 3114) was included as a fixed effect to capture any unobserved variation in living conditions or social norms not accounted for by SEP and year. In sum, endogeneity bias stemming from phenotypic correlation can be resolved by using multivariate, multilevel models with correlated residuals and random effects, thus increasing the chance of detecting tradeoffs.

The general formula for our multivariate model is given below. For simplicity, we illustrate here only one measure of maternal and offspring condition, but the full set of measures included nutritional status and body shape for maternal condition, and live birth, birth weight, body length, gestational age, head circumference and placental weight for offspring condition (see Table 1). All of these outcome variables were modelled using Gaussian distributions except for nutritional status and body shape (both cumulative ordinal distributions), the probability of a live birth (Bernoulli) and parity (Poisson). Gaussian variables were standardised prior to analysis to facilitate comparison and interpretation of model parameters.

Our model structure for the ordinal response of nutritional status for birth record i is given by

$$\eqalign{& {\rm Nutritional\;statu}{\rm s}_i\approx {\rm Ordered}-{\rm logit}( {{\rm \phi }_i, \;\kappa } ) \cr & {\rm \phi }_i = \beta _1{\rm parit}{\rm y}_i\;\ldots + \beta _x{\rm offspring\;se}{\rm x}_i + \mu _{{\rm mother}[ i ] }^{{\rm NS}} + \varepsilon _i^{{\rm NS}} } $$

Here ϕ is the linear predictor for the probability of observing a mother across the κ categories, given the global intercepts μ _k = {μ _bad, μ _normal, μ _good} for each nutritional status. The log-cumulative-odds for observing a mother in any particular category k are therefore given by μ _k − ϕ. The β are regression coefficients for each of x fixed effects, which include the main predictor parity as well as other covariates including maternal age, place of residence, year, socio-economic position (SEP), marital status, maternal height, age of menarche, twin birth and offspring sex (see Tables 1 and 2). Here, height and age of menarche reflect nutrition and general conditions during a woman's development (taller and earlier menarche = better conditions; Bogin, Reference Bogin1999; Stearns & Koella, Reference Stearns and Koella1986), while SEP, marital status and year reflect current conditions (higher SEP and married, later years = better); twin and male births (Galbarczyk, Klimek, Nenko, & Jasienska, Reference Galbarczyk, Klimek, Nenko and Jasienska2019) are expected to impose additional costs on mothers. Note that SEP was modelled as a monotonic effect (Bürkner & Charpentier, Reference Bürkner and Charpentier2020) and year using a spline to allow for nonlinearity. A priori, we expected year to be nonlinear because, while conditions generally improved over the study period, World War I (1914–1918) and the Spanish flu (1918) presented major shocks that could have disrupted a linear increase. Random intercepts for nutritional status μ ^NS are specified for the mother, along with an observation-level random effect ɛ capturing residual variation (also known as overdispersion; Nakagawa, Johnson, & Schielzeth, Reference Nakagawa, Johnson and Schielzeth2017), which can then be correlated across outcomes to address the endogeneity bias (Steele et al., Reference Steele, Vignoles and Jenkins2007) as specified below.

Table 2. List of predictions and relevant model parameters that test them. The column on the right indicates qualitative support for each prediction. Note that for Predictions 1 and 2, the main effects of parity (β ₁) are the most pertinent tests; negative correlations (ρ) between residuals or mother-level intercepts of parity and outcomes would indicate energetic tradeoffs, but caused by some unobserved variable rather than by parity itself.

Similarly, our basic model structure for birth weight is given by

$$\eqalign{& {\rm Birth\;weigh}{\rm t}_i\approx {\rm Normal}( {{\rm \mu }_i, \;1} ) \cr & {\mu }_i = \mu _0 + \beta _1{\rm parit}{\rm y}_i\;\ldots + \beta _x{\rm offspring\;se}{\rm x}_i + \mu _{{\rm mother}[ i ] }^{{\rm BW}} + \varepsilon _i^{{\rm BW}} } $$

Here the global intercept is represented by μ ₀. We fix the scale of the residual standard deviation (σ = 1) because this parameter is redundant with the observation-level random effect ɛ. While mathematically equivalent, this parameterisation facilitates correlating residual error across the Gaussian and non-Gaussian responses.

Finally, for our count measure of parity the basic model structure is

$$\eqalign{& {\rm Parit}{\rm y}_i\approx {\rm Poisson}( {{\rm \lambda }_i} ) \cr & {\rm log}( {{\rm \lambda }_i} ) = \mu _0 + \ldots + \beta _x{\rm age\;at\;menarch}{\rm e}_i + \mu _{{\rm mother}[ i ] }^{\rm P} + \varepsilon _i^{\rm P} } $$

This model included the same fixed effects except twin birth and offspring sex as these are not relevant predictors of parity. To account for associations between the response measures, we correlate maternal and observation-level random effects. For example, we estimate correlations between the maternal random effects such that:

$$\eqalign{& \left[{\matrix{ {\mu_{{\rm mother}}^{{\rm NS}} } \cr {\mu_{{\rm mother}}^{{\rm BW}} } \cr {\mu_{{\rm mother}}^{\rm P} } \cr } } \right]\approx {\rm MVNormal}\left({\left[{\matrix{ 0 \cr 0 \cr 0 \cr } } \right], \;{\bf S}} \right)\cr & {\bf S} = \left({\matrix{ {\sigma_{\mu^{{\rm NS}}}} & 0 & 0 \cr 0 & {\sigma_{\mu^{{\rm BW}}}} & 0 \cr 0 & 0 & {\sigma_{\mu^{\rm P}}} \cr } } \right)\;{\bf R}\;\left({\matrix{ {\sigma_{\mu^{{\rm NS}}}} & 0 & 0 \cr 0 & {\sigma_{\mu^{{\rm BW}}}} & 0 \cr 0 & 0 & {\sigma_{\mu^{\rm P}}} \cr } } \right)\;\cr & {\bf R} = \left({\matrix{ 1 & {\rho_{\mu^{{\rm NS}}, \;\mu^{{\rm BW}}\;}} & {\rho_{\mu^{{\rm NS}}, \;\mu^{\rm P}\;}} \cr {\rho_{\mu^{{\rm BW}}, \;\mu^{{\rm NS}}\;}} & 1 & {\rho_{\mu^{{\rm BW}}, \;\mu^{\rm P}\;}} \cr {\rho_{\mu^{\rm P}, \;\mu^{{\rm NS}}\;}} & {\rho_{\mu^{\rm P}, \;\mu^{{\rm BW}}\;}} & 1 \cr } } \right)\;} $$

The covariance matrix S is parameterised such that inference is made directly on the correlation matrix R using the standard deviations (σ) of random effects (McElreath, Reference McElreath2020). The same parameterisation is used for the correlations of observation-level random effects. These correlation matrices can also provide evidence for hypothesised energetic tradeoffs between traits, e.g. a negative correlation between residuals of parity and birth weight (see Table 2), although these tradeoffs would be caused by an unobserved variable rather than by parity per se.

Common priors were placed on all model parameters to conservatively regularise our estimates and thus enhance the robustness of our inferences (Gelman et al., Reference Gelman, Carlin, Stern, Dunson, Vehtari and Rubin2013; McElreath, Reference McElreath2020). In particular, we specified

$$\eqalign{& \mu _k, \;\mu _0, \;\beta \approx {\rm Normal}( {0, \;1} ) \cr & \varepsilon \;, \;\mu _{{\rm mother}}, \;\mu _{{\rm residence}}\approx {\rm Exponential}( 3 ) \cr & {\bf R}\;\approx \;{\rm LKJcorr}( 2 ) } $$

Table 2 lists all predictions and specifies which model parameters constitute direct tests of these predictions. We further provide expected associations between control variables and outcomes that do not constitute direct tests of our predictions.

To assess support for predictions, we use fully Bayesian, quantitative inference (Gelman et al., Reference Gelman, Carlin, Stern, Dunson, Vehtari and Rubin2013; McElreath, Reference McElreath2020; McShane, Gal, Gelman, Robert, & Tackett, Reference McShane, Gal, Gelman, Robert and Tackett2019). For each estimated parameter that tests a prediction, we plot the posterior distribution and report the proportion of the posterior that supports the prediction (e.g. denoted as P _<0 for predicted negative associations), thus directly quantifying the degree of support for the prediction. Note that in contrast to frequentist p-values, which indicate the probability of observing the data under the null hypothesis, p(data|H₀), the reported posterior probabilities estimate the support of hypothesised positive or negative associations given the data, i.e. p(H₁|data). In short, larger values of P _<0 etc. indicate greater support for our predictions, which allows readers to quantitatively evaluate the results rather than relying on dichotomous inference.

Data issues