```
source("0__helpers.R")
library(brms)
opts_chunk$set(warning=FALSE, cache=F,cache.lazy=F,tidy=FALSE,autodep=TRUE,dev=c('png','pdf'),fig.width=20,fig.height=12.5,out.width='1440px',out.height='900px')
```

Here we show the effect size estimates for paternal age in different robustness analyses.

```
if(!file.exists("coefs/robustness_comparison.rdata")) {
main_mods = c("coefs/krmh/m3_children_linear.rds",
"coefs/rpqa/m3_children_linear.rds",
"coefs/ddb/m3_children_linear.rds",
"coefs/swed/m3_children_linear.rds")
paths = c(main_mods,
list.files("coefs", full.names = TRUE, pattern = "^r\\d.+rds$", recursive = T))
filenames = c(main_mods,
list.files("coefs", full.names = TRUE, pattern = "^r\\d.+rds$", recursive = T))
i=1
effect_estimates = data.frame()
for (i in seq_along(paths)) {
filename = filenames[i]
try({
model = readRDS(paths[i])
if (class(model) == "brmsfit") {
chg = paternal_age_10y_effect(model)[3,]
chg$model = filename
chg$robustness_analysis = as.numeric(str_match(filename, "/r(\\d+)")[,2])
if (is.na(chg$robustness_analysis)) {
chg$robustness_analysis = 0
}
chg$population = str_match(paths[i], "\\w+/(\\w+)/")[,2]
effect_estimates = rbind(chg, effect_estimates)
}
rm(model)
}, silent = T)
}
save(effect_estimates, file = "coefs/robustness_comparison.rdata")
} else {
load("coefs/robustness_comparison.rdata")
}
effect_estimates$median_estimate = as.numeric(effect_estimates$median_estimate)
effect_estimates = effect_estimates %>% arrange(robustness_analysis)
pops = c("krmh", "rpqa", "ddb", "swed")
effect_estimates$upper95 = effect_estimates$lower95 = NA
for(i in seq_along(pops)) {
pop = pops[i]
effs = effect_estimates[which(effect_estimates$population == pop), ]
effs$lower95 = effs$median_estimate[1]
effs$upper95 = effs$median_estimate
effs$lower95[1] = as.numeric(str_match(effs$ci_95[1], "\\[(-?[0-9.]+);")[,2])
effs$upper95[1] = as.numeric(str_match(effs$ci_95[1], ";(-?[0-9.]+)]")[,2])
effect_estimates[effect_estimates$population == pop, ] = effs
}
pops = c("krmh"='Krummhörn', "rpqa" = 'Québec', "ddb" = 'Historical Sweden', "swed" = '20th-century Sweden')
effect_estimates$population = pops[effect_estimates$population]
```

**Legend**: In reference to *m3*, the main reported model, the robustness models were implemented as follows: *r1* relaxed exclusion criteria (not in 20th-century Sweden), *r2* had only birth cohort as a covariate, *r3* adjusted for birth order as a continuous variable, *r4* adjusted for number of dependent siblings instead of birth order, *r5* interacted birth order with number of siblings, *r6* did not adjust for birth order, *r7* adjusted only for parental loss in the first 5 years, *r8* adjusted for being the first-/last-born adult son, *r9* adjusted for a continuous nonlinear thin-splate spline for birth year instead of 5-year bins, *r10* added a group-level slope for paternal age, *r11* included separate group-level effects for each parent instead of one per marriage, *r12* added a moderation by anchor sex, *r13* adjusted for paternal age at first birth, *r14* compared a model with linear group fixed effects, *r15* added a moderator by region and group-level effects by church parish (not in 20th-century Sweden), *r16* was restricted to Skellefteå (only in historical Sweden), *r17* tested whether hypothetical cases of Down’s syndrome could explain the effects, *r18* reversed hurdle Poisson and Poisson distribution for the respective populations, *r19* used a normal distribution, *r20* did not adjust for maternal age, *r21* adjusted for maternal age as a continuous variable, *r22* relaxed exclusion criteria and included 30 more years of birth cohorts, allowing for more potential censoring, *r23* used Student’s t distributions for population-level priors and half-Cauchy priors for the family variance component, *r24* used noninformative priors, which should lead to results comparable with maximum likelihood, *r25* controlled for migration status (not in 20th-century Sweden), *r26* separate parental age contributions (only in 20th-century Sweden).

```
robustness_comparison = ggplot(effect_estimates, aes(x = factor(robustness_analysis), y = median_estimate, ymin = lower95, ymax = upper95)) +
geom_hline(aes(yintercept = ifelse(robustness_analysis == 0, median_estimate, NA)), linetype = 'dashed') +
geom_linerange(aes(linetype = robustness_analysis == 0)) +
geom_text(aes(label = ifelse(robustness_analysis == 0, NA, robustness_analysis), group = effect), vjust = -0.3) +
scale_linetype_manual(values = c("FALSE" = 'dashed', "TRUE" = 'solid'), guide = F) +
facet_wrap(~ population, scales = "free", nrow = 2, ncol = 2) +
scale_x_discrete("", expand = c(0.07,0.07)) +
scale_y_continuous("Fewer children per 10 years of paternal age (in % change)") +
theme(axis.ticks.y = element_blank(), axis.text.y = element_blank()) +
coord_flip()
robustness_comparison
```

Here we use the same axes to make effect sizes comparable across populations.

```
robustness_comparison_same_axes = ggplot(effect_estimates, aes(x = factor(robustness_analysis), y = median_estimate, ymin = lower95, ymax = upper95)) +
geom_hline(aes(yintercept = ifelse(robustness_analysis == 0, median_estimate, NA)), linetype = 'dashed') +
geom_linerange(aes(linetype = robustness_analysis == 0)) +
geom_text(aes(label = ifelse(robustness_analysis == 0, NA, robustness_analysis), group = effect), vjust = -0.3) +
scale_linetype_manual(values = c("FALSE" = 'dashed', "TRUE" = 'solid'), guide = F) +
facet_wrap(~ population, nrow = 2, ncol = 2) +
scale_x_discrete("", expand = c(0.07,0.07)) +
scale_y_continuous("Fewer children per 10 years of paternal age (in % change)") +
theme(axis.ticks.y = element_blank(), axis.text.y = element_blank()) +
coord_flip()
robustness_comparison_same_axes
```

In this interactive plot, move your mouse below the number to make a tooltip with the name of the robustness model appear. Clicking will bring up details for this robustness analysis in this population.

```
library(rbokeh)
effect_estimates = effect_estimates %>%
mutate(
r_label = str_match(model, "/(r\\d+_.+?)\\.rds$")[,2],
r_label = if_else(is.na(r_label),"m3 model with 95% CI", r_label),
r = paste0(r_label, " ", median_estimate),
lty = if_else(robustness_analysis == 0, 1, 2),
link_to_analysis = if_else(robustness_analysis == 0,
str_replace(
str_replace(
str_replace(model, "coefs/", "2_"),
"\\.rds", ""),
"/m", "_main_models.html#m"),
str_replace(
str_replace(
str_replace(model, "coefs/", "5_"),
"\\.rds", ""),
"/r", "_robustness.html#r"))
)
bokeh_coef_plot = function(., population) {
link_to_analysis = .[["link_to_analysis"]]
figure(., title = population, tools = NULL, toolbar_location = NULL) %>%
ly_points(median_estimate, robustness_analysis,
color = 'black',
alpha = 0.05,
glyph = 'cross',
hover = "<strong>@r_label</strong> @median_estimate",
url = link_to_analysis) %>%
ly_text(median_estimate, robustness_analysis,
text = robustness_analysis,
font_size = "10pt"
) %>%
ly_segments(y0 = robustness_analysis,
x0 = lower95,
y1 = robustness_analysis,
x1 = upper95,
type = 2
) %>%
ly_segments(x0 = ifelse(robustness_analysis == 0, median_estimate, NA),
x1 = ifelse(robustness_analysis == 0, median_estimate, NA),
y0 = 0,
y1 = 26,
type = 2) %>%
y_axis(label = "", visible = F) %>%
x_axis(label = "Fewer children per 10 years of paternal age (in % change)")
}
grid_plot(list(
list(
"20th-century Sweden" = effect_estimates %>% filter(population == "20th-century Sweden") %>%
bokeh_coef_plot("20th-century Sweden"),
"Historical Sweden" = effect_estimates %>% filter(population == "Historical Sweden") %>%
bokeh_coef_plot("Historical Sweden")
),
list(
"Krummhörn" = effect_estimates %>% filter(population == "Krummhörn") %>%
bokeh_coef_plot("Krummhörn"),
"Québec" = effect_estimates %>% filter(population == "Québec") %>%
bokeh_coef_plot("Québec")
)
),
same_axes = c(F, T),width = 1000,height=800
)
```

In this comparison table, the model names are clickable and lead to the model details. The list is dynamic and can be subset and ordered by population and model.

```
effect_estimates = effect_estimates %>% mutate(
link = paste0("<a title='View model details' href='",link_to_analysis,"'>",r_label,"</a>"),
population = factor(population)
)
saveRDS(effect_estimates, file = "coefs/robustness_comparison_table.rdata")
DT::datatable(effect_estimates %>% select(population, link, median_estimate, ci_95),
filter = 'top', escape = c(1,3,4), options = list(
pageLength = 50, autoWidth = TRUE
), rownames = F)
```

```
ggsave(plot = robustness_comparison, "../paternal_age_fitness/library/robustness_comparison.png", width = 17.8, height = 17.8*0.625)
ggsave(plot = robustness_comparison, "../paternal_age_fitness/library/robustness_comparison.pdf", width = 17.8, height = 17.8*0.625, scale = 3, device = pdf, units = "cm", dpi = 600)
```