Demographic data base (Sweden 1760-1880), main models

Loading details

source("0__helpers.R")
opts_chunk$set(warning=TRUE, cache=F,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

We adjusted for average paternal age within families to isolate the effect of paternal age differences between siblings. We further adjusted for birth cohort in five-year groupings (small groupings at the edge of the range were lumped) to account for secular changes in mortality and fertility, as well as residual censoring. We adjusted for parental deaths in the first 45 years of life to remove effects related to orphanhood and parental senescence (in categories of 0-1, 2-5, 6-10, …, 45+, unknown) for both parents separately. Parental loss at 45+ served as the reference category. We adjusted for maternal age (up to 20, 21-34, 35+), which we binned to reduce multicollinearity with paternal age and to capture nonlinear effects. A maternal age of 21-34 served as the reference category. We also adjusted for number of siblings continuously, number of older siblings (0-5, 5+), and being born last. Being first-born served as the reference category.

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.

Priors

We used weakly informative priors that are documented in each model summary.

m1: No sibling comparison

Here, we ignore the pedigree structure of the data to see whether it matters for the estimation of the paternal age effect.

Model summary

Full summary

model_summary = summary(model, use_cache = FALSE, priors = TRUE)
print(model_summary)

##  Family: hurdle_poisson(log) 
## Formula: children ~ paternalage + birth_cohort + male + maternalage.factor + paternal_loss + maternal_loss + older_siblings + nr.siblings + last_born 
##          hu ~ paternalage + birth_cohort + male + maternalage.factor + paternal_loss + maternal_loss + older_siblings + nr.siblings + last_born
##    Data: model_data (Number of observations: 56663) 
## Samples: 6 chains, each with iter = 800; warmup = 300; thin = 1; 
##          total post-warmup samples = 3000
##     ICs: LOO = Not computed; WAIC = Not computed
##  
## Priors: 
## b ~ normal(0,5)
## b_hu ~ normal(0,5)
## 
## Population-Level Effects: 
##                           Estimate Est.Error l-95% CI u-95% CI Eff.Sample
## Intercept                     1.34      0.06     1.23     1.45        140
## paternalage                   0.01      0.00     0.00     0.01       3000
## birth_cohort1750M1755        -0.09      0.10    -0.29     0.10        445
## birth_cohort1755M1760         0.11      0.08    -0.04     0.25        247
## birth_cohort1760M1765         0.19      0.06     0.06     0.32        187
## birth_cohort1765M1770         0.14      0.06     0.01     0.26        180
## birth_cohort1770M1775         0.12      0.07    -0.01     0.25        188
## birth_cohort1775M1780         0.09      0.06    -0.03     0.21        186
## birth_cohort1780M1785         0.18      0.06     0.06     0.31        168
## birth_cohort1785M1790         0.17      0.06     0.05     0.29        156
## birth_cohort1790M1795         0.08      0.06    -0.03     0.20        154
## birth_cohort1795M1800         0.10      0.06     0.00     0.22        144
## birth_cohort1800M1805         0.05      0.06    -0.06     0.16        139
## birth_cohort1805M1810         0.07      0.06    -0.04     0.18        137
## birth_cohort1810M1815         0.08      0.06    -0.03     0.19        132
## birth_cohort1815M1820         0.15      0.05     0.04     0.26        129
## birth_cohort1820M1825         0.15      0.05     0.05     0.25        131
## birth_cohort1825M1830         0.11      0.05     0.01     0.22        131
## birth_cohort1830M1835         0.13      0.05     0.03     0.24        132
## birth_cohort1835M1840         0.12      0.05     0.02     0.23        133
## birth_cohort1840M1845         0.10      0.05     0.00     0.21        131
## birth_cohort1845M1850         0.11      0.05     0.01     0.22        130
## male1                         0.04      0.01     0.02     0.05       3000
## maternalage.factor1020        0.04      0.03    -0.01     0.09       3000
## maternalage.factor3559        0.06      0.01     0.05     0.08       3000
## paternal_loss01               0.08      0.03     0.03     0.14       3000
## paternal_loss15               0.02      0.02    -0.02     0.05       3000
## paternal_loss510             -0.02      0.02    -0.06     0.01       3000
## paternal_loss1015            -0.03      0.01    -0.06     0.00       3000
## paternal_loss1520            -0.07      0.01    -0.10    -0.05       3000
## paternal_loss2025            -0.03      0.01    -0.06    -0.01       3000
## paternal_loss2530            -0.02      0.01    -0.04     0.00       3000
## paternal_loss3035            -0.02      0.01    -0.04     0.00       3000
## paternal_loss3540             0.02      0.01     0.00     0.04       3000
## paternal_loss4045             0.02      0.01    -0.01     0.04       3000
## maternal_loss01               0.08      0.04    -0.01     0.16       3000
## maternal_loss15               0.03      0.02    -0.02     0.07       3000
## maternal_loss510             -0.02      0.02    -0.05     0.02       3000
## maternal_loss1015            -0.05      0.02    -0.08    -0.02       3000
## maternal_loss1520            -0.04      0.02    -0.07    -0.01       3000
## maternal_loss2025            -0.08      0.01    -0.10    -0.05       3000
## maternal_loss2530            -0.04      0.01    -0.07    -0.02       3000
## maternal_loss3035            -0.03      0.01    -0.05    -0.01       3000
## maternal_loss3540            -0.01      0.01    -0.03     0.01       3000
## maternal_loss4045            -0.02      0.01    -0.04     0.00       3000
## older_siblings1               0.00      0.01    -0.02     0.02       3000
## older_siblings2               0.00      0.01    -0.02     0.02       3000
## older_siblings3              -0.01      0.01    -0.03     0.01       2405
## older_siblings4              -0.04      0.01    -0.06    -0.01       3000
## older_siblings5P             -0.05      0.01    -0.08    -0.02       1975
## nr.siblings                   0.03      0.00     0.03     0.03       2567
## last_born1                   -0.01      0.01    -0.03     0.00       3000
## hu_Intercept                 -0.07      0.14    -0.35     0.22         77
## hu_paternalage               -0.04      0.01    -0.07    -0.01       2125
## hu_birth_cohort1750M1755      0.33      0.24    -0.14     0.82        675
## hu_birth_cohort1755M1760      0.11      0.19    -0.27     0.50        141
## hu_birth_cohort1760M1765      0.02      0.18    -0.33     0.37        132
## hu_birth_cohort1765M1770     -0.25      0.18    -0.59     0.10        124
## hu_birth_cohort1770M1775     -0.04      0.18    -0.38     0.31        131
## hu_birth_cohort1775M1780      0.11      0.16    -0.20     0.44        173
## hu_birth_cohort1780M1785      0.05      0.17    -0.28     0.37        115
## hu_birth_cohort1785M1790      0.16      0.15    -0.14     0.48        109
## hu_birth_cohort1790M1795      0.34      0.14     0.06     0.63         86
## hu_birth_cohort1795M1800      0.15      0.14    -0.12     0.44         83
## hu_birth_cohort1800M1805      0.11      0.14    -0.16     0.39         77
## hu_birth_cohort1805M1810      0.06      0.14    -0.22     0.34         79
## hu_birth_cohort1810M1815      0.15      0.14    -0.12     0.42         80
## hu_birth_cohort1815M1820     -0.02      0.13    -0.29     0.24         73
## hu_birth_cohort1820M1825     -0.12      0.13    -0.39     0.15         75
## hu_birth_cohort1825M1830     -0.12      0.13    -0.39     0.14         75
## hu_birth_cohort1830M1835     -0.10      0.13    -0.37     0.17         73
## hu_birth_cohort1835M1840     -0.11      0.13    -0.37     0.15         75
## hu_birth_cohort1840M1845     -0.09      0.13    -0.36     0.17         72
## hu_birth_cohort1845M1850     -0.07      0.13    -0.34     0.19         74
## hu_male1                      0.04      0.02     0.01     0.08       3000
## hu_maternalage.factor1020     0.04      0.08    -0.11     0.19       3000
## hu_maternalage.factor3559     0.07      0.02     0.02     0.12       2199
## hu_paternal_loss01            0.75      0.08     0.59     0.91       3000
## hu_paternal_loss15            0.59      0.05     0.49     0.69       3000
## hu_paternal_loss510           0.62      0.04     0.54     0.71       3000
## hu_paternal_loss1015          0.48      0.04     0.40     0.56       3000
## hu_paternal_loss1520          0.39      0.04     0.32     0.46       3000
## hu_paternal_loss2025          0.30      0.03     0.24     0.37       3000
## hu_paternal_loss2530          0.22      0.03     0.15     0.28       1644
## hu_paternal_loss3035          0.16      0.03     0.10     0.23       1748
## hu_paternal_loss3540          0.12      0.03     0.05     0.18       1731
## hu_paternal_loss4045          0.04      0.03    -0.02     0.11       3000
## hu_maternal_loss01            1.65      0.11     1.44     1.86       3000
## hu_maternal_loss15            0.88      0.06     0.77     1.00       3000
## hu_maternal_loss510           0.75      0.05     0.66     0.84       3000
## hu_maternal_loss1015          0.71      0.05     0.62     0.80       3000
## hu_maternal_loss1520          0.58      0.04     0.50     0.67       3000
## hu_maternal_loss2025          0.41      0.04     0.33     0.48       3000
## hu_maternal_loss2530          0.29      0.03     0.22     0.35       3000
## hu_maternal_loss3035          0.24      0.03     0.18     0.30       3000
## hu_maternal_loss3540          0.15      0.03     0.09     0.21       3000
## hu_maternal_loss4045          0.09      0.03     0.03     0.15       3000
## hu_older_siblings1           -0.01      0.03    -0.06     0.04       1773
## hu_older_siblings2            0.02      0.03    -0.03     0.08       1596
## hu_older_siblings3            0.04      0.03    -0.02     0.11       1261
## hu_older_siblings4            0.03      0.04    -0.05     0.11       1594
## hu_older_siblings5P          -0.02      0.04    -0.10     0.06       1229
## hu_nr.siblings                0.03      0.00     0.02     0.03       1794
## hu_last_born1                 0.00      0.02    -0.05     0.05       3000
##                           Rhat
## Intercept                 1.03
## paternalage               1.00
## birth_cohort1750M1755     1.01
## birth_cohort1755M1760     1.01
## birth_cohort1760M1765     1.02
## birth_cohort1765M1770     1.02
## birth_cohort1770M1775     1.02
## birth_cohort1775M1780     1.02
## birth_cohort1780M1785     1.02
## birth_cohort1785M1790     1.02
## birth_cohort1790M1795     1.02
## 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
## 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
## hu_Intercept              1.07
## hu_paternalage            1.00
## hu_birth_cohort1750M1755  1.02
## hu_birth_cohort1755M1760  1.03
## hu_birth_cohort1760M1765  1.04
## hu_birth_cohort1765M1770  1.04
## hu_birth_cohort1770M1775  1.04
## hu_birth_cohort1775M1780  1.04
## hu_birth_cohort1780M1785  1.04
## hu_birth_cohort1785M1790  1.04
## hu_birth_cohort1790M1795  1.06
## hu_birth_cohort1795M1800  1.06
## hu_birth_cohort1800M1805  1.06
## hu_birth_cohort1805M1810  1.06
## hu_birth_cohort1810M1815  1.06
## hu_birth_cohort1815M1820  1.07
## hu_birth_cohort1820M1825  1.07
## hu_birth_cohort1825M1830  1.07
## hu_birth_cohort1830M1835  1.07
## hu_birth_cohort1835M1840  1.07
## hu_birth_cohort1840M1845  1.07
## hu_birth_cohort1845M1850  1.07
## hu_male1                  1.00
## hu_maternalage.factor1020 1.00
## hu_maternalage.factor3559 1.00
## hu_paternal_loss01        1.00
## hu_paternal_loss15        1.00
## hu_paternal_loss510       1.00
## hu_paternal_loss1015      1.00
## hu_paternal_loss1520      1.00
## hu_paternal_loss2025      1.00
## hu_paternal_loss2530      1.00
## hu_paternal_loss3035      1.00
## hu_paternal_loss3540      1.00
## hu_paternal_loss4045      1.00
## hu_maternal_loss01        1.00
## hu_maternal_loss15        1.00
## hu_maternal_loss510       1.00
## hu_maternal_loss1015      1.00
## hu_maternal_loss1520      1.00
## hu_maternal_loss2025      1.00
## hu_maternal_loss2530      1.00
## hu_maternal_loss3035      1.00
## hu_maternal_loss3540      1.00
## hu_maternal_loss4045      1.00
## hu_older_siblings1        1.00
## hu_older_siblings2        1.00
## hu_older_siblings3        1.00
## hu_older_siblings4        1.00
## hu_older_siblings5P       1.00
## hu_nr.siblings            1.00
## hu_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	3.827	3.422	4.28
paternalage	1.005	0.9956	1.015
birth_cohort1750M1755	0.9132	0.7479	1.104
birth_cohort1755M1760	1.112	0.9607	1.288
birth_cohort1760M1765	1.212	1.067	1.375
birth_cohort1765M1770	1.147	1.011	1.303
birth_cohort1770M1775	1.126	0.9921	1.284
birth_cohort1775M1780	1.095	0.9671	1.238
birth_cohort1780M1785	1.201	1.061	1.36
birth_cohort1785M1790	1.186	1.056	1.338
birth_cohort1790M1795	1.083	0.9675	1.218
birth_cohort1795M1800	1.108	0.9951	1.241
birth_cohort1800M1805	1.052	0.9426	1.178
birth_cohort1805M1810	1.073	0.9638	1.197
birth_cohort1810M1815	1.086	0.9746	1.212
birth_cohort1815M1820	1.156	1.042	1.29
birth_cohort1820M1825	1.162	1.048	1.289
birth_cohort1825M1830	1.117	1.006	1.243
birth_cohort1830M1835	1.14	1.028	1.272
birth_cohort1835M1840	1.132	1.018	1.263
birth_cohort1840M1845	1.108	0.9965	1.235
birth_cohort1845M1850	1.118	1.01	1.243
male1	1.037	1.025	1.049
maternalage.factor1020	1.036	0.9857	1.09
maternalage.factor3559	1.067	1.05	1.084
paternal_loss01	1.087	1.026	1.151
paternal_loss15	1.018	0.9819	1.056
paternal_loss510	0.9772	0.9452	1.009
paternal_loss1015	0.9714	0.9454	0.9985
paternal_loss1520	0.9283	0.9051	0.9528
paternal_loss2025	0.9686	0.9463	0.9907
paternal_loss2530	0.9803	0.9589	1.003
paternal_loss3035	0.9781	0.9575	0.999
paternal_loss3540	1.02	0.999	1.04
paternal_loss4045	1.017	0.9947	1.04
maternal_loss01	1.079	0.9872	1.177
maternal_loss15	1.027	0.9822	1.074
maternal_loss510	0.9819	0.9477	1.016
maternal_loss1015	0.9515	0.921	0.9827
maternal_loss1520	0.9603	0.9309	0.9909
maternal_loss2025	0.9266	0.9018	0.953
maternal_loss2530	0.958	0.936	0.98
maternal_loss3035	0.971	0.9499	0.9921
maternal_loss3540	0.9894	0.97	1.008
maternal_loss4045	0.9772	0.9585	0.9973
older_siblings1	0.9992	0.982	1.017
older_siblings2	0.9988	0.9801	1.018
older_siblings3	0.9928	0.9709	1.014
older_siblings4	0.9632	0.9378	0.9883
older_siblings5P	0.951	0.9241	0.977
nr.siblings	1.032	1.03	1.035
last_born1	0.9881	0.9721	1.004
hu_Intercept	0.9295	0.7052	1.243
hu_paternalage	0.9624	0.9364	0.9887
hu_birth_cohort1750M1755	1.396	0.8689	2.281
hu_birth_cohort1755M1760	1.117	0.766	1.649
hu_birth_cohort1760M1765	1.017	0.7168	1.443
hu_birth_cohort1765M1770	0.7816	0.5562	1.106
hu_birth_cohort1770M1775	0.9619	0.6842	1.357
hu_birth_cohort1775M1780	1.12	0.8196	1.551
hu_birth_cohort1780M1785	1.051	0.7578	1.452
hu_birth_cohort1785M1790	1.173	0.8691	1.611
hu_birth_cohort1790M1795	1.409	1.06	1.872
hu_birth_cohort1795M1800	1.168	0.8859	1.55
hu_birth_cohort1800M1805	1.116	0.8479	1.481
hu_birth_cohort1805M1810	1.06	0.8016	1.406
hu_birth_cohort1810M1815	1.165	0.8889	1.529
hu_birth_cohort1815M1820	0.9755	0.7478	1.277
hu_birth_cohort1820M1825	0.8851	0.6764	1.159
hu_birth_cohort1825M1830	0.8844	0.676	1.153
hu_birth_cohort1830M1835	0.9035	0.6925	1.18
hu_birth_cohort1835M1840	0.8979	0.6882	1.167
hu_birth_cohort1840M1845	0.9094	0.6986	1.19
hu_birth_cohort1845M1850	0.9293	0.7132	1.208
hu_male1	1.045	1.008	1.083
hu_maternalage.factor1020	1.043	0.8984	1.214
hu_maternalage.factor3559	1.073	1.024	1.124
hu_paternal_loss01	2.123	1.811	2.49
hu_paternal_loss15	1.806	1.636	1.997
hu_paternal_loss510	1.868	1.72	2.026
hu_paternal_loss1015	1.619	1.496	1.751
hu_paternal_loss1520	1.473	1.376	1.58
hu_paternal_loss2025	1.352	1.268	1.442
hu_paternal_loss2530	1.241	1.165	1.326
hu_paternal_loss3035	1.178	1.105	1.253
hu_paternal_loss3540	1.126	1.056	1.195
hu_paternal_loss4045	1.045	0.9753	1.12
hu_maternal_loss01	5.182	4.224	6.427
hu_maternal_loss15	2.423	2.167	2.715
hu_maternal_loss510	2.117	1.931	2.319
hu_maternal_loss1015	2.039	1.864	2.227
hu_maternal_loss1520	1.79	1.651	1.946
hu_maternal_loss2025	1.503	1.393	1.623
hu_maternal_loss2530	1.333	1.248	1.425
hu_maternal_loss3035	1.267	1.194	1.345
hu_maternal_loss3540	1.161	1.097	1.228
hu_maternal_loss4045	1.095	1.031	1.164
hu_older_siblings1	0.9909	0.9407	1.045
hu_older_siblings2	1.025	0.9696	1.087
hu_older_siblings3	1.041	0.9762	1.114
hu_older_siblings4	1.03	0.9534	1.113
hu_older_siblings5P	0.9808	0.9049	1.065
hu_nr.siblings	1.027	1.019	1.036
hu_last_born1	0.9956	0.9484	1.046

Paternal age effect

pander::pander(paternal_age_10y_effect(model))

effect	median_estimate	ci_95	ci_80
estimate father 25y	2.34	[1.98;2.72]	[2.10;2.57]
estimate father 35y	2.39	[2.02;2.78]	[2.16;2.63]
percentage change	2.37	[0.69;4.08]	[1.27;3.49]
OR/IRR	1.01	[1.00;1.01]	[1.00;1.01]
OR hurdle	0.96	[0.94;0.99]	[0.95;0.98]

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)

plot of chunk unnamed-chunk-12

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

plot of chunk unnamed-chunk-13

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

plot of chunk unnamed-chunk-14

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

plot of chunk unnamed-chunk-14

Rhat

Did the 6 chains converge?

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

plot of chunk unnamed-chunk-15

Effective sample size over average sample size

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

plot of chunk unnamed-chunk-16

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

m2: Sibling comparison, no paternal age effect

Here, we compared siblings by including a random intercept for the family, but we modelled no effect for paternal age differences among siblings.

Model summary

Full summary

model_summary = summary(model, use_cache = FALSE, priors = TRUE)
print(model_summary)

##  Family: hurdle_poisson(log) 
## Formula: children ~ birth_cohort + male + maternalage.factor + paternalage.mean + paternal_loss + maternal_loss + older_siblings + nr.siblings + last_born + (1 | idParents) 
##          hu ~ 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: 56663) 
## Samples: 6 chains, each with iter = 1500; warmup = 1000; thin = 1; 
##          total post-warmup samples = 3000
##     ICs: LOO = Not computed; WAIC = Not computed
##  
## Priors: 
## b ~ normal(0,5)
## sd ~ student_t(3, 0, 5)
## b_hu ~ normal(0,5)
## sd_hu ~ student_t(3, 0, 10)
## 
## Group-Level Effects: 
## ~idParents (Number of levels: 14746) 
##                  Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
## sd(Intercept)        0.36      0.01     0.34     0.37        895    1
## sd(hu_Intercept)     0.82      0.02     0.79     0.85       1092    1
## 
## Population-Level Effects: 
##                           Estimate Est.Error l-95% CI u-95% CI Eff.Sample
## Intercept                     1.32      0.08     1.17     1.49        191
## birth_cohort1750M1755        -0.16      0.12    -0.40     0.07        638
## birth_cohort1755M1760         0.09      0.10    -0.10     0.28        281
## birth_cohort1760M1765         0.15      0.09    -0.02     0.32        217
## birth_cohort1765M1770         0.10      0.09    -0.07     0.27        224
## birth_cohort1770M1775         0.06      0.09    -0.12     0.23        212
## birth_cohort1775M1780         0.04      0.09    -0.13     0.21        193
## birth_cohort1780M1785         0.15      0.09    -0.03     0.32        203
## birth_cohort1785M1790         0.13      0.08    -0.04     0.29        196
## birth_cohort1790M1795         0.03      0.08    -0.13     0.19        186
## birth_cohort1795M1800         0.03      0.08    -0.14     0.18        180
## birth_cohort1800M1805        -0.02      0.08    -0.18     0.12        174
## birth_cohort1805M1810        -0.01      0.08    -0.17     0.14        175
## birth_cohort1810M1815         0.02      0.08    -0.14     0.17        176
## birth_cohort1815M1820         0.07      0.08    -0.09     0.22        169
## birth_cohort1820M1825         0.09      0.08    -0.07     0.24        165
## birth_cohort1825M1830         0.04      0.08    -0.11     0.19        169
## birth_cohort1830M1835         0.07      0.08    -0.09     0.22        168
## birth_cohort1835M1840         0.07      0.08    -0.09     0.22        164
## birth_cohort1840M1845         0.04      0.08    -0.11     0.19        164
## birth_cohort1845M1850         0.05      0.08    -0.10     0.20        164
## male1                         0.04      0.01     0.03     0.05       3000
## maternalage.factor1020        0.04      0.03    -0.01     0.10       3000
## maternalage.factor3559        0.05      0.01     0.03     0.07       1547
## paternalage.mean              0.01      0.01    -0.01     0.02       1388
## paternal_loss01               0.05      0.04    -0.02     0.12       1860
## paternal_loss15               0.02      0.02    -0.03     0.06       1155
## paternal_loss510             -0.03      0.02    -0.07     0.01       1207
## paternal_loss1015            -0.03      0.02    -0.06     0.01       1132
## paternal_loss1520            -0.08      0.02    -0.11    -0.04       1205
## paternal_loss2025            -0.03      0.02    -0.07     0.00       1167
## paternal_loss2530            -0.03      0.01    -0.05     0.00       1055
## paternal_loss3035            -0.02      0.01    -0.05     0.01       1118
## paternal_loss3540             0.02      0.01    -0.01     0.04        777
## paternal_loss4045             0.03      0.01     0.00     0.05       1556
## maternal_loss01               0.06      0.05    -0.05     0.16       1844
## maternal_loss15               0.00      0.03    -0.06     0.05       1656
## maternal_loss510             -0.03      0.02    -0.07     0.02       1604
## maternal_loss1015            -0.04      0.02    -0.09     0.00       1694
## maternal_loss1520            -0.04      0.02    -0.08     0.00       1822
## maternal_loss2025            -0.07      0.02    -0.11    -0.04       1440
## maternal_loss2530            -0.03      0.02    -0.06     0.00       1448
## maternal_loss3035            -0.02      0.01    -0.05     0.00       1200
## maternal_loss3540             0.00      0.01    -0.03     0.02       1504
## maternal_loss4045            -0.02      0.01    -0.05     0.00       1809
## older_siblings1               0.01      0.01    -0.01     0.02       1466
## older_siblings2               0.01      0.01    -0.02     0.03       1291
## older_siblings3               0.00      0.01    -0.02     0.03       1065
## older_siblings4              -0.03      0.02    -0.06     0.00       1210
## older_siblings5P             -0.03      0.02    -0.07     0.00        940
## nr.siblings                   0.03      0.00     0.02     0.03        775
## last_born1                   -0.02      0.01    -0.04     0.00       3000
## hu_Intercept                 -0.06      0.21    -0.48     0.36        170
## hu_birth_cohort1750M1755      0.36      0.30    -0.22     0.94        347
## hu_birth_cohort1755M1760      0.08      0.27    -0.44     0.59        254
## hu_birth_cohort1760M1765     -0.05      0.24    -0.51     0.41        210
## hu_birth_cohort1765M1770     -0.34      0.24    -0.83     0.12        199
## hu_birth_cohort1770M1775     -0.11      0.24    -0.58     0.38        198
## hu_birth_cohort1775M1780      0.05      0.23    -0.39     0.52        190
## hu_birth_cohort1780M1785     -0.01      0.23    -0.48     0.45        197
## hu_birth_cohort1785M1790      0.13      0.22    -0.30     0.55        183
## hu_birth_cohort1790M1795      0.31      0.21    -0.11     0.72        173
## hu_birth_cohort1795M1800      0.11      0.21    -0.29     0.51        165
## hu_birth_cohort1800M1805      0.03      0.21    -0.38     0.43        162
## hu_birth_cohort1805M1810     -0.01      0.21    -0.42     0.39        158
## hu_birth_cohort1810M1815      0.06      0.21    -0.35     0.45        162
## hu_birth_cohort1815M1820     -0.10      0.20    -0.51     0.29        153
## hu_birth_cohort1820M1825     -0.21      0.20    -0.61     0.17        151
## hu_birth_cohort1825M1830     -0.21      0.20    -0.60     0.19        150
## hu_birth_cohort1830M1835     -0.19      0.20    -0.58     0.20        153
## hu_birth_cohort1835M1840     -0.21      0.20    -0.61     0.18        155
## hu_birth_cohort1840M1845     -0.20      0.20    -0.60     0.20        153
## hu_birth_cohort1845M1850     -0.16      0.20    -0.55     0.24        151
## hu_male1                      0.04      0.02     0.01     0.08       3000
## hu_maternalage.factor1020     0.05      0.09    -0.12     0.23       3000
## hu_maternalage.factor3559     0.08      0.03     0.03     0.13       3000
## hu_paternalage.mean          -0.03      0.02    -0.07     0.01       2052
## hu_paternal_loss01            0.87      0.09     0.69     1.05       3000
## hu_paternal_loss15            0.68      0.06     0.55     0.80       3000
## hu_paternal_loss510           0.69      0.05     0.59     0.80       1207
## hu_paternal_loss1015          0.53      0.05     0.44     0.62       1189
## hu_paternal_loss1520          0.43      0.04     0.34     0.51       1223
## hu_paternal_loss2025          0.32      0.04     0.24     0.40       1164
## hu_paternal_loss2530          0.23      0.04     0.15     0.31       1196
## hu_paternal_loss3035          0.16      0.04     0.08     0.24       1270
## hu_paternal_loss3540          0.12      0.04     0.05     0.19       1627
## hu_paternal_loss4045          0.04      0.04    -0.04     0.12       1949
## hu_maternal_loss01            1.86      0.12     1.62     2.09       3000
## hu_maternal_loss15            1.03      0.07     0.90     1.16       3000
## hu_maternal_loss510           0.87      0.06     0.76     0.98       3000
## hu_maternal_loss1015          0.81      0.05     0.70     0.91       3000
## hu_maternal_loss1520          0.67      0.05     0.58     0.77       2112
## hu_maternal_loss2025          0.46      0.04     0.37     0.54       3000
## hu_maternal_loss2530          0.31      0.04     0.24     0.39       3000
## hu_maternal_loss3035          0.24      0.04     0.17     0.32       2011
## hu_maternal_loss3540          0.14      0.03     0.08     0.21       3000
## hu_maternal_loss4045          0.08      0.03     0.01     0.14       3000
## hu_older_siblings1           -0.03      0.03    -0.09     0.02       3000
## hu_older_siblings2           -0.01      0.03    -0.08     0.05       1774
## hu_older_siblings3           -0.02      0.04    -0.09     0.05       1643
## hu_older_siblings4           -0.04      0.04    -0.13     0.04       1505
## hu_older_siblings5P          -0.09      0.05    -0.18     0.00       1182
## hu_nr.siblings                0.04      0.01     0.03     0.05       1349
## hu_last_born1                 0.01      0.03    -0.04     0.06       3000
##                           Rhat
## Intercept                 1.03
## birth_cohort1750M1755     1.01
## birth_cohort1755M1760     1.01
## birth_cohort1760M1765     1.02
## birth_cohort1765M1770     1.02
## birth_cohort1770M1775     1.02
## birth_cohort1775M1780     1.03
## 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.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.00
## paternal_loss01           1.00
## paternal_loss15           1.00
## paternal_loss510          1.00
## paternal_loss1015         1.01
## paternal_loss1520         1.00
## paternal_loss2025         1.01
## paternal_loss2530         1.01
## paternal_loss3035         1.01
## 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.01
## last_born1                1.00
## hu_Intercept              1.02
## hu_birth_cohort1750M1755  1.01
## hu_birth_cohort1755M1760  1.02
## hu_birth_cohort1760M1765  1.02
## hu_birth_cohort1765M1770  1.02
## hu_birth_cohort1770M1775  1.02
## hu_birth_cohort1775M1780  1.02
## hu_birth_cohort1780M1785  1.02
## hu_birth_cohort1785M1790  1.02
## hu_birth_cohort1790M1795  1.02
## hu_birth_cohort1795M1800  1.02
## hu_birth_cohort1800M1805  1.03
## hu_birth_cohort1805M1810  1.03
## hu_birth_cohort1810M1815  1.03
## hu_birth_cohort1815M1820  1.03
## hu_birth_cohort1820M1825  1.03
## hu_birth_cohort1825M1830  1.03
## hu_birth_cohort1830M1835  1.03
## hu_birth_cohort1835M1840  1.03
## hu_birth_cohort1840M1845  1.03
## hu_birth_cohort1845M1850  1.03
## hu_male1                  1.00
## hu_maternalage.factor1020 1.00
## hu_maternalage.factor3559 1.00
## hu_paternalage.mean       1.01
## hu_paternal_loss01        1.00
## hu_paternal_loss15        1.00
## hu_paternal_loss510       1.00
## hu_paternal_loss1015      1.01
## hu_paternal_loss1520      1.01
## hu_paternal_loss2025      1.01
## hu_paternal_loss2530      1.01
## hu_paternal_loss3035      1.01
## hu_paternal_loss3540      1.00
## hu_paternal_loss4045      1.00
## hu_maternal_loss01        1.00
## hu_maternal_loss15        1.00
## hu_maternal_loss510       1.00
## hu_maternal_loss1015      1.00
## hu_maternal_loss1520      1.00
## hu_maternal_loss2025      1.00
## hu_maternal_loss2530      1.00
## hu_maternal_loss3035      1.00
## hu_maternal_loss3540      1.00
## hu_maternal_loss4045      1.00
## hu_older_siblings1        1.00
## hu_older_siblings2        1.00
## hu_older_siblings3        1.00
## hu_older_siblings4        1.00
## hu_older_siblings5P       1.00
## hu_nr.siblings            1.00
## hu_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	3.757	3.212	4.44
birth_cohort1750M1755	0.8533	0.6727	1.072
birth_cohort1755M1760	1.099	0.9066	1.321
birth_cohort1760M1765	1.165	0.9775	1.377
birth_cohort1765M1770	1.107	0.9317	1.31
birth_cohort1770M1775	1.061	0.8894	1.256
birth_cohort1775M1780	1.045	0.876	1.234
birth_cohort1780M1785	1.164	0.9709	1.374
birth_cohort1785M1790	1.134	0.9616	1.331
birth_cohort1790M1795	1.029	0.8751	1.209
birth_cohort1795M1800	1.026	0.8705	1.192
birth_cohort1800M1805	0.9775	0.8325	1.133
birth_cohort1805M1810	0.9924	0.8444	1.152
birth_cohort1810M1815	1.017	0.867	1.181
birth_cohort1815M1820	1.071	0.915	1.242
birth_cohort1820M1825	1.091	0.9348	1.265
birth_cohort1825M1830	1.046	0.8955	1.21
birth_cohort1830M1835	1.071	0.9172	1.242
birth_cohort1835M1840	1.069	0.9144	1.242
birth_cohort1840M1845	1.045	0.8954	1.214
birth_cohort1845M1850	1.054	0.9009	1.224
male1	1.041	1.027	1.055
maternalage.factor1020	1.044	0.9859	1.109
maternalage.factor3559	1.053	1.032	1.073
paternalage.mean	1.009	0.9942	1.025
paternal_loss01	1.049	0.9766	1.128
paternal_loss15	1.015	0.9677	1.066
paternal_loss510	0.971	0.9307	1.011
paternal_loss1015	0.9725	0.9387	1.009
paternal_loss1520	0.9229	0.8934	0.9564
paternal_loss2025	0.9669	0.9369	0.9981
paternal_loss2530	0.9744	0.9472	1.002
paternal_loss3035	0.9805	0.9545	1.009
paternal_loss3540	1.019	0.9922	1.046
paternal_loss4045	1.027	1.002	1.054
maternal_loss01	1.063	0.9554	1.175
maternal_loss15	0.9953	0.9396	1.053
maternal_loss510	0.9741	0.9314	1.021
maternal_loss1015	0.96	0.9179	1.001
maternal_loss1520	0.9609	0.923	1.001
maternal_loss2025	0.9289	0.8981	0.9603
maternal_loss2530	0.9697	0.9418	0.9987
maternal_loss3035	0.9766	0.9506	1.004
maternal_loss3540	0.9994	0.9748	1.023
maternal_loss4045	0.9789	0.9559	1.004
older_siblings1	1.005	0.9862	1.025
older_siblings2	1.006	0.985	1.029
older_siblings3	1.004	0.9783	1.031
older_siblings4	0.9715	0.9423	1.001
older_siblings5P	0.967	0.9366	0.9988
nr.siblings	1.028	1.024	1.032
last_born1	0.979	0.961	0.9976
hu_Intercept	0.9399	0.6214	1.437
hu_birth_cohort1750M1755	1.433	0.8021	2.572
hu_birth_cohort1755M1760	1.081	0.6429	1.797
hu_birth_cohort1760M1765	0.95	0.6029	1.507
hu_birth_cohort1765M1770	0.7125	0.4376	1.126
hu_birth_cohort1770M1775	0.9001	0.5608	1.458
hu_birth_cohort1775M1780	1.055	0.6746	1.676
hu_birth_cohort1780M1785	0.9892	0.6186	1.575
hu_birth_cohort1785M1790	1.14	0.7437	1.729
hu_birth_cohort1790M1795	1.36	0.894	2.056
hu_birth_cohort1795M1800	1.12	0.7449	1.667
hu_birth_cohort1800M1805	1.034	0.6839	1.536
hu_birth_cohort1805M1810	0.9887	0.6584	1.471
hu_birth_cohort1810M1815	1.065	0.7065	1.574
hu_birth_cohort1815M1820	0.9022	0.6017	1.342
hu_birth_cohort1820M1825	0.812	0.5457	1.188
hu_birth_cohort1825M1830	0.8129	0.5502	1.208
hu_birth_cohort1830M1835	0.828	0.5582	1.226
hu_birth_cohort1835M1840	0.8115	0.5447	1.199
hu_birth_cohort1840M1845	0.8221	0.55	1.219
hu_birth_cohort1845M1850	0.8525	0.5742	1.271
hu_male1	1.046	1.007	1.086
hu_maternalage.factor1020	1.055	0.8888	1.256
hu_maternalage.factor3559	1.085	1.029	1.142
hu_paternalage.mean	0.9714	0.9357	1.009
hu_paternal_loss01	2.393	1.986	2.87
hu_paternal_loss15	1.968	1.736	2.22
hu_paternal_loss510	2.003	1.803	2.22
hu_paternal_loss1015	1.693	1.546	1.859
hu_paternal_loss1520	1.533	1.402	1.672
hu_paternal_loss2025	1.384	1.271	1.499
hu_paternal_loss2530	1.256	1.158	1.364
hu_paternal_loss3035	1.174	1.086	1.267
hu_paternal_loss3540	1.128	1.047	1.209
hu_paternal_loss4045	1.04	0.9622	1.122
hu_maternal_loss01	6.408	5.059	8.108
hu_maternal_loss15	2.797	2.462	3.189
hu_maternal_loss510	2.378	2.137	2.661
hu_maternal_loss1015	2.239	2.019	2.492
hu_maternal_loss1520	1.96	1.778	2.153
hu_maternal_loss2025	1.581	1.443	1.717
hu_maternal_loss2530	1.369	1.266	1.482
hu_maternal_loss3035	1.277	1.191	1.371
hu_maternal_loss3540	1.155	1.081	1.239
hu_maternal_loss4045	1.079	1.008	1.151
hu_older_siblings1	0.9692	0.9174	1.023
hu_older_siblings2	0.9857	0.9242	1.053
hu_older_siblings3	0.9792	0.9131	1.055
hu_older_siblings4	0.9562	0.8802	1.046
hu_older_siblings5P	0.9121	0.8341	0.9974
hu_nr.siblings	1.044	1.033	1.056
hu_last_born1	1.007	0.9568	1.062

Paternal age effect

pander::pander(paternal_age_10y_effect(model))

Marginal effect plots

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

plot of chunk unnamed-chunk-23

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

plot of chunk unnamed-chunk-24

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

plot of chunk unnamed-chunk-25

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

plot of chunk unnamed-chunk-25

Rhat

Did the 6 chains converge?

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

plot of chunk unnamed-chunk-26

Effective sample size over average sample size

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

## Error in factor(x, labels = names(sort(x))): invalid 'labels'; length 29599 should be 1 or 29591

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

m3: Sibling comparison, linear paternal age effect

Here, we compared siblings by including a random intercept for the family, and we modelled a linear effect for paternal age differences among siblings.

Model summary

Full summary

model_summary = summary(model, use_cache = FALSE, priors = TRUE)
print(model_summary)

##  Family: hurdle_poisson(log) 
## Formula: children ~ paternalage + birth_cohort + male + maternalage.factor + paternalage.mean + paternal_loss + maternal_loss + older_siblings + nr.siblings + last_born + (1 | idParents) 
##          hu ~ 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: 56663) 
## Samples: 6 chains, each with iter = 800; warmup = 300; thin = 1; 
##          total post-warmup samples = 3000
##     ICs: LOO = Not computed; WAIC = Not computed
##  
## Priors: 
## b ~ normal(0,5)
## sd ~ student_t(3, 0, 5)
## b_hu ~ normal(0,5)
## sd_hu ~ student_t(3, 0, 10)
## 
## Group-Level Effects: 
## ~idParents (Number of levels: 14746) 
##                  Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
## sd(Intercept)        0.36      0.01     0.34     0.37       1049    1
## sd(hu_Intercept)     0.82      0.02     0.78     0.85       1292    1
## 
## Population-Level Effects: 
##                           Estimate Est.Error l-95% CI u-95% CI Eff.Sample
## Intercept                     1.32      0.08     1.16     1.48        247
## paternalage                  -0.05      0.02    -0.09    -0.02        896
## birth_cohort1750M1755        -0.17      0.12    -0.40     0.06        607
## birth_cohort1755M1760         0.09      0.10    -0.10     0.27        380
## birth_cohort1760M1765         0.15      0.09    -0.02     0.31        285
## birth_cohort1765M1770         0.10      0.09    -0.07     0.27        284
## birth_cohort1770M1775         0.06      0.09    -0.12     0.23        279
## birth_cohort1775M1780         0.04      0.09    -0.13     0.22        263
## birth_cohort1780M1785         0.15      0.09    -0.01     0.32        287
## birth_cohort1785M1790         0.12      0.08    -0.04     0.29        271
## birth_cohort1790M1795         0.03      0.08    -0.13     0.18        247
## birth_cohort1795M1800         0.02      0.08    -0.13     0.18        234
## birth_cohort1800M1805        -0.03      0.08    -0.18     0.13        230
## birth_cohort1805M1810        -0.01      0.08    -0.16     0.14        228
## birth_cohort1810M1815         0.01      0.08    -0.14     0.17        223
## birth_cohort1815M1820         0.07      0.08    -0.08     0.22        218
## birth_cohort1820M1825         0.08      0.08    -0.07     0.24        221
## birth_cohort1825M1830         0.04      0.08    -0.11     0.19        223
## birth_cohort1830M1835         0.07      0.08    -0.08     0.22        217
## birth_cohort1835M1840         0.07      0.08    -0.08     0.22        219
## birth_cohort1840M1845         0.04      0.08    -0.10     0.19        219
## birth_cohort1845M1850         0.05      0.08    -0.10     0.21        219
## male1                         0.04      0.01     0.03     0.05       3000
## maternalage.factor1020        0.04      0.03    -0.02     0.10       3000
## maternalage.factor3559        0.06      0.01     0.04     0.09       3000
## paternalage.mean              0.06      0.02     0.02     0.10        957
## paternal_loss01               0.06      0.04    -0.01     0.13       2274
## paternal_loss15               0.03      0.02    -0.02     0.07       1751
## paternal_loss510             -0.02      0.02    -0.06     0.02       1638
## paternal_loss1015            -0.02      0.02    -0.06     0.02       1590
## paternal_loss1520            -0.07      0.02    -0.11    -0.04       1533
## paternal_loss2025            -0.03      0.02    -0.06     0.00       1561
## paternal_loss2530            -0.02      0.01    -0.05     0.01       1443
## paternal_loss3035            -0.02      0.01    -0.04     0.01       1691
## paternal_loss3540             0.02      0.01     0.00     0.05       1855
## paternal_loss4045             0.03      0.01     0.00     0.05       3000
## maternal_loss01               0.07      0.05    -0.04     0.17       3000
## maternal_loss15               0.00      0.03    -0.05     0.06       3000
## maternal_loss510             -0.02      0.02    -0.07     0.03       1929
## maternal_loss1015            -0.04      0.02    -0.08     0.01       2085
## maternal_loss1520            -0.04      0.02    -0.08     0.00       2052
## maternal_loss2025            -0.07      0.02    -0.11    -0.03       1705
## maternal_loss2530            -0.03      0.02    -0.06     0.00       1788
## maternal_loss3035            -0.02      0.01    -0.05     0.01       1966
## maternal_loss3540             0.00      0.01    -0.02     0.03       1975
## maternal_loss4045            -0.02      0.01    -0.04     0.00       3000
## older_siblings1               0.02      0.01    -0.01     0.04       1769
## older_siblings2               0.03      0.01     0.00     0.06       1239
## older_siblings3               0.04      0.02     0.00     0.07       1042
## older_siblings4               0.01      0.02    -0.03     0.06       1025
## older_siblings5P              0.03      0.03    -0.03     0.08        996
## nr.siblings                   0.02      0.00     0.02     0.03       1097
## last_born1                   -0.02      0.01    -0.04     0.00       3000
## hu_Intercept                 -0.07      0.20    -0.44     0.32        184
## hu_paternalage                0.05      0.06    -0.06     0.16        948
## hu_birth_cohort1750M1755      0.37      0.28    -0.18     0.93        658
## hu_birth_cohort1755M1760      0.09      0.24    -0.38     0.59        358
## hu_birth_cohort1760M1765     -0.04      0.22    -0.46     0.40        298
## hu_birth_cohort1765M1770     -0.33      0.22    -0.76     0.09        279
## hu_birth_cohort1770M1775     -0.09      0.22    -0.55     0.34        286
## hu_birth_cohort1775M1780      0.06      0.22    -0.37     0.48        262
## hu_birth_cohort1780M1785      0.00      0.21    -0.42     0.43        271
## hu_birth_cohort1785M1790      0.14      0.20    -0.26     0.52        246
## hu_birth_cohort1790M1795      0.33      0.19    -0.06     0.68        188
## hu_birth_cohort1795M1800      0.13      0.19    -0.25     0.50        181
## hu_birth_cohort1800M1805      0.05      0.19    -0.34     0.41        210
## hu_birth_cohort1805M1810      0.00      0.19    -0.37     0.36        169
## hu_birth_cohort1810M1815      0.08      0.19    -0.30     0.43        181
## hu_birth_cohort1815M1820     -0.09      0.19    -0.46     0.27        176
## hu_birth_cohort1820M1825     -0.19      0.19    -0.57     0.16        175
## hu_birth_cohort1825M1830     -0.19      0.18    -0.55     0.16        201
## hu_birth_cohort1830M1835     -0.18      0.19    -0.55     0.17        171
## hu_birth_cohort1835M1840     -0.20      0.18    -0.57     0.15        171
## hu_birth_cohort1840M1845     -0.19      0.18    -0.55     0.16        175
## hu_birth_cohort1845M1850     -0.15      0.18    -0.51     0.20        173
## hu_male1                      0.04      0.02     0.01     0.08       3000
## hu_maternalage.factor1020     0.06      0.09    -0.12     0.24       3000
## hu_maternalage.factor3559     0.07      0.03     0.01     0.13       3000
## hu_paternalage.mean          -0.08      0.06    -0.19     0.03       1009
## hu_paternal_loss01            0.87      0.09     0.68     1.05       3000
## hu_paternal_loss15            0.67      0.06     0.55     0.79       1867
## hu_paternal_loss510           0.69      0.05     0.58     0.79       3000
## hu_paternal_loss1015          0.52      0.05     0.43     0.62       1712
## hu_paternal_loss1520          0.42      0.05     0.33     0.51       1607
## hu_paternal_loss2025          0.32      0.04     0.24     0.40       1494
## hu_paternal_loss2530          0.22      0.04     0.14     0.30       1463
## hu_paternal_loss3035          0.16      0.04     0.08     0.23       1857
## hu_paternal_loss3540          0.12      0.04     0.05     0.19       2071
## hu_paternal_loss4045          0.04      0.04    -0.04     0.12       3000
## hu_maternal_loss01            1.86      0.12     1.63     2.10       3000
## hu_maternal_loss15            1.02      0.07     0.88     1.16       3000
## hu_maternal_loss510           0.86      0.06     0.75     0.98       3000
## hu_maternal_loss1015          0.80      0.06     0.69     0.92       3000
## hu_maternal_loss1520          0.67      0.05     0.57     0.77       3000
## hu_maternal_loss2025          0.46      0.05     0.37     0.55       3000
## hu_maternal_loss2530          0.31      0.04     0.23     0.40       3000
## hu_maternal_loss3035          0.24      0.04     0.17     0.32       1995
## hu_maternal_loss3540          0.14      0.03     0.08     0.21       3000
## hu_maternal_loss4045          0.08      0.04     0.01     0.15       3000
## hu_older_siblings1           -0.04      0.03    -0.10     0.02       3000
## hu_older_siblings2           -0.04      0.04    -0.12     0.05       1069
## hu_older_siblings3           -0.05      0.05    -0.15     0.05       1098
## hu_older_siblings4           -0.09      0.06    -0.21     0.04       1014
## hu_older_siblings5P          -0.15      0.08    -0.31     0.01        918
## hu_nr.siblings                0.05      0.01     0.03     0.06       1241
## hu_last_born1                 0.01      0.03    -0.04     0.06       3000
##                           Rhat
## Intercept                 1.04
## paternalage               1.00
## birth_cohort1750M1755     1.01
## birth_cohort1755M1760     1.02
## birth_cohort1760M1765     1.03
## birth_cohort1765M1770     1.03
## birth_cohort1770M1775     1.03
## birth_cohort1775M1780     1.04
## birth_cohort1780M1785     1.03
## birth_cohort1785M1790     1.04
## birth_cohort1790M1795     1.04
## birth_cohort1795M1800     1.04
## birth_cohort1800M1805     1.04
## birth_cohort1805M1810     1.04
## birth_cohort1810M1815     1.05
## birth_cohort1815M1820     1.05
## birth_cohort1820M1825     1.05
## birth_cohort1825M1830     1.04
## birth_cohort1830M1835     1.05
## birth_cohort1835M1840     1.05
## birth_cohort1840M1845     1.04
## birth_cohort1845M1850     1.04
## 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.00
## nr.siblings               1.00
## last_born1                1.00
## hu_Intercept              1.02
## hu_paternalage            1.01
## hu_birth_cohort1750M1755  1.00
## hu_birth_cohort1755M1760  1.01
## hu_birth_cohort1760M1765  1.01
## hu_birth_cohort1765M1770  1.01
## hu_birth_cohort1770M1775  1.01
## hu_birth_cohort1775M1780  1.01
## hu_birth_cohort1780M1785  1.01
## hu_birth_cohort1785M1790  1.01
## hu_birth_cohort1790M1795  1.01
## hu_birth_cohort1795M1800  1.02
## hu_birth_cohort1800M1805  1.02
## hu_birth_cohort1805M1810  1.02
## hu_birth_cohort1810M1815  1.02
## hu_birth_cohort1815M1820  1.02
## hu_birth_cohort1820M1825  1.02
## hu_birth_cohort1825M1830  1.02
## hu_birth_cohort1830M1835  1.02
## hu_birth_cohort1835M1840  1.02
## hu_birth_cohort1840M1845  1.02
## hu_birth_cohort1845M1850  1.02
## hu_male1                  1.00
## hu_maternalage.factor1020 1.00
## hu_maternalage.factor3559 1.00
## hu_paternalage.mean       1.01
## hu_paternal_loss01        1.00
## hu_paternal_loss15        1.00
## hu_paternal_loss510       1.00
## hu_paternal_loss1015      1.00
## hu_paternal_loss1520      1.00
## hu_paternal_loss2025      1.00
## hu_paternal_loss2530      1.00
## hu_paternal_loss3035      1.00
## hu_paternal_loss3540      1.00
## hu_paternal_loss4045      1.00
## hu_maternal_loss01        1.00
## hu_maternal_loss15        1.00
## hu_maternal_loss510       1.00
## hu_maternal_loss1015      1.00
## hu_maternal_loss1520      1.00
## hu_maternal_loss2025      1.00
## hu_maternal_loss2530      1.00
## hu_maternal_loss3035      1.00
## hu_maternal_loss3540      1.00
## hu_maternal_loss4045      1.00
## hu_older_siblings1        1.00
## hu_older_siblings2        1.01
## hu_older_siblings3        1.01
## hu_older_siblings4        1.01
## hu_older_siblings5P       1.01
## hu_nr.siblings            1.01
## hu_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	3.74	3.194	4.392
paternalage	0.9483	0.9128	0.9849
birth_cohort1750M1755	0.843	0.6693	1.057
birth_cohort1755M1760	1.097	0.9086	1.316
birth_cohort1760M1765	1.16	0.9757	1.37
birth_cohort1765M1770	1.103	0.9314	1.306
birth_cohort1770M1775	1.061	0.8882	1.26
birth_cohort1775M1780	1.046	0.8768	1.243
birth_cohort1780M1785	1.163	0.9866	1.373
birth_cohort1785M1790	1.131	0.9612	1.336
birth_cohort1790M1795	1.025	0.8771	1.2
birth_cohort1795M1800	1.021	0.8766	1.197
birth_cohort1800M1805	0.9744	0.8359	1.137
birth_cohort1805M1810	0.9878	0.8502	1.151
birth_cohort1810M1815	1.013	0.8731	1.182
birth_cohort1815M1820	1.068	0.9225	1.24
birth_cohort1820M1825	1.087	0.9363	1.266
birth_cohort1825M1830	1.042	0.8992	1.212
birth_cohort1830M1835	1.069	0.9226	1.245
birth_cohort1835M1840	1.069	0.9231	1.243
birth_cohort1840M1845	1.045	0.9045	1.213
birth_cohort1845M1850	1.054	0.9078	1.228
male1	1.041	1.027	1.054
maternalage.factor1020	1.037	0.9778	1.103
maternalage.factor3559	1.066	1.044	1.089
paternalage.mean	1.061	1.02	1.105
paternal_loss01	1.063	0.9862	1.144
paternal_loss15	1.027	0.9797	1.077
paternal_loss510	0.9817	0.9414	1.025
paternal_loss1015	0.9815	0.945	1.018
paternal_loss1520	0.9303	0.8992	0.9623
paternal_loss2025	0.973	0.943	1.004
paternal_loss2530	0.9801	0.9521	1.008
paternal_loss3035	0.9851	0.959	1.013
paternal_loss3540	1.022	0.996	1.048
paternal_loss4045	1.028	1.004	1.054
maternal_loss01	1.068	0.9613	1.184
maternal_loss15	1.003	0.9486	1.059
maternal_loss510	0.9801	0.9353	1.027
maternal_loss1015	0.9645	0.9252	1.008
maternal_loss1520	0.9648	0.9273	1.003
maternal_loss2025	0.9327	0.9001	0.9663
maternal_loss2530	0.9732	0.9433	1.003
maternal_loss3035	0.9793	0.953	1.006
maternal_loss3540	1.002	0.9778	1.027
maternal_loss4045	0.9805	0.9568	1.004
older_siblings1	1.016	0.9949	1.037
older_siblings2	1.029	0.9999	1.059
older_siblings3	1.037	1	1.073
older_siblings4	1.012	0.9684	1.057
older_siblings5P	1.028	0.9711	1.086
nr.siblings	1.024	1.019	1.029
last_born1	0.9788	0.9602	0.9974
hu_Intercept	0.9314	0.6461	1.376
hu_paternalage	1.051	0.9455	1.171
hu_birth_cohort1750M1755	1.446	0.8342	2.536
hu_birth_cohort1755M1760	1.09	0.6843	1.803
hu_birth_cohort1760M1765	0.9626	0.6318	1.494
hu_birth_cohort1765M1770	0.7218	0.467	1.1
hu_birth_cohort1770M1775	0.9125	0.5758	1.399
hu_birth_cohort1775M1780	1.067	0.6925	1.614
hu_birth_cohort1780M1785	1.004	0.6566	1.531
hu_birth_cohort1785M1790	1.151	0.7711	1.684
hu_birth_cohort1790M1795	1.385	0.9371	1.964
hu_birth_cohort1795M1800	1.136	0.7771	1.648
hu_birth_cohort1800M1805	1.047	0.7134	1.511
hu_birth_cohort1805M1810	1.001	0.6917	1.438
hu_birth_cohort1810M1815	1.081	0.7441	1.541
hu_birth_cohort1815M1820	0.9139	0.632	1.304
hu_birth_cohort1820M1825	0.824	0.5646	1.169
hu_birth_cohort1825M1830	0.8243	0.5752	1.177
hu_birth_cohort1830M1835	0.8385	0.5775	1.185
hu_birth_cohort1835M1840	0.8204	0.5667	1.162
hu_birth_cohort1840M1845	0.8305	0.5771	1.176
hu_birth_cohort1845M1850	0.8625	0.5999	1.22
hu_male1	1.045	1.006	1.088
hu_maternalage.factor1020	1.057	0.8897	1.267
hu_maternalage.factor3559	1.073	1.011	1.138
hu_paternalage.mean	0.9262	0.829	1.03
hu_paternal_loss01	2.376	1.98	2.855
hu_paternal_loss15	1.955	1.732	2.207
hu_paternal_loss510	1.991	1.791	2.213
hu_paternal_loss1015	1.685	1.533	1.857
hu_paternal_loss1520	1.526	1.393	1.667
hu_paternal_loss2025	1.378	1.269	1.498
hu_paternal_loss2530	1.252	1.155	1.353
hu_paternal_loss3035	1.171	1.087	1.264
hu_paternal_loss3540	1.126	1.047	1.209
hu_paternal_loss4045	1.04	0.9654	1.124
hu_maternal_loss01	6.399	5.082	8.148
hu_maternal_loss15	2.779	2.415	3.204
hu_maternal_loss510	2.37	2.112	2.665
hu_maternal_loss1015	2.232	1.991	2.499
hu_maternal_loss1520	1.958	1.776	2.162
hu_maternal_loss2025	1.578	1.447	1.725
hu_maternal_loss2530	1.368	1.26	1.494
hu_maternal_loss3035	1.277	1.187	1.38
hu_maternal_loss3540	1.154	1.079	1.233
hu_maternal_loss4045	1.079	1.007	1.157
hu_older_siblings1	0.9591	0.9037	1.018
hu_older_siblings2	0.9648	0.888	1.048
hu_older_siblings3	0.9487	0.8608	1.05
hu_older_siblings4	0.9175	0.8107	1.04
hu_older_siblings5P	0.859	0.7305	1.006
hu_nr.siblings	1.049	1.035	1.064
hu_last_born1	1.008	0.9586	1.061

Paternal age effect

pander::pander(paternal_age_10y_effect(model))

effect	median_estimate	ci_95	ci_80
estimate father 25y	2.26	[1.78;2.79]	[1.95;2.61]
estimate father 35y	2.10	[1.64;2.62]	[1.79;2.43]
percentage change	-7.29	[-13.40; -1.07]	[-11.15; -3.33]
OR/IRR	0.95	[0.91;0.98]	[0.92;0.97]
OR hurdle	1.05	[0.95;1.17]	[0.98;1.13]

Marginal effect plots

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

plot of chunk unnamed-chunk-34

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

plot of chunk unnamed-chunk-35

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

plot of chunk unnamed-chunk-36

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

plot of chunk unnamed-chunk-36

Rhat

Did the 6 chains converge?

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

plot of chunk unnamed-chunk-37

Effective sample size over average sample size

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

## Error in factor(x, labels = names(sort(x))): invalid 'labels'; length 29601 should be 1 or 29598

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

m4: Sibling comparison, nonlinear paternal age effect

Here, we compared siblings by including a random intercept for the family, and we modelled a possibly nonlinear effect for paternal age differences among siblings.

Model summary

Full summary

model_summary = summary(model, use_cache = FALSE, priors = TRUE)

## Warning: There were 22 divergent transitions after warmup. Increasing adapt_delta above 0.8 may help.
## See http://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup

print(model_summary)

##  Family: hurdle_poisson(log) 
## Formula: children ~ s(paternalage) + birth_cohort + male + maternalage.factor + paternalage.mean + paternal_loss + maternal_loss + older_siblings + nr.siblings + last_born + (1 | idParents) 
##          hu ~ s(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: 56663) 
## Samples: 6 chains, each with iter = 800; warmup = 300; thin = 1; 
##          total post-warmup samples = 3000
##     ICs: LOO = Not computed; WAIC = Not computed
##  
## Priors: 
## b ~ normal(0,5)
## sd ~ student_t(3, 0, 5)
## sds ~ student_t(3, 0, 10)
## b_hu ~ normal(0,5)
## sd_hu ~ student_t(3, 0, 10)
## sds_hu ~ student_t(3, 0, 10)
## 
## Smooth Terms: 
##                        Estimate Est.Error l-95% CI u-95% CI Eff.Sample
## sds(spaternalage_1)        0.28      0.19     0.03     0.75        904
## sds(hu_spaternalage_1)     0.36      0.33     0.01     1.30       1069
##                        Rhat
## sds(spaternalage_1)    1.01
## sds(hu_spaternalage_1) 1.00
## 
## Group-Level Effects: 
## ~idParents (Number of levels: 14746) 
##                  Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
## sd(Intercept)        0.36      0.01     0.35     0.37       1239 1.01
## sd(hu_Intercept)     0.82      0.02     0.79     0.85        818 1.00
## 
## Population-Level Effects: 
##                           Estimate Est.Error l-95% CI u-95% CI Eff.Sample
## Intercept                     1.15      0.11     0.93     1.35        220
## birth_cohort1750M1755        -0.17      0.12    -0.39     0.06        484
## birth_cohort1755M1760         0.09      0.10    -0.09     0.28        226
## birth_cohort1760M1765         0.15      0.09    -0.02     0.33        177
## birth_cohort1765M1770         0.10      0.09    -0.07     0.28        166
## birth_cohort1770M1775         0.06      0.09    -0.10     0.25        166
## birth_cohort1775M1780         0.05      0.09    -0.11     0.22        165
## birth_cohort1780M1785         0.16      0.09    -0.01     0.32        156
## birth_cohort1785M1790         0.13      0.08    -0.03     0.30        151
## birth_cohort1790M1795         0.03      0.08    -0.13     0.19        158
## birth_cohort1795M1800         0.02      0.08    -0.13     0.18        142
## birth_cohort1800M1805        -0.02      0.08    -0.18     0.14        141
## birth_cohort1805M1810        -0.01      0.08    -0.16     0.14        139
## birth_cohort1810M1815         0.02      0.08    -0.13     0.17        138
## birth_cohort1815M1820         0.07      0.08    -0.08     0.22        136
## birth_cohort1820M1825         0.09      0.08    -0.06     0.24        133
## birth_cohort1825M1830         0.04      0.08    -0.10     0.20        134
## birth_cohort1830M1835         0.07      0.08    -0.08     0.23        133
## birth_cohort1835M1840         0.07      0.08    -0.08     0.23        133
## birth_cohort1840M1845         0.05      0.08    -0.10     0.21        133
## birth_cohort1845M1850         0.06      0.08    -0.09     0.21        132
## male1                         0.04      0.01     0.03     0.05       3000
## maternalage.factor1020        0.04      0.03    -0.02     0.10       3000
## maternalage.factor3559        0.06      0.01     0.04     0.08       3000
## paternalage.mean              0.06      0.02     0.01     0.10        622
## paternal_loss01               0.06      0.04    -0.01     0.13       1762
## paternal_loss15               0.03      0.02    -0.02     0.08       1283
## paternal_loss510             -0.02      0.02    -0.06     0.02       1318
## paternal_loss1015            -0.02      0.02    -0.06     0.02       1142
## paternal_loss1520            -0.08      0.02    -0.11    -0.04       1168
## paternal_loss2025            -0.03      0.02    -0.06     0.00       1199
## paternal_loss2530            -0.02      0.02    -0.05     0.01       1113
## paternal_loss3035            -0.02      0.01    -0.05     0.01       1199
## paternal_loss3540             0.02      0.01    -0.01     0.04       1202
## paternal_loss4045             0.02      0.01     0.00     0.05       1669
## maternal_loss01               0.07      0.06    -0.04     0.17       2209
## maternal_loss15               0.00      0.03    -0.05     0.06       1979
## maternal_loss510             -0.02      0.02    -0.06     0.03       1580
## maternal_loss1015            -0.03      0.02    -0.08     0.01       1776
## maternal_loss1520            -0.03      0.02    -0.07     0.00       1668
## maternal_loss2025            -0.07      0.02    -0.10    -0.04       1857
## maternal_loss2530            -0.03      0.02    -0.06     0.00       1739
## maternal_loss3035            -0.02      0.01    -0.05     0.01       1842
## maternal_loss3540             0.00      0.01    -0.02     0.03       1771
## maternal_loss4045            -0.02      0.01    -0.04     0.00       2209
## older_siblings1               0.01      0.01    -0.01     0.04        919
## older_siblings2               0.02      0.01     0.00     0.05        740
## older_siblings3               0.03      0.02    -0.01     0.07        649
## older_siblings4               0.00      0.02    -0.04     0.05        628
## older_siblings5P              0.02      0.03    -0.04     0.07        587
## nr.siblings                   0.02      0.00     0.02     0.03        659
## last_born1                   -0.02      0.01    -0.04     0.00       3000
## spaternalage_1               -0.05      0.04    -0.14     0.03       1110
## hu_Intercept                  0.09      0.28    -0.45     0.64        364
## hu_birth_cohort1750M1755      0.36      0.28    -0.16     0.92        770
## hu_birth_cohort1755M1760      0.08      0.24    -0.39     0.55        465
## hu_birth_cohort1760M1765     -0.04      0.22    -0.46     0.38        369
## hu_birth_cohort1765M1770     -0.34      0.22    -0.78     0.08        347
## hu_birth_cohort1770M1775     -0.10      0.22    -0.55     0.33        347
## hu_birth_cohort1775M1780      0.06      0.21    -0.36     0.45        314
## hu_birth_cohort1780M1785     -0.01      0.21    -0.42     0.40        340
## hu_birth_cohort1785M1790      0.14      0.20    -0.25     0.52        292
## hu_birth_cohort1790M1795      0.32      0.19    -0.05     0.69        276
## hu_birth_cohort1795M1800      0.12      0.18    -0.24     0.47        268
## hu_birth_cohort1800M1805      0.04      0.19    -0.32     0.40        268
## hu_birth_cohort1805M1810      0.00      0.18    -0.36     0.34        258
## hu_birth_cohort1810M1815      0.07      0.18    -0.29     0.40        256
## hu_birth_cohort1815M1820     -0.10      0.18    -0.45     0.23        250
## hu_birth_cohort1820M1825     -0.20      0.18    -0.56     0.13        249
## hu_birth_cohort1825M1830     -0.20      0.18    -0.55     0.13        249
## hu_birth_cohort1830M1835     -0.18      0.18    -0.54     0.15        251
## hu_birth_cohort1835M1840     -0.20      0.18    -0.55     0.12        254
## hu_birth_cohort1840M1845     -0.19      0.18    -0.54     0.14        253
## hu_birth_cohort1845M1850     -0.16      0.18    -0.51     0.18        251
## hu_male1                      0.04      0.02     0.01     0.08       3000
## hu_maternalage.factor1020     0.06      0.09    -0.12     0.24       3000
## hu_maternalage.factor3559     0.07      0.03     0.01     0.13       3000
## hu_paternalage.mean          -0.07      0.06    -0.18     0.04        670
## hu_paternal_loss01            0.87      0.09     0.68     1.05       3000
## hu_paternal_loss15            0.67      0.06     0.55     0.79       1731
## hu_paternal_loss510           0.69      0.05     0.59     0.79       1450
## hu_paternal_loss1015          0.52      0.05     0.43     0.61       1503
## hu_paternal_loss1520          0.42      0.04     0.33     0.51       1377
## hu_paternal_loss2025          0.32      0.04     0.24     0.40       1219
## hu_paternal_loss2530          0.23      0.04     0.15     0.31       1260
## hu_paternal_loss3035          0.16      0.04     0.08     0.24       1440
## hu_paternal_loss3540          0.12      0.04     0.05     0.19       1382
## hu_paternal_loss4045          0.04      0.04    -0.03     0.12       3000
## hu_maternal_loss01            1.85      0.12     1.61     2.08       3000
## hu_maternal_loss15            1.02      0.07     0.89     1.15       2186
## hu_maternal_loss510           0.87      0.06     0.75     0.98       1803
## hu_maternal_loss1015          0.81      0.06     0.69     0.91       3000
## hu_maternal_loss1520          0.67      0.05     0.58     0.77       3000
## hu_maternal_loss2025          0.46      0.05     0.37     0.55       1666
## hu_maternal_loss2530          0.32      0.04     0.24     0.40       2060
## hu_maternal_loss3035          0.24      0.04     0.17     0.32       1829
## hu_maternal_loss3540          0.14      0.04     0.07     0.21       3000
## hu_maternal_loss4045          0.08      0.03     0.01     0.14       3000
## hu_older_siblings1           -0.04      0.03    -0.10     0.03       1202
## hu_older_siblings2           -0.03      0.04    -0.11     0.05        787
## hu_older_siblings3           -0.04      0.06    -0.15     0.07        728
## hu_older_siblings4           -0.07      0.07    -0.20     0.06        703
## hu_older_siblings5P          -0.13      0.09    -0.31     0.03        635
## hu_nr.siblings                0.05      0.01     0.03     0.06        682
## hu_last_born1                 0.01      0.03    -0.04     0.06       3000
## hu_spaternalage_1             0.03      0.08    -0.13     0.21        656
##                           Rhat
## Intercept                 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.01
## birth_cohort1795M1800     1.01
## birth_cohort1800M1805     1.01
## birth_cohort1805M1810     1.01
## birth_cohort1810M1815     1.01
## birth_cohort1815M1820     1.01
## birth_cohort1820M1825     1.01
## birth_cohort1825M1830     1.01
## birth_cohort1830M1835     1.01
## birth_cohort1835M1840     1.01
## birth_cohort1840M1845     1.01
## birth_cohort1845M1850     1.01
## 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.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.01
## older_siblings3           1.01
## older_siblings4           1.01
## older_siblings5P          1.01
## nr.siblings               1.01
## last_born1                1.00
## spaternalage_1            1.00
## hu_Intercept              1.01
## hu_birth_cohort1750M1755  1.01
## hu_birth_cohort1755M1760  1.01
## hu_birth_cohort1760M1765  1.02
## hu_birth_cohort1765M1770  1.02
## hu_birth_cohort1770M1775  1.02
## hu_birth_cohort1775M1780  1.02
## hu_birth_cohort1780M1785  1.02
## hu_birth_cohort1785M1790  1.02
## hu_birth_cohort1790M1795  1.02
## hu_birth_cohort1795M1800  1.02
## hu_birth_cohort1800M1805  1.02
## hu_birth_cohort1805M1810  1.03
## hu_birth_cohort1810M1815  1.02
## hu_birth_cohort1815M1820  1.03
## hu_birth_cohort1820M1825  1.02
## hu_birth_cohort1825M1830  1.03
## hu_birth_cohort1830M1835  1.03
## hu_birth_cohort1835M1840  1.02
## hu_birth_cohort1840M1845  1.02
## hu_birth_cohort1845M1850  1.02
## hu_male1                  1.00
## hu_maternalage.factor1020 1.00
## hu_maternalage.factor3559 1.00
## hu_paternalage.mean       1.01
## hu_paternal_loss01        1.00
## hu_paternal_loss15        1.00
## hu_paternal_loss510       1.00
## hu_paternal_loss1015      1.00
## hu_paternal_loss1520      1.00
## hu_paternal_loss2025      1.00
## hu_paternal_loss2530      1.00
## hu_paternal_loss3035      1.00
## hu_paternal_loss3540      1.00
## hu_paternal_loss4045      1.00
## hu_maternal_loss01        1.00
## hu_maternal_loss15        1.00
## hu_maternal_loss510       1.00
## hu_maternal_loss1015      1.00
## hu_maternal_loss1520      1.00
## hu_maternal_loss2025      1.00
## hu_maternal_loss2530      1.00
## hu_maternal_loss3035      1.00
## hu_maternal_loss3540      1.00
## hu_maternal_loss4045      1.00
## hu_older_siblings1        1.00
## hu_older_siblings2        1.00
## hu_older_siblings3        1.00
## hu_older_siblings4        1.00
## hu_older_siblings5P       1.00
## hu_nr.siblings            1.00
## hu_last_born1             1.00
## hu_spaternalage_1         1.01
## 
## Samples were drawn using sampling(NUTS). For each parameter, Eff.Sample 
## is a crude measure of effective sample size, and Rhat is the potential 
## scale reduction factor on split chains (at convergence, Rhat = 1).

Table of fixed effects

fixed_eff = data.frame(model_summary$fixed, check.names = F)
fixed_eff$Est.Error = fixed_eff$Eff.Sample = fixed_eff$Rhat = NULL
fixed_eff$`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	3.145	2.543	3.853
birth_cohort1750M1755	0.8428	0.6777	1.058
birth_cohort1755M1760	1.098	0.9141	1.327
birth_cohort1760M1765	1.162	0.9812	1.389
birth_cohort1765M1770	1.105	0.9339	1.327
birth_cohort1770M1775	1.065	0.9004	1.278
birth_cohort1775M1780	1.052	0.8949	1.252
birth_cohort1780M1785	1.168	0.9939	1.382
birth_cohort1785M1790	1.135	0.9703	1.343
birth_cohort1790M1795	1.027	0.8786	1.21
birth_cohort1795M1800	1.023	0.8749	1.201
birth_cohort1800M1805	0.9768	0.8376	1.15
birth_cohort1805M1810	0.9906	0.8503	1.155
birth_cohort1810M1815	1.015	0.8757	1.191
birth_cohort1815M1820	1.07	0.9215	1.249
birth_cohort1820M1825	1.09	0.9381	1.273
birth_cohort1825M1830	1.045	0.901	1.22
birth_cohort1830M1835	1.072	0.9213	1.254
birth_cohort1835M1840	1.072	0.9245	1.255
birth_cohort1840M1845	1.048	0.9006	1.228
birth_cohort1845M1850	1.057	0.9119	1.235
male1	1.041	1.027	1.055
maternalage.factor1020	1.037	0.9783	1.102
maternalage.factor3559	1.065	1.043	1.087
paternalage.mean	1.057	1.014	1.101
paternal_loss01	1.061	0.9857	1.14
paternal_loss15	1.026	0.9788	1.078
paternal_loss510	0.9805	0.9394	1.024
paternal_loss1015	0.9796	0.9439	1.016
paternal_loss1520	0.9276	0.8969	0.9602
paternal_loss2025	0.9699	0.9401	1.001
paternal_loss2530	0.9768	0.948	1.006
paternal_loss3035	0.9813	0.9536	1.01
paternal_loss3540	1.018	0.9918	1.045
paternal_loss4045	1.025	0.9993	1.053
maternal_loss01	1.069	0.9605	1.19
maternal_loss15	1.004	0.9487	1.06
maternal_loss510	0.9816	0.9397	1.028
maternal_loss1015	0.9662	0.9258	1.009
maternal_loss1520	0.9668	0.93	1.005
maternal_loss2025	0.933	0.902	0.9654
maternal_loss2530	0.9744	0.9452	1.003
maternal_loss3035	0.9804	0.9537	1.007
maternal_loss3540	1.003	0.9782	1.027
maternal_loss4045	0.9806	0.9571	1.005
older_siblings1	1.014	0.9936	1.036
older_siblings2	1.024	0.9954	1.054
older_siblings3	1.029	0.9917	1.068
older_siblings4	1.003	0.9591	1.048
older_siblings5P	1.016	0.9586	1.075
nr.siblings	1.024	1.019	1.03
last_born1	0.9789	0.9613	0.997
spaternalage_1	0.9514	0.8717	1.032
hu_Intercept	1.092	0.6354	1.902
hu_birth_cohort1750M1755	1.436	0.8551	2.5
hu_birth_cohort1755M1760	1.082	0.6798	1.727
hu_birth_cohort1760M1765	0.9567	0.6331	1.461
hu_birth_cohort1765M1770	0.7139	0.4603	1.086
hu_birth_cohort1770M1775	0.9023	0.5794	1.39
hu_birth_cohort1775M1780	1.057	0.7001	1.564
hu_birth_cohort1780M1785	0.9943	0.6593	1.499
hu_birth_cohort1785M1790	1.149	0.7768	1.68
hu_birth_cohort1790M1795	1.378	0.9465	1.994
hu_birth_cohort1795M1800	1.129	0.7899	1.593
hu_birth_cohort1800M1805	1.041	0.7295	1.493
hu_birth_cohort1805M1810	0.9956	0.698	1.405
hu_birth_cohort1810M1815	1.071	0.7447	1.499
hu_birth_cohort1815M1820	0.9081	0.6383	1.259
hu_birth_cohort1820M1825	0.8175	0.5721	1.14
hu_birth_cohort1825M1830	0.8179	0.5776	1.144
hu_birth_cohort1830M1835	0.8314	0.5843	1.161
hu_birth_cohort1835M1840	0.8147	0.5758	1.129
hu_birth_cohort1840M1845	0.8256	0.5838	1.154
hu_birth_cohort1845M1850	0.8556	0.6021	1.197
hu_male1	1.046	1.007	1.086
hu_maternalage.factor1020	1.059	0.8866	1.266
hu_maternalage.factor3559	1.072	1.014	1.138
hu_paternalage.mean	0.9311	0.8317	1.038
hu_paternal_loss01	2.375	1.976	2.849
hu_paternal_loss15	1.948	1.725	2.2
hu_paternal_loss510	1.992	1.805	2.209
hu_paternal_loss1015	1.683	1.531	1.849
hu_paternal_loss1520	1.527	1.394	1.662
hu_paternal_loss2025	1.379	1.27	1.496
hu_paternal_loss2530	1.255	1.158	1.361
hu_paternal_loss3035	1.175	1.088	1.267
hu_paternal_loss3540	1.13	1.049	1.212
hu_paternal_loss4045	1.043	0.9663	1.127
hu_maternal_loss01	6.365	5.026	8
hu_maternal_loss15	2.781	2.433	3.171
hu_maternal_loss510	2.376	2.118	2.663
hu_maternal_loss1015	2.238	1.997	2.495
hu_maternal_loss1520	1.956	1.778	2.158
hu_maternal_loss2025	1.58	1.441	1.727
hu_maternal_loss2530	1.37	1.266	1.486
hu_maternal_loss3035	1.277	1.183	1.377
hu_maternal_loss3540	1.155	1.076	1.24
hu_maternal_loss4045	1.078	1.009	1.152
hu_older_siblings1	0.9643	0.9054	1.026
hu_older_siblings2	0.9746	0.8945	1.056
hu_older_siblings3	0.962	0.8616	1.071
hu_older_siblings4	0.9334	0.8154	1.067
hu_older_siblings5P	0.876	0.7332	1.03
hu_nr.siblings	1.048	1.033	1.064
hu_last_born1	1.008	0.9574	1.058
hu_spaternalage_1	1.035	0.8753	1.228

Paternal age effect

pander::pander(paternal_age_10y_effect(model))

effect	median_estimate	ci_95	ci_80
estimate father 25y	2.23	[1.76;2.78]	[1.90;2.57]
estimate father 35y	2.13	[1.67;2.69]	[1.81;2.47]
percentage change	-4.43	[-11.65; 2.93]	[ -9.10; 0.47]
OR/IRR	0.95	[0.87;1.03]	[0.90;1.00]
OR hurdle	1.03	[0.88;1.23]	[0.94;1.13]

Marginal effect plots

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

plot of chunk unnamed-chunk-45

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

plot of chunk unnamed-chunk-46

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

plot of chunk unnamed-chunk-47

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

plot of chunk unnamed-chunk-47

Rhat

Did the 6 chains converge?

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

plot of chunk unnamed-chunk-48

Effective sample size over average sample size

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

## Error in factor(x, labels = names(sort(x))): invalid 'labels'; length 29619 should be 1 or 29610

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

Model comparison

We compare the four models using an approximate leave-one-out cross-validation information criterion as implemented in brms and loo and the Watanabe-Akaike information criterion.

Approximate leave-one-out (LOO) cross-validation

	LOOIC	SE
m1	191798	758.7
m2	185251	723
m3	185226	722.9
m4	185248	723
m1 - m2	6547	217.8
m1 - m3	6572	218.1
m1 - m4	6550	218.3
m2 - m3	25.39	15.36
m2 - m4	3.12	16.38
m3 - m4	-22.27	14.58

Watanabe-Akaike information criterion

	WAIC	SE
m1	191798	758.7
m2	184258	715.3
m3	184250	715.4
m4	184259	715.4
m1 - m2	7540	219.5
m1 - m3	7548	219.8
m1 - m4	7539	219.9
m2 - m3	7.65	11.75
m2 - m4	-1.11	12.89
m3 - m4	-8.76	10.85