All of the following models have the following in common:

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.

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.

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.

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
```

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
```

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
```

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
```

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
}
```

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"))
```