20th century Sweden sensitivity analyses

Loading details

source("0__helpers.R")
opts_chunk$set(warning=TRUE, 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',cache.extra=file.info('swed1.rdata')[, 'mtime'])

make_path = function(file) {
    get_coefficient_path(file, "swed")
} 
# options for each chunk calling knit_child
opts_chunk$set(warning=FALSE, message = FALSE, echo = FALSE)

Analysis description

Data subset

The swed_subset_children.1 dataset is based on the full dataset of all participants where paternal age is known and birth years are from 1947 to 1959. The subset contains 88707 randomly drawn participants from 56679 families.

Model description

All of the following models are the same as our main model m3, except for the noted changes to test robustness.

s1: Mediation via age

Here, we tested whether the effect on reproductive success is mediated by age (mortality). Unfortunately for the cohort we use for computing an effect on number of children, we have both left and right censoring (because the death register was only established in 1962, so some people may have died before then, and because many people born in 1958 are still alive). So we can only include people who died after their 11th-15th birthday, and have to enter an unknown age for those who have not died before their 51st-63rd birthday.

Model summary

Full summary

model_summary = summary(model, use_cache = FALSE, priors = TRUE)
print(model_summary)
##  Family: poisson (log) 
## Formula: children ~ paternalage + age + birth_cohort + male + maternalage.factor + paternalage.mean + paternal_loss + maternal_loss + older_siblings + nr.siblings + last_born + (1 | idParents) 
##    Data: model_data (Number of observations: 117726) 
## Samples: 6 chains, each with iter = 1000; warmup = 500; thin = 1; 
##          total post-warmup samples = 3000
##    WAIC: Not computed
##  
## Priors: 
## b ~ normal(0,5)
## sd ~ student_t(3, 0, 5)
## 
## Group-Level Effects: 
## ~idParents (Number of levels: 75475) 
##               Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
## sd(Intercept)     0.01      0.01        0     0.02        405 1.01
## 
## Population-Level Effects: 
##                        Estimate Est.Error l-95% CI u-95% CI Eff.Sample
## Intercept                  0.40      0.02     0.35     0.44       3000
## paternalage               -0.05      0.01    -0.08    -0.03        840
## age5063                    0.20      0.02     0.16     0.24       3000
## age025                    -2.87      0.14    -3.17    -2.61       3000
## ageunknown                 0.35      0.02     0.31     0.38       3000
## birth_cohort1950M1955     -0.01      0.01    -0.02     0.00       3000
## birth_cohort1955M1960      0.00      0.01    -0.01     0.01       3000
## male1                     -0.05      0.00    -0.06    -0.04       3000
## maternalage.factor1420     0.07      0.01     0.06     0.09       3000
## maternalage.factor3561    -0.01      0.01    -0.03     0.01       3000
## paternalage.mean           0.01      0.01    -0.02     0.04        869
## paternal_loss01            0.17      0.28    -0.41     0.67       3000
## paternal_loss15           -0.05      0.08    -0.20     0.10       3000
## paternal_loss510          -0.06      0.03    -0.13     0.00       3000
## paternal_loss1015         -0.02      0.02    -0.06     0.01       3000
## paternal_loss1520          0.00      0.01    -0.03     0.03       3000
## paternal_loss2025          0.00      0.01    -0.02     0.02       3000
## paternal_loss2530         -0.02      0.01    -0.04     0.00       3000
## paternal_loss3035         -0.01      0.01    -0.02     0.01       3000
## paternal_loss3540         -0.01      0.01    -0.02     0.01       3000
## paternal_loss4045         -0.01      0.01    -0.03     0.00       3000
## paternal_lossunclear      -0.06      0.01    -0.07    -0.05       3000
## maternal_loss01           -0.22      0.23    -0.70     0.20       3000
## maternal_loss15           -0.01      0.10    -0.21     0.20       3000
## maternal_loss510           0.02      0.05    -0.07     0.11       3000
## maternal_loss1015         -0.02      0.03    -0.07     0.04       3000
## maternal_loss1520         -0.01      0.02    -0.05     0.03       3000
## maternal_loss2025          0.02      0.02    -0.01     0.06       3000
## maternal_loss2530          0.00      0.01    -0.03     0.03       3000
## maternal_loss3035          0.00      0.01    -0.03     0.02       3000
## maternal_loss3540          0.00      0.01    -0.02     0.02       3000
## maternal_loss4045         -0.01      0.01    -0.03     0.01       3000
## maternal_lossunclear      -0.02      0.01    -0.03    -0.01       3000
## older_siblings1            0.02      0.01     0.00     0.03       1123
## older_siblings2            0.03      0.01     0.01     0.05        843
## older_siblings3            0.04      0.02     0.01     0.08        989
## older_siblings4            0.01      0.02    -0.03     0.06       1065
## older_siblings5P          -0.04      0.03    -0.10     0.02        813
## nr.siblings                0.04      0.00     0.03     0.04       1075
## last_born1                 0.01      0.00     0.00     0.02       3000
##                        Rhat
## Intercept                 1
## paternalage               1
## age5063                   1
## age025                    1
## ageunknown                1
## birth_cohort1950M1955     1
## birth_cohort1955M1960     1
## male1                     1
## maternalage.factor1420    1
## maternalage.factor3561    1
## paternalage.mean          1
## paternal_loss01           1
## paternal_loss15           1
## paternal_loss510          1
## paternal_loss1015         1
## paternal_loss1520         1
## paternal_loss2025         1
## paternal_loss2530         1
## paternal_loss3035         1
## paternal_loss3540         1
## paternal_loss4045         1
## paternal_lossunclear      1
## maternal_loss01           1
## maternal_loss15           1
## maternal_loss510          1
## maternal_loss1015         1
## maternal_loss1520         1
## maternal_loss2025         1
## maternal_loss2530         1
## maternal_loss3035         1
## maternal_loss3540         1
## maternal_loss4045         1
## maternal_lossunclear      1
## older_siblings1           1
## older_siblings2           1
## older_siblings3           1
## older_siblings4           1
## older_siblings5P          1
## nr.siblings               1
## last_born1                1
## 
## Samples were drawn using sampling(NUTS). For each parameter, Eff.Sample 
## is a crude measure of effective sample size, and Rhat is the potential 
## scale reduction factor on split chains (at convergence, Rhat = 1).

Table of fixed effects

fixed_eff = data.frame(model_summary$fixed, check.names = F)
fixed_eff$Est.Error = fixed_eff$Eff.Sample = fixed_eff$Rhat = NULL
fixed_eff$OR = exp(fixed_eff$Estimate)
fixed_eff$OR_low = exp(fixed_eff$`l-95% CI`)
fixed_eff$OR_high = exp(fixed_eff$`u-95% CI`)
pander::pander(fixed_eff)
  Estimate l-95% CI u-95% CI OR OR_low OR_high
Intercept 0.3954 0.3526 0.4395 1.485 1.423 1.552
paternalage -0.0527 -0.08027 -0.02564 0.9487 0.9229 0.9747
age5063 0.1968 0.1558 0.2382 1.217 1.169 1.269
age025 -2.874 -3.166 -2.607 0.05648 0.04219 0.07376
ageunknown 0.3459 0.3141 0.3783 1.413 1.369 1.46
birth_cohort1950M1955 -0.005947 -0.01622 0.004462 0.9941 0.9839 1.004
birth_cohort1955M1960 0.000356 -0.01135 0.01194 1 0.9887 1.012
male1 -0.05134 -0.05977 -0.04311 0.95 0.942 0.9578
maternalage.factor1420 0.07475 0.05981 0.09033 1.078 1.062 1.095
maternalage.factor3561 -0.00977 -0.02523 0.006405 0.9903 0.9751 1.006
paternalage.mean 0.007943 -0.01958 0.03602 1.008 0.9806 1.037
paternal_loss01 0.1668 -0.4096 0.6711 1.182 0.6639 1.956
paternal_loss15 -0.04718 -0.2035 0.1017 0.9539 0.8159 1.107
paternal_loss510 -0.06122 -0.1255 0.004951 0.9406 0.8821 1.005
paternal_loss1015 -0.02264 -0.06019 0.01294 0.9776 0.9416 1.013
paternal_loss1520 -0.0004874 -0.02704 0.026 0.9995 0.9733 1.026
paternal_loss2025 -0.002718 -0.02476 0.02012 0.9973 0.9755 1.02
paternal_loss2530 -0.01904 -0.03791 -0.0002248 0.9811 0.9628 0.9998
paternal_loss3035 -0.005255 -0.02298 0.01145 0.9948 0.9773 1.012
paternal_loss3540 -0.005259 -0.02075 0.01064 0.9948 0.9795 1.011
paternal_loss4045 -0.01196 -0.02746 0.00324 0.9881 0.9729 1.003
paternal_lossunclear -0.05787 -0.06959 -0.046 0.9438 0.9328 0.955
maternal_loss01 -0.2177 -0.7006 0.1999 0.8044 0.4963 1.221
maternal_loss15 -0.007657 -0.2109 0.1984 0.9924 0.8098 1.219
maternal_loss510 0.0186 -0.07007 0.1069 1.019 0.9323 1.113
maternal_loss1015 -0.01713 -0.07272 0.03708 0.983 0.9299 1.038
maternal_loss1520 -0.006065 -0.04531 0.03378 0.994 0.9557 1.034
maternal_loss2025 0.02439 -0.008897 0.05763 1.025 0.9911 1.059
maternal_loss2530 0.0001224 -0.02803 0.02921 1 0.9724 1.03
maternal_loss3035 -0.001922 -0.026 0.02045 0.9981 0.9743 1.021
maternal_loss3540 0.001399 -0.01868 0.02075 1.001 0.9815 1.021
maternal_loss4045 -0.008895 -0.02769 0.0106 0.9911 0.9727 1.011
maternal_lossunclear -0.02213 -0.0325 -0.01206 0.9781 0.968 0.988
older_siblings1 0.0187 0.004737 0.03233 1.019 1.005 1.033
older_siblings2 0.03017 0.007255 0.0542 1.031 1.007 1.056
older_siblings3 0.0426 0.009681 0.07687 1.044 1.01 1.08
older_siblings4 0.01366 -0.03312 0.06071 1.014 0.9674 1.063
older_siblings5P -0.04236 -0.1021 0.01854 0.9585 0.9029 1.019
nr.siblings 0.03724 0.0316 0.04254 1.038 1.032 1.043
last_born1 0.00826 -0.001796 0.01796 1.008 0.9982 1.018

Paternal age effect

pander::pander(paternal_age_10y_effect(model))
effect median_estimate ci_95 ci_80
estimate father 25y 1.39 [1.33;1.44] [1.35;1.42]
estimate father 35y 1.31 [1.26;1.37] [1.28;1.35]
percentage change -5.1 [-7.71;-2.53] [-6.91;-3.42]
OR/IRR -0.05 [-0.08;-0.03] [-0.07;-0.03]

Marginal effect plots

plot.brmsMarginalEffects_shades(
    x = marginal_effects(model, re_formula = NA, probs = c(0.025,0.975)),
    y = marginal_effects(model, re_formula = NA, probs = c(0.1,0.9)), 
    ask = FALSE)

Coefficient plot

Coefficient estimates (95% and 80% credibility).

mcmc_intervals(as.matrix(model$fit), regex_pars = "b_[^I]")

mcmc_areas(as.matrix(model$fit), regex_pars = "b")

Diagnostics

These plots were made to diagnose misfit and nonconvergence.

Posterior predictive checks

brms::pp_check(model, re_formula = NA, type = "dens_overlay")

brms::pp_check(model, re_formula = NA, type = "hist")

Rhat

Did the 6 chains converge?

stanplot(model, pars = "^b_[^I]", type = 'rhat')

Effective sample size by average sample size

stanplot(model, pars = "^b", type = 'ess')

Monte Carlo SE

stanplot(model, pars = "^b", type = 'mcse')

Trace plots

if(any( summary(model)$fixed[,"Rhat"] > 1.1)) { # only do traceplots if not converged
    plot(model, N = 3, ask = FALSE)
}

Further plots

stanplot(model, pars = "^b", type = 'diag')

File name

coefs/swed/s1_control_age.rds

Cluster script

opts_chunk$set(echo = FALSE)
clusterscript = str_replace(basename(model_filename), "\\.rds",".html")
cat("[Cluster script](" , clusterscript, ")", sep = "")

Cluster script

s2: Mediation via reproductive timing

Here, we tested whether the paternal age effect on reproductive succes is mediated by reproductive timing (as indexed by anchors’ ages at first and last birth). Because age at first and last birth are by definition only available for anchors who had at least one child, this analysis has to be restricted to such anchors. Hence, paternal age effects on mortality until age 1 and 15 cannot, in principle, be mediated by reproductive timing of the anchors.

Model summary

Full summary

model_summary = summary(model, use_cache = FALSE, priors = TRUE)
print(model_summary)
##  Family: poisson (log) 
## Formula: children ~ paternalage + age_at_1st_child + age_at_last_child + birth_cohort + male + maternalage.factor + paternalage.mean + paternal_loss + maternal_loss + older_siblings + nr.siblings + last_born + (1 | idParents) 
##    Data: model_data (Number of observations: 94472) 
## Samples: 6 chains, each with iter = 800; warmup = 300; thin = 1; 
##          total post-warmup samples = 3000
##    WAIC: Not computed
##  
## Priors: 
## b ~ normal(0,5)
## sd ~ student_t(3, 0, 5)
## 
## Group-Level Effects: 
## ~idParents (Number of levels: 64717) 
##               Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
## sd(Intercept)        0         0        0     0.01       3000    1
## 
## Population-Level Effects: 
##                        Estimate Est.Error l-95% CI u-95% CI Eff.Sample
## Intercept                  0.47      0.02     0.43     0.50       3000
## paternalage                0.00      0.01    -0.03     0.02        884
## age_at_1st_child          -0.53      0.00    -0.54    -0.52       3000
## age_at_last_child          0.51      0.00     0.50     0.52       3000
## birth_cohort1950M1955      0.02      0.01     0.01     0.03       3000
## birth_cohort1955M1960      0.06      0.01     0.05     0.07       3000
## male1                      0.00      0.00    -0.01     0.01       3000
## maternalage.factor1420    -0.01      0.01    -0.03     0.00       3000
## maternalage.factor3561     0.01      0.01    -0.01     0.02       3000
## paternalage.mean           0.02      0.01     0.00     0.05        901
## paternal_loss01           -0.14      0.27    -0.71     0.37       3000
## paternal_loss15           -0.05      0.08    -0.19     0.10       3000
## paternal_loss510          -0.02      0.03    -0.09     0.04       3000
## paternal_loss1015         -0.02      0.02    -0.05     0.02       3000
## paternal_loss1520         -0.02      0.01    -0.05     0.01       3000
## paternal_loss2025         -0.02      0.01    -0.05     0.00       3000
## paternal_loss2530         -0.02      0.01    -0.04     0.00       3000
## paternal_loss3035         -0.01      0.01    -0.03     0.01       3000
## paternal_loss3540         -0.01      0.01    -0.03     0.00       3000
## paternal_loss4045         -0.01      0.01    -0.02     0.01       3000
## paternal_lossunclear       0.01      0.01     0.00     0.02       3000
## maternal_loss01           -0.04      0.21    -0.47     0.35       3000
## maternal_loss15            0.03      0.10    -0.18     0.22       3000
## maternal_loss510           0.01      0.05    -0.08     0.11       3000
## maternal_loss1015          0.01      0.03    -0.05     0.06       3000
## maternal_loss1520         -0.01      0.02    -0.05     0.03       3000
## maternal_loss2025         -0.01      0.02    -0.04     0.03       3000
## maternal_loss2530         -0.01      0.01    -0.04     0.02       3000
## maternal_loss3035          0.00      0.01    -0.02     0.03       3000
## maternal_loss3540         -0.01      0.01    -0.03     0.01       3000
## maternal_loss4045         -0.01      0.01    -0.03     0.01       3000
## maternal_lossunclear       0.01      0.01     0.00     0.02       3000
## older_siblings1            0.00      0.01    -0.01     0.01       1107
## older_siblings2            0.00      0.01    -0.02     0.02        957
## older_siblings3            0.00      0.02    -0.03     0.03       1031
## older_siblings4           -0.01      0.02    -0.06     0.03       1045
## older_siblings5P          -0.03      0.03    -0.09     0.03        968
## nr.siblings                0.01      0.00     0.00     0.02        988
## last_born1                 0.00      0.00    -0.01     0.01       3000
##                        Rhat
## Intercept              1.00
## paternalage            1.01
## age_at_1st_child       1.00
## age_at_last_child      1.00
## birth_cohort1950M1955  1.00
## birth_cohort1955M1960  1.00
## male1                  1.00
## maternalage.factor1420 1.00
## maternalage.factor3561 1.00
## paternalage.mean       1.01
## paternal_loss01        1.00
## paternal_loss15        1.00
## paternal_loss510       1.00
## paternal_loss1015      1.00
## paternal_loss1520      1.00
## paternal_loss2025      1.00
## paternal_loss2530      1.00
## paternal_loss3035      1.00
## paternal_loss3540      1.00
## paternal_loss4045      1.00
## paternal_lossunclear   1.00
## maternal_loss01        1.00
## maternal_loss15        1.00
## maternal_loss510       1.00
## maternal_loss1015      1.00
## maternal_loss1520      1.00
## maternal_loss2025      1.00
## maternal_loss2530      1.00
## maternal_loss3035      1.00
## maternal_loss3540      1.00
## maternal_loss4045      1.00
## maternal_lossunclear   1.00
## older_siblings1        1.00
## older_siblings2        1.00
## older_siblings3        1.00
## older_siblings4        1.01
## older_siblings5P       1.00
## nr.siblings            1.00
## last_born1             1.00
## 
## Samples were drawn using sampling(NUTS). For each parameter, Eff.Sample 
## is a crude measure of effective sample size, and Rhat is the potential 
## scale reduction factor on split chains (at convergence, Rhat = 1).

Table of fixed effects

fixed_eff = data.frame(model_summary$fixed, check.names = F)
fixed_eff$Est.Error = fixed_eff$Eff.Sample = fixed_eff$Rhat = NULL
fixed_eff$OR = exp(fixed_eff$Estimate)
fixed_eff$OR_low = exp(fixed_eff$`l-95% CI`)
fixed_eff$OR_high = exp(fixed_eff$`u-95% CI`)
pander::pander(fixed_eff)
  Estimate l-95% CI u-95% CI OR OR_low OR_high
Intercept 0.4684 0.4319 0.5047 1.597 1.54 1.657
paternalage -0.002744 -0.03026 0.02366 0.9973 0.9702 1.024
age_at_1st_child -0.5327 -0.5419 -0.5235 0.587 0.5816 0.5924
age_at_last_child 0.51 0.5024 0.5178 1.665 1.653 1.678
birth_cohort1950M1955 0.02073 0.01043 0.0315 1.021 1.01 1.032
birth_cohort1955M1960 0.05997 0.04831 0.07201 1.062 1.049 1.075
male1 0.001377 -0.007947 0.009968 1.001 0.9921 1.01
maternalage.factor1420 -0.01213 -0.0274 0.003278 0.9879 0.973 1.003
maternalage.factor3561 0.008294 -0.007585 0.02376 1.008 0.9924 1.024
paternalage.mean 0.02292 -0.003996 0.05078 1.023 0.996 1.052
paternal_loss01 -0.1402 -0.7126 0.3712 0.8692 0.4904 1.45
paternal_loss15 -0.04511 -0.1936 0.102 0.9559 0.824 1.107
paternal_loss510 -0.02468 -0.08548 0.03869 0.9756 0.9181 1.039
paternal_loss1015 -0.01679 -0.05465 0.02032 0.9834 0.9468 1.021
paternal_loss1520 -0.02211 -0.04956 0.005697 0.9781 0.9517 1.006
paternal_loss2025 -0.02252 -0.04529 -5.073e-05 0.9777 0.9557 0.9999
paternal_loss2530 -0.01728 -0.03638 0.0008725 0.9829 0.9643 1.001
paternal_loss3035 -0.01052 -0.02776 0.007093 0.9895 0.9726 1.007
paternal_loss3540 -0.0107 -0.02583 0.00456 0.9894 0.9745 1.005
paternal_loss4045 -0.008815 -0.02455 0.006862 0.9912 0.9758 1.007
paternal_lossunclear 0.008101 -0.003614 0.01981 1.008 0.9964 1.02
maternal_loss01 -0.03884 -0.4745 0.3494 0.9619 0.6222 1.418
maternal_loss15 0.02523 -0.1842 0.2202 1.026 0.8318 1.246
maternal_loss510 0.01205 -0.08096 0.1068 1.012 0.9222 1.113
maternal_loss1015 0.006282 -0.04911 0.06107 1.006 0.9521 1.063
maternal_loss1520 -0.009696 -0.05085 0.02993 0.9904 0.9504 1.03
maternal_loss2025 -0.005388 -0.03894 0.0274 0.9946 0.9618 1.028
maternal_loss2530 -0.007119 -0.03613 0.02056 0.9929 0.9645 1.021
maternal_loss3035 0.003013 -0.02171 0.02717 1.003 0.9785 1.028
maternal_loss3540 -0.009724 -0.03047 0.01084 0.9903 0.97 1.011
maternal_loss4045 -0.006065 -0.02515 0.01287 0.994 0.9752 1.013
maternal_lossunclear 0.006305 -0.00393 0.01676 1.006 0.9961 1.017
older_siblings1 0.001041 -0.01301 0.01462 1.001 0.9871 1.015
older_siblings2 -0.0003881 -0.02334 0.02218 0.9996 0.9769 1.022
older_siblings3 -0.0002695 -0.033 0.03429 0.9997 0.9675 1.035
older_siblings4 -0.01194 -0.05613 0.03398 0.9881 0.9454 1.035
older_siblings5P -0.0291 -0.0853 0.02746 0.9713 0.9182 1.028
nr.siblings 0.009865 0.004498 0.01542 1.01 1.005 1.016
last_born1 0.002988 -0.006336 0.01178 1.003 0.9937 1.012

Paternal age effect

pander::pander(paternal_age_10y_effect(model))
effect median_estimate ci_95 ci_80
estimate father 25y 2.22 [2.19;2.26] [2.2;2.25]
estimate father 35y 2.22 [2.16;2.27] [2.18;2.25]
percentage change -0.26 [-2.98;2.39] [-2.06;1.5]
OR/IRR 0 [-0.03;0.02] [-0.02;0.01]

Marginal effect plots

plot.brmsMarginalEffects_shades(
    x = marginal_effects(model, re_formula = NA, probs = c(0.025,0.975)),
    y = marginal_effects(model, re_formula = NA, probs = c(0.1,0.9)), 
    ask = FALSE)

Coefficient plot

Coefficient estimates (95% and 80% credibility).

mcmc_intervals(as.matrix(model$fit), regex_pars = "b_[^I]")

mcmc_areas(as.matrix(model$fit), regex_pars = "b")

Diagnostics

These plots were made to diagnose misfit and nonconvergence.

Posterior predictive checks

brms::pp_check(model, re_formula = NA, type = "dens_overlay")

brms::pp_check(model, re_formula = NA, type = "hist")

Rhat

Did the 6 chains converge?

stanplot(model, pars = "^b_[^I]", type = 'rhat')

Effective sample size by average sample size

stanplot(model, pars = "^b", type = 'ess')

Monte Carlo SE

stanplot(model, pars = "^b", type = 'mcse')

Trace plots

if(any( summary(model)$fixed[,"Rhat"] > 1.1)) { # only do traceplots if not converged
    plot(model, N = 3, ask = FALSE)
}

Further plots

stanplot(model, pars = "^b", type = 'diag')

File name

coefs/swed/s2_reproductive_timing.rds

Cluster script

opts_chunk$set(echo = FALSE)
clusterscript = str_replace(basename(model_filename), "\\.rds",".html")
cat("[Cluster script](" , clusterscript, ")", sep = "")

Cluster script