Demographic data base (Sweden 1760-1880), 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, "ddb")
}
# options for each chunk calling knit_child
opts_chunk$set(warning=FALSE, message = FALSE, echo = FALSE)

Analysis description

Data subset

The ddb.1 dataset contains only those participants where paternal age is known and the birthdate is between 1760 and 1880.

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: 56010) 
## 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: 14708) 
##               Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
## sd(Intercept)     1.41      0.04     1.33     1.49        989    1
## 
## Population-Level Effects: 
##                        Estimate Est.Error l-95% CI u-95% CI Eff.Sample
## Intercept                  5.11      0.62     3.99     6.40        336
## paternalage               -1.23      0.11    -1.44    -1.01       1110
## birth_cohort1750M1755     -0.18      1.21    -2.19     2.55        783
## birth_cohort1755M1760     -1.04      0.76    -2.50     0.47        510
## birth_cohort1760M1765     -1.40      0.68    -2.79    -0.14        400
## birth_cohort1765M1770     -0.27      0.79    -1.78     1.26        503
## birth_cohort1770M1775     -0.40      0.78    -1.91     1.11        493
## birth_cohort1775M1780     -1.94      0.62    -3.27    -0.81        339
## birth_cohort1780M1785     -1.98      0.62    -3.28    -0.87        338
## birth_cohort1785M1790     -1.68      0.62    -2.99    -0.52        338
## birth_cohort1790M1795     -1.83      0.61    -3.15    -0.74        319
## birth_cohort1795M1800     -1.43      0.61    -2.75    -0.35        317
## birth_cohort1800M1805     -1.48      0.60    -2.77    -0.36        319
## birth_cohort1805M1810     -1.67      0.60    -2.98    -0.56        316
## birth_cohort1810M1815     -1.77      0.60    -3.06    -0.66        320
## birth_cohort1815M1820     -1.74      0.59    -2.99    -0.65        317
## birth_cohort1820M1825     -1.31      0.60    -2.58    -0.24        312
## birth_cohort1825M1830     -1.18      0.60    -2.47    -0.09        315
## birth_cohort1830M1835     -1.33      0.60    -2.61    -0.23        314
## birth_cohort1835M1840     -1.25      0.59    -2.54    -0.17        318
## birth_cohort1840M1845     -0.91      0.60    -2.18     0.19        316
## birth_cohort1845M1850     -0.80      0.60    -2.05     0.29        317
## male1                     -0.37      0.05    -0.46    -0.29       3000
## maternalage.factor1020    -0.09      0.20    -0.46     0.30       3000
## maternalage.factor3559    -0.21      0.07    -0.34    -0.08       3000
## paternalage.mean           1.26      0.11     1.04     1.49       1171
## paternal_loss01           -0.99      0.16    -1.31    -0.67       2106
## paternal_loss15           -0.46      0.13    -0.70    -0.21       1576
## paternal_loss510          -0.40      0.11    -0.62    -0.18       1539
## paternal_loss1015         -0.30      0.11    -0.51    -0.09       1469
## paternal_loss1520         -0.20      0.10    -0.40     0.00       1433
## paternal_loss2025         -0.08      0.10    -0.27     0.11       1646
## paternal_loss2530         -0.07      0.10    -0.26     0.12       1521
## paternal_loss3035         -0.14      0.09    -0.31     0.04       1609
## paternal_loss3540          0.00      0.09    -0.18     0.19       1686
## paternal_loss4045          0.09      0.10    -0.10     0.29       3000
## maternal_loss01           -3.45      0.14    -3.74    -3.18       2162
## maternal_loss15           -1.21      0.12    -1.43    -0.98       1565
## maternal_loss510          -0.94      0.11    -1.15    -0.73       1779
## maternal_loss1015         -0.55      0.11    -0.76    -0.34       2147
## maternal_loss1520         -0.44      0.11    -0.65    -0.23       2268
## maternal_loss2025         -0.44      0.10    -0.63    -0.24       2175
## maternal_loss2530         -0.45      0.09    -0.63    -0.27       1952
## maternal_loss3035         -0.31      0.09    -0.48    -0.14       2124
## maternal_loss3540         -0.17      0.08    -0.33     0.00       2218
## maternal_loss4045         -0.08      0.09    -0.26     0.10       3000
## older_siblings1            0.60      0.07     0.46     0.74       3000
## older_siblings2            0.92      0.09     0.74     1.09       1599
## older_siblings3            1.05      0.10     0.85     1.26       1383
## older_siblings4            1.39      0.13     1.15     1.65       1249
## older_siblings5P           1.99      0.17     1.66     2.33       1097
## nr.siblings               -0.25      0.01    -0.28    -0.22       1300
## last_born1                -0.08      0.06    -0.19     0.04       3000
##                        Rhat
## Intercept              1.02
## paternalage            1.00
## birth_cohort1750M1755  1.01
## birth_cohort1755M1760  1.01
## birth_cohort1760M1765  1.01
## birth_cohort1765M1770  1.01
## birth_cohort1770M1775  1.01
## birth_cohort1775M1780  1.02
## birth_cohort1780M1785  1.02
## birth_cohort1785M1790  1.02
## birth_cohort1790M1795  1.02
## birth_cohort1795M1800  1.02
## birth_cohort1800M1805  1.02
## birth_cohort1805M1810  1.02
## birth_cohort1810M1815  1.02
## birth_cohort1815M1820  1.02
## birth_cohort1820M1825  1.02
## birth_cohort1825M1830  1.02
## birth_cohort1830M1835  1.02
## birth_cohort1835M1840  1.02
## birth_cohort1840M1845  1.02
## birth_cohort1845M1850  1.02
## male1                  1.00
## maternalage.factor1020 1.00
## maternalage.factor3559 1.00
## paternalage.mean       1.00
## 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
## 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
## older_siblings1        1.00
## older_siblings2        1.00
## older_siblings3        1.00
## older_siblings4        1.00
## 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

Estimates are exp(b). When they are referring to the hurdle (hu) component, or a dichotomous outcome, they are odds ratios, when they are referring to a Poisson component, they are hazard ratios. In both cases, they are presented with 95% credibility intervals. To see the effects on the response scale (probability or number of children), consult the marginal effect plots.

fixed_eff = data.frame(model_summary$fixed, check.names = F)
fixed_eff$Est.Error = fixed_eff$Eff.Sample = fixed_eff$Rhat = NULL
fixed_eff$`Odds/hazard ratio` = exp(fixed_eff$Estimate)
fixed_eff$`OR/HR low 95%` = exp(fixed_eff$`l-95% CI`)
fixed_eff$`OR/HR high 95%` = exp(fixed_eff$`u-95% CI`)
fixed_eff = fixed_eff %>% select(`Odds/hazard ratio`, `OR/HR low 95%`, `OR/HR high 95%`)
pander::pander(fixed_eff)

	Odds/hazard ratio	OR/HR low 95%	OR/HR high 95%
Intercept	165.5	54.03	604.2
paternalage	0.2925	0.2365	0.363
birth_cohort1750M1755	0.8326	0.1115	12.83
birth_cohort1755M1760	0.3518	0.08236	1.597
birth_cohort1760M1765	0.2463	0.06114	0.8664
birth_cohort1765M1770	0.766	0.1688	3.509
birth_cohort1770M1775	0.6694	0.1481	3.031
birth_cohort1775M1780	0.1433	0.03789	0.4443
birth_cohort1780M1785	0.1378	0.0376	0.4195
birth_cohort1785M1790	0.1864	0.05028	0.5929
birth_cohort1790M1795	0.1598	0.04271	0.4777
birth_cohort1795M1800	0.24	0.06378	0.708
birth_cohort1800M1805	0.2266	0.06297	0.6947
birth_cohort1805M1810	0.1888	0.05069	0.5739
birth_cohort1810M1815	0.1708	0.04702	0.5163
birth_cohort1815M1820	0.1749	0.05021	0.5219
birth_cohort1820M1825	0.2688	0.07543	0.7895
birth_cohort1825M1830	0.3065	0.08468	0.9121
birth_cohort1830M1835	0.2642	0.07365	0.7955
birth_cohort1835M1840	0.2863	0.07922	0.8417
birth_cohort1840M1845	0.403	0.1136	1.208
birth_cohort1845M1850	0.4515	0.1292	1.332
male1	0.6889	0.6293	0.7517
maternalage.factor1020	0.9131	0.6281	1.354
maternalage.factor3559	0.8082	0.7114	0.9186
paternalage.mean	3.534	2.823	4.447
paternal_loss01	0.3711	0.2701	0.5097
paternal_loss15	0.633	0.495	0.8143
paternal_loss510	0.6719	0.537	0.8384
paternal_loss1015	0.7412	0.601	0.9113
paternal_loss1520	0.8209	0.6737	0.9966
paternal_loss2025	0.9232	0.7631	1.119
paternal_loss2530	0.9343	0.7735	1.127
paternal_loss3035	0.8733	0.7318	1.045
paternal_loss3540	0.9994	0.8389	1.204
paternal_loss4045	1.096	0.905	1.34
maternal_loss01	0.03174	0.02387	0.04157
maternal_loss15	0.2994	0.2389	0.3762
maternal_loss510	0.3901	0.316	0.481
maternal_loss1015	0.5754	0.4688	0.7126
maternal_loss1520	0.6413	0.5208	0.7947
maternal_loss2025	0.6464	0.5344	0.7838
maternal_loss2530	0.637	0.53	0.7657
maternal_loss3035	0.73	0.6164	0.8724
maternal_loss3540	0.843	0.7158	0.9971
maternal_loss4045	0.9251	0.7731	1.103
older_siblings1	1.828	1.588	2.1
older_siblings2	2.497	2.091	2.96
older_siblings3	2.863	2.348	3.522
older_siblings4	4.024	3.148	5.188
older_siblings5P	7.311	5.282	10.31
nr.siblings	0.7765	0.7548	0.7988
last_born1	0.9269	0.8283	1.039

Paternal age effect

pander::pander(paternal_age_10y_effect(model))

effect	median_estimate	ci_95	ci_80
estimate father 25y	0.94	[0.92;0.95]	[0.93;0.95]
estimate father 35y	0.92	[0.89;0.94]	[0.9;0.93]
percentage change	-1.82	[-3.14;-1.08]	[-2.63;-1.28]
OR/IRR	0.29	[0.24;0.36]	[0.25;0.34]

Marginal effect plots

In these marginal effect plots, we set all predictors except the one shown on the X axis to their mean and in the case of factors to their reference level. We then plot the estimated association between the X axis predictor and the outcome on the response scale (e.g. probability of survival/marriage or number of children).

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

Here, we plotted the 95% posterior densities for the unexponentiated model coefficients (b_). The darkly shaded area represents the 50% credibility interval, the dark line represent the posterior mean estimate.

mcmc_areas(as.matrix(model$fit), regex_pars = "b_[^I]", point_est = "mean", prob = 0.50, prob_outer = 0.95) + ggtitle("Posterior densities with means and 50% intervals") + analysis_theme + theme(axis.text = element_text(size = 12), panel.grid = element_blank()) + xlab("Coefficient size")

Diagnostics

These plots were made to diagnose misfit and nonconvergence.

Posterior predictive checks

In posterior predictive checks, we test whether we can approximately reproduce the real data distribution from our model.

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 over average sample size

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

Trace plots

Trace plots are only shown in the case of nonconvergence.

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

File/cluster script name

This model was stored in the file: coefs/ddb/e1_survive1y.rds.

Click the following link to see the script used to generate this model:

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: 45540) 
## 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: 13947) 
##               Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
## sd(Intercept)     1.26      0.06     1.14     1.38        572    1
## 
## Population-Level Effects: 
##                        Estimate Est.Error l-95% CI u-95% CI Eff.Sample
## Intercept                  4.88      0.71     3.73     6.67         89
## paternalage               -0.12      0.16    -0.43     0.20        682
## birth_cohort1750M1755     -0.23      1.02    -2.24     1.82        180
## birth_cohort1755M1760      0.52      1.03    -1.47     2.68        179
## birth_cohort1760M1765      0.04      0.83    -1.77     1.59        125
## birth_cohort1765M1770      0.84      0.98    -1.08     2.83        157
## birth_cohort1770M1775     -0.61      0.74    -2.30     0.65        103
## birth_cohort1775M1780     -0.74      0.74    -2.50     0.48         93
## birth_cohort1780M1785     -0.43      0.74    -2.17     0.85         99
## birth_cohort1785M1790     -0.81      0.71    -2.57     0.32         88
## birth_cohort1790M1795     -1.15      0.68    -2.87    -0.09         82
## birth_cohort1795M1800     -0.61      0.68    -2.39     0.48         80
## birth_cohort1800M1805     -0.60      0.68    -2.34     0.44         80
## birth_cohort1805M1810     -0.82      0.68    -2.54     0.23         80
## birth_cohort1810M1815     -0.51      0.67    -2.23     0.54         79
## birth_cohort1815M1820      0.00      0.67    -1.72     1.05         80
## birth_cohort1820M1825      0.35      0.67    -1.34     1.40         80
## birth_cohort1825M1830     -0.13      0.67    -1.79     0.89         79
## birth_cohort1830M1835     -0.01      0.67    -1.73     1.03         80
## birth_cohort1835M1840      0.16      0.67    -1.58     1.19         80
## birth_cohort1840M1845      0.41      0.67    -1.34     1.46         80
## birth_cohort1845M1850      0.01      0.67    -1.70     1.07         80
## male1                     -0.35      0.07    -0.48    -0.23       3000
## maternalage.factor1020    -0.27      0.34    -0.89     0.45       3000
## maternalage.factor3559    -0.01      0.09    -0.19     0.16       2292
## paternalage.mean           0.20      0.17    -0.13     0.53        739
## paternal_loss01           -1.87      0.21    -2.28    -1.46       1202
## paternal_loss15           -1.34      0.17    -1.66    -1.00        959
## paternal_loss510          -1.39      0.15    -1.70    -1.11        905
## paternal_loss1015         -1.05      0.15    -1.35    -0.75        924
## paternal_loss1520         -0.73      0.15    -1.03    -0.43        933
## paternal_loss2025         -0.68      0.15    -0.98    -0.39        924
## paternal_loss2530         -0.69      0.14    -0.98    -0.40        850
## paternal_loss3035         -0.44      0.15    -0.73    -0.14        908
## paternal_loss3540         -0.22      0.15    -0.53     0.09       1001
## paternal_loss4045         -0.35      0.17    -0.67    -0.02       1110
## maternal_loss01           -2.01      0.24    -2.46    -1.51       1737
## maternal_loss15           -1.57      0.16    -1.87    -1.26       1682
## maternal_loss510          -1.19      0.15    -1.47    -0.90       1336
## maternal_loss1015         -1.00      0.14    -1.27    -0.71       1336
## maternal_loss1520         -0.84      0.14    -1.12    -0.57       1311
## maternal_loss2025         -0.64      0.14    -0.91    -0.36       1450
## maternal_loss2530         -0.62      0.14    -0.87    -0.34       1412
## maternal_loss3035         -0.53      0.13    -0.77    -0.28       1376
## maternal_loss3540         -0.42      0.12    -0.66    -0.17       1592
## maternal_loss4045         -0.23      0.14    -0.50     0.04       1638
## older_siblings1           -0.23      0.11    -0.45     0.00       1301
## older_siblings2           -0.26      0.13    -0.53     0.00        852
## older_siblings3           -0.25      0.16    -0.57     0.06        719
## older_siblings4           -0.31      0.20    -0.70     0.07        696
## older_siblings5P          -0.10      0.25    -0.59     0.40        617
## nr.siblings               -0.08      0.02    -0.12    -0.04        875
## last_born1                -0.03      0.08    -0.19     0.13       3000
##                        Rhat
## Intercept              1.05
## paternalage            1.00
## birth_cohort1750M1755  1.03
## birth_cohort1755M1760  1.03
## birth_cohort1760M1765  1.04
## birth_cohort1765M1770  1.03
## birth_cohort1770M1775  1.04
## birth_cohort1775M1780  1.05
## birth_cohort1780M1785  1.05
## birth_cohort1785M1790  1.05
## birth_cohort1790M1795  1.06
## birth_cohort1795M1800  1.06
## birth_cohort1800M1805  1.06
## birth_cohort1805M1810  1.06
## birth_cohort1810M1815  1.06
## birth_cohort1815M1820  1.06
## birth_cohort1820M1825  1.06
## birth_cohort1825M1830  1.06
## birth_cohort1830M1835  1.06
## birth_cohort1835M1840  1.06
## birth_cohort1840M1845  1.06
## birth_cohort1845M1850  1.06
## male1                  1.00
## maternalage.factor1020 1.00
## maternalage.factor3559 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.01
## paternal_loss3540      1.01
## paternal_loss4045      1.01
## 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
## older_siblings1        1.00
## older_siblings2        1.00
## older_siblings3        1.00
## older_siblings4        1.01
## older_siblings5P       1.01
## 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$`Odds/hazard ratio` = exp(fixed_eff$Estimate)
fixed_eff$`OR/HR low 95%` = exp(fixed_eff$`l-95% CI`)
fixed_eff$`OR/HR high 95%` = exp(fixed_eff$`u-95% CI`)
fixed_eff = fixed_eff %>% select(`Odds/hazard ratio`, `OR/HR low 95%`, `OR/HR high 95%`)
pander::pander(fixed_eff)

	Odds/hazard ratio	OR/HR low 95%	OR/HR high 95%
Intercept	131.4	41.64	785.2
paternalage	0.8868	0.649	1.227
birth_cohort1750M1755	0.7912	0.1066	6.172
birth_cohort1755M1760	1.684	0.2301	14.61
birth_cohort1760M1765	1.046	0.1704	4.918
birth_cohort1765M1770	2.32	0.3406	16.93
birth_cohort1770M1775	0.5407	0.09981	1.909
birth_cohort1775M1780	0.4752	0.08211	1.623
birth_cohort1780M1785	0.6485	0.1145	2.35
birth_cohort1785M1790	0.4436	0.07638	1.379
birth_cohort1790M1795	0.3169	0.05651	0.9105
birth_cohort1795M1800	0.5411	0.09169	1.615
birth_cohort1800M1805	0.549	0.09602	1.558
birth_cohort1805M1810	0.4386	0.07876	1.264
birth_cohort1810M1815	0.5993	0.1072	1.722
birth_cohort1815M1820	1	0.1796	2.853
birth_cohort1820M1825	1.425	0.2611	4.059
birth_cohort1825M1830	0.8802	0.1667	2.437
birth_cohort1830M1835	0.9892	0.1771	2.807
birth_cohort1835M1840	1.179	0.2063	3.29
birth_cohort1840M1845	1.514	0.262	4.292
birth_cohort1845M1850	1.013	0.1831	2.905
male1	0.7025	0.6178	0.7983
maternalage.factor1020	0.7601	0.409	1.571
maternalage.factor3559	0.9926	0.828	1.179
paternalage.mean	1.223	0.8804	1.701
paternal_loss01	0.1538	0.102	0.2324
paternal_loss15	0.2621	0.1895	0.3686
paternal_loss510	0.2483	0.1832	0.33
paternal_loss1015	0.3505	0.2588	0.4712
paternal_loss1520	0.4837	0.3567	0.6522
paternal_loss2025	0.5051	0.3742	0.6782
paternal_loss2530	0.5009	0.3752	0.6704
paternal_loss3035	0.6465	0.4834	0.8715
paternal_loss3540	0.8002	0.5887	1.09
paternal_loss4045	0.7043	0.5123	0.9795
maternal_loss01	0.1345	0.08518	0.2198
maternal_loss15	0.2091	0.1544	0.2832
maternal_loss510	0.3041	0.2303	0.4049
maternal_loss1015	0.3697	0.2821	0.4934
maternal_loss1520	0.4303	0.3279	0.5664
maternal_loss2025	0.5288	0.4034	0.6959
maternal_loss2530	0.538	0.4172	0.7084
maternal_loss3035	0.59	0.4632	0.7546
maternal_loss3540	0.6539	0.5145	0.8413
maternal_loss4045	0.7922	0.6072	1.04
older_siblings1	0.7963	0.6377	1.001
older_siblings2	0.7733	0.591	1.001
older_siblings3	0.7775	0.563	1.065
older_siblings4	0.7315	0.4957	1.072
older_siblings5P	0.9027	0.554	1.485
nr.siblings	0.9234	0.885	0.9644
last_born1	0.9683	0.8251	1.141

Paternal age effect

pander::pander(paternal_age_10y_effect(model))

effect	median_estimate	ci_95	ci_80
estimate father 25y	0.93	[0.92;0.95]	[0.92;0.95]
estimate father 35y	0.93	[0.91;0.95]	[0.92;0.94]
percentage change	-0.16	[-0.68;0.27]	[-0.47;0.12]
OR/IRR	0.89	[0.65;1.23]	[0.71;1.09]

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

mcmc_areas(as.matrix(model$fit), regex_pars = "b_[^I]", point_est = "mean", prob = 0.50, prob_outer = 0.95) + ggtitle("Posterior densities with means and 50% intervals") + analysis_theme + theme(axis.text = element_text(size = 12), panel.grid = element_blank()) + xlab("Coefficient size")

Diagnostics

These plots were made to diagnose misfit and nonconvergence.

Posterior predictive checks

In posterior predictive checks, we test whether we can approximately reproduce the real data distribution from our model.

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 over average sample size

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

Trace plots

Trace plots are only shown in the case of nonconvergence.

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

File/cluster script name

This model was stored in the file: coefs/ddb/e2_surviveR.rds.

Click the following link to see the script used to generate this model:

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: 43048) 
## 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: 13901) 
##               Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
## sd(Intercept)     1.11      0.03     1.06     1.17        922    1
## 
## Population-Level Effects: 
##                        Estimate Est.Error l-95% CI u-95% CI Eff.Sample
## Intercept                  0.96      0.24     0.50     1.43        212
## paternalage                0.42      0.07     0.28     0.56        593
## birth_cohort1750M1755     -0.35      0.32    -1.01     0.27        485
## birth_cohort1755M1760      0.05      0.30    -0.54     0.64        328
## birth_cohort1760M1765      0.21      0.28    -0.32     0.77        283
## birth_cohort1765M1770      0.32      0.28    -0.21     0.87        251
## birth_cohort1770M1775      0.31      0.28    -0.23     0.84        245
## birth_cohort1775M1780      0.39      0.27    -0.14     0.91        254
## birth_cohort1780M1785      0.41      0.28    -0.12     0.98        256
## birth_cohort1785M1790      0.19      0.26    -0.31     0.70        215
## birth_cohort1790M1795     -0.11      0.24    -0.58     0.36        192
## birth_cohort1795M1800      0.18      0.23    -0.28     0.64        188
## birth_cohort1800M1805      0.32      0.23    -0.13     0.78        182
## birth_cohort1805M1810      0.44      0.23     0.01     0.88        181
## birth_cohort1810M1815      0.29      0.23    -0.15     0.73        170
## birth_cohort1815M1820      0.35      0.23    -0.08     0.79        167
## birth_cohort1820M1825      0.41      0.23    -0.02     0.84        161
## birth_cohort1825M1830      0.38      0.23    -0.05     0.81        166
## birth_cohort1830M1835      0.43      0.23     0.00     0.86        173
## birth_cohort1835M1840      0.39      0.23    -0.04     0.83        167
## birth_cohort1840M1845      0.30      0.22    -0.13     0.73        163
## birth_cohort1845M1850      0.28      0.23    -0.14     0.73        172
## male1                     -0.16      0.02    -0.20    -0.11       3000
## maternalage.factor1020     0.05      0.12    -0.18     0.29       2490
## maternalage.factor3559    -0.07      0.04    -0.15     0.00       1723
## paternalage.mean          -0.41      0.07    -0.55    -0.26        628
## paternal_loss01           -0.94      0.12    -1.17    -0.72       1839
## paternal_loss15           -0.73      0.08    -0.88    -0.58       1167
## paternal_loss510          -0.85      0.07    -0.98    -0.71       1158
## paternal_loss1015         -0.64      0.06    -0.76    -0.52       1064
## paternal_loss1520         -0.56      0.06    -0.67    -0.45        879
## paternal_loss2025         -0.46      0.05    -0.57    -0.36       1015
## paternal_loss2530         -0.30      0.05    -0.40    -0.20       1029
## paternal_loss3035         -0.23      0.05    -0.33    -0.13       1121
## paternal_loss3540         -0.20      0.05    -0.29    -0.11       1257
## paternal_loss4045         -0.04      0.05    -0.14     0.06       1515
## maternal_loss01           -1.24      0.15    -1.54    -0.93       3000
## maternal_loss15           -1.07      0.09    -1.24    -0.89       1711
## maternal_loss510          -0.99      0.07    -1.14    -0.85       1511
## maternal_loss1015         -1.00      0.07    -1.12    -0.87       1038
## maternal_loss1520         -0.85      0.06    -0.97    -0.73       1455
## maternal_loss2025         -0.59      0.06    -0.70    -0.48        833
## maternal_loss2530         -0.33      0.05    -0.43    -0.23       1207
## maternal_loss3035         -0.28      0.05    -0.38    -0.19       1218
## maternal_loss3540         -0.18      0.05    -0.27    -0.09        932
## maternal_loss4045         -0.08      0.05    -0.17     0.01       1409
## older_siblings1           -0.06      0.04    -0.14     0.02        936
## older_siblings2           -0.18      0.05    -0.28    -0.08        693
## older_siblings3           -0.25      0.07    -0.39    -0.13        617
## older_siblings4           -0.31      0.08    -0.47    -0.15        625
## older_siblings5P          -0.40      0.11    -0.61    -0.19        558
## nr.siblings                0.01      0.01    -0.01     0.03        596
## last_born1                -0.01      0.03    -0.07     0.06       3000
##                        Rhat
## Intercept              1.02
## paternalage            1.01
## birth_cohort1750M1755  1.01
## birth_cohort1755M1760  1.01
## birth_cohort1760M1765  1.01
## birth_cohort1765M1770  1.02
## birth_cohort1770M1775  1.02
## birth_cohort1775M1780  1.02
## birth_cohort1780M1785  1.02
## birth_cohort1785M1790  1.02
## birth_cohort1790M1795  1.03
## birth_cohort1795M1800  1.03
## birth_cohort1800M1805  1.03
## birth_cohort1805M1810  1.03
## birth_cohort1810M1815  1.03
## birth_cohort1815M1820  1.03
## birth_cohort1820M1825  1.03
## birth_cohort1825M1830  1.03
## birth_cohort1830M1835  1.03
## birth_cohort1835M1840  1.03
## birth_cohort1840M1845  1.03
## birth_cohort1845M1850  1.03
## male1                  1.00
## maternalage.factor1020 1.00
## maternalage.factor3559 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.01
## paternal_loss2025      1.00
## paternal_loss2530      1.01
## paternal_loss3035      1.00
## paternal_loss3540      1.00
## 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.00
## maternal_loss2530      1.00
## maternal_loss3035      1.00
## maternal_loss3540      1.00
## maternal_loss4045      1.00
## older_siblings1        1.00
## older_siblings2        1.00
## older_siblings3        1.01
## older_siblings4        1.00
## 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$`Odds/hazard ratio` = exp(fixed_eff$Estimate)
fixed_eff$`OR/HR low 95%` = exp(fixed_eff$`l-95% CI`)
fixed_eff$`OR/HR high 95%` = exp(fixed_eff$`u-95% CI`)
fixed_eff = fixed_eff %>% select(`Odds/hazard ratio`, `OR/HR low 95%`, `OR/HR high 95%`)
pander::pander(fixed_eff)

	Odds/hazard ratio	OR/HR low 95%	OR/HR high 95%
Intercept	2.614	1.652	4.195
paternalage	1.524	1.323	1.745
birth_cohort1750M1755	0.7012	0.3633	1.306
birth_cohort1755M1760	1.056	0.582	1.903
birth_cohort1760M1765	1.236	0.7284	2.154
birth_cohort1765M1770	1.383	0.8113	2.39
birth_cohort1770M1775	1.369	0.7927	2.315
birth_cohort1775M1780	1.476	0.8734	2.478
birth_cohort1780M1785	1.512	0.8911	2.653
birth_cohort1785M1790	1.208	0.7336	2.022
birth_cohort1790M1795	0.8986	0.5587	1.437
birth_cohort1795M1800	1.196	0.7575	1.896
birth_cohort1800M1805	1.384	0.8755	2.176
birth_cohort1805M1810	1.554	1.013	2.419
birth_cohort1810M1815	1.336	0.8585	2.08
birth_cohort1815M1820	1.422	0.9233	2.202
birth_cohort1820M1825	1.507	0.9776	2.324
birth_cohort1825M1830	1.458	0.9519	2.247
birth_cohort1830M1835	1.539	1	2.374
birth_cohort1835M1840	1.484	0.9626	2.293
birth_cohort1840M1845	1.346	0.8738	2.075
birth_cohort1845M1850	1.326	0.8699	2.068
male1	0.8553	0.8156	0.8989
maternalage.factor1020	1.051	0.8339	1.335
maternalage.factor3559	0.9316	0.8637	1.003
paternalage.mean	0.6648	0.5774	0.769
paternal_loss01	0.3895	0.311	0.4879
paternal_loss15	0.4817	0.4129	0.5621
paternal_loss510	0.427	0.3746	0.4913
paternal_loss1015	0.5264	0.4666	0.5928
paternal_loss1520	0.5691	0.5101	0.6366
paternal_loss2025	0.6319	0.5675	0.7007
paternal_loss2530	0.7413	0.6701	0.8216
paternal_loss3035	0.7952	0.7219	0.8759
paternal_loss3540	0.8168	0.7449	0.8982
paternal_loss4045	0.9583	0.8662	1.059
maternal_loss01	0.2905	0.2137	0.3937
maternal_loss15	0.3443	0.2892	0.4105
maternal_loss510	0.3706	0.3208	0.4289
maternal_loss1015	0.3687	0.3249	0.4185
maternal_loss1520	0.4273	0.3776	0.4832
maternal_loss2025	0.5553	0.4986	0.6176
maternal_loss2530	0.7187	0.6478	0.7938
maternal_loss3035	0.7523	0.6843	0.8307
maternal_loss3540	0.8368	0.7645	0.9145
maternal_loss4045	0.9249	0.8458	1.014
older_siblings1	0.9421	0.8712	1.017
older_siblings2	0.835	0.7549	0.9221
older_siblings3	0.7751	0.6773	0.8818
older_siblings4	0.7332	0.6274	0.8646
older_siblings5P	0.6692	0.5453	0.8234
nr.siblings	1.008	0.9895	1.027
last_born1	0.9924	0.9298	1.059

Paternal age effect

pander::pander(paternal_age_10y_effect(model))

effect	median_estimate	ci_95	ci_80
estimate father 25y	0.68	[0.56;0.76]	[0.6;0.74]
estimate father 35y	0.76	[0.67;0.81]	[0.7;0.8]
percentage change	7.94	[4.47;12.38]	[5.48;10.83]
OR/IRR	1.53	[1.32;1.75]	[1.38;1.67]

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

mcmc_areas(as.matrix(model$fit), regex_pars = "b_[^I]", point_est = "mean", prob = 0.50, prob_outer = 0.95) + ggtitle("Posterior densities with means and 50% intervals") + analysis_theme + theme(axis.text = element_text(size = 12), panel.grid = element_blank()) + xlab("Coefficient size")

Diagnostics

These plots were made to diagnose misfit and nonconvergence.

Posterior predictive checks

In posterior predictive checks, we test whether we can approximately reproduce the real data distribution from our model.

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 over average sample size

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

Trace plots

Trace plots are only shown in the case of nonconvergence.

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

File/cluster script name

This model was stored in the file: coefs/ddb/e3_ever_married.rds.

Click the following link to see the script used to generate this model:

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: 25518) 
## 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: 11129) 
##               Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
## sd(Intercept)     0.45      0.01     0.44     0.46       1596    1
## 
## Population-Level Effects: 
##                        Estimate Est.Error l-95% CI u-95% CI Eff.Sample
## Intercept                  0.83      0.09     0.66     1.00        237
## paternalage               -0.06      0.02    -0.09    -0.02       1114
## birth_cohort1750M1755     -0.12      0.11    -0.34     0.11        980
## birth_cohort1755M1760      0.12      0.09    -0.06     0.30        388
## birth_cohort1760M1765      0.20      0.09     0.03     0.37        275
## birth_cohort1765M1770      0.17      0.09     0.00     0.33        264
## birth_cohort1770M1775      0.11      0.09    -0.07     0.28        274
## birth_cohort1775M1780      0.06      0.09    -0.11     0.24        253
## birth_cohort1780M1785      0.17      0.09    -0.01     0.34        262
## birth_cohort1785M1790      0.18      0.08     0.02     0.34        237
## birth_cohort1790M1795      0.06      0.08    -0.10     0.22        230
## birth_cohort1795M1800      0.02      0.08    -0.13     0.18        217
## birth_cohort1800M1805     -0.01      0.08    -0.16     0.15        211
## birth_cohort1805M1810      0.03      0.08    -0.12     0.19        216
## birth_cohort1810M1815      0.04      0.08    -0.11     0.19        209
## birth_cohort1815M1820      0.13      0.08    -0.02     0.28        204
## birth_cohort1820M1825      0.15      0.08     0.00     0.30        201
## birth_cohort1825M1830      0.12      0.08    -0.02     0.28        201
## birth_cohort1830M1835      0.14      0.08    -0.01     0.29        202
## birth_cohort1835M1840      0.14      0.08    -0.01     0.30        201
## birth_cohort1840M1845      0.12      0.08    -0.03     0.27        200
## birth_cohort1845M1850      0.12      0.08    -0.03     0.27        200
## spouses                    0.24      0.02     0.20     0.28       1821
## male1                      0.03      0.03    -0.02     0.08       1763
## maternalage.factor1020     0.03      0.03    -0.03     0.09       4312
## maternalage.factor3559     0.06      0.01     0.04     0.08       3560
## paternalage.mean           0.06      0.02     0.02     0.10       1176
## paternal_loss01            0.09      0.04     0.02     0.16       2426
## paternal_loss15            0.00      0.03    -0.05     0.05       1592
## paternal_loss510          -0.01      0.02    -0.05     0.04       1627
## paternal_loss1015          0.00      0.02    -0.05     0.03       1546
## paternal_loss1520         -0.08      0.02    -0.11    -0.04       1474
## paternal_loss2025         -0.01      0.02    -0.05     0.02       1421
## paternal_loss2530         -0.01      0.02    -0.04     0.02       1474
## paternal_loss3035          0.00      0.01    -0.03     0.03       1586
## paternal_loss3540          0.03      0.01     0.00     0.06       1643
## paternal_loss4045          0.03      0.01     0.00     0.05       2436
## maternal_loss01            0.06      0.06    -0.05     0.17       2899
## maternal_loss15           -0.01      0.03    -0.07     0.05       2111
## maternal_loss510          -0.03      0.02    -0.08     0.02       2071
## maternal_loss1015         -0.04      0.02    -0.08     0.01       1938
## maternal_loss1520         -0.04      0.02    -0.09    -0.01       2331
## maternal_loss2025         -0.07      0.02    -0.10    -0.03       1750
## maternal_loss2530         -0.02      0.02    -0.05     0.01       1546
## maternal_loss3035         -0.01      0.01    -0.04     0.02       1602
## maternal_loss3540          0.00      0.01    -0.02     0.03       1729
## maternal_loss4045         -0.02      0.01    -0.05     0.00       2808
## older_siblings1            0.01      0.01    -0.02     0.03       2446
## older_siblings2            0.02      0.01    -0.01     0.05       1702
## older_siblings3            0.02      0.02    -0.01     0.06       1334
## older_siblings4            0.01      0.02    -0.03     0.06       1506
## older_siblings5P           0.02      0.03    -0.04     0.07       1381
## nr.siblings                0.02      0.00     0.02     0.03       1114
## last_born1                -0.01      0.01    -0.03     0.01       6000
## spouses:male1              0.07      0.02     0.02     0.11       1743
##                        Rhat
## Intercept              1.01
## paternalage            1.00
## birth_cohort1750M1755  1.00
## birth_cohort1755M1760  1.01
## birth_cohort1760M1765  1.01
## birth_cohort1765M1770  1.01
## birth_cohort1770M1775  1.01
## birth_cohort1775M1780  1.01
## birth_cohort1780M1785  1.01
## birth_cohort1785M1790  1.01
## birth_cohort1790M1795  1.02
## birth_cohort1795M1800  1.02
## birth_cohort1800M1805  1.02
## birth_cohort1805M1810  1.02
## birth_cohort1810M1815  1.02
## birth_cohort1815M1820  1.02
## birth_cohort1820M1825  1.02
## birth_cohort1825M1830  1.02
## birth_cohort1830M1835  1.02
## birth_cohort1835M1840  1.02
## birth_cohort1840M1845  1.02
## birth_cohort1845M1850  1.02
## spouses                1.00
## male1                  1.00
## maternalage.factor1020 1.00
## maternalage.factor3559 1.00
## paternalage.mean       1.00
## 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
## 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
## older_siblings1        1.00
## older_siblings2        1.00
## older_siblings3        1.01
## older_siblings4        1.00
## older_siblings5P       1.01
## nr.siblings            1.01
## last_born1             1.00
## spouses:male1          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$`Odds/hazard ratio` = exp(fixed_eff$Estimate)
fixed_eff$`OR/HR low 95%` = exp(fixed_eff$`l-95% CI`)
fixed_eff$`OR/HR high 95%` = exp(fixed_eff$`u-95% CI`)
fixed_eff = fixed_eff %>% select(`Odds/hazard ratio`, `OR/HR low 95%`, `OR/HR high 95%`)
pander::pander(fixed_eff)

	Odds/hazard ratio	OR/HR low 95%	OR/HR high 95%
Intercept	2.299	1.939	2.722
paternalage	0.9463	0.9108	0.9841
birth_cohort1750M1755	0.8912	0.7116	1.112
birth_cohort1755M1760	1.128	0.9439	1.35
birth_cohort1760M1765	1.227	1.033	1.452
birth_cohort1765M1770	1.179	0.9965	1.395
birth_cohort1770M1775	1.114	0.9361	1.327
birth_cohort1775M1780	1.065	0.8946	1.27
birth_cohort1780M1785	1.184	0.9918	1.408
birth_cohort1785M1790	1.198	1.018	1.412
birth_cohort1790M1795	1.064	0.9055	1.252
birth_cohort1795M1800	1.023	0.8741	1.201
birth_cohort1800M1805	0.9923	0.8524	1.161
birth_cohort1805M1810	1.03	0.8832	1.205
birth_cohort1810M1815	1.041	0.8951	1.213
birth_cohort1815M1820	1.134	0.9789	1.323
birth_cohort1820M1825	1.156	0.9991	1.346
birth_cohort1825M1830	1.132	0.9754	1.32
birth_cohort1830M1835	1.15	0.9912	1.341
birth_cohort1835M1840	1.156	0.9939	1.346
birth_cohort1840M1845	1.122	0.9671	1.308
birth_cohort1845M1850	1.124	0.9673	1.311
spouses	1.269	1.217	1.323
male1	1.031	0.9768	1.087
maternalage.factor1020	1.031	0.9735	1.093
maternalage.factor3559	1.06	1.038	1.084
paternalage.mean	1.062	1.02	1.106
paternal_loss01	1.094	1.016	1.178
paternal_loss15	1	0.9489	1.053
paternal_loss510	0.9934	0.9506	1.039
paternal_loss1015	0.9951	0.9555	1.035
paternal_loss1520	0.9277	0.8934	0.9621
paternal_loss2025	0.9878	0.9547	1.021
paternal_loss2530	0.9907	0.9613	1.022
paternal_loss3035	0.9995	0.9713	1.029
paternal_loss3540	1.03	1.003	1.058
paternal_loss4045	1.026	0.9995	1.053
maternal_loss01	1.063	0.9519	1.187
maternal_loss15	0.9875	0.9302	1.049
maternal_loss510	0.9684	0.9216	1.019
maternal_loss1015	0.9636	0.9209	1.009
maternal_loss1520	0.9563	0.9175	0.9949
maternal_loss2025	0.9346	0.9025	0.9673
maternal_loss2530	0.983	0.953	1.014
maternal_loss3035	0.9872	0.9596	1.016
maternal_loss3540	1.001	0.9758	1.026
maternal_loss4045	0.9789	0.9556	1.003
older_siblings1	1.005	0.9841	1.026
older_siblings2	1.02	0.9923	1.048
older_siblings3	1.025	0.9886	1.062
older_siblings4	1.015	0.9717	1.058
older_siblings5P	1.015	0.9579	1.072
nr.siblings	1.024	1.018	1.03
last_born1	0.9882	0.9701	1.006
spouses:male1	1.068	1.017	1.122

Paternal age effect

pander::pander(paternal_age_10y_effect(model))

effect	median_estimate	ci_95	ci_80
estimate father 25y	3.53	[3.04;4.11]	[3.21;3.91]
estimate father 35y	3.34	[2.87;3.9]	[3.02;3.71]
percentage change	-5.36	[-8.92;-1.59]	[-7.71;-2.96]
OR/IRR	0.95	[0.91;0.98]	[0.92;0.97]

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

mcmc_areas(as.matrix(model$fit), regex_pars = "b_[^I]", point_est = "mean", prob = 0.50, prob_outer = 0.95) + ggtitle("Posterior densities with means and 50% intervals") + analysis_theme + theme(axis.text = element_text(size = 12), panel.grid = element_blank()) + xlab("Coefficient size")

Diagnostics

These plots were made to diagnose misfit and nonconvergence.

Posterior predictive checks

In posterior predictive checks, we test whether we can approximately reproduce the real data distribution from our model.

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 over average sample size

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

Trace plots

Trace plots are only shown in the case of nonconvergence.

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

File/cluster script name

This model was stored in the file: coefs/ddb/e4_children.rds.

Click the following link to see the script used to generate this model:

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.

All episodes

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