Québec (Registre de la Population du Québec Ancien 1630-1750), 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, "rpqa")
}
# options for each chunk calling knit_child
opts_chunk$set(warning=FALSE, message = FALSE, echo = FALSE)

Analysis description

Data subset

The rpqa.1 dataset contains only those participants where paternal age is known and the birthdate is between 1630 and 1750.

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: 61493) 
## 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: 11940) 
##               Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
## sd(Intercept)     1.06      0.03        1     1.12        623    1
## 
## Population-Level Effects: 
##                        Estimate Est.Error l-95% CI u-95% CI Eff.Sample
## Intercept                  3.74      0.25     3.26     4.25        398
## paternalage               -0.44      0.07    -0.58    -0.31        789
## birth_cohort1675M1680      0.20      0.30    -0.37     0.77        627
## birth_cohort1680M1685     -0.57      0.24    -1.07    -0.08        435
## birth_cohort1685M1690     -1.46      0.22    -1.90    -1.04        349
## birth_cohort1690M1695     -0.84      0.22    -1.29    -0.40        366
## birth_cohort1695M1700     -0.83      0.22    -1.27    -0.42        349
## birth_cohort1700M1705     -1.30      0.21    -1.75    -0.89        320
## birth_cohort1705M1710     -0.78      0.22    -1.22    -0.36        349
## birth_cohort1710M1715     -1.43      0.21    -1.87    -1.03        317
## birth_cohort1715M1720     -1.29      0.21    -1.71    -0.89        315
## birth_cohort1720M1725     -1.28      0.21    -1.71    -0.90        314
## birth_cohort1725M1730     -1.65      0.20    -2.07    -1.26        309
## birth_cohort1730M1735     -1.87      0.20    -2.28    -1.47        201
## birth_cohort1735M1740     -1.67      0.20    -2.09    -1.28        312
## male1                     -0.42      0.03    -0.49    -0.36       3000
## maternalage.factor1420    -0.18      0.07    -0.32    -0.04       2048
## maternalage.factor3550    -0.09      0.06    -0.19     0.02       1510
## paternalage.mean           0.54      0.07     0.40     0.69        856
## paternal_loss01           -1.02      0.11    -1.24    -0.79       1413
## paternal_loss15           -0.47      0.10    -0.67    -0.28       1287
## paternal_loss510          -0.35      0.09    -0.53    -0.18       1164
## paternal_loss1015         -0.23      0.08    -0.39    -0.06       1142
## paternal_loss1520         -0.33      0.08    -0.47    -0.18       1045
## paternal_loss2025         -0.16      0.07    -0.30    -0.01       1055
## paternal_loss2530         -0.25      0.07    -0.38    -0.11       1119
## paternal_loss3035         -0.18      0.07    -0.32    -0.05       1131
## paternal_loss3540         -0.16      0.07    -0.29    -0.03       1195
## paternal_loss4045         -0.03      0.08    -0.18     0.12       1565
## paternal_lossunclear      -0.36      0.08    -0.50    -0.21       1140
## maternal_loss01           -2.20      0.10    -2.41    -2.00        996
## maternal_loss15           -0.70      0.09    -0.87    -0.53        994
## maternal_loss510          -0.44      0.08    -0.60    -0.28       1063
## maternal_loss1015         -0.41      0.08    -0.56    -0.25       1092
## maternal_loss1520         -0.39      0.08    -0.53    -0.23       1326
## maternal_loss2025         -0.27      0.08    -0.42    -0.12       1145
## maternal_loss2530         -0.17      0.07    -0.31    -0.03       1188
## maternal_loss3035         -0.21      0.07    -0.34    -0.08       1007
## maternal_loss3540         -0.24      0.06    -0.36    -0.12       1219
## maternal_loss4045         -0.03      0.07    -0.16     0.11       1549
## maternal_lossunclear      -0.27      0.08    -0.42    -0.13       1348
## older_siblings1            0.45      0.06     0.33     0.56       1665
## older_siblings2            0.54      0.06     0.42     0.66       1354
## older_siblings3            0.79      0.07     0.65     0.94       1186
## older_siblings4            0.83      0.08     0.68     0.98       1071
## older_siblings5P           1.02      0.09     0.84     1.20        912
## nr.siblings               -0.09      0.01    -0.11    -0.08       1267
## last_born1                 0.00      0.05    -0.10     0.10       3000
##                        Rhat
## Intercept              1.02
## paternalage            1.00
## birth_cohort1675M1680  1.01
## birth_cohort1680M1685  1.02
## birth_cohort1685M1690  1.02
## birth_cohort1690M1695  1.02
## birth_cohort1695M1700  1.02
## birth_cohort1700M1705  1.02
## birth_cohort1705M1710  1.02
## birth_cohort1710M1715  1.02
## birth_cohort1715M1720  1.02
## birth_cohort1720M1725  1.03
## birth_cohort1725M1730  1.03
## birth_cohort1730M1735  1.03
## birth_cohort1735M1740  1.03
## male1                  1.00
## maternalage.factor1420 1.00
## maternalage.factor3550 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
## paternal_lossunclear   1.00
## maternal_loss01        1.00
## maternal_loss15        1.01
## maternal_loss510       1.00
## maternal_loss1015      1.01
## maternal_loss1520      1.00
## maternal_loss2025      1.01
## maternal_loss2530      1.00
## maternal_loss3035      1.01
## maternal_loss3540      1.00
## maternal_loss4045      1.00
## maternal_lossunclear   1.00
## older_siblings1        1.00
## older_siblings2        1.00
## older_siblings3        1.00
## older_siblings4        1.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	42.3	25.95	69.83
paternalage	0.643	0.5615	0.7354
birth_cohort1675M1680	1.216	0.6914	2.169
birth_cohort1680M1685	0.5651	0.3431	0.9193
birth_cohort1685M1690	0.2332	0.1501	0.3541
birth_cohort1690M1695	0.4322	0.2765	0.6711
birth_cohort1695M1700	0.435	0.2812	0.6595
birth_cohort1700M1705	0.2719	0.1746	0.4095
birth_cohort1705M1710	0.4602	0.2951	0.6991
birth_cohort1710M1715	0.2394	0.1543	0.3556
birth_cohort1715M1720	0.2765	0.1808	0.4088
birth_cohort1720M1725	0.2774	0.1807	0.4085
birth_cohort1725M1730	0.1919	0.1263	0.2825
birth_cohort1730M1735	0.1548	0.1023	0.2302
birth_cohort1735M1740	0.189	0.1238	0.2783
male1	0.6559	0.6155	0.696
maternalage.factor1420	0.837	0.7272	0.9602
maternalage.factor3550	0.9167	0.824	1.024
paternalage.mean	1.721	1.491	1.985
paternal_loss01	0.3594	0.2895	0.4522
paternal_loss15	0.6225	0.5097	0.7536
paternal_loss510	0.7041	0.589	0.8371
paternal_loss1015	0.7978	0.6802	0.9426
paternal_loss1520	0.7197	0.6219	0.8363
paternal_loss2025	0.8561	0.7422	0.9897
paternal_loss2530	0.7824	0.6809	0.8994
paternal_loss3035	0.8313	0.7246	0.952
paternal_loss3540	0.8524	0.7461	0.9733
paternal_loss4045	0.9682	0.8331	1.124
paternal_lossunclear	0.7006	0.6048	0.8104
maternal_loss01	0.1109	0.09014	0.1348
maternal_loss15	0.4959	0.4178	0.5909
maternal_loss510	0.6449	0.5486	0.758
maternal_loss1015	0.6664	0.5726	0.7758
maternal_loss1520	0.6804	0.5876	0.7912
maternal_loss2025	0.7636	0.6577	0.891
maternal_loss2530	0.8454	0.7305	0.9682
maternal_loss3035	0.8086	0.7136	0.9186
maternal_loss3540	0.7885	0.7008	0.8899
maternal_loss4045	0.9697	0.8537	1.111
maternal_lossunclear	0.761	0.656	0.8802
older_siblings1	1.562	1.392	1.75
older_siblings2	1.718	1.518	1.944
older_siblings3	2.211	1.924	2.561
older_siblings4	2.29	1.964	2.67
older_siblings5P	2.78	2.326	3.321
nr.siblings	0.9103	0.8977	0.9226
last_born1	0.9954	0.9039	1.102

Paternal age effect

pander::pander(paternal_age_10y_effect(model))

effect	median_estimate	ci_95	ci_80
estimate father 25y	0.92	[0.91;0.93]	[0.91;0.92]
estimate father 35y	0.91	[0.9;0.92]	[0.9;0.91]
percentage change	-1.03	[-1.51;-0.67]	[-1.32;-0.78]
OR/IRR	0.64	[0.56;0.74]	[0.59;0.7]

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/rpqa/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: 49246) 
## 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: 11383) 
##               Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
## sd(Intercept)      1.1      0.04     1.01     1.19        957    1
## 
## Population-Level Effects: 
##                        Estimate Est.Error l-95% CI u-95% CI Eff.Sample
## Intercept                  5.29      0.47     4.44     6.26        430
## paternalage                0.09      0.12    -0.15     0.31       1143
## birth_cohort1675M1680     -1.23      0.48    -2.27    -0.36        518
## birth_cohort1680M1685     -1.80      0.45    -2.77    -0.99        443
## birth_cohort1685M1690     -1.91      0.45    -2.88    -1.10        396
## birth_cohort1690M1695     -1.26      0.46    -2.25    -0.43        483
## birth_cohort1695M1700     -1.77      0.44    -2.73    -0.98        392
## birth_cohort1700M1705     -2.00      0.44    -2.98    -1.23        409
## birth_cohort1705M1710     -1.82      0.44    -2.78    -1.05        391
## birth_cohort1710M1715     -2.51      0.43    -3.46    -1.76        381
## birth_cohort1715M1720     -1.81      0.44    -2.77    -1.05        431
## birth_cohort1720M1725     -2.36      0.43    -3.31    -1.59        376
## birth_cohort1725M1730     -3.01      0.43    -3.96    -2.27        400
## birth_cohort1730M1735     -2.66      0.43    -3.62    -1.91        379
## birth_cohort1735M1740     -2.16      0.43    -3.11    -1.40        378
## male1                     -0.22      0.05    -0.32    -0.13       3000
## maternalage.factor1420     0.08      0.12    -0.15     0.33       3000
## maternalage.factor3550     0.10      0.08    -0.07     0.26       3000
## paternalage.mean          -0.06      0.12    -0.30     0.19       1071
## paternal_loss01           -0.90      0.18    -1.24    -0.56       3000
## paternal_loss15           -0.55      0.15    -0.83    -0.27       1955
## paternal_loss510          -0.71      0.13    -0.96    -0.45       1660
## paternal_loss1015         -0.40      0.13    -0.65    -0.15       1733
## paternal_loss1520         -0.49      0.12    -0.72    -0.26       1792
## paternal_loss2025         -0.41      0.11    -0.64    -0.18       1698
## paternal_loss2530         -0.32      0.11    -0.54    -0.10       1669
## paternal_loss3035         -0.15      0.12    -0.38     0.07       1628
## paternal_loss3540         -0.03      0.12    -0.26     0.20       1936
## paternal_loss4045          0.05      0.13    -0.20     0.30       2013
## paternal_lossunclear      -0.48      0.11    -0.70    -0.27       1691
## maternal_loss01           -0.60      0.19    -0.96    -0.22       3000
## maternal_loss15           -0.85      0.12    -1.08    -0.61       3000
## maternal_loss510          -0.52      0.11    -0.74    -0.29       2516
## maternal_loss1015         -0.52      0.11    -0.72    -0.30       3000
## maternal_loss1520         -0.22      0.12    -0.45     0.02       3000
## maternal_loss2025         -0.02      0.12    -0.25     0.22       3000
## maternal_loss2530         -0.09      0.11    -0.30     0.13       3000
## maternal_loss3035         -0.14      0.10    -0.34     0.06       3000
## maternal_loss3540         -0.12      0.10    -0.30     0.08       3000
## maternal_loss4045         -0.21      0.10    -0.40    -0.01       3000
## maternal_lossunclear      -0.11      0.11    -0.33     0.11       3000
## older_siblings1           -0.06      0.09    -0.25     0.12       2082
## older_siblings2           -0.09      0.10    -0.29     0.11       1437
## older_siblings3           -0.09      0.11    -0.31     0.14       1371
## older_siblings4           -0.07      0.12    -0.31     0.18       1104
## older_siblings5P          -0.07      0.15    -0.37     0.24        804
## nr.siblings                0.01      0.01    -0.02     0.03       1260
## last_born1                -0.16      0.08    -0.31    -0.01       3000
##                        Rhat
## Intercept              1.01
## paternalage            1.00
## birth_cohort1675M1680  1.01
## birth_cohort1680M1685  1.01
## birth_cohort1685M1690  1.01
## birth_cohort1690M1695  1.01
## birth_cohort1695M1700  1.01
## birth_cohort1700M1705  1.01
## birth_cohort1705M1710  1.01
## birth_cohort1710M1715  1.01
## birth_cohort1715M1720  1.01
## birth_cohort1720M1725  1.01
## birth_cohort1725M1730  1.01
## birth_cohort1730M1735  1.01
## birth_cohort1735M1740  1.01
## male1                  1.00
## maternalage.factor1420 1.00
## maternalage.factor3550 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
## paternal_lossunclear   1.00
## maternal_loss01        1.00
## maternal_loss15        1.00
## maternal_loss510       1.00
## maternal_loss1015      1.00
## maternal_loss1520      1.00
## maternal_loss2025      1.00
## maternal_loss2530      1.00
## maternal_loss3035      1.00
## maternal_loss3540      1.00
## maternal_loss4045      1.00
## maternal_lossunclear   1.00
## older_siblings1        1.00
## older_siblings2        1.00
## older_siblings3        1.00
## older_siblings4        1.00
## 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	197.6	84.56	523.3
paternalage	1.092	0.858	1.37
birth_cohort1675M1680	0.2936	0.1038	0.7002
birth_cohort1680M1685	0.1657	0.0624	0.3715
birth_cohort1685M1690	0.1483	0.05632	0.3334
birth_cohort1690M1695	0.2837	0.1051	0.6515
birth_cohort1695M1700	0.1699	0.06522	0.375
birth_cohort1700M1705	0.1357	0.05104	0.291
birth_cohort1705M1710	0.1618	0.06179	0.3492
birth_cohort1710M1715	0.08141	0.03152	0.172
birth_cohort1715M1720	0.1632	0.06272	0.3483
birth_cohort1720M1725	0.09458	0.03658	0.2036
birth_cohort1725M1730	0.04918	0.01898	0.1029
birth_cohort1730M1735	0.07026	0.02684	0.1475
birth_cohort1735M1740	0.1158	0.04447	0.2454
male1	0.7998	0.7291	0.8792
maternalage.factor1420	1.086	0.8583	1.392
maternalage.factor3550	1.103	0.9343	1.298
paternalage.mean	0.9464	0.7442	1.215
paternal_loss01	0.4059	0.2891	0.5736
paternal_loss15	0.5786	0.4351	0.7667
paternal_loss510	0.4923	0.3848	0.6375
paternal_loss1015	0.6685	0.5218	0.8586
paternal_loss1520	0.6127	0.4861	0.7736
paternal_loss2025	0.6608	0.5275	0.8324
paternal_loss2530	0.7239	0.5841	0.9043
paternal_loss3035	0.8594	0.6829	1.075
paternal_loss3540	0.9694	0.7681	1.216
paternal_loss4045	1.051	0.816	1.348
paternal_lossunclear	0.6159	0.4953	0.7626
maternal_loss01	0.5502	0.3839	0.8058
maternal_loss15	0.4291	0.3406	0.5413
maternal_loss510	0.5973	0.479	0.7448
maternal_loss1015	0.5972	0.4857	0.7439
maternal_loss1520	0.8056	0.6399	1.016
maternal_loss2025	0.9781	0.7755	1.25
maternal_loss2530	0.9136	0.7386	1.137
maternal_loss3035	0.8669	0.7134	1.067
maternal_loss3540	0.8909	0.7383	1.082
maternal_loss4045	0.8129	0.6697	0.987
maternal_lossunclear	0.8946	0.7188	1.112
older_siblings1	0.9374	0.7794	1.127
older_siblings2	0.9105	0.75	1.112
older_siblings3	0.9153	0.7357	1.147
older_siblings4	0.9299	0.7318	1.193
older_siblings5P	0.9351	0.6875	1.276
nr.siblings	1.006	0.9837	1.028
last_born1	0.8514	0.7318	0.9902

Paternal age effect

pander::pander(paternal_age_10y_effect(model))

effect	median_estimate	ci_95	ci_80
estimate father 25y	0.94	[0.93;0.95]	[0.93;0.95]
estimate father 35y	0.94	[0.93;0.95]	[0.94;0.95]
percentage change	0.09	[-0.17;0.36]	[-0.07;0.26]
OR/IRR	1.09	[0.86;1.37]	[0.94;1.27]

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/rpqa/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: 49479) 
## 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: 11467) 
##               Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
## sd(Intercept)     0.78      0.02     0.73     0.82        917    1
## 
## Population-Level Effects: 
##                        Estimate Est.Error l-95% CI u-95% CI Eff.Sample
## Intercept                  1.84      0.13     1.59     2.09       1046
## paternalage                0.00      0.06    -0.11     0.11       1325
## birth_cohort1675M1680     -0.06      0.09    -0.24     0.10        845
## birth_cohort1680M1685     -0.03      0.08    -0.20     0.14        764
## birth_cohort1685M1690      0.05      0.09    -0.13     0.23        831
## birth_cohort1690M1695      0.21      0.09     0.03     0.38        666
## birth_cohort1695M1700      0.28      0.09     0.12     0.46        639
## birth_cohort1700M1705      0.31      0.09     0.14     0.48        629
## birth_cohort1705M1710      0.25      0.08     0.09     0.41        622
## birth_cohort1710M1715      0.19      0.08     0.03     0.36        611
## birth_cohort1715M1720      0.20      0.08     0.04     0.35        521
## birth_cohort1720M1725      0.27      0.08     0.12     0.43        555
## birth_cohort1725M1730     -0.04      0.08    -0.18     0.11        511
## birth_cohort1730M1735     -0.07      0.08    -0.21     0.08        508
## birth_cohort1735M1740     -0.05      0.07    -0.19     0.10        522
## male1                     -0.86      0.03    -0.91    -0.80       3000
## maternalage.factor1420    -0.04      0.06    -0.15     0.08       3000
## maternalage.factor3550    -0.01      0.04    -0.10     0.07       3000
## paternalage.mean          -0.03      0.06    -0.15     0.08       1388
## paternal_loss01           -0.39      0.11    -0.60    -0.18       3000
## paternal_loss15           -0.32      0.08    -0.47    -0.17       1513
## paternal_loss510          -0.26      0.07    -0.40    -0.12       1434
## paternal_loss1015         -0.17      0.07    -0.30    -0.03       1453
## paternal_loss1520         -0.35      0.06    -0.47    -0.23       1195
## paternal_loss2025         -0.21      0.06    -0.33    -0.09       1246
## paternal_loss2530         -0.18      0.06    -0.29    -0.07       1215
## paternal_loss3035         -0.22      0.06    -0.33    -0.11       1253
## paternal_loss3540         -0.20      0.05    -0.31    -0.09       1382
## paternal_loss4045         -0.12      0.06    -0.24    -0.01       1954
## paternal_lossunclear      -0.54      0.06    -0.65    -0.43       1185
## maternal_loss01           -0.63      0.11    -0.84    -0.42       3000
## maternal_loss15           -0.15      0.07    -0.30     0.00       2116
## maternal_loss510          -0.23      0.06    -0.36    -0.11       2011
## maternal_loss1015         -0.26      0.06    -0.37    -0.14       2009
## maternal_loss1520         -0.23      0.06    -0.35    -0.12       1990
## maternal_loss2025         -0.17      0.06    -0.28    -0.05       2028
## maternal_loss2530         -0.10      0.06    -0.21     0.02       2127
## maternal_loss3035         -0.05      0.05    -0.15     0.06       2073
## maternal_loss3540         -0.05      0.05    -0.15     0.04       3000
## maternal_loss4045          0.06      0.05    -0.05     0.16       3000
## maternal_lossunclear      -0.36      0.05    -0.46    -0.26       1623
## older_siblings1           -0.08      0.05    -0.17     0.01       2051
## older_siblings2           -0.12      0.05    -0.22    -0.02       1669
## older_siblings3           -0.07      0.06    -0.19     0.04       1398
## older_siblings4           -0.09      0.06    -0.21     0.04       1232
## older_siblings5P          -0.08      0.08    -0.23     0.07        821
## nr.siblings                0.02      0.01     0.00     0.03       1923
## last_born1                 0.04      0.04    -0.05     0.12       3000
##                        Rhat
## Intercept              1.01
## paternalage            1.01
## birth_cohort1675M1680  1.01
## birth_cohort1680M1685  1.01
## birth_cohort1685M1690  1.01
## birth_cohort1690M1695  1.01
## birth_cohort1695M1700  1.01
## birth_cohort1700M1705  1.01
## birth_cohort1705M1710  1.01
## birth_cohort1710M1715  1.01
## birth_cohort1715M1720  1.02
## birth_cohort1720M1725  1.01
## birth_cohort1725M1730  1.01
## birth_cohort1730M1735  1.02
## birth_cohort1735M1740  1.01
## male1                  1.00
## maternalage.factor1420 1.00
## maternalage.factor3550 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.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	6.309	4.914	8.084
paternalage	0.9961	0.8938	1.114
birth_cohort1675M1680	0.938	0.79	1.104
birth_cohort1680M1685	0.9697	0.8197	1.145
birth_cohort1685M1690	1.049	0.8811	1.254
birth_cohort1690M1695	1.235	1.032	1.458
birth_cohort1695M1700	1.329	1.124	1.577
birth_cohort1700M1705	1.357	1.153	1.61
birth_cohort1705M1710	1.287	1.095	1.514
birth_cohort1710M1715	1.212	1.029	1.433
birth_cohort1715M1720	1.217	1.036	1.423
birth_cohort1720M1725	1.311	1.123	1.533
birth_cohort1725M1730	0.9642	0.8322	1.117
birth_cohort1730M1735	0.931	0.8068	1.081
birth_cohort1735M1740	0.9537	0.8262	1.105
male1	0.425	0.4018	0.4478
maternalage.factor1420	0.9631	0.8607	1.08
maternalage.factor3550	0.987	0.9079	1.074
paternalage.mean	0.9692	0.8624	1.086
paternal_loss01	0.6778	0.5491	0.8364
paternal_loss15	0.726	0.6235	0.8466
paternal_loss510	0.7745	0.6724	0.8902
paternal_loss1015	0.8466	0.7407	0.9672
paternal_loss1520	0.7057	0.627	0.7934
paternal_loss2025	0.8128	0.7223	0.9176
paternal_loss2530	0.8373	0.7454	0.9354
paternal_loss3035	0.7991	0.7164	0.893
paternal_loss3540	0.8198	0.7344	0.9152
paternal_loss4045	0.8837	0.7861	0.9935
paternal_lossunclear	0.5809	0.5195	0.6474
maternal_loss01	0.5318	0.4296	0.656
maternal_loss15	0.8624	0.7437	1.003
maternal_loss510	0.793	0.6996	0.8974
maternal_loss1015	0.7709	0.689	0.8658
maternal_loss1520	0.7936	0.706	0.889
maternal_loss2025	0.8464	0.7576	0.9521
maternal_loss2530	0.9085	0.812	1.015
maternal_loss3035	0.9506	0.8574	1.061
maternal_loss3540	0.9479	0.8627	1.045
maternal_loss4045	1.058	0.9536	1.171
maternal_lossunclear	0.6982	0.6339	0.7737
older_siblings1	0.9209	0.8409	1.005
older_siblings2	0.888	0.801	0.9824
older_siblings3	0.9308	0.8298	1.038
older_siblings4	0.9155	0.8079	1.037
older_siblings5P	0.9269	0.7943	1.074
nr.siblings	1.016	1.005	1.028
last_born1	1.038	0.9559	1.132

Paternal age effect

pander::pander(paternal_age_10y_effect(model))

effect	median_estimate	ci_95	ci_80
estimate father 25y	0.84	[0.83;0.85]	[0.83;0.85]
estimate father 35y	0.84	[0.83;0.86]	[0.83;0.85]
percentage change	-0.02	[-0.83;0.74]	[-0.55;0.46]
OR/IRR	1	[0.89;1.11]	[0.93;1.07]

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/rpqa/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: 35351) 
## 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: 10511) 
##               Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
## sd(Intercept)     0.36         0     0.35     0.37        845 1.01
## 
## Population-Level Effects: 
##                        Estimate Est.Error l-95% CI u-95% CI Eff.Sample
## Intercept                  1.96      0.03     1.89     2.02        933
## paternalage                0.01      0.01    -0.01     0.03       1552
## birth_cohort1675M1680      0.04      0.02     0.01     0.07       1408
## birth_cohort1680M1685      0.05      0.02     0.02     0.08       1172
## birth_cohort1685M1690      0.05      0.02     0.02     0.08       1019
## birth_cohort1690M1695      0.05      0.02     0.01     0.08        842
## birth_cohort1695M1700      0.03      0.02     0.00     0.07        780
## birth_cohort1700M1705      0.01      0.02    -0.02     0.04        779
## birth_cohort1705M1710     -0.03      0.02    -0.06     0.01        728
## birth_cohort1710M1715     -0.04      0.02    -0.08    -0.01        752
## birth_cohort1715M1720     -0.07      0.02    -0.11    -0.04        672
## birth_cohort1720M1725     -0.08      0.02    -0.11    -0.05        684
## birth_cohort1725M1730     -0.10      0.02    -0.14    -0.07        695
## birth_cohort1730M1735     -0.11      0.02    -0.14    -0.08        677
## birth_cohort1735M1740     -0.11      0.02    -0.15    -0.08        704
## spouses                    0.05      0.01     0.04     0.06       3000
## male1                     -0.12      0.01    -0.14    -0.10       3000
## maternalage.factor1420     0.00      0.01    -0.02     0.02       3000
## maternalage.factor3550     0.01      0.01    -0.01     0.02       3000
## paternalage.mean          -0.03      0.01    -0.05     0.00       1430
## paternal_loss01           -0.04      0.02    -0.09     0.00       1375
## paternal_loss15           -0.02      0.02    -0.06     0.01        990
## paternal_loss510          -0.02      0.02    -0.05     0.01        817
## paternal_loss1015         -0.02      0.01    -0.04     0.01        816
## paternal_loss1520         -0.01      0.01    -0.04     0.01        848
## paternal_loss2025         -0.02      0.01    -0.05     0.00        867
## paternal_loss2530         -0.01      0.01    -0.04     0.01        975
## paternal_loss3035         -0.02      0.01    -0.04     0.00       1125
## paternal_loss3540         -0.02      0.01    -0.04     0.00       1241
## paternal_loss4045         -0.01      0.01    -0.04     0.01       1750
## paternal_lossunclear      -0.07      0.01    -0.10    -0.05       1055
## maternal_loss01           -0.02      0.02    -0.07     0.03       2250
## maternal_loss15           -0.02      0.02    -0.05     0.01       1388
## maternal_loss510           0.03      0.01     0.00     0.05       1266
## maternal_loss1015          0.00      0.01    -0.02     0.03       1328
## maternal_loss1520          0.01      0.01    -0.01     0.04       1222
## maternal_loss2025         -0.03      0.01    -0.05    -0.01       1378
## maternal_loss2530          0.01      0.01    -0.01     0.03       1534
## maternal_loss3035          0.00      0.01    -0.02     0.02       1628
## maternal_loss3540          0.01      0.01    -0.01     0.03       1883
## maternal_loss4045          0.02      0.01     0.00     0.03       3000
## maternal_lossunclear      -0.02      0.02    -0.05     0.01       1158
## older_siblings1           -0.03      0.01    -0.04    -0.01       3000
## older_siblings2           -0.02      0.01    -0.03     0.00       2275
## older_siblings3           -0.03      0.01    -0.05    -0.01       2090
## older_siblings4           -0.03      0.01    -0.05     0.00       2041
## older_siblings5P          -0.04      0.01    -0.07    -0.01       1801
## nr.siblings                0.01      0.00     0.01     0.01       1348
## last_born1                 0.00      0.01    -0.01     0.02       3000
## spouses:male1              0.18      0.01     0.16     0.20       3000
##                        Rhat
## Intercept              1.01
## paternalage            1.00
## birth_cohort1675M1680  1.00
## birth_cohort1680M1685  1.00
## birth_cohort1685M1690  1.00
## birth_cohort1690M1695  1.00
## birth_cohort1695M1700  1.00
## birth_cohort1700M1705  1.00
## birth_cohort1705M1710  1.00
## birth_cohort1710M1715  1.00
## birth_cohort1715M1720  1.00
## birth_cohort1720M1725  1.00
## birth_cohort1725M1730  1.00
## birth_cohort1730M1735  1.00
## birth_cohort1735M1740  1.00
## spouses                1.00
## male1                  1.00
## maternalage.factor1420 1.00
## maternalage.factor3550 1.00
## paternalage.mean       1.00
## paternal_loss01        1.00
## paternal_loss15        1.00
## paternal_loss510       1.01
## paternal_loss1015      1.00
## paternal_loss1520      1.00
## paternal_loss2025      1.00
## paternal_loss2530      1.00
## paternal_loss3035      1.00
## paternal_loss3540      1.00
## paternal_loss4045      1.00
## paternal_lossunclear   1.00
## maternal_loss01        1.00
## maternal_loss15        1.00
## maternal_loss510       1.00
## maternal_loss1015      1.00
## maternal_loss1520      1.00
## maternal_loss2025      1.00
## maternal_loss2530      1.00
## maternal_loss3035      1.00
## maternal_loss3540      1.00
## maternal_loss4045      1.00
## maternal_lossunclear   1.00
## older_siblings1        1.00
## older_siblings2        1.00
## older_siblings3        1.00
## older_siblings4        1.00
## older_siblings5P       1.00
## nr.siblings            1.00
## 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	7.09	6.648	7.535
paternalage	1.01	0.9868	1.032
birth_cohort1675M1680	1.037	1.006	1.069
birth_cohort1680M1685	1.054	1.02	1.087
birth_cohort1685M1690	1.051	1.016	1.086
birth_cohort1690M1695	1.048	1.014	1.083
birth_cohort1695M1700	1.035	1.001	1.069
birth_cohort1700M1705	1.009	0.9765	1.042
birth_cohort1705M1710	0.9751	0.9437	1.008
birth_cohort1710M1715	0.959	0.9275	0.9915
birth_cohort1715M1720	0.9309	0.8998	0.9628
birth_cohort1720M1725	0.9229	0.893	0.9534
birth_cohort1725M1730	0.9007	0.871	0.9309
birth_cohort1730M1735	0.8955	0.8655	0.9262
birth_cohort1735M1740	0.8915	0.8619	0.9224
spouses	1.051	1.037	1.065
male1	0.8866	0.8665	0.9077
maternalage.factor1420	0.9961	0.9765	1.016
maternalage.factor3550	1.007	0.9904	1.024
paternalage.mean	0.9726	0.9491	0.9977
paternal_loss01	0.9591	0.9167	1.003
paternal_loss15	0.9791	0.9423	1.015
paternal_loss510	0.9814	0.9508	1.013
paternal_loss1015	0.9848	0.9573	1.014
paternal_loss1520	0.9855	0.9583	1.013
paternal_loss2025	0.9797	0.955	1.004
paternal_loss2530	0.988	0.9648	1.011
paternal_loss3035	0.979	0.9583	0.9995
paternal_loss3540	0.984	0.9653	1.003
paternal_loss4045	0.986	0.9656	1.006
paternal_lossunclear	0.9286	0.9017	0.9555
maternal_loss01	0.9811	0.9357	1.029
maternal_loss15	0.9798	0.9484	1.013
maternal_loss510	1.026	0.9958	1.055
maternal_loss1015	1.002	0.9756	1.029
maternal_loss1520	1.012	0.9856	1.038
maternal_loss2025	0.9696	0.9471	0.9941
maternal_loss2530	1.01	0.9885	1.034
maternal_loss3035	1	0.9806	1.02
maternal_loss3540	1.012	0.994	1.03
maternal_loss4045	1.017	0.9992	1.034
maternal_lossunclear	0.9773	0.949	1.007
older_siblings1	0.9733	0.9583	0.9878
older_siblings2	0.9828	0.9662	1.001
older_siblings3	0.9707	0.9514	0.9896
older_siblings4	0.9748	0.9529	0.9962
older_siblings5P	0.9582	0.9323	0.9854
nr.siblings	1.009	1.007	1.012
last_born1	1.003	0.9879	1.019
spouses:male1	1.2	1.179	1.22

Paternal age effect

pander::pander(paternal_age_10y_effect(model))

effect	median_estimate	ci_95	ci_80
estimate father 25y	7.5	[7.25;7.76]	[7.33;7.67]
estimate father 35y	7.57	[7.28;7.88]	[7.38;7.78]
percentage change	0.94	[-1.32;3.23]	[-0.5;2.47]
OR/IRR	1.01	[0.99;1.03]	[1;1.02]

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/rpqa/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.