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

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.

*m1*: No sibling comparison

Here, we ignore the pedigree structure of the data to see whether it matters for the estimation of the paternal age effect.

```
model_filename = make_path("m1_children_linear_noranef")
summarise_model()
m1 = model
```

*m2*: Sibling comparison, no paternal age effect

Here, we compared siblings by including a random intercept for the family, but we modelled no effect for paternal age differences among siblings.

```
model_filename = make_path("m2_children_linear_no_diff")
summarise_model()
m2 = model
```

*m3*: Sibling comparison, linear paternal age effect

Here, we compared siblings by including a random intercept for the family, and we modelled a linear effect for paternal age differences among siblings.

```
model_filename = make_path("m3_children_linear")
summarise_model()
m3 = model
```

*m4*: Sibling comparison, nonlinear paternal age effect

Here, we compared siblings by including a random intercept for the family, and we modelled a possibly nonlinear effect for paternal age differences among siblings.

```
model_filename = make_path("m4_children_nonlinear")
summarise_model()
m4 = model
```

## Model comparison

We compare the four models using an approximative leave-one-out cross-validation information criterion as implemented in `brms`

and `loo`

and the Watanabe-Akaike information criterion.

```
model_comparison = LOO(m1, m2, m3, m4, pointwise = FALSE)
print(model_comparison)
plot(model_comparison$m1)
plot(model_comparison$m2)
plot(model_comparison$m3)
plot(model_comparison$m4)
WAIC(m1, m2, m3, m4, compare = TRUE, pointwise = FALSE)
```