Unpacking mating success and testing Bateman’s principles in a human population

Human marriage systems, characterized by long-term partnerships and extended windows of parental care, differ from the mating systems of pulsed or seasonally breeding non-human animals in which Bateman’s principles were originally tested. These features, paradigmatic of but not unique to humans, complicate the accurate measurement of mating success in evaluating Bateman’s three principles. Here, we unpack the concept of mating success into distinct components: number of partners, number of years partnered, the timing of partnerships, and the quality of partners. Drawing on longitudinal records of marriage and reproduction collected in a natural-fertility East African population over a 20-year period, we test and compare various models of the relationship between mating success and reproductive success (RS), and show that an accurate assessment of male and female reproductive behaviour requires consideration of all major components of mating success. Furthermore, we demonstrate that while Bateman’s third principle holds when mating success is defined in terms of years married, women’s fitness increases whereas men’s fitness decreases from an increase in the number of marriage partners, holding constant the total effective duration of marriages. We discuss these findings in terms of the distinct, sex-specific pathways through which RS can be optimized, and comment on the contribution of this approach to the broader study of sexual selection.

marriage histories for each individual. All individuals are censored at last household visit, and appear only once, as a single row, in the completed database. Therefore, our models do not require adjustments for repeated measurements.
In any given year, there were sometimes inconsistencies in the collected datafor example, conflicting records as to the year in which a wife and ex-husband claimed their relationship began-these conflicts were typically adjudicated at 10 subsequent visits or with relatives or neighbors.
The database constitutes an almost complete census of a single village over the time-period of investigation. If households were missed (due to temporary travel absence, sickness, or a witchcraft case-the latter leading individuals to withdraw from social interactions) data were collected from neighbors and/or 15 in subsequent visits. Village census (excluding members of the migrant Sukuma community who live outside of the main village cluster but on village land) varied from 1,001 (in 1995) to 1,540 (in 2010), with an almost constant sex ratio of about 98 males per 100 females.
RS is measured as the number of offspring reaching their fifth birthday, after 20 which mortality is low. To measure the effects of marriage success on RS, we consider: i) number of spouses ever married (this ranges between 0-5 for both women and men); and, ii) the number of years a focal individual has been married. For men, who can be married to >1 spouse at a time in the case of concurrent polygyny, this value can exceed the number of years for which they 25 have wives. We draw on recent methodological advances using Gaussian random fields [1] to explicitly account for the effects of: iii) the timing of marriages, and iv) the age of each spouse in each year. The dataset and code supporting this article are included the Supplementary Materials and will be maintained at https://github.com/ctross/batemanpimbwe. 30

Paternity ascription and naming customs
It is possible that there might be a fraction of misattributed paternity due to 'hidden copulations' in this population, and how these cases are coded could differentially influence the apparent variance of male and female fitness. Genetic measures of paternity are not available for this population, but repeated 35 household censuses over 20 years reveal remarkable stability to the claimed paternity of children; children are seen as wealth in much of Africa ("wealth in people" [2]) and paternity is rarely denied by men who have sired offspring. Although women might infrequently manipulate attributions of paternity, in such cases there is normally much public discussion on the topic; typically a consensus is reached before the child reaches 3 years of age, often on the basis of physical resemblance. Pimbwe name their children at birth with a "traditional" name, typically the traditional name of a recently deceased ancestor or a nickname depicting some characteristic of the child (e.g., "Tabu," Troublesome). Very soon after birth, the child is also given the name of his/her father; if 45 the mother fails to divulge the name of the child's father, the child takes the name of his mother's father. As the child grows she or he is usually given a religious name (for schooling or other administrative purposes). While religious names can change over the years (in approximately 3% of cases, e.g., Maria to Marietta, Marisella or Meri, or to an entirely new name by choice), Pimbwe 50 names generally remain stable in our demographic records. A child's father's name changed only in 11 of 1,587 cases. In each of these cases, an initially recorded maternal grandfather's name was replaced by the name of the child's presumed father, because the man in question claimed the child, typically marrying the child's mother. In two cases there was no resolution on a father's 55 name, and the child stayed with his maternal grandfather's name.

Marriage and mating
Mating in humans is often regulated through the concept of legitimacy conferred through marriage. We view this as an opportunity for, rather than a challenge to, theorizing sexual selection. Humans are unique in institutionalizing formal 60 norms that constrain mating access. These constraints are variable, both between populations [3], and even within populations, where there may be different kinds of "marriages" [4,5]. Furthermore, marriage can both constrain extra-pair mating and result from extra-pair mating. Accordingly, evolutionary perspectives on sexual selection in humans prefer not to treat marital institutions as entirely exogenous, but rather focus on the norms that pertain to reproduction as coevolving with strategies dictated by individual fitness interests (e.g., [6,7]). In short, while marriage and mating are not equivalent, marriage is a particularly human contract that can shed light on the operation of sexual selection in our species. 70

Fertility rates
Initial analysis of the data suggest that there is no evidence of demographic transition in this highly rural population. There is no evidence of a decline in the rate at which women or men (at or over age 30) produce surviving children as a function of birth year (See Table S1)-in fact, the estimated coefficient 75 of year of birth for women is significantly positive. The R 2 coefficient is very small in both males and females, which itself indicates that birth year plays little role in explaining fertility rate. See also the scatter plots in Figure S1.
While simple analyses such as this can be biased by differential survivorship of individuals contingent on their fertility, and differential reporting of births 80 contingent on respondents' age, the lack of any evidence of a negative association between birth year and fertility suggests that demographic transition in Mpimbwe is minimal.
Moreover, and quite contrary to demographic transition, we find that the frequency of polygyny is actually increasing in recent years, suggesting that 85 males are still aiming to increase reproductive output by acquiring multiple wives. See Figure S2.
Accordingly, we refer to the study population as being largely characterized by natural fertility. [

Scatter plots and correlations of the data
To gain a sense of the data distributions, we plot bivariate scatters and correlation coefficients for all individual-level data used in the analysis. Figures 95 S3-S6 present plots for various subsets of data: 1) those ever married, 2) those at least 45 years of age, 3) those who have married more than once, and 4) the full dataset.
[ Figure S3  Note, however, that high-dimensional relationships can be difficult to see in raw scatters. Most variables included here are influenced by age. Given this consideration, we suggest that Fig. S4, limited to individuals of 45 years and 105 over, is most informative.

Notes about the statistical analyses 2.1 Elasticities
In the main text, we use elasticity parameters, rather than standard slope coefficients, to estimate the relationship between an input variable (marriage 110 success) and an output variable (reproductive success). This is required to ensure a properly specified statistical model that reflects the constraints on the plausible values of the data.
An elasticity parameter, ω, for example, represents the percent change in one variable, y, as a function of the percent change in another, x. It provides 115 an estimation of the relationship between variables that allows for diminishing marginal returns. It is defined as: A slope parameter, β, by contrast, assumes that there are never diminishing marginal returns between variables: Notice, that in a standard linear regression model like: the right-hand side implies that RS can be real-valued or even negative. Also, the coefficients controlling the effects of age, spouse number, and marriage years are constant for all values of the inputs. These are assumptions that we know are discordant with the data generating process a priori. There is 125 both theoretical and empirical [8] evidence that there are typically diminishing marginal returns to each of these input variables in humans-for example, we don't expect the effect of a marginal decade of age on RS to be the same for a 20 year old woman as for an 80 year old woman. Therefore, we do not use a model which assumes this to be the case, because doing so can easily lead to incorrect 130 estimates of the values we hope to measure. The main model we use-essentially a Cobb-Douglas production function-is standard in economics for modeling production [9] and reproduction [10] where outcomes are non-negative and diminishing marginal returns to inputs are possible. While theoretical models of reproduction often consider the linear relationship between variables, such approximations are normally justified in light of simplifying assumptions that don't hold in empirical cases. For data on human reproduction, where diminishing marginal returns are commonly found, accurate estimation of the relationship between variables suggests the use of elasticities, so that the slope ∂y ∂x depends upon the magnitude of the independent variable 140 being evaluated. We note here that when we consider the AIC of the basic linear model: and the AIC of a simple model that instead uses elasticities: for males, the linear model has an AIC of 218.91 points more than the model which uses elasticities. For females, the difference is starker at 251.86. In other 145 words, by using the assumption of a constant slope, we severely reduce our predictive accuracy relative to a better-specified generative model. Our use of elasticities, however, does not mean that our methods cannot speak to the simple slope of RS on marriage or mating success. Should a slope estimate at a given point be desired, it is recoverable by taking the partial 150 derivative of the reproduction function with respect to a given input. This involves nothing more than a transformation of our parameter estimates.
For example, using a simple function for reproduction: r ∼ αx γ y ω , where r is reproductive output, x > 0 is age, y > 0 is years married, α > 0 is an intercept and ω and γ are elasticities, the slope of years married on reproduction, for 155 individuals of a given age class, is: As an example of this transformation, we present the corresponding male and female values of the implied slope coefficients from our model of marital years described in the main paper ( Figure S7). We also present the corresponding plot from a robustness check using only individuals of age 45 or older, dropping 160 the two stage model ( Figure S8, and see discussion in section 4.2)-finding essentially identical inferences using either methodological approach. Here we find that because of diminishing marginal returns to years married for females, and approximately constant marginal returns to years married for males, the slope of reproductive success on marriage success starts out approximately equal 165 for males and females of a given age, but a sex difference becomes pronounced as number of years married increases.

Main model
To accommodate the fact that some of our empirical data show strong signs of non-linearity (specifically in the form of diminishing marginal fitness returns to marriage success for females), we measure the Bateman gradient using an elasticity parameter, which indicates the percent change in RS with respect to the percent change in marriage success. This is a standard approach in economics to modeling reproduction and marriage [8], where diminishing marginal fitness returns to marriage success occur and would lead to inaccurate estimates of standard slope coefficients. Standard slope coefficients, however, can be calculated from our parameter estimates assuming a fixed level of inputs (as demonstrated above in section 2.1). We note that our reproductive success data 180 also show strong signs of zero-inflation and over-dispersion ( Figure S9). This leads standard linear regression models to fail quality checks-as we show at the end of section 4.3.

Main equations
To appropriately model our zero-inflated reproductive success data (see Figure   185 S9), we use a two-stage modeling framework [1] and a modeling structure that lets us decompose the paths through which spouse number affects RS (see Figure S10).
[ Figure S9 where E [i] is the exposure time to the possibility of reproduction-i.e., years lived in the interval between age 11 and death/censoring-and S(i) is a function returning the sex of individual i. See the results of this submodel in section 3.2. In cases where M [i] = 1, we fit our main model linking marriage success and reproductive success. Specifically, we model the reproductive success, R, of individual, i, using a negative binomial outcome distribution [8]: where the term µ [i] B [S(i)] defines the shape parameter of a Gamma distribution, 210 and B [S(i)] defines the inverse scale parameter. This is equivalent to using a Gamma-Poisson mixture model, which has been recommended for modeling over-dispersed fertility-related outcomes-i.e., where the variance exceeds the mean-which are commonly found in polygynyous societies [11]. We can then define a model of the mean of RS, µ [i] , using a standard log link function: where the new variables in the regression are: spouse number, N [i] , and effective marital years, Y [i] -this variable is constructed as described below. Note that the regression parameters, β, are unique to the sex of individual i, and that the elasticity on Y [i] is constrained to be positive, to ensure identification given the reflection invariance that could otherwise arise from the multiplication of 220 parameters.
The variable Y [i] (number of "effective marital years") is obtained for an individual by calculating a weighted sum of the number of years in which he or she has been married to each of his or her spouses. More formally, Y [i] is the sum: Value to a focal individual of a marital year at a given age Value to a focal individual of marriage to spouse of a given age , if n > 0 (10) where the first factor gives the estimated value to the focal individual of having a spouse at a given age (timing weight) and the second factor gives the estimated value to the focal individual of having a spouse of a given age (spousal quality weight). We note that effective marital years and spouse number in the main model are not strongly correlated (for men ρ = 0.3; for women ρ = 0.05).
In the first factor, θ is an intercept and ν is an age-specific random effect. In the second factor, φ is an intercept and ψ is an age-specific random effect. The function A(n, i, e) gives the age index of the n-th spouse of individual i, in individual i's e-th year of life. Eq. 10 imposes the following: 1) if an individual is not married in a given year, the number of quality-weighted spouse 235 years is not incremented in that year; 2) the maximum possible value of a quality-weighted spouse year is 1, and is approached when an individual of a given sex is married during their peak reproductive period to a spouse also at their peak reproductive period; and 3) if an individual is married, the number of quality-weighted spouse years is incremented, but the weights on spousal 240 quality and marriage timing may reduce this value to something less than 1, but greater than 0.
To test our predictions concerning how one's definition of mating success affects the apparent Bateman gradient, we can refit the model, varying our definition of Y [i] . For example, we can replace logistic(θ [S(i)] + ν [e,S(i)] ) and/or with the constant value 1, and investigate how the relationship between reproductive success, R, and marriage success, Y , changes as we account for the timing of marriages and the quality of spouses, alone and in combination. All random effects in this model are estimated using Gaussian random fields (see details in section 2.2.2).

250
In cases where M [i] = 0, we fit a separate model linking marriage success and age (also referred to as exposure time). Because marriage does not exactly correspond to mating, we have to make a special allowance that never-married individuals can reproduce. So, as in the main model, we model the reproductive success, R, of individual, i, again using a negative binomial outcome distribution 255 [8].
However, we now define an independent, simplified model of mean RS in this class of individuals: where the only predictor is E [i] , the age of individual i minus 11 years-i.e., the number of years in which individual i has been old enough to theoretically 260 reproduce. See the results of this submodel in section 3.3.

Random effects
We model the random effects vectors ν and ψ using Gaussian random field submodels, which allow for the effects marriage timing and spousal quality to take on arbitrary functional forms, while still partially sharing information across neighboring parameters. More specifically, we model: where η ν , η ψ ∈ (0, ∞) serve to scale variance, and L ν and L ψ are factors given by the Cholesky decomposition of the correlation matrices ρ ν and ρ ψ . The correlation matrices are defined using a distance decay function: where κ ν , κ ψ ∈ (0, 1) control the maximum correlations, τ ν , τ ψ ∈ (0, ∞) control the decay rates, and C is a constant which normalizes the maximum distance between age categories to 1. The use of a multivariate normal distribution to generate correlated random 275 effects carries several benefits: 1) it allows for partially-pooled estimation of marriage timing and spousal quality weights, 2) it imposes no a priori functional form (e.g., linear, quadratic, etc.) on these random effects, and 3) it allows for a reduction in the effective parameter complexity of the model by reducing the extent to which neighboring random effects parameters can vary independently.

Priors
We define weak priors over the top-level parameters. These priors reflect the fact that we have little information about the values of model parameters a priori, and allow the data and likelihood to dominate the posterior. Using subscripts of m and f for the male and female parameters, in the main regression equations, 285 we use: The parameters on the random effects have weakly informative priors: The correlation parameters have Beta priors: The decay parameters have half-normal priors: The variance parameters have half-normal priors:

Software implementation
Models are fit using Hamiltonian Monte Carlo [12] and the Stan 2.16.0 C++ library [13]. Markov Chain Conte Carlo methods are needed to accurately sample to posterior distribution of this complex, multi-level model. Models are 310 implemented using an R [14] workflow, through the rstan interface. All code used in modeling is included in the Supplementary Materials folder. The code will also be maintained on https://github.com/ctross/batemanpimbwe.

Model fit
315 Figure S11 gives the traceplots for each of the models fit in the main paper (ordered as in Table 1 of the main text). We see that in each model both chains appear to have converged to the same posterior region and to have mixed thoroughly.

Probability of remaining never-married
The first submodel of the main analysis investigates the sex-specific probability of remaining never-married as a function of age. We find that males have higher probability of remaining never-married until later in life than females. However, both males and females reach high levels of marriage probability by 325 their mid-thirties ( Figure S12).

RS of never-married individuals
From theory, we would expect reproductive success to follow a zero-inflated negative binomial generative model, where, for example, RS=0 in all cases 330 where an individual has never acquired even a single mate, and RS ∈ N 0 if an individual has ever acquired a mate. Because marriage is not equivalent to mating, however, we test this expectation, using a model that allows for never-married individuals to reproduce. Specifically, we model sex-and agespecific reproductive success among never-married individuals using a standard 335 negative binomial model. We find that never-married females have reliably higher average levels of reproductive success than never-married males ( Figure  S13). For both males and females, never-married individuals-especially the young-are predicted to have near zero RS. The predicted mean number of offspring does not reach the value of one for females until age thirty, and 340 for males it does not reach the value of one until after age 55 ( Figure S13)essentially confirming our above-stated expectation, and supporting our use of a zero-inflated negative binomial generative model. Figure S14 provides direct density estimates of RS among never-married individuals, and leads to exactly the same inference-although marriage is not equivalent to mating, our RS 345 data differ by only a small margin from those that would be produced by a zero-inflated negative binomial model; nearly all density for RS in never-married males and females is concentrated on the value RS=0.
[ Figure S13 about here.] [ Figure S14 about here.] 350 We note that estimates of male reproductive success, and particularly the reproductive success of never-married males, may be underestimated if such males have sired offspring outside of the village where this study was conducted. Children born in non-marital relationships (either prior to marriage or with an extra-pair individual) are more likely to appear in a mother's than a father's 355 record, for two reasons: first, a young child has a higher likelihood of living with its mother than its father, and therefore being observed as a household resident; second, on account of the above bias and the fact that some nonmarital relations occur with individuals outside the focal village, illegitimate children born to the men of the focal village with outside women are more 360 likely to be lost from the record. While such issues are unavoidable properties of studies of human fertility, the difference in average reproductive success for reproductively complete (age 45 or older) females and males in our sample is only -0.18, which is small in proportion to the average number of kids, 6.09.
Our main model-a zero-inflated negative binomial model-serves mostly 365 to automatically adjust the sample to exclude individuals who are too young to reproduce. As a robustness check, we also analyze the relationship between RS and marriage success using only those individuals of age 45 or over. This eliminates the zero-inflation problem, and negates the need to use a two-stage model, or to restrict inclusion into the sample to only those individuals with at

Timing and quality curves
In order to interpret our measures of marriage success in the main text, we present the estimated sex-specific weights for the value of having a spouse at a 375 given age (marriage timing) and the value of having a spouse of a given age (spousal quality). With regards to the marriage timing weights, we note two key empirical results: 1) marriages occurring when individuals of either sex are younger (< 40 years old) have larger implications for reproductive success than marriages which occur at older ages; and, 2) the reproductive value of marriages 380 for males diminishes less with age than for females, as in indicated by the male curve ( Figure S15.a) bottoming after age 60, while the same curve for females ( Figure S15.b) bottoms earlier. These curves effectively show sex differences in age-specific reproduction for this population, reflecting, for example, the effects of menopause on the reproduction of females. These curves were estimated endogenously in the model, not hard-coded a priori ; the fact that they reflect known patterns in humans give us confidence in the appropriateness of our modeling technique.
[ Figure S15 about here.] We note from the spousal quality curves that: 1) for males, acquisition of spouses who are younger (< 40 years old) has a larger impact on reproductive success than acquisition of older spouses, as indicated by the high density near value equals 1 for younger spouses, and high density near value equals 0 for older spouses ( Figure S16.a); and, 2) for females also, acquisition of younger spouses (< 60 years old) has a larger impact on reproductive success than acquisition of 395 older spouses ( Figure S16.b).
[ Figure S16 about here.] The similar patterning of the estimated curves between Figures S15 and S16 are notable-and expected-given that the spousal quality curve for a male, for example, should theoretically follow a similar form as the marriage timing 400 curve for a female.

Bateman's second and third principles (Full presentation of results)
In the main text, we commented on spouse number and the fully weighted model. Considering only spouse number, we found no evidence of a difference in 405 the opportunity for sexual selection between men and women of age 45 or older. Nor, in the corresponding regression model using only spouse number in the full sample of individuals, was there evidence of a reliably positive relationship between spouse number and RS for either males or females. Likewise, we failed to find any indication of sex differences in the effects of spouse number on RS.

410
Defining marriage success using marital years rather than spouse number, however, we find marginally reliable evidence of a difference in variance in marriage success between males and females of age 45 or older: log(I sm /I sf ) = 0.26 (0.00, 0.53) ( Table 2 in main text; n m = 171, n f = 176). Also, we now find reliably positive effects of marriage success on RS, for males, β m = 415 0.88 (0.77, 1.00), and females, β f = 0.46 (0.37, 0.55) ( Table 1 in main text); n m = 447, n f = 627), and evidence that the effect is stronger in males than females, β m − β f = 0.42 (0.27, 0.57) ( Table 2 in main text), as is predicted under the Bateman model.
Further inclusion of weighting functions on martial years leads to better predictive accuracy, as measured by WAIC, and yields stronger elasticity estimates, as we describe in the main text.

Age threshold of 55 years for variance measures
Taking an age cut-off at age 45 might be considered too early as an indicator of 425 overall fitness, especially for men. Accordingly, we calculate the opportunity for selection, and the opportunity for sexual selection, for a smaller sample of individuals older than 55 years of age.

Opportunity for selection measures
Sex-specific levels of variation in RS for men (n m = 97) and women (n f = 104) that this reflects a more recent increase in the frequency of concurrent polygyny (see, for example, Figure S2).

Opportunity for sexual selection measures
Here, using the same older sample, we examine the opportunity for sexual selection, specifically for spousal numbers and (weighted) marital years. Like 440 the younger sample (≥45 years), the older sample shows no sex difference in variation in spouse number. However, unlike in the ≥ age 45 sample (where we see a marginally reliable effect), the older sample also shows no sex difference in variation in marital years. But, as we show in Table S2, we see that as long as we include weights on the number of marital years, we detect increased 445 male-to-female variance in mating success even in the subset of individuals over the age of 55-n m = 97 and n f = 104-replicating the result from the main analysis. [

450
In the main text, we attempt to fully mine the information content of the database by using the full sample of individuals (those over age 11). This gives us a sample size of 447 males and 627 females, but introduces some added complexity in the modeling, by requiring a two-stage model design to deal with zero-inflation [1]. Here, we present the results of a comparable analysis using only adults of age 45 and over. In this case, our sample size is reduced to 171 males and 176 females, but the issue of zero inflation is eliminated. In this case, we can simplify the statistics, omit the two-stage model, and avoid using the measure 'ever married' to permit inclusion into the main regression analysis. In this case, there are only four males and a single female with zero 460 years married-possibly due to missing data or simply concealment of specific relationships by respondents-we therefore treat the marriage success measures of these individuals as if they were missing data, constrained to fall within the observed population-level distribution.
[ The results are presented in Tables S3 and S4, which correspond directly to Tables 1 and 2 of the main text. The qualitative results of the main text are largely unchanged. We do note that in this subset of data, the negative effect of spouse number on RS in the full model remains for males, but the 470 positive effect for females disappears. The magnitude of the sex difference, however, remains mostly constant and this contrast remains reliable. The fact that the ≥45 sample shows a slightly different patterning than the full sample reported in the main paper suggests the possibility that the advantage of multiple spouses to Pimbwe women has become more prominent among younger 475 and more contemporary women than it was in the past.

Simple linear regression models
We also replicate the main analysis using simpler tools. First, we use the same subset of individuals-those with at least one year of marriage-as we used in the main paper. The results are presented in Table S5. As we found in 480 our main analysis, in the univariate models spouse number is not a reliable predictor of RS in men or women. Moreover, in the univariate models, marital years are a reliable predictor of RS in men and women-and the effect is roughly twice as large among men, just as we found in the main analysis using elasticities. Assuming that we use the methods described in section 2.1 to 485 calculate slopes from elasticities-and evaluate these calculations for individuals of age 55 with 40 years of marriage-then the slopes derived from the elasticites estimated in the main text are: 0.163 (0.135, 0.190) for males, and 0.075 (0.059, 0.095) for females. These values are in near numerical correspondence with the slopes estimated directly in the linear models. Finally, in the model with both predictors, we find that the effect of spouse number on RS holding constant marital years is smaller in males than in females, just as we found in the main analysis. There is a difference on this point between the OLS and Bayesian models in that the negative effect of spouse number on male RS in model (3) does not reach significance in the OLS model. The difference between this 495 OLS model and the Bayesian model presented in the main paper is that the Bayesian model includes timing and spousal quality weights on the number of marital years, increasing the predictive accuracy of the model. The OLS model by contrast uses raw marital years.
[ Next, we use the same subset of individuals as in section 4.2-those individuals of age 45 or over. The results of these models are presented in Table S6. As we found in our Bayesian analyses, in the univariate models spouse number is not a reliable predictor of RS in men or women. Moreover, in the univariate models, marital years are a reliable predictor of RS in men and women-and 505 the effect is a little more than twice as large among men, just as we found in the Bayesian analysis using elasticities. Finally, in the model with both predictors, we find that males suffer from increasing spouse number, holding constant years married, while females benefit. There is a difference on this point between the OLS and Bayesian models in that the negative effect of spouse number on male 510 RS in model (3) does not reach significance in the OLS model. As noted before, the Bayesian model includes timing and spousal quality weights on the number of marital years, increasing the predictive accuracy of the model. The OLS model by contrast uses raw marital years.
[ Finally, we fit the OLS models using the whole dataset, including children as young as 11 years old with RS = 0 and martial years = 0. As we found in our Bayesian analyses, in the univariate models, marital years are a reliable positive predictor of RS in men and women-and the effect is a little more than twice as large among men. In the univariate models here, however, spouse 520 number is now a reliable positive predictor of RS in both women and men. The effect is even slightly larger in males than females. As we expand on later, this result is likely driven by the model failing to correctly deal with zero-inflation. In the model with both predictors, we find that males still benefit more from increasing marital years, holding constant spouse number, than do females.

525
Females continue to benefit more than males from increasing spouse number holding constant years married, but, in this case, both effects are on the positive side of zero. Again, the Bayesian model includes timing and spousal quality weights on the number of marital years, increasing the predictive accuracy of the model. The OLS model by contrast uses raw marital years.

530
This last model would lead us to somewhat different conclusions than all other models presented in the main text and this supplement. In this set of models, for both sexes, the strong effects of spouse number on RS primarily reflect the marginal change in RS between those-frequently very young-individuals with zero spouses and those with at least single spouse, an effect seen in many 535 other populations, e.g., [15,16,17]. In the model presented in the main paper, however, the effect of spouse number on RS primarily reflects the marginal change in RS between those individuals with a single spouse and those with more than a single spouse.
Note, however, that this last OLS model deals poorly with the zero-inflated 540 outcome data, and the QQ-plots [18] of the residuals suggest that the assumptions of the linear regression analysis are violated-making the results of the model hard to evaluate. In Figure S17, note how the QQ-plots of the residuals for both male and female outcomes show a dramatic divergence of the distribution of model residuals from that which is assumed by the regression model.

545
By contrast, if we use only those individuals age 45 and older, the QQ-plots in Figure S18 demonstrate that the residuals are distributed as assumed by the linear regression model. In unison, these figures suggest that our statistical inferences with OLS are valid in the case where only post-reproductive individuals are included in the analysis, but not in the case where standard OLS tools 550 are applied to the full dataset. We suspect that some of the OLS results differ due to deviations from the assumption of normality.
[    Figure S7: Sex-specific slope estimates of RS on years married calculated from regression parameters in the 'marital years model' using the full dataset as presented in the main text. Estimates, nevertheless, are predicted value for males and females at age 45. Note that the slope depends on years married if the elasticity on years married is not equal to one. The dashed vertical bars show the 90% density intervals of the data for years married for the subset of individuals between ages 44 and 46. The model predictions are best considered only in this range. We see that the male slope estimate is higher than the corresponding female slope estimate over almost the entire data range-a result concordant with our findings as presented in the main paper.  Figure S8: Sex-specific slope estimates of RS on years married calculated from regression parameters in the 'marital years model' as presented in section 4.2 of this supplement. Estimates are predicted values for males and females at age 45. Note that the slope depends on years married if the elasticity on years married is not equal to one. The dashed vertical bars show the 90% density intervals of the data for years married for the subset of individuals between ages 45 and 46. The model predictions are best considered only in this range. We see that the male slope estimate is higher than the corresponding female slope estimate over almost the entire data range-a result concordant with our findings as presented in the main paper.  Figure S9: Histogram of male (red) and female (blue) RS, matched to the counter-factual histogram if the same number of observations were drawn from a Poisson distribution with the same sex-specific mean (dark grey). Note that our data are highly non-Gaussian. We observe a large peak at RS = 0, due mostly to young individuals who have never been married. We deal with this density at RS = 0, by using a variant of a standard zero-inflated model. Note also that RS is over-dispersed relative to a Poisson distribution (this means that the data distribution has higher density in both tail regions than the standard Poisson distribution). We deal with this by using a negative binomial regression framework.