Sweden (1947-2000), selective episodes

Loading details

source("0__helpers.R")
opts_chunk$set(warning=TRUE, cache = FALSE,tidy=FALSE,dev=c('png'),fig.width=20,fig.height=12.5)

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_survival.1 dataset contains only those participants where paternal age is known and birth years are from 1969 to 2000. We use it for the analyses pertaining to survival. We randomly sampled 100,000 families from the whole dataset.

The swed_subset_children.1 dataset contains only those participants where paternal age is known and birth years are from 1947 to 1959. We use it for the analyses pertaining to marriage, divorce and reproductive success. We randomly sampled 100,000 families from the whole dataset.

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

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.

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.

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

Model summary

Full summary

model_summary = summary(model, use_cache = FALSE, priors = TRUE)
print(model_summary)
##  Family: bernoulli (cauchit) 
## Formula: survive1y ~ paternalage + 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: 363744) 
## 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: 200000) 
##               Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
## sd(Intercept)    21.41       0.8    19.76    22.96        205 1.03
## 
## Population-Level Effects: 
##                        Estimate Est.Error l-95% CI u-95% CI Eff.Sample
## Intercept                 96.06      4.72    87.17   105.59        463
## paternalage              -10.72      1.94   -14.50    -6.88        488
## birth_cohort1970M1975     -5.41      1.94    -9.34    -1.77        251
## birth_cohort1975M1980     -3.12      2.00    -7.06     0.75        317
## birth_cohort1980M1985      4.66      2.20     0.36     8.89        405
## birth_cohort1985M1990      3.90      2.14    -0.41     8.23        327
## birth_cohort1990M1995      5.97      2.20     1.72    10.42        361
## birth_cohort1995M2000     14.75      3.05     8.92    20.72        368
## male1                     -4.69      1.00    -6.73    -2.73        337
## maternalage.factor1420     0.42      2.18    -3.58     4.93        578
## maternalage.factor3561    -4.91      1.75    -8.19    -1.47        753
## paternalage.mean          10.23      2.04     6.20    14.17        665
## paternal_loss01            0.33      4.99    -9.73    10.01       3000
## paternal_loss15           -5.30      4.50   -13.41     4.04       1098
## paternal_loss510           3.95      4.33    -4.09    12.77       1661
## paternal_loss1015          3.54      3.74    -3.41    11.17       1277
## paternal_loss1520         -3.51      2.94    -8.89     2.43        708
## paternal_loss2025         -2.70      2.99    -8.26     3.67        464
## paternal_loss2530          0.87      2.68    -4.04     6.38        907
## paternal_loss3035          2.02      2.93    -3.34     7.97        514
## paternal_loss3540         -4.42      2.66    -9.17     1.01        819
## paternal_loss4045         -1.32      4.98   -10.96     8.96       3000
## maternal_loss01           -5.21      5.62   -16.28     5.58        765
## maternal_loss15            1.90      4.72    -7.17    11.19       3000
## maternal_loss510           1.54      4.54    -6.87    10.56       3000
## maternal_loss1015         -0.35      4.18    -8.03     8.20       1203
## maternal_loss1520         -6.18      3.25   -12.29     0.52        941
## maternal_loss2025         -6.09      3.05   -11.88     0.33        951
## maternal_loss2530         -4.50      3.09   -10.11     1.93        978
## maternal_loss3035         -2.91      3.27    -9.07     3.99        792
## maternal_loss3540         -1.60      3.96    -8.63     6.56        848
## maternal_loss4045          0.97      4.70    -8.54    10.14       3000
## older_siblings1           -6.34      1.37    -9.06    -3.75        249
## older_siblings2            0.01      2.01    -3.87     3.98        147
## older_siblings3           14.49      2.50     9.72    19.44        326
## older_siblings4           16.49      2.71    11.00    21.83        602
## older_siblings5P          25.63      3.53    18.62    32.67        625
## nr.siblings               -7.96      0.40    -8.76    -7.19        420
## last_born1               -12.79      1.02   -14.73   -10.83        245
##                        Rhat
## Intercept              1.01
## paternalage            1.01
## birth_cohort1970M1975  1.01
## birth_cohort1975M1980  1.01
## birth_cohort1980M1985  1.00
## birth_cohort1985M1990  1.01
## birth_cohort1990M1995  1.01
## birth_cohort1995M2000  1.01
## male1                  1.02
## maternalage.factor1420 1.01
## maternalage.factor3561 1.01
## paternalage.mean       1.01
## paternal_loss01        1.00
## paternal_loss15        1.00
## paternal_loss510       1.00
## paternal_loss1015      1.00
## paternal_loss1520      1.01
## paternal_loss2025      1.01
## paternal_loss2530      1.01
## paternal_loss3035      1.01
## paternal_loss3540      1.01
## paternal_loss4045      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.01
## maternal_loss2530      1.00
## maternal_loss3035      1.00
## maternal_loss3540      1.01
## maternal_loss4045      1.00
## older_siblings1        1.02
## older_siblings2        1.06
## older_siblings3        1.02
## older_siblings4        1.01
## older_siblings5P       1.01
## nr.siblings            1.01
## last_born1             1.03
## 
## 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 96.06 87.17 105.6 5.236e+41 7.169e+37 7.219e+45
paternalage -10.72 -14.5 -6.876 2.215e-05 5.064e-07 0.001033
birth_cohort1970M1975 -5.409 -9.344 -1.769 0.004475 8.75e-05 0.1705
birth_cohort1975M1980 -3.117 -7.06 0.7462 0.04427 0.0008587 2.109
birth_cohort1980M1985 4.659 0.3597 8.886 105.5 1.433 7230
birth_cohort1985M1990 3.904 -0.4149 8.228 49.61 0.6604 3745
birth_cohort1990M1995 5.968 1.715 10.42 390.9 5.559 33363
birth_cohort1995M2000 14.75 8.923 20.72 2556476 7504 996516331
male1 -4.693 -6.735 -2.734 0.009161 0.001189 0.06493
maternalage.factor1420 0.4207 -3.583 4.931 1.523 0.0278 138.5
maternalage.factor3561 -4.912 -8.194 -1.472 0.007358 0.0002762 0.2295
paternalage.mean 10.23 6.197 14.17 27741 491.5 1420899
paternal_loss01 0.3262 -9.732 10.01 1.386 5.937e-05 22327
paternal_loss15 -5.304 -13.41 4.044 0.004974 1.502e-06 57.06
paternal_loss510 3.951 -4.088 12.77 51.97 0.01677 352889
paternal_loss1015 3.542 -3.406 11.17 34.55 0.03319 70655
paternal_loss1520 -3.514 -8.889 2.433 0.02978 0.0001379 11.39
paternal_loss2025 -2.7 -8.257 3.67 0.06721 0.0002594 39.27
paternal_loss2530 0.8731 -4.045 6.378 2.394 0.01751 588.8
paternal_loss3035 2.018 -3.34 7.971 7.52 0.03542 2895
paternal_loss3540 -4.425 -9.17 1.006 0.01198 0.0001041 2.735
paternal_loss4045 -1.32 -10.96 8.964 0.2672 1.741e-05 7818
maternal_loss01 -5.208 -16.28 5.575 0.005475 8.498e-08 263.8
maternal_loss15 1.898 -7.174 11.19 6.672 0.0007666 72578
maternal_loss510 1.542 -6.873 10.56 4.676 0.001035 38393
maternal_loss1015 -0.3531 -8.034 8.2 0.7025 0.0003243 3641
maternal_loss1520 -6.179 -12.29 0.5236 0.002073 4.595e-06 1.688
maternal_loss2025 -6.093 -11.88 0.3274 0.002258 6.936e-06 1.387
maternal_loss2530 -4.496 -10.11 1.926 0.01115 4.053e-05 6.861
maternal_loss3035 -2.908 -9.073 3.989 0.05461 0.0001147 54.01
maternal_loss3540 -1.598 -8.632 6.563 0.2023 0.0001782 708.5
maternal_loss4045 0.9686 -8.537 10.14 2.634 0.0001961 25350
older_siblings1 -6.344 -9.063 -3.752 0.001758 0.0001159 0.02347
older_siblings2 0.01168 -3.867 3.979 1.012 0.02092 53.47
older_siblings3 14.49 9.718 19.44 1966735 16620 276670549
older_siblings4 16.49 11 21.83 14519530 59968 3.018e+09
older_siblings5P 25.63 18.62 32.67 1.348e+11 121468561 1.54e+14
nr.siblings -7.964 -8.763 -7.185 0.0003476 0.0001564 0.0007578
last_born1 -12.79 -14.73 -10.83 2.8e-06 3.996e-07 1.978e-05

Paternal age effect

pander::pander(paternal_age_10y_effect(model))
effect median_estimate ci_95 ci_80
estimate father 25y 0.9965 [0.9963;0.9968] [0.9964;0.9967]
estimate father 35y 0.9961 [0.9958;0.9964] [0.9959;0.9963]
percentage change -0.0456 [-0.0645;-0.0288] [-0.0578;-0.0336]
OR/IRR 0 [0;0.001] [0;3e-04]

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)

plot of chunk unnamed-chunk-12plot of chunk unnamed-chunk-12plot of chunk unnamed-chunk-12plot of chunk unnamed-chunk-12plot of chunk unnamed-chunk-12plot of chunk unnamed-chunk-12plot of chunk unnamed-chunk-12plot of chunk unnamed-chunk-12plot of chunk unnamed-chunk-12plot of chunk unnamed-chunk-12

Coefficient plot

Coefficient estimates (95% and 80% credibility).

mcmc_intervals(as.matrix(model$fit), regex_pars = "b_[^I]")
plot of chunk unnamed-chunk-13

plot of chunk unnamed-chunk-13

mcmc_areas(as.matrix(model$fit), regex_pars = "b")
plot of chunk unnamed-chunk-13

plot of chunk unnamed-chunk-13

Diagnostics

These plots were made to diagnose misfit and nonconvergence.

Posterior predictive checks

brms::pp_check(model, re_formula = NA, type = "dens_overlay")
plot of chunk unnamed-chunk-14

plot of chunk unnamed-chunk-14

brms::pp_check(model, re_formula = NA, type = "hist")
plot of chunk unnamed-chunk-14

plot of chunk unnamed-chunk-14

Rhat

Did the 6 chains converge?

stanplot(model, pars = "^b_[^I]", type = 'rhat')
plot of chunk unnamed-chunk-15

plot of chunk unnamed-chunk-15

Effective sample size by average sample size

stanplot(model, pars = "^b", type = 'ess')
plot of chunk unnamed-chunk-16

plot of chunk unnamed-chunk-16

Monte Carlo SE

stanplot(model, pars = "^b", type = 'mcse')
plot of chunk unnamed-chunk-17

plot of chunk unnamed-chunk-17

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')
plot of chunk unnamed-chunk-19

plot of chunk unnamed-chunk-19

File name

coefs/swed/e1_survive1y.rds

Cluster script

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

Cluster script

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

Model summary

Full summary

model_summary = summary(model, use_cache = FALSE, priors = TRUE)
print(model_summary)
##  Family: bernoulli (cauchit) 
## Formula: surviveR ~ paternalage + 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: 242223) 
## Samples: 6 chains, each with iter = 1500; warmup = 1000; 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: 142552) 
##               Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
## sd(Intercept)    12.03     14.82     0.18    44.82          3 3.16
## 
## Population-Level Effects: 
##                        Estimate Est.Error l-95% CI u-95% CI Eff.Sample
## Intercept                164.48     18.24   128.07   198.74         13
## paternalage               -1.45      4.31    -9.43     6.83         16
## birth_cohort1970M1975     -5.91      4.58   -14.91     2.93         38
## birth_cohort1975M1980      0.40      4.61    -8.68     9.73        136
## birth_cohort1980M1985      4.10      4.23    -4.18    12.52       3000
## birth_cohort1985M1990      5.10      4.65    -3.67    14.13         23
## male1                     -8.81      4.07   -17.11    -0.69        449
## maternalage.factor1420     0.37      4.77    -9.02     9.63        112
## maternalage.factor3561     2.78      4.38    -6.17    11.42       1518
## paternalage.mean          -6.09      4.58   -14.76     2.81         18
## paternal_loss01            0.15      4.98    -9.57    10.14       3000
## paternal_loss15           -2.03      5.33   -11.89     8.41         61
## paternal_loss510          -3.53      5.25   -13.97     6.59         23
## paternal_loss1015         -0.60      4.79   -10.00     8.31       3000
## paternal_loss1520         -2.03      4.79   -11.49     7.18       3000
## paternal_loss2025         -1.02      5.09   -10.78     9.00        244
## paternal_loss2530         -0.96      4.78    -9.99     8.84       3000
## paternal_loss3035         -1.09      4.84   -10.80     8.69       3000
## paternal_loss3540         -0.07      5.18    -9.68     9.69        145
## paternal_loss4045          0.10      4.96    -9.51     9.71       3000
## maternal_loss01            0.20      5.22   -10.28    10.09       3000
## maternal_loss15           -1.93      5.28   -12.22     8.06       3000
## maternal_loss510          -2.65      5.26   -12.38     7.57        143
## maternal_loss1015         -1.65      4.98   -11.48     7.88       3000
## maternal_loss1520         -0.03      5.10    -9.89     9.85       3000
## maternal_loss2025         -1.24      5.36   -12.54     9.32         99
## maternal_loss2530         -1.75      4.85   -11.41     7.64       3000
## maternal_loss3035         -1.14      4.89   -10.40     8.51       3000
## maternal_loss3540         -0.76      5.12   -10.77     8.67        216
## maternal_loss4045         -0.67      4.92   -10.22     8.94       3000
## older_siblings1            4.39      4.31    -3.89    12.88        835
## older_siblings2           -1.03      5.32   -12.24     8.81         19
## older_siblings3           -0.28      5.02   -10.10     9.45         59
## older_siblings4            0.04      4.72    -9.62     9.52       3000
## older_siblings5P           1.55      4.91    -7.83    11.28       3000
## nr.siblings               -7.16      1.07    -8.87    -4.60         13
## last_born1                 7.96      4.37    -0.74    16.56       1761
##                        Rhat
## Intercept              1.14
## paternalage            1.13
## birth_cohort1970M1975  1.05
## birth_cohort1975M1980  1.03
## birth_cohort1980M1985  1.02
## birth_cohort1985M1990  1.08
## male1                  1.02
## maternalage.factor1420 1.03
## maternalage.factor3561 1.01
## paternalage.mean       1.10
## paternal_loss01        1.00
## paternal_loss15        1.06
## paternal_loss510       1.07
## paternal_loss1015      1.01
## paternal_loss1520      1.01
## paternal_loss2025      1.02
## paternal_loss2530      1.01
## paternal_loss3035      1.00
## paternal_loss3540      1.03
## paternal_loss4045      1.00
## maternal_loss01        1.02
## maternal_loss15        1.04
## maternal_loss510       1.04
## maternal_loss1015      1.01
## maternal_loss1520      1.01
## maternal_loss2025      1.05
## maternal_loss2530      1.01
## maternal_loss3035      1.01
## maternal_loss3540      1.03
## maternal_loss4045      1.01
## older_siblings1        1.00
## older_siblings2        1.15
## older_siblings3        1.03
## older_siblings4        1.01
## older_siblings5P       1.01
## nr.siblings            1.23
## last_born1             1.01
## 
## 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 164.5 128.1 198.7 2.716e+71 4.16e+55 2.05e+86
paternalage -1.447 -9.427 6.831 0.2353 8.053e-05 926
birth_cohort1970M1975 -5.914 -14.91 2.929 0.002701 3.358e-07 18.7
birth_cohort1975M1980 0.3976 -8.68 9.731 1.488 0.0001699 16831
birth_cohort1980M1985 4.095 -4.183 12.52 60.07 0.01525 273667
birth_cohort1985M1990 5.096 -3.673 14.13 163.3 0.02539 1365784
male1 -8.813 -17.11 -0.6904 0.0001488 3.712e-08 0.5014
maternalage.factor1420 0.3747 -9.023 9.628 1.455 0.0001206 15183
maternalage.factor3561 2.78 -6.174 11.42 16.12 0.002082 90896
paternalage.mean -6.094 -14.76 2.811 0.002257 3.904e-07 16.63
paternal_loss01 0.1468 -9.571 10.14 1.158 6.971e-05 25324
paternal_loss15 -2.027 -11.89 8.408 0.1317 6.885e-06 4485
paternal_loss510 -3.528 -13.97 6.591 0.02937 8.6e-07 728.4
paternal_loss1015 -0.6017 -9.999 8.31 0.5479 4.546e-05 4062
paternal_loss1520 -2.032 -11.49 7.175 0.1311 1.023e-05 1307
paternal_loss2025 -1.018 -10.78 9.003 0.3613 2.081e-05 8126
paternal_loss2530 -0.9557 -9.986 8.838 0.3845 4.604e-05 6894
paternal_loss3035 -1.087 -10.8 8.688 0.3373 2.049e-05 5931
paternal_loss3540 -0.07058 -9.683 9.691 0.9319 6.235e-05 16168
paternal_loss4045 0.1044 -9.508 9.711 1.11 7.428e-05 16500
maternal_loss01 0.1976 -10.28 10.09 1.218 3.431e-05 24157
maternal_loss15 -1.929 -12.22 8.062 0.1454 4.919e-06 3171
maternal_loss510 -2.655 -12.38 7.571 0.07033 4.196e-06 1940
maternal_loss1015 -1.651 -11.48 7.879 0.1918 1.029e-05 2642
maternal_loss1520 -0.03434 -9.891 9.852 0.9662 5.065e-05 18993
maternal_loss2025 -1.238 -12.54 9.316 0.2899 3.571e-06 11114
maternal_loss2530 -1.749 -11.41 7.641 0.174 1.111e-05 2082
maternal_loss3035 -1.136 -10.4 8.511 0.3212 3.057e-05 4968
maternal_loss3540 -0.7648 -10.77 8.667 0.4654 2.105e-05 5810
maternal_loss4045 -0.6703 -10.22 8.935 0.5115 3.659e-05 7595
older_siblings1 4.389 -3.892 12.88 80.54 0.02041 390959
older_siblings2 -1.029 -12.24 8.811 0.3574 4.841e-06 6704
older_siblings3 -0.2772 -10.1 9.447 0.7579 4.112e-05 12664
older_siblings4 0.04359 -9.618 9.52 1.045 6.655e-05 13633
older_siblings5P 1.552 -7.827 11.28 4.723 0.0003987 79410
nr.siblings -7.165 -8.872 -4.604 0.0007735 0.0001403 0.01001
last_born1 7.96 -0.7404 16.56 2864 0.4769 15587840

Paternal age effect

pander::pander(paternal_age_10y_effect(model))
effect median_estimate ci_95 ci_80
estimate father 25y 0.9976 [0.9973;0.998] [0.9974;0.9979]
estimate father 35y 0.9976 [0.9973;0.9981] [0.9974;0.9979]
percentage change -0.0025 [-0.0176;0.0104] [-0.0132;0.0067]
OR/IRR 0.251 [1e-04;925.9538] [8e-04;65.491]

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)

plot of chunk unnamed-chunk-25plot of chunk unnamed-chunk-25plot of chunk unnamed-chunk-25plot of chunk unnamed-chunk-25plot of chunk unnamed-chunk-25plot of chunk unnamed-chunk-25plot of chunk unnamed-chunk-25plot of chunk unnamed-chunk-25plot of chunk unnamed-chunk-25plot of chunk unnamed-chunk-25

Coefficient plot

Coefficient estimates (95% and 80% credibility).

mcmc_intervals(as.matrix(model$fit), regex_pars = "b_[^I]")
plot of chunk unnamed-chunk-26

plot of chunk unnamed-chunk-26

mcmc_areas(as.matrix(model$fit), regex_pars = "b")
plot of chunk unnamed-chunk-26

plot of chunk unnamed-chunk-26

Diagnostics

These plots were made to diagnose misfit and nonconvergence.

Posterior predictive checks

brms::pp_check(model, re_formula = NA, type = "dens_overlay")
plot of chunk unnamed-chunk-27

plot of chunk unnamed-chunk-27

brms::pp_check(model, re_formula = NA, type = "hist")
plot of chunk unnamed-chunk-27

plot of chunk unnamed-chunk-27

Rhat

Did the 6 chains converge?

stanplot(model, pars = "^b_[^I]", type = 'rhat')
plot of chunk unnamed-chunk-28

plot of chunk unnamed-chunk-28

Effective sample size by average sample size

stanplot(model, pars = "^b", type = 'ess')
plot of chunk unnamed-chunk-29

plot of chunk unnamed-chunk-29

Monte Carlo SE

stanplot(model, pars = "^b", type = 'mcse')
plot of chunk unnamed-chunk-30

plot of chunk unnamed-chunk-30

Trace plots

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

plot of chunk unnamed-chunk-31plot of chunk unnamed-chunk-31plot of chunk unnamed-chunk-31plot of chunk unnamed-chunk-31plot of chunk unnamed-chunk-31plot of chunk unnamed-chunk-31plot of chunk unnamed-chunk-31plot of chunk unnamed-chunk-31plot of chunk unnamed-chunk-31plot of chunk unnamed-chunk-31plot of chunk unnamed-chunk-31plot of chunk unnamed-chunk-31plot of chunk unnamed-chunk-31

Further plots

stanplot(model, pars = "^b", type = 'diag')
plot of chunk unnamed-chunk-32

plot of chunk unnamed-chunk-32

File name

coefs/swed/e2_surviveR.rds

Cluster script

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

Cluster script

e3: Probability of ever marrying

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

Model summary

Full summary

model_summary = summary(model, use_cache = FALSE, priors = TRUE)
print(model_summary)
##  Family: bernoulli (cauchit) 
## Formula: ever_married ~ paternalage + 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: 1392421) 
## 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: 874603) 
##               Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
## sd(Intercept)     0.93      0.01     0.91     0.94        572 1.01
## 
## Population-Level Effects: 
##                        Estimate Est.Error l-95% CI u-95% CI Eff.Sample
## Intercept                  2.35      0.02     2.31     2.39       3000
## paternalage                0.14      0.02     0.11     0.18        940
## birth_cohort1950M1955     -0.46      0.01    -0.48    -0.45       3000
## birth_cohort1955M1960     -0.89      0.01    -0.90    -0.87       3000
## male1                     -0.47      0.01    -0.48    -0.46       3000
## maternalage.factor1420    -0.05      0.01    -0.07    -0.03       3000
## maternalage.factor3561    -0.01      0.01    -0.03     0.01       3000
## paternalage.mean          -0.22      0.02    -0.26    -0.19        991
## paternal_loss01           -0.18      0.22    -0.59     0.27       3000
## paternal_loss15           -0.15      0.07    -0.28    -0.01       3000
## paternal_loss510          -0.19      0.03    -0.25    -0.12       3000
## paternal_loss1015         -0.20      0.02    -0.24    -0.15       3000
## paternal_loss1520         -0.17      0.02    -0.20    -0.14       3000
## paternal_loss2025         -0.13      0.01    -0.15    -0.10       3000
## paternal_loss2530         -0.10      0.01    -0.12    -0.08       3000
## paternal_loss3035         -0.10      0.01    -0.12    -0.08       3000
## paternal_loss3540         -0.10      0.01    -0.11    -0.08       3000
## paternal_loss4045         -0.08      0.01    -0.10    -0.06       3000
## paternal_lossunclear       0.00      0.01    -0.02     0.01       3000
## maternal_loss01           -0.51      0.22    -0.93    -0.05       3000
## maternal_loss15           -0.32      0.10    -0.50    -0.12       3000
## maternal_loss510          -0.22      0.05    -0.31    -0.12       3000
## maternal_loss1015         -0.20      0.03    -0.26    -0.14       3000
## maternal_loss1520         -0.18      0.02    -0.23    -0.14       3000
## maternal_loss2025         -0.07      0.02    -0.11    -0.03       3000
## maternal_loss2530         -0.09      0.02    -0.13    -0.06       3000
## maternal_loss3035         -0.07      0.01    -0.10    -0.04       3000
## maternal_loss3540         -0.08      0.01    -0.10    -0.06       3000
## maternal_loss4045         -0.06      0.01    -0.08    -0.03       3000
## maternal_lossunclear       0.04      0.01     0.03     0.06       3000
## older_siblings1           -0.01      0.01    -0.03     0.00       1268
## older_siblings2           -0.03      0.01    -0.06     0.00        949
## older_siblings3           -0.07      0.02    -0.11    -0.03       1082
## older_siblings4           -0.15      0.03    -0.21    -0.10       1066
## older_siblings5P          -0.24      0.04    -0.31    -0.17        880
## nr.siblings                0.02      0.00     0.01     0.03        999
## last_born1                 0.05      0.01     0.04     0.06       3000
##                        Rhat
## Intercept              1.00
## paternalage            1.01
## 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.01
## older_siblings3        1.01
## older_siblings4        1.01
## older_siblings5P       1.01
## nr.siblings            1.01
## 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 2.353 2.312 2.394 10.52 10.1 10.95
paternalage 0.1418 0.1089 0.1765 1.152 1.115 1.193
birth_cohort1950M1955 -0.4612 -0.4766 -0.4461 0.6305 0.6209 0.6401
birth_cohort1955M1960 -0.8871 -0.9025 -0.8707 0.4118 0.4055 0.4186
male1 -0.473 -0.4837 -0.4625 0.6232 0.6165 0.6297
maternalage.factor1420 -0.04575 -0.06579 -0.02505 0.9553 0.9363 0.9753
maternalage.factor3561 -0.01041 -0.02866 0.00771 0.9896 0.9717 1.008
paternalage.mean -0.2246 -0.259 -0.1912 0.7989 0.7718 0.826
paternal_loss01 -0.1785 -0.5904 0.2664 0.8365 0.5541 1.305
paternal_loss15 -0.1461 -0.2842 -0.01102 0.8641 0.7526 0.989
paternal_loss510 -0.1853 -0.2495 -0.1205 0.8309 0.7792 0.8865
paternal_loss1015 -0.1955 -0.2361 -0.1542 0.8224 0.7897 0.8571
paternal_loss1520 -0.1681 -0.2003 -0.1355 0.8452 0.8185 0.8733
paternal_loss2025 -0.1259 -0.1518 -0.09911 0.8817 0.8592 0.9056
paternal_loss2530 -0.09939 -0.1236 -0.0763 0.9054 0.8837 0.9265
paternal_loss3035 -0.1026 -0.1237 -0.08148 0.9025 0.8836 0.9218
paternal_loss3540 -0.09553 -0.1141 -0.07698 0.9089 0.8921 0.9259
paternal_loss4045 -0.07848 -0.09886 -0.05857 0.9245 0.9059 0.9431
paternal_lossunclear -0.002007 -0.01649 0.0132 0.998 0.9836 1.013
maternal_loss01 -0.5136 -0.9337 -0.05031 0.5983 0.3931 0.9509
maternal_loss15 -0.3155 -0.5008 -0.1205 0.7294 0.606 0.8865
maternal_loss510 -0.2187 -0.314 -0.1223 0.8035 0.7305 0.8848
maternal_loss1015 -0.1975 -0.2617 -0.1357 0.8208 0.7698 0.8731
maternal_loss1520 -0.1839 -0.2295 -0.136 0.832 0.7949 0.8729
maternal_loss2025 -0.06734 -0.1071 -0.02816 0.9349 0.8985 0.9722
maternal_loss2530 -0.09459 -0.1276 -0.05932 0.9097 0.8802 0.9424
maternal_loss3035 -0.07304 -0.1025 -0.04397 0.9296 0.9025 0.957
maternal_loss3540 -0.08144 -0.1046 -0.05818 0.9218 0.9006 0.9435
maternal_loss4045 -0.0554 -0.07805 -0.03289 0.9461 0.9249 0.9676
maternal_lossunclear 0.04495 0.03151 0.05854 1.046 1.032 1.06
older_siblings1 -0.01366 -0.03059 0.002694 0.9864 0.9699 1.003
older_siblings2 -0.03255 -0.062 -0.004635 0.968 0.9399 0.9954
older_siblings3 -0.0705 -0.1123 -0.02879 0.9319 0.8938 0.9716
older_siblings4 -0.1522 -0.2111 -0.09716 0.8588 0.8097 0.9074
older_siblings5P -0.2402 -0.3137 -0.1673 0.7865 0.7307 0.8459
nr.siblings 0.01831 0.01107 0.02555 1.018 1.011 1.026
last_born1 0.04864 0.03728 0.06039 1.05 1.038 1.062

Paternal age effect

pander::pander(paternal_age_10y_effect(model))
effect median_estimate ci_95 ci_80
estimate father 25y 0.8554 [0.8538;0.8569] [0.8544;0.8564]
estimate father 35y 0.8636 [0.8617;0.8656] [0.8624;0.8649]
percentage change 0.8217 [0.6339;1.0141] [0.6923;0.9505]
OR/IRR 1.1521 [1.115;1.193] [1.1265;1.1783]

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)

plot of chunk unnamed-chunk-38plot of chunk unnamed-chunk-38plot of chunk unnamed-chunk-38plot of chunk unnamed-chunk-38plot of chunk unnamed-chunk-38plot of chunk unnamed-chunk-38plot of chunk unnamed-chunk-38plot of chunk unnamed-chunk-38plot of chunk unnamed-chunk-38plot of chunk unnamed-chunk-38

Coefficient plot

Coefficient estimates (95% and 80% credibility).

mcmc_intervals(as.matrix(model$fit), regex_pars = "b_[^I]")
plot of chunk unnamed-chunk-39

plot of chunk unnamed-chunk-39

mcmc_areas(as.matrix(model$fit), regex_pars = "b")
plot of chunk unnamed-chunk-39

plot of chunk unnamed-chunk-39

Diagnostics

These plots were made to diagnose misfit and nonconvergence.

Posterior predictive checks

brms::pp_check(model, re_formula = NA, type = "dens_overlay")
plot of chunk unnamed-chunk-40

plot of chunk unnamed-chunk-40

brms::pp_check(model, re_formula = NA, type = "hist")
plot of chunk unnamed-chunk-40

plot of chunk unnamed-chunk-40

Rhat

Did the 6 chains converge?

stanplot(model, pars = "^b_[^I]", type = 'rhat')
plot of chunk unnamed-chunk-41

plot of chunk unnamed-chunk-41

Effective sample size by average sample size

stanplot(model, pars = "^b", type = 'ess')
plot of chunk unnamed-chunk-42

plot of chunk unnamed-chunk-42

Monte Carlo SE

stanplot(model, pars = "^b", type = 'mcse')
plot of chunk unnamed-chunk-43

plot of chunk unnamed-chunk-43

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')
plot of chunk unnamed-chunk-45

plot of chunk unnamed-chunk-45

File name

coefs/swed/e3_ever_married.rds

Cluster script

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

Cluster script

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

Model summary

Full summary

model_summary = summary(model, use_cache = FALSE, priors = TRUE)
print(model_summary)
##  Family: poisson (log) 
## Formula: children ~ paternalage + birth_cohort + spouses * male + maternalage.factor + paternalage.mean + paternal_loss + maternal_loss + older_siblings + nr.siblings + last_born + (1 | idParents) 
##    Data: model_data (Number of observations: 1032368) 
## Samples: 6 chains, each with iter = 1500; warmup = 500; thin = 1; 
##          total post-warmup samples = 6000
##    WAIC: Not computed
##  
## Priors: 
## b ~ normal(0,5)
## sd ~ student_t(3, 0, 5)
## 
## Group-Level Effects: 
## ~idParents (Number of levels: 719113) 
##               Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
## sd(Intercept)        0         0        0        0       2848    1
## 
## Population-Level Effects: 
##                        Estimate Est.Error l-95% CI u-95% CI Eff.Sample
## Intercept                  0.57      0.01     0.56     0.58       6000
## paternalage               -0.04      0.00    -0.05    -0.03       2828
## birth_cohort1950M1955      0.06      0.00     0.05     0.06       6000
## birth_cohort1955M1960      0.08      0.00     0.08     0.08       6000
## spouses                    0.10      0.00     0.09     0.10       6000
## male1                     -0.04      0.00    -0.05    -0.03       6000
## maternalage.factor1420     0.03      0.00     0.02     0.03       6000
## maternalage.factor3561     0.01      0.00     0.00     0.01       6000
## paternalage.mean           0.03      0.00     0.02     0.04       2939
## paternal_loss01            0.10      0.06    -0.02     0.22       6000
## paternal_loss15            0.02      0.02    -0.02     0.06       6000
## paternal_loss510           0.00      0.01    -0.02     0.01       6000
## paternal_loss1015          0.01      0.01    -0.01     0.02       6000
## paternal_loss1520          0.00      0.00     0.00     0.01       6000
## paternal_loss2025          0.00      0.00    -0.01     0.01       6000
## paternal_loss2530          0.00      0.00    -0.01     0.00       6000
## paternal_loss3035          0.00      0.00    -0.01     0.00       6000
## paternal_loss3540          0.00      0.00    -0.01     0.00       6000
## paternal_loss4045          0.00      0.00    -0.01     0.00       6000
## paternal_lossunclear      -0.01      0.00    -0.02    -0.01       6000
## maternal_loss01           -0.07      0.08    -0.23     0.09       6000
## maternal_loss15           -0.01      0.03    -0.07     0.05       6000
## maternal_loss510           0.02      0.01    -0.01     0.05       6000
## maternal_loss1015          0.00      0.01    -0.02     0.02       6000
## maternal_loss1520          0.00      0.01    -0.01     0.02       6000
## maternal_loss2025          0.01      0.01     0.00     0.02       6000
## maternal_loss2530          0.01      0.00     0.00     0.01       6000
## maternal_loss3035          0.00      0.00     0.00     0.01       6000
## maternal_loss3540          0.00      0.00    -0.01     0.01       6000
## maternal_loss4045          0.00      0.00    -0.01     0.00       6000
## maternal_lossunclear       0.00      0.00     0.00     0.00       6000
## older_siblings1            0.01      0.00     0.00     0.01       3530
## older_siblings2            0.01      0.00     0.00     0.02       2955
## older_siblings3            0.01      0.01     0.00     0.02       2952
## older_siblings4            0.00      0.01    -0.01     0.01       3123
## older_siblings5P          -0.04      0.01    -0.06    -0.02       2899
## nr.siblings                0.03      0.00     0.03     0.03       2892
## last_born1                 0.00      0.00     0.00     0.00       6000
## spouses:male1              0.05      0.00     0.04     0.05       6000
##                        Rhat
## Intercept                 1
## paternalage               1
## birth_cohort1950M1955     1
## birth_cohort1955M1960     1
## spouses                   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
## spouses:male1             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.5711 0.5604 0.5819 1.77 1.751 1.79
paternalage -0.03901 -0.04754 -0.03024 0.9617 0.9536 0.9702
birth_cohort1950M1955 0.05592 0.05242 0.05938 1.058 1.054 1.061
birth_cohort1955M1960 0.0809 0.0771 0.08458 1.084 1.08 1.088
spouses 0.09545 0.09188 0.09909 1.1 1.096 1.104
male1 -0.03903 -0.04604 -0.03197 0.9617 0.955 0.9685
maternalage.factor1420 0.02926 0.0242 0.03429 1.03 1.024 1.035
maternalage.factor3561 0.006941 0.002121 0.01153 1.007 1.002 1.012
paternalage.mean 0.0287 0.01978 0.03729 1.029 1.02 1.038
paternal_loss01 0.09964 -0.01695 0.2171 1.105 0.9832 1.242
paternal_loss15 0.02306 -0.01833 0.06327 1.023 0.9818 1.065
paternal_loss510 -0.003866 -0.02293 0.01481 0.9961 0.9773 1.015
paternal_loss1015 0.006088 -0.005324 0.01744 1.006 0.9947 1.018
paternal_loss1520 0.004145 -0.004451 0.01254 1.004 0.9956 1.013
paternal_loss2025 -0.001725 -0.008904 0.00526 0.9983 0.9911 1.005
paternal_loss2530 -0.001666 -0.007554 0.004273 0.9983 0.9925 1.004
paternal_loss3035 -0.0004285 -0.005792 0.00496 0.9996 0.9942 1.005
paternal_loss3540 -0.00177 -0.006531 0.003036 0.9982 0.9935 1.003
paternal_loss4045 -0.004514 -0.009496 0.0004539 0.9955 0.9905 1
paternal_lossunclear -0.01315 -0.01687 -0.009354 0.9869 0.9833 0.9907
maternal_loss01 -0.06821 -0.2295 0.08862 0.9341 0.795 1.093
maternal_loss15 -0.009866 -0.07458 0.05302 0.9902 0.9281 1.054
maternal_loss510 0.01917 -0.01029 0.0481 1.019 0.9898 1.049
maternal_loss1015 -0.0003591 -0.01748 0.01665 0.9996 0.9827 1.017
maternal_loss1520 0.003117 -0.01029 0.01631 1.003 0.9898 1.016
maternal_loss2025 0.00693 -0.004044 0.01781 1.007 0.996 1.018
maternal_loss2530 0.005332 -0.003586 0.01443 1.005 0.9964 1.015
maternal_loss3035 0.004706 -0.002797 0.01232 1.005 0.9972 1.012
maternal_loss3540 -0.0009199 -0.007219 0.005389 0.9991 0.9928 1.005
maternal_loss4045 -0.001168 -0.00732 0.004776 0.9988 0.9927 1.005
maternal_lossunclear -0.001524 -0.004782 0.001752 0.9985 0.9952 1.002
older_siblings1 0.008765 0.004494 0.01314 1.009 1.005 1.013
older_siblings2 0.01224 0.004909 0.01972 1.012 1.005 1.02
older_siblings3 0.008292 -0.002362 0.01937 1.008 0.9976 1.02
older_siblings4 0.0002728 -0.01413 0.01469 1 0.986 1.015
older_siblings5P -0.04201 -0.06013 -0.02359 0.9589 0.9416 0.9767
nr.siblings 0.03043 0.02868 0.03217 1.031 1.029 1.033
last_born1 0.001151 -0.00187 0.004116 1.001 0.9981 1.004
spouses:male1 0.04699 0.0417 0.05223 1.048 1.043 1.054

Paternal age effect

pander::pander(paternal_age_10y_effect(model))
effect median_estimate ci_95 ci_80
estimate father 25y 2.0521 [2.0416;2.0629] [2.0452;2.0591]
estimate father 35y 1.9735 [1.958;1.9892] [1.9635;1.9837]
percentage change -3.8239 [-4.6432;-2.9784] [-4.3831;-3.2768]
OR/IRR 0.9618 [0.9536;0.9702] [0.9562;0.9672]

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)

plot of chunk unnamed-chunk-51plot of chunk unnamed-chunk-51plot of chunk unnamed-chunk-51plot of chunk unnamed-chunk-51plot of chunk unnamed-chunk-51plot of chunk unnamed-chunk-51plot of chunk unnamed-chunk-51plot of chunk unnamed-chunk-51plot of chunk unnamed-chunk-51plot of chunk unnamed-chunk-51plot of chunk unnamed-chunk-51plot of chunk unnamed-chunk-51

Coefficient plot

Coefficient estimates (95% and 80% credibility).

mcmc_intervals(as.matrix(model$fit), regex_pars = "b_[^I]")
plot of chunk unnamed-chunk-52

plot of chunk unnamed-chunk-52

mcmc_areas(as.matrix(model$fit), regex_pars = "b")
plot of chunk unnamed-chunk-52

plot of chunk unnamed-chunk-52

Diagnostics

These plots were made to diagnose misfit and nonconvergence.

Posterior predictive checks

brms::pp_check(model, re_formula = NA, type = "dens_overlay")
plot of chunk unnamed-chunk-53

plot of chunk unnamed-chunk-53

brms::pp_check(model, re_formula = NA, type = "hist")
plot of chunk unnamed-chunk-53

plot of chunk unnamed-chunk-53

Rhat

Did the 6 chains converge?

stanplot(model, pars = "^b_[^I]", type = 'rhat')
plot of chunk unnamed-chunk-54

plot of chunk unnamed-chunk-54

Effective sample size by average sample size

stanplot(model, pars = "^b", type = 'ess')
plot of chunk unnamed-chunk-55

plot of chunk unnamed-chunk-55

Monte Carlo SE

stanplot(model, pars = "^b", type = 'mcse')
plot of chunk unnamed-chunk-56

plot of chunk unnamed-chunk-56

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')
plot of chunk unnamed-chunk-58

plot of chunk unnamed-chunk-58

File name

coefs/swed/e4_children.rds

Cluster script

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

Cluster script

e5: Probability of ever divorcing

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 summary

Full summary

model_summary = summary(model, use_cache = FALSE, priors = TRUE)
print(model_summary)
##  Family: bernoulli (cauchit) 
## Formula: ever_divorced ~ paternalage + spouses + 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: 1032368) 
## 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: 719113) 
##               Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
## sd(Intercept)     0.81      0.01     0.79     0.83        427 1.01
## 
## Population-Level Effects: 
##                        Estimate Est.Error l-95% CI u-95% CI Eff.Sample
## Intercept                 -0.79      0.03    -0.84    -0.74       3000
## paternalage               -0.03      0.02    -0.07     0.02       1573
## spouses                   14.53      0.20    14.15    14.91       3000
## birth_cohort1950M1955      0.04      0.01     0.02     0.05       3000
## birth_cohort1955M1960      0.14      0.01     0.12     0.16       3000
## male1                     -0.03      0.01    -0.05    -0.02       3000
## maternalage.factor1420     0.16      0.01     0.13     0.18       3000
## maternalage.factor3561    -0.11      0.01    -0.14    -0.09       3000
## paternalage.mean          -0.17      0.02    -0.21    -0.12       1566
## paternal_loss01           -0.21      0.34    -0.94     0.37       3000
## paternal_loss15            0.35      0.09     0.17     0.52       3000
## paternal_loss510           0.37      0.04     0.29     0.45       3000
## paternal_loss1015          0.28      0.03     0.23     0.33       3000
## paternal_loss1520          0.24      0.02     0.20     0.28       3000
## paternal_loss2025          0.18      0.02     0.15     0.22       3000
## paternal_loss2530          0.15      0.02     0.12     0.18       3000
## paternal_loss3035          0.13      0.01     0.10     0.16       3000
## paternal_loss3540          0.08      0.01     0.06     0.10       3000
## paternal_loss4045          0.06      0.01     0.04     0.09       3000
## paternal_lossunclear       0.08      0.01     0.06     0.09       3000
## maternal_loss01            0.55      0.27    -0.02     1.07       3000
## maternal_loss15            0.06      0.14    -0.23     0.32       3000
## maternal_loss510           0.22      0.07     0.08     0.34       3000
## maternal_loss1015          0.16      0.04     0.08     0.25       3000
## maternal_loss1520          0.19      0.03     0.13     0.25       3000
## maternal_loss2025          0.15      0.03     0.09     0.20       3000
## maternal_loss2530          0.08      0.02     0.03     0.12       3000
## maternal_loss3035          0.13      0.02     0.09     0.16       3000
## maternal_loss3540          0.07      0.02     0.04     0.10       3000
## maternal_loss4045          0.04      0.02     0.02     0.07       3000
## maternal_lossunclear       0.01      0.01     0.00     0.03       3000
## older_siblings1           -0.03      0.01    -0.05    -0.01       1818
## older_siblings2           -0.03      0.02    -0.06     0.01       1650
## older_siblings3            0.02      0.03    -0.04     0.07       1571
## older_siblings4            0.04      0.04    -0.04     0.11       1608
## older_siblings5P           0.06      0.05    -0.03     0.15       1696
## nr.siblings                0.00      0.00     0.00     0.01       1599
## last_born1                -0.02      0.01    -0.04    -0.01       3000
##                        Rhat
## Intercept                 1
## paternalage               1
## spouses                   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.792 -0.8434 -0.7425 0.4529 0.4302 0.4759
paternalage -0.0293 -0.0727 0.01541 0.9711 0.9299 1.016
spouses 14.53 14.15 14.91 2043904 1402661 2995616
birth_cohort1950M1955 0.03555 0.01715 0.05409 1.036 1.017 1.056
birth_cohort1955M1960 0.1374 0.1184 0.1564 1.147 1.126 1.169
male1 -0.03175 -0.04514 -0.01885 0.9688 0.9559 0.9813
maternalage.factor1420 0.1576 0.1343 0.1806 1.171 1.144 1.198
maternalage.factor3561 -0.1132 -0.138 -0.08789 0.893 0.8711 0.9159
paternalage.mean -0.1681 -0.2128 -0.1241 0.8453 0.8083 0.8833
paternal_loss01 -0.2132 -0.9382 0.3722 0.808 0.3913 1.451
paternal_loss15 0.3472 0.1689 0.5188 1.415 1.184 1.68
paternal_loss510 0.3734 0.2929 0.4523 1.453 1.34 1.572
paternal_loss1015 0.2788 0.2262 0.332 1.322 1.254 1.394
paternal_loss1520 0.2358 0.1975 0.2771 1.266 1.218 1.319
paternal_loss2025 0.1842 0.1487 0.2186 1.202 1.16 1.244
paternal_loss2530 0.1463 0.1162 0.1771 1.158 1.123 1.194
paternal_loss3035 0.1325 0.1044 0.1591 1.142 1.11 1.173
paternal_loss3540 0.07995 0.05651 0.1044 1.083 1.058 1.11
paternal_loss4045 0.06469 0.03976 0.09017 1.067 1.041 1.094
paternal_lossunclear 0.07539 0.05729 0.09315 1.078 1.059 1.098
maternal_loss01 0.5481 -0.02177 1.073 1.73 0.9785 2.923
maternal_loss15 0.05638 -0.2298 0.322 1.058 0.7947 1.38
maternal_loss510 0.2155 0.08332 0.3425 1.24 1.087 1.408
maternal_loss1015 0.1645 0.08161 0.2455 1.179 1.085 1.278
maternal_loss1520 0.1937 0.1308 0.2543 1.214 1.14 1.29
maternal_loss2025 0.1461 0.09409 0.1977 1.157 1.099 1.219
maternal_loss2530 0.07823 0.03424 0.1211 1.081 1.035 1.129
maternal_loss3035 0.1254 0.08864 0.1612 1.134 1.093 1.175
maternal_loss3540 0.06587 0.03505 0.09715 1.068 1.036 1.102
maternal_loss4045 0.04474 0.01503 0.07471 1.046 1.015 1.078
maternal_lossunclear 0.01157 -0.004306 0.0279 1.012 0.9957 1.028
older_siblings1 -0.02717 -0.04861 -0.005686 0.9732 0.9526 0.9943
older_siblings2 -0.02658 -0.06438 0.009751 0.9738 0.9376 1.01
older_siblings3 0.01572 -0.03778 0.06871 1.016 0.9629 1.071
older_siblings4 0.03694 -0.0357 0.109 1.038 0.9649 1.115
older_siblings5P 0.05653 -0.0337 0.1486 1.058 0.9669 1.16
nr.siblings 0.004679 -0.004383 0.0137 1.005 0.9956 1.014
last_born1 -0.02484 -0.03866 -0.01092 0.9755 0.9621 0.9891

Paternal age effect

pander::pander(paternal_age_10y_effect(model))
effect median_estimate ci_95 ci_80
estimate father 25y 0.1999 [0.1968;0.2032] [0.1978;0.202]
estimate father 35y 0.1967 [0.1924;0.2014] [0.1938;0.1998]
percentage change -0.3226 [-0.7798;0.1685] [-0.619;-0.0029]
OR/IRR 0.9705 [0.9299;1.0155] [0.9442;0.9997]

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)

plot of chunk unnamed-chunk-64plot of chunk unnamed-chunk-64plot of chunk unnamed-chunk-64plot of chunk unnamed-chunk-64plot of chunk unnamed-chunk-64plot of chunk unnamed-chunk-64plot of chunk unnamed-chunk-64plot of chunk unnamed-chunk-64plot of chunk unnamed-chunk-64plot of chunk unnamed-chunk-64plot of chunk unnamed-chunk-64

Coefficient plot

Coefficient estimates (95% and 80% credibility).

mcmc_intervals(as.matrix(model$fit), regex_pars = "b_[^I]")
plot of chunk unnamed-chunk-65

plot of chunk unnamed-chunk-65

mcmc_areas(as.matrix(model$fit), regex_pars = "b")
plot of chunk unnamed-chunk-65

plot of chunk unnamed-chunk-65

Diagnostics

These plots were made to diagnose misfit and nonconvergence.

Posterior predictive checks

brms::pp_check(model, re_formula = NA, type = "dens_overlay")
plot of chunk unnamed-chunk-66

plot of chunk unnamed-chunk-66

brms::pp_check(model, re_formula = NA, type = "hist")
plot of chunk unnamed-chunk-66

plot of chunk unnamed-chunk-66

Rhat

Did the 6 chains converge?

stanplot(model, pars = "^b_[^I]", type = 'rhat')
plot of chunk unnamed-chunk-67

plot of chunk unnamed-chunk-67

Effective sample size by average sample size

stanplot(model, pars = "^b", type = 'ess')
plot of chunk unnamed-chunk-68

plot of chunk unnamed-chunk-68

Monte Carlo SE

stanplot(model, pars = "^b", type = 'mcse')
plot of chunk unnamed-chunk-69

plot of chunk unnamed-chunk-69

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')
plot of chunk unnamed-chunk-71

plot of chunk unnamed-chunk-71

File name

coefs/swed/e5_divorce.rds

Cluster script

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

Cluster script

All episodes

Here we show the effect of paternal age for each episode.