Comparing paternal age effects in our selective episode models

Loading details

source("0__helpers.R")
opts_chunk$set(warning=FALSE, cache=F,cache.lazy=F,tidy=FALSE,autodep=TRUE,dev=c('png','pdf'), fig.width=17.8,fig.height = 17.8*0.625,out.width='1440px',out.height='900px')
 
get_paternalage_effect = function(model, population, var = "paternalage") {
    brms::marginal_effects(readRDS(get_coefficient_path(model, population)), var)
}

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).

Marginal effect plot

episode = "e1_survive1y"
e1 = bind_rows(
    "Krummhörn" = get_paternalage_effect(episode, "krmh"),
    "Québec" = get_paternalage_effect(episode, "rpqa"),
    "20th-century Sweden" =  get_paternalage_effect(episode, "swed"),
    "Historical Sweden" = get_paternalage_effect(episode, "ddb"),
.id = "Population"
) %>% mutate(Population = str_sub(Population,1, -13))

e1_plot = ggplot(e1,  aes(x = paternalage*10, y = Estimate, colour = Population, fill = Population, ymin = lowerCI, ymax = upperCI)) + 
    geom_smooth(stat = 'identity', position = position_identity()) +
    ggtitle("Selective episode - e1") + 
  scale_y_continuous("Probability of survival to first year", limits = c(0,1))  + 
    scale_x_continuous("Paternal age", breaks = c(20,30,40,50,60,70)) + 
    analysis_theme +
    pop_colour + pop_fill +
    theme(legend.position = c(-1,0),
                legend.justification = c(0,0))
e1_plot

Effect size comparison

comp_plot(episode)

##            Population median_estimate         ci_95         ci_80
## 1           Krummhörn           -2.15 [-5.35;-0.21] [-4.03;-0.82]
## 2              Québec           -1.03 [-1.51;-0.67] [-1.32;-0.78]
## 3 20th-century Sweden           -0.05 [-0.06;-0.03] [-0.06;-0.03]
## 4   Historical Sweden           -1.82 [-3.14;-1.08] [-2.63;-1.28]

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).

Marginal effect plot

episode = "e2_surviveR"
e2 = bind_rows(
    "Krummhörn" = get_paternalage_effect(episode, "krmh"),
    "Québec" = get_paternalage_effect(episode, "rpqa"),
    "20th-century Sweden" =  get_paternalage_effect(episode, "swed"),
    "Historical Sweden" = get_paternalage_effect(episode, "ddb"),
.id = "Population"
) %>% mutate(Population = str_sub(Population,1, -13))

e2_plot = ggplot(e2,  aes(x = paternalage*10, y = Estimate, colour = Population, fill = Population, ymin = lowerCI, ymax = upperCI)) + 
    geom_smooth(stat = 'identity', position = position_identity()) +
    ggtitle("Selective episode - e2") + 
  scale_y_continuous("Probability of surviving to age 15", limits = c(0,1))  + 
    scale_x_continuous("Paternal age", breaks = c(20,30,40,50,60,70)) + 
    analysis_theme +
    pop_colour + pop_fill +
    theme(legend.position = c(-1,0),
                legend.justification = c(0,0))
e2_plot

Effect size comparison

comp_plot(episode)

##            Population median_estimate        ci_95        ci_80
## 1           Krummhörn            0.06 [-2.69;2.73] [-1.73;1.75]
## 2              Québec            0.09 [-0.17;0.36] [-0.07;0.26]
## 3 20th-century Sweden            0.03     [0;0.06]  [0.01;0.05]
## 4   Historical Sweden           -0.16 [-0.68;0.27] [-0.47;0.12]

e3: Probability of ever marrying

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

Marginal effect plot

episode = "e3_ever_married"
e3 = bind_rows(
    "Krummhörn" = get_paternalage_effect(episode, "krmh"),
    "Québec" = get_paternalage_effect(episode, "rpqa"),
    "20th-century Sweden" =  get_paternalage_effect(episode, "swed"),
    "Historical Sweden" = get_paternalage_effect(episode, "ddb"),
.id = "Population"
) %>% mutate(Population = str_sub(Population,1, -13))

e3_plot = ggplot(e3,  aes(x = paternalage*10, y = Estimate, colour = Population, fill = Population, ymin = lowerCI, ymax = upperCI)) + 
    geom_smooth(stat = 'identity', position = position_identity()) +
    ggtitle("Selective episode - e3") + 
  scale_y_continuous("Probability of ever marrying", limits = c(0,1))  + 
    scale_x_continuous("Paternal age", breaks = c(20,30,40,50,60,70)) + 
    analysis_theme +
    pop_colour + pop_fill +
    theme(legend.position = c(-1,0),
                legend.justification = c(0,0))
e3_plot

Effect size comparison

comp_plot(episode)

##            Population median_estimate         ci_95          ci_80
## 1           Krummhörn           -5.19 [-14.79;2.01] [-11.08;-0.37]
## 2              Québec           -0.02  [-0.83;0.74]   [-0.55;0.46]
## 3 20th-century Sweden            0.82   [0.63;1.01]    [0.69;0.95]
## 4   Historical Sweden            7.94  [4.47;12.38]   [5.48;10.83]

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).

Marginal effect plot

episode = "e4_children"
e4 = bind_rows(
    "Krummhörn" = get_paternalage_effect(episode, "krmh"),
    "Québec" = get_paternalage_effect(episode, "rpqa"),
    "20th-century Sweden" =  get_paternalage_effect(episode, "swed"),
    "Historical Sweden" = get_paternalage_effect(episode, "ddb"),
.id = "Population"
) %>% mutate(Population = str_sub(Population,1, -13))

e4_plot = ggplot(e4,  aes(x = paternalage*10, y = Estimate, colour = Population, fill = Population, ymin = lowerCI, ymax = upperCI)) + 
    geom_smooth(stat = 'identity', position = position_identity()) +
    ggtitle("Selective episode - e4") + 
  scale_y_continuous("Probability of surviving to age 15", limits = c(0,NA))  + 
    scale_x_continuous("Paternal age", breaks = c(20,30,40,50,60,70)) + 
    analysis_theme +
    pop_colour + pop_fill +
    theme(legend.position = c(-1,0),
                legend.justification = c(0,0))
e4_plot

Effect size comparison

comp_plot(episode)

##            Population median_estimate         ci_95         ci_80
## 1           Krummhörn           15.62  [3.92;29.31]  [7.67;24.68]
## 2              Québec            0.94  [-1.32;3.23]   [-0.5;2.47]
## 3 20th-century Sweden           -3.82 [-4.64;-2.98] [-4.38;-3.28]
## 4   Historical Sweden           -5.36 [-8.92;-1.59] [-7.71;-2.96]

All episodes

filenames = c("m3_children_linear.rds","m3_children_linear.rds","m3_children_linear.rds","m3_children_linear.rds", 
        list.files("coefs", full.names = TRUE, pattern = "^e._.+rds$", recursive = T))
paths = 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", 
                    list.files("coefs", full.names = TRUE, pattern = "^e._.+rds$", recursive = T))

i=1
effect_estimates = data.frame()
models = list()
for (i in seq_along(paths)) {
    filename = filenames[i]
    model = readRDS(paths[i])
    if (class(model) == "brmsfit") {
        models[[filename]] = model
        chg = paternal_age_10y_effect(model)[3,]
        chg$model = filename
        chg$population = str_match(paths[i], "\\w+/(\\w+)/")[,2]
        effect_estimates = rbind(chg, effect_estimates)
    }
}
effect_estimates$median_estimate = as.numeric(effect_estimates$median_estimate)
effect_estimates = effect_estimates %>% arrange(episode)
effect_estimates$episode = str_match(effect_estimates$model, "([em]\\d)_")[,2]
effect_estimates$lower95 = as.numeric(str_match(effect_estimates$ci_95, "\\[(-?[0-9.]+);")[,2])
effect_estimates$upper95 = as.numeric(str_match(effect_estimates$ci_95, ";(-?[0-9.]+)]")[,2])
effect_estimates$lower80 = as.numeric(str_match(effect_estimates$ci_80, "\\[(-?[0-9.]+);")[,2])
effect_estimates$upper80 = as.numeric(str_match(effect_estimates$ci_80, ";(-?[0-9.]+)]")[,2])

pops = c("krmh"='Krummhörn', "rpqa" = 'Québec', "ddb" = 'Historical Sweden', "swed" = '20th-century Sweden')
eps = c("m3"='0. Number of children\n(all)', "e1" = '1. Infant survival', "e2" = '2. Survival to 15', "e3" = '3. Ever married', "e4" = "4. Number of children\n(conditional on ever marrying)")
effect_estimates$Population = pops[effect_estimates$population]
effect_estimates$Population = factor(effect_estimates$Population, levels = pops)
effect_estimates$episode = factor(eps[effect_estimates$episode])
effect_estimates$episode = factor(effect_estimates$episode, levels = rev(levels(effect_estimates$episode)))

episode_comparison = ggplot(effect_estimates %>% filter(!is.na(episode)), aes(x = factor(episode), y = median_estimate, ymin = lower95, ymax = upper95, color = Population)) + 
    geom_linerange(position = position_dodge(width = 0.4)) +
    geom_pointrange(aes(ymin = lower80, ymax = upper80), size = 1, position = position_dodge(width = 0.4)) +
    geom_hline(yintercept = 0, linetype = "dashed") +
    xlab("") +
    pop_colour +
    scale_y_continuous("Paternal age effect") +
    coord_flip()    + 
    theme(legend.position = c(1,1),
    legend.justification = c(1,1))
episode_comparison