Selective Episodes

Model description

All of the following models have the following in common:

Estimation

We fit all models using brms v. 1.2.0, a Bayesian regression analysis statistical package. brms uses Stan, a probabilistic programming langugage to fit models using Hamiltonian Monte Carlo.

Covariates

As in our main models we adjust for average paternal age in the family, birth cohort (birth years in five equally large bins), for male sex, for age at paternal and maternal loss (0-1, 2-5, 6-10, …, 41-45, 45+, unknown), for maternal age (bins of 14-20, 20-35 and 35-50), for the number of siblings, for the number of older siblings (0-5, 5+) and for being last born.

Model stratification

We added random intercepts for each family (father-mother dyad). We then controlled for the average paternal age in the family. Hence, the paternal age effects in the plot are split into those between families and those within families or between siblings. We are interested in the effect of paternal age between siblings, as this effect cannot be explained by e.g. genetic propensities of the father to reproduce late.

e1: Survival to first year

Here, we predict the probability that the anchor survives the first year of life. All children born to this father are compared, if their death date is known or their survival can be inferred (from later marriage or children).

model_filename = make_path("e1_survive1y")
summarise_model()
e1 = model

e2: Probability of surviving the first 15 years of life

Here, we predict the probability that the anchor survives the first fifteen of life. All children born to this father who lived at least one year are compared, if their death date is known or their survival can be inferred (from later marriage or children).

model_filename = make_path("e2_surviveR")
summarise_model()
e2 = model

e3: Probability of ever marrying

Here, we predict the probability that the anchor ever marries. All anchors who reached reproductive age (15) are included.

model_filename = make_path("e3_ever_married")
summarise_model()
e3 = model

e4: Number of children

Here, we predict the number of children that the anchor had. To separate this effect from previous selective episodes, we include only ever-married anchors and control for their number of spouses (interacted with sex, because men tend to have more additional children from further spouses).

model_filename = make_path("e4_children")
summarise_model()
e4 = model

e5: Probability of ever divorcing

Here, we predict the probability of ever divorcing. To separate this effect from previous selective episodes, we include only ever-married anchors. Divorce data was only analysed in modern Sweden.

model_filename = make_path("e5_divorce")
if(file.exists(model_filename)) {
    cat(summarise_model())
    e5 = model
}

All episodes

Here we show the effect of paternal age for each episode.

episodes = bind_rows(
    e1 = paternal_age_10y_effect(e1),
    e2 = paternal_age_10y_effect(e2),
    e3 = paternal_age_10y_effect(e3),
    e4 = paternal_age_10y_effect(e4), 
    e5 = if( exists("e5") ) { paternal_age_10y_effect(e5) } else { data.frame() },
    .id = "episode"
)
episodes$median_estimate = as.numeric(episodes$median_estimate)
episodes$lower95 = as.numeric(str_match(episodes$ci_95, "\\[(-?[0-9.]+);")[,2])
episodes$upper95 = as.numeric(str_match(episodes$ci_95, ";(-?[0-9.]+)]")[,2])

ggplot(episodes %>% filter(effect == "percentage change") %>% mutate(episode = factor(episode,levels =  rev(unique(episode)))), aes(x = episode, y = median_estimate, ymin = lower95, ymax = upper95)) + 
    geom_hline(yintercept = 0, linetype = 'dashed') +
    geom_pointrange() + 
    xlab("Selective episode") +
    ylab("Percentage change in outcome by paternal age") +
    coord_flip()

saveRDS(episodes, file = make_path("episodes"))