Québec (Registre de la Population du Québec Ancien 1630-1750), main effects

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

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: 68724) 
## 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                     2.14      0.02     2.10     2.17        874
## paternalage                  -0.02      0.00    -0.02    -0.01       3000
## birth_cohort1675M1680         0.01      0.01    -0.02     0.04        690
## birth_cohort1680M1685         0.03      0.01     0.01     0.06        679
## birth_cohort1685M1690         0.04      0.01     0.01     0.07        650
## birth_cohort1690M1695         0.03      0.01     0.01     0.06        595
## birth_cohort1695M1700         0.03      0.01     0.00     0.05        589
## birth_cohort1700M1705         0.01      0.01    -0.01     0.03        562
## birth_cohort1705M1710        -0.01      0.01    -0.03     0.02        557
## birth_cohort1710M1715        -0.01      0.01    -0.04     0.01        542
## birth_cohort1715M1720        -0.04      0.01    -0.07    -0.02        535
## birth_cohort1720M1725        -0.04      0.01    -0.06    -0.02        510
## birth_cohort1725M1730        -0.07      0.01    -0.09    -0.04        528
## birth_cohort1730M1735        -0.05      0.01    -0.07    -0.02        531
## birth_cohort1735M1740        -0.06      0.01    -0.08    -0.03        513
## male1                         0.11      0.00     0.10     0.12       3000
## maternalage.factor1420       -0.01      0.01    -0.02     0.01       3000
## maternalage.factor3550        0.01      0.01     0.00     0.02       3000
## paternal_loss01              -0.01      0.02    -0.04     0.03       2154
## paternal_loss15               0.00      0.01    -0.02     0.02       1736
## paternal_loss510              0.00      0.01    -0.02     0.02       1931
## paternal_loss1015            -0.01      0.01    -0.03     0.00       1663
## paternal_loss1520            -0.02      0.01    -0.04    -0.01       1913
## paternal_loss2025            -0.02      0.01    -0.04    -0.01       1685
## paternal_loss2530            -0.01      0.01    -0.03     0.00       1732
## paternal_loss3035            -0.02      0.01    -0.04    -0.01       1608
## paternal_loss3540            -0.02      0.01    -0.04    -0.01       1696
## paternal_loss4045             0.00      0.01    -0.02     0.01       1926
## paternal_lossunclear         -0.04      0.01    -0.05    -0.03       1557
## maternal_loss01              -0.03      0.02    -0.07     0.01       2724
## maternal_loss15              -0.02      0.01    -0.04     0.01       3000
## maternal_loss510              0.01      0.01    -0.01     0.02       3000
## maternal_loss1015            -0.01      0.01    -0.03     0.01       3000
## maternal_loss1520             0.00      0.01    -0.02     0.02       3000
## maternal_loss2025            -0.03      0.01    -0.04    -0.01       3000
## maternal_loss2530            -0.02      0.01    -0.04    -0.01       3000
## maternal_loss3035            -0.02      0.01    -0.03     0.00       3000
## maternal_loss3540            -0.01      0.01    -0.02     0.00       3000
## maternal_loss4045             0.00      0.01    -0.02     0.01       3000
## maternal_lossunclear         -0.02      0.01    -0.04    -0.01       3000
## older_siblings1              -0.02      0.01    -0.03     0.00       3000
## older_siblings2              -0.01      0.01    -0.03     0.00       2020
## older_siblings3              -0.02      0.01    -0.04    -0.01       2067
## older_siblings4               0.00      0.01    -0.02     0.01       1965
## older_siblings5P             -0.01      0.01    -0.03     0.00       1731
## nr.siblings                   0.01      0.00     0.01     0.01       2459
## last_born1                    0.01      0.01    -0.01     0.02       3000
## hu_Intercept                 -0.87      0.07    -1.00    -0.73        431
## hu_paternalage                0.00      0.01    -0.03     0.02       1759
## hu_birth_cohort1675M1680      0.06      0.06    -0.06     0.19        561
## hu_birth_cohort1680M1685      0.20      0.06     0.08     0.32        540
## hu_birth_cohort1685M1690      0.30      0.06     0.17     0.43        573
## hu_birth_cohort1690M1695      0.08      0.06    -0.04     0.19        481
## hu_birth_cohort1695M1700      0.10      0.06    -0.02     0.21        465
## hu_birth_cohort1700M1705      0.23      0.06     0.12     0.33        461
## hu_birth_cohort1705M1710      0.14      0.06     0.03     0.25        427
## hu_birth_cohort1710M1715      0.39      0.05     0.28     0.50        367
## hu_birth_cohort1715M1720      0.28      0.05     0.18     0.38        437
## hu_birth_cohort1720M1725      0.30      0.05     0.20     0.40        420
## hu_birth_cohort1725M1730      0.63      0.05     0.53     0.74        326
## hu_birth_cohort1730M1735      0.68      0.05     0.59     0.78        332
## hu_birth_cohort1735M1740      0.53      0.05     0.43     0.62        397
## hu_male1                      0.40      0.02     0.37     0.43       2748
## hu_maternalage.factor1420     0.00      0.03    -0.07     0.06       1348
## hu_maternalage.factor3550     0.07      0.02     0.02     0.11       1745
## hu_paternal_loss01            0.59      0.07     0.47     0.72       1651
## hu_paternal_loss15            0.33      0.05     0.24     0.43       1211
## hu_paternal_loss510           0.32      0.04     0.24     0.40       1147
## hu_paternal_loss1015          0.19      0.04     0.12     0.27       1164
## hu_paternal_loss1520          0.28      0.04     0.21     0.35       1080
## hu_paternal_loss2025          0.18      0.03     0.12     0.25       1120
## hu_paternal_loss2530          0.18      0.03     0.12     0.25       1137
## hu_paternal_loss3035          0.14      0.03     0.08     0.20       1128
## hu_paternal_loss3540          0.10      0.03     0.04     0.16       1143
## hu_paternal_loss4045          0.06      0.04    -0.01     0.13       1273
## hu_paternal_lossunclear       0.35      0.03     0.29     0.41        919
## hu_maternal_loss01            1.04      0.07     0.91     1.18       2375
## hu_maternal_loss15            0.42      0.04     0.34     0.50       1830
## hu_maternal_loss510           0.29      0.04     0.22     0.36       2149
## hu_maternal_loss1015          0.26      0.04     0.19     0.34       1952
## hu_maternal_loss1520          0.20      0.04     0.12     0.27       1960
## hu_maternal_loss2025          0.16      0.03     0.10     0.23       1892
## hu_maternal_loss2530          0.07      0.03     0.00     0.13       1760
## hu_maternal_loss3035          0.09      0.03     0.03     0.14       1797
## hu_maternal_loss3540          0.08      0.03     0.02     0.14       2160
## hu_maternal_loss4045          0.00      0.03    -0.05     0.06       1979
## hu_maternal_lossunclear       0.19      0.03     0.14     0.25       1672
## hu_older_siblings1           -0.05      0.03    -0.10     0.01       1246
## hu_older_siblings2           -0.08      0.03    -0.13    -0.01       1212
## hu_older_siblings3           -0.11      0.03    -0.18    -0.05       1200
## hu_older_siblings4           -0.12      0.03    -0.18    -0.05       1196
## hu_older_siblings5P          -0.15      0.03    -0.21    -0.09       1014
## hu_nr.siblings                0.02      0.00     0.01     0.02       2021
## hu_last_born1                 0.02      0.03    -0.03     0.08       2115
##                           Rhat
## Intercept                 1.00
## paternalage               1.00
## birth_cohort1675M1680     1.00
## birth_cohort1680M1685     1.00
## birth_cohort1685M1690     1.00
## birth_cohort1690M1695     1.01
## birth_cohort1695M1700     1.00
## birth_cohort1700M1705     1.00
## birth_cohort1705M1710     1.01
## birth_cohort1710M1715     1.01
## birth_cohort1715M1720     1.00
## birth_cohort1720M1725     1.00
## birth_cohort1725M1730     1.01
## birth_cohort1730M1735     1.00
## birth_cohort1735M1740     1.00
## male1                     1.00
## maternalage.factor1420    1.00
## maternalage.factor3550    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.00
## nr.siblings               1.00
## last_born1                1.00
## hu_Intercept              1.01
## hu_paternalage            1.00
## hu_birth_cohort1675M1680  1.01
## hu_birth_cohort1680M1685  1.02
## hu_birth_cohort1685M1690  1.01
## hu_birth_cohort1690M1695  1.02
## hu_birth_cohort1695M1700  1.01
## hu_birth_cohort1700M1705  1.01
## hu_birth_cohort1705M1710  1.02
## hu_birth_cohort1710M1715  1.02
## hu_birth_cohort1715M1720  1.02
## hu_birth_cohort1720M1725  1.02
## hu_birth_cohort1725M1730  1.02
## hu_birth_cohort1730M1735  1.02
## hu_birth_cohort1735M1740  1.02
## hu_male1                  1.00
## hu_maternalage.factor1420 1.00
## hu_maternalage.factor3550 1.00
## hu_paternal_loss01        1.00
## hu_paternal_loss15        1.00
## hu_paternal_loss510       1.01
## 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_paternal_lossunclear   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.01
## 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_maternal_lossunclear   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	8.468	8.199	8.748
paternalage	0.9841	0.9784	0.9902
birth_cohort1675M1680	1.013	0.9846	1.042
birth_cohort1680M1685	1.035	1.006	1.063
birth_cohort1685M1690	1.041	1.012	1.07
birth_cohort1690M1695	1.034	1.007	1.061
birth_cohort1695M1700	1.026	0.9987	1.052
birth_cohort1700M1705	1.011	0.9858	1.035
birth_cohort1705M1710	0.9936	0.9689	1.018
birth_cohort1710M1715	0.9899	0.9648	1.014
birth_cohort1715M1720	0.9584	0.9349	0.9823
birth_cohort1720M1725	0.961	0.9377	0.9837
birth_cohort1725M1730	0.9352	0.9133	0.9568
birth_cohort1730M1735	0.9552	0.9326	0.9777
birth_cohort1735M1740	0.9459	0.9242	0.968
male1	1.117	1.109	1.126
maternalage.factor1420	0.995	0.9792	1.011
maternalage.factor3550	1.013	1.001	1.025
paternal_loss01	0.9942	0.9606	1.028
paternal_loss15	1	0.9772	1.023
paternal_loss510	0.9974	0.9775	1.016
paternal_loss1015	0.9867	0.9698	1.004
paternal_loss1520	0.9763	0.9593	0.9929
paternal_loss2025	0.9785	0.9626	0.994
paternal_loss2530	0.9879	0.9725	1.003
paternal_loss3035	0.9764	0.9622	0.9917
paternal_loss3540	0.9784	0.9637	0.993
paternal_loss4045	0.9973	0.9812	1.014
paternal_lossunclear	0.9609	0.9471	0.9751
maternal_loss01	0.9695	0.9342	1.007
maternal_loss15	0.9831	0.9618	1.005
maternal_loss510	1.006	0.9873	1.024
maternal_loss1015	0.9899	0.9716	1.008
maternal_loss1520	0.9982	0.9809	1.016
maternal_loss2025	0.9723	0.9563	0.9887
maternal_loss2530	0.9759	0.9607	0.9911
maternal_loss3035	0.9826	0.9678	0.997
maternal_loss3540	0.9892	0.9764	1.003
maternal_loss4045	0.9987	0.9843	1.014
maternal_lossunclear	0.9758	0.9625	0.9894
older_siblings1	0.9828	0.9698	0.9962
older_siblings2	0.9855	0.9708	0.9995
older_siblings3	0.9796	0.9642	0.9938
older_siblings4	0.9964	0.9801	1.013
older_siblings5P	0.9886	0.9738	1.004
nr.siblings	1.007	1.005	1.008
last_born1	1.007	0.9933	1.021
hu_Intercept	0.4195	0.3673	0.4818
hu_paternalage	0.999	0.974	1.024
hu_birth_cohort1675M1680	1.067	0.9402	1.21
hu_birth_cohort1680M1685	1.224	1.085	1.381
hu_birth_cohort1685M1690	1.349	1.189	1.531
hu_birth_cohort1690M1695	1.083	0.9587	1.212
hu_birth_cohort1695M1700	1.102	0.9846	1.231
hu_birth_cohort1700M1705	1.256	1.129	1.397
hu_birth_cohort1705M1710	1.153	1.034	1.282
hu_birth_cohort1710M1715	1.478	1.329	1.643
hu_birth_cohort1715M1720	1.324	1.193	1.468
hu_birth_cohort1720M1725	1.352	1.222	1.496
hu_birth_cohort1725M1730	1.886	1.704	2.091
hu_birth_cohort1730M1735	1.979	1.798	2.183
hu_birth_cohort1735M1740	1.695	1.538	1.865
hu_male1	1.494	1.451	1.539
hu_maternalage.factor1420	0.9983	0.9346	1.064
hu_maternalage.factor3550	1.07	1.023	1.118
hu_paternal_loss01	1.81	1.593	2.055
hu_paternal_loss15	1.396	1.271	1.532
hu_paternal_loss510	1.374	1.274	1.489
hu_paternal_loss1015	1.211	1.125	1.307
hu_paternal_loss1520	1.317	1.233	1.415
hu_paternal_loss2025	1.199	1.124	1.28
hu_paternal_loss2530	1.201	1.126	1.28
hu_paternal_loss3035	1.149	1.081	1.224
hu_paternal_loss3540	1.108	1.043	1.177
hu_paternal_loss4045	1.065	0.9938	1.144
hu_paternal_lossunclear	1.42	1.337	1.514
hu_maternal_loss01	2.842	2.496	3.242
hu_maternal_loss15	1.521	1.408	1.648
hu_maternal_loss510	1.339	1.247	1.44
hu_maternal_loss1015	1.303	1.213	1.403
hu_maternal_loss1520	1.218	1.132	1.305
hu_maternal_loss2025	1.178	1.103	1.256
hu_maternal_loss2530	1.068	1.002	1.139
hu_maternal_loss3035	1.089	1.03	1.154
hu_maternal_loss3540	1.085	1.024	1.146
hu_maternal_loss4045	1.004	0.9468	1.064
hu_maternal_lossunclear	1.215	1.148	1.282
hu_older_siblings1	0.9548	0.9025	1.009
hu_older_siblings2	0.9275	0.8747	0.9856
hu_older_siblings3	0.8935	0.839	0.9497
hu_older_siblings4	0.89	0.8342	0.9541
hu_older_siblings5P	0.8628	0.8104	0.9183
hu_nr.siblings	1.017	1.011	1.022
hu_last_born1	1.025	0.9708	1.085

Paternal age effect

pander::pander(paternal_age_10y_effect(model))

effect	median_estimate	ci_95	ci_80
estimate father 25y	5.82	[5.57;6.06]	[5.65;5.98]
estimate father 35y	5.73	[5.48;5.97]	[5.56;5.89]
percentage change	-1.57	[-2.51;-0.59]	[-2.17;-0.90]
OR/IRR	0.98	[0.98;0.99]	[0.98;0.99]
OR hurdle	1.00	[0.97;1.02]	[0.98;1.01]

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/rpqa/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: 68724) 
## 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: 12205) 
##                  Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
## sd(Intercept)        0.27      0.00     0.27     0.28       1172    1
## sd(hu_Intercept)     0.63      0.01     0.60     0.66       1530    1
## 
## Population-Level Effects: 
##                           Estimate Est.Error l-95% CI u-95% CI Eff.Sample
## Intercept                     2.11      0.03     2.06     2.16       1126
## birth_cohort1675M1680         0.02      0.02    -0.01     0.05       1141
## birth_cohort1680M1685         0.05      0.02     0.02     0.08       1010
## birth_cohort1685M1690         0.05      0.02     0.02     0.08       1046
## birth_cohort1690M1695         0.05      0.02     0.02     0.08        700
## birth_cohort1695M1700         0.04      0.02     0.00     0.07        677
## birth_cohort1700M1705         0.02      0.02    -0.01     0.05        717
## birth_cohort1705M1710         0.00      0.02    -0.03     0.03        583
## birth_cohort1710M1715        -0.01      0.02    -0.04     0.02        667
## birth_cohort1715M1720        -0.04      0.02    -0.07     0.00        682
## birth_cohort1720M1725        -0.04      0.02    -0.07    -0.01        655
## birth_cohort1725M1730        -0.06      0.02    -0.10    -0.03        660
## birth_cohort1730M1735        -0.05      0.02    -0.08    -0.02        561
## birth_cohort1735M1740        -0.06      0.02    -0.09    -0.03        646
## male1                         0.11      0.00     0.10     0.12       3000
## maternalage.factor1420       -0.01      0.01    -0.03     0.01       3000
## maternalage.factor3550        0.01      0.01    -0.01     0.02       3000
## paternalage.mean             -0.02      0.01    -0.03    -0.01       1783
## paternal_loss01              -0.03      0.02    -0.07     0.01       1767
## paternal_loss15              -0.01      0.02    -0.04     0.02       1359
## paternal_loss510             -0.01      0.01    -0.03     0.02       1115
## paternal_loss1015            -0.01      0.01    -0.04     0.01       1116
## paternal_loss1520            -0.02      0.01    -0.04     0.00       1088
## paternal_loss2025            -0.02      0.01    -0.04     0.00        955
## paternal_loss2530            -0.01      0.01    -0.03     0.01       1085
## paternal_loss3035            -0.02      0.01    -0.04     0.00       1163
## paternal_loss3540            -0.02      0.01    -0.04     0.00       1374
## paternal_loss4045             0.00      0.01    -0.02     0.02       3000
## paternal_lossunclear         -0.05      0.01    -0.08    -0.03       1260
## maternal_loss01              -0.05      0.02    -0.10     0.00       3000
## maternal_loss15              -0.03      0.01    -0.05     0.00       1866
## maternal_loss510              0.01      0.01    -0.01     0.04       1666
## maternal_loss1015            -0.01      0.01    -0.03     0.02       2000
## maternal_loss1520             0.01      0.01    -0.02     0.03       2001
## maternal_loss2025            -0.02      0.01    -0.04     0.00       2038
## maternal_loss2530            -0.01      0.01    -0.03     0.01       1983
## maternal_loss3035            -0.01      0.01    -0.03     0.01       2140
## maternal_loss3540             0.00      0.01    -0.02     0.02       3000
## maternal_loss4045             0.00      0.01    -0.01     0.02       3000
## maternal_lossunclear         -0.02      0.01    -0.04     0.01       1415
## older_siblings1              -0.02      0.01    -0.04    -0.01       3000
## older_siblings2              -0.02      0.01    -0.04    -0.01       2357
## older_siblings3              -0.03      0.01    -0.05    -0.01       2221
## older_siblings4              -0.02      0.01    -0.03     0.00       2218
## older_siblings5P             -0.03      0.01    -0.04    -0.01       1580
## nr.siblings                   0.01      0.00     0.00     0.01       1876
## last_born1                    0.00      0.01    -0.01     0.02       3000
## hu_Intercept                 -0.85      0.09    -1.03    -0.68        805
## hu_birth_cohort1675M1680      0.08      0.07    -0.06     0.23        860
## hu_birth_cohort1680M1685      0.23      0.07     0.10     0.37        870
## hu_birth_cohort1685M1690      0.34      0.07     0.21     0.48        805
## hu_birth_cohort1690M1695      0.09      0.07    -0.04     0.23        766
## hu_birth_cohort1695M1700      0.11      0.07    -0.02     0.24        672
## hu_birth_cohort1700M1705      0.26      0.07     0.13     0.39        654
## hu_birth_cohort1705M1710      0.16      0.06     0.03     0.29        615
## hu_birth_cohort1710M1715      0.44      0.06     0.32     0.57        628
## hu_birth_cohort1715M1720      0.29      0.06     0.17     0.42        600
## hu_birth_cohort1720M1725      0.32      0.06     0.19     0.44        599
## hu_birth_cohort1725M1730      0.68      0.06     0.56     0.80        566
## hu_birth_cohort1730M1735      0.73      0.06     0.62     0.85        564
## hu_birth_cohort1735M1740      0.57      0.06     0.46     0.70        567
## hu_male1                      0.44      0.02     0.40     0.47       3000
## hu_maternalage.factor1420     0.04      0.04    -0.03     0.11       3000
## hu_maternalage.factor3550     0.08      0.02     0.04     0.13       3000
## hu_paternalage.mean          -0.02      0.02    -0.05     0.01       3000
## hu_paternal_loss01            0.61      0.07     0.48     0.76       3000
## hu_paternal_loss15            0.35      0.05     0.24     0.45       1779
## hu_paternal_loss510           0.32      0.05     0.23     0.42       1681
## hu_paternal_loss1015          0.20      0.04     0.11     0.28       1584
## hu_paternal_loss1520          0.30      0.04     0.22     0.39       1467
## hu_paternal_loss2025          0.19      0.04     0.11     0.26       1516
## hu_paternal_loss2530          0.19      0.04     0.11     0.27       1325
## hu_paternal_loss3035          0.15      0.04     0.08     0.22       1651
## hu_paternal_loss3540          0.11      0.04     0.03     0.18       1567
## hu_paternal_loss4045          0.07      0.04    -0.01     0.15       3000
## hu_paternal_lossunclear       0.39      0.04     0.31     0.46       1583
## hu_maternal_loss01            1.09      0.07     0.95     1.24       3000
## hu_maternal_loss15            0.42      0.05     0.32     0.51       3000
## hu_maternal_loss510           0.28      0.04     0.19     0.36       3000
## hu_maternal_loss1015          0.27      0.04     0.19     0.36       3000
## hu_maternal_loss1520          0.21      0.04     0.13     0.29       3000
## hu_maternal_loss2025          0.18      0.04     0.10     0.26       3000
## hu_maternal_loss2530          0.07      0.04    -0.01     0.14       3000
## hu_maternal_loss3035          0.10      0.04     0.03     0.17       3000
## hu_maternal_loss3540          0.09      0.03     0.02     0.15       3000
## hu_maternal_loss4045         -0.01      0.04    -0.07     0.06       3000
## hu_maternal_lossunclear       0.23      0.04     0.16     0.31       2095
## hu_older_siblings1           -0.04      0.03    -0.10     0.02       3000
## hu_older_siblings2           -0.07      0.03    -0.13     0.00       2531
## hu_older_siblings3           -0.11      0.03    -0.17    -0.04       2264
## hu_older_siblings4           -0.11      0.04    -0.18    -0.04       2335
## hu_older_siblings5P          -0.14      0.03    -0.20    -0.07       1561
## hu_nr.siblings                0.01      0.00     0.01     0.02       2142
## hu_last_born1                 0.01      0.03    -0.05     0.07       3000
##                           Rhat
## Intercept                 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.01
## paternal_loss2530         1.01
## 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
## hu_Intercept              1.00
## hu_birth_cohort1675M1680  1.00
## hu_birth_cohort1680M1685  1.00
## hu_birth_cohort1685M1690  1.00
## hu_birth_cohort1690M1695  1.00
## hu_birth_cohort1695M1700  1.00
## hu_birth_cohort1700M1705  1.00
## hu_birth_cohort1705M1710  1.00
## hu_birth_cohort1710M1715  1.00
## hu_birth_cohort1715M1720  1.00
## hu_birth_cohort1720M1725  1.00
## hu_birth_cohort1725M1730  1.00
## hu_birth_cohort1730M1735  1.00
## hu_birth_cohort1735M1740  1.00
## hu_male1                  1.00
## hu_maternalage.factor1420 1.00
## hu_maternalage.factor3550 1.00
## hu_paternalage.mean       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_paternal_lossunclear   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_maternal_lossunclear   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	8.256	7.818	8.708
birth_cohort1675M1680	1.021	0.9913	1.052
birth_cohort1680M1685	1.049	1.018	1.081
birth_cohort1685M1690	1.053	1.019	1.088
birth_cohort1690M1695	1.05	1.016	1.085
birth_cohort1695M1700	1.037	1.005	1.072
birth_cohort1700M1705	1.023	0.9919	1.056
birth_cohort1705M1710	0.9997	0.9694	1.033
birth_cohort1710M1715	0.9894	0.9586	1.023
birth_cohort1715M1720	0.9623	0.9337	0.9952
birth_cohort1720M1725	0.9603	0.9318	0.9903
birth_cohort1725M1730	0.9378	0.9093	0.9695
birth_cohort1730M1735	0.9536	0.9252	0.9848
birth_cohort1735M1740	0.9431	0.9154	0.9735
male1	1.118	1.108	1.127
maternalage.factor1420	0.9911	0.9721	1.01
maternalage.factor3550	1.006	0.9932	1.02
paternalage.mean	0.9808	0.9704	0.9914
paternal_loss01	0.9675	0.928	1.007
paternal_loss15	0.9867	0.9562	1.018
paternal_loss510	0.994	0.9668	1.021
paternal_loss1015	0.9877	0.963	1.013
paternal_loss1520	0.979	0.9565	1.002
paternal_loss2025	0.9796	0.9585	1.002
paternal_loss2530	0.9891	0.968	1.01
paternal_loss3035	0.9786	0.9594	0.9982
paternal_loss3540	0.9789	0.9602	0.9975
paternal_loss4045	0.9975	0.9776	1.018
paternal_lossunclear	0.9496	0.9258	0.974
maternal_loss01	0.9515	0.9085	0.9972
maternal_loss15	0.9742	0.9471	1.004
maternal_loss510	1.014	0.989	1.039
maternal_loss1015	0.993	0.9699	1.017
maternal_loss1520	1.007	0.9837	1.03
maternal_loss2025	0.9816	0.9598	1.004
maternal_loss2530	0.9877	0.9685	1.007
maternal_loss3035	0.9924	0.9741	1.011
maternal_loss3540	1.001	0.9846	1.018
maternal_loss4045	1.004	0.9873	1.02
maternal_lossunclear	0.9814	0.9578	1.006
older_siblings1	0.9802	0.9654	0.9949
older_siblings2	0.9781	0.9623	0.9938
older_siblings3	0.972	0.9554	0.9886
older_siblings4	0.9848	0.967	1.003
older_siblings5P	0.9736	0.9563	0.9914
nr.siblings	1.007	1.005	1.009
last_born1	1.002	0.9876	1.018
hu_Intercept	0.426	0.3579	0.5089
hu_birth_cohort1675M1680	1.084	0.9442	1.256
hu_birth_cohort1680M1685	1.263	1.102	1.452
hu_birth_cohort1685M1690	1.408	1.228	1.621
hu_birth_cohort1690M1695	1.094	0.9605	1.261
hu_birth_cohort1695M1700	1.114	0.9759	1.273
hu_birth_cohort1700M1705	1.298	1.139	1.483
hu_birth_cohort1705M1710	1.174	1.033	1.335
hu_birth_cohort1710M1715	1.553	1.371	1.767
hu_birth_cohort1715M1720	1.343	1.19	1.527
hu_birth_cohort1720M1725	1.371	1.212	1.548
hu_birth_cohort1725M1730	1.968	1.749	2.236
hu_birth_cohort1730M1735	2.084	1.854	2.348
hu_birth_cohort1735M1740	1.774	1.579	2.004
hu_male1	1.546	1.497	1.598
hu_maternalage.factor1420	1.038	0.9663	1.114
hu_maternalage.factor3550	1.086	1.036	1.14
hu_paternalage.mean	0.9816	0.9506	1.013
hu_paternal_loss01	1.849	1.61	2.129
hu_paternal_loss15	1.414	1.267	1.569
hu_paternal_loss510	1.38	1.261	1.52
hu_paternal_loss1015	1.216	1.119	1.328
hu_paternal_loss1520	1.352	1.244	1.473
hu_paternal_loss2025	1.21	1.12	1.303
hu_paternal_loss2530	1.209	1.117	1.305
hu_paternal_loss3035	1.16	1.078	1.247
hu_paternal_loss3540	1.114	1.034	1.195
hu_paternal_loss4045	1.075	0.9935	1.163
hu_paternal_lossunclear	1.473	1.366	1.592
hu_maternal_loss01	2.986	2.589	3.454
hu_maternal_loss15	1.52	1.377	1.667
hu_maternal_loss510	1.324	1.214	1.437
hu_maternal_loss1015	1.314	1.214	1.427
hu_maternal_loss1520	1.234	1.139	1.339
hu_maternal_loss2025	1.2	1.109	1.298
hu_maternal_loss2530	1.068	0.9912	1.149
hu_maternal_loss3035	1.101	1.026	1.18
hu_maternal_loss3540	1.093	1.025	1.165
hu_maternal_loss4045	0.9944	0.9284	1.065
hu_maternal_lossunclear	1.262	1.171	1.359
hu_older_siblings1	0.9577	0.9026	1.017
hu_older_siblings2	0.9331	0.8749	0.9967
hu_older_siblings3	0.899	0.8401	0.9622
hu_older_siblings4	0.8994	0.8389	0.9654
hu_older_siblings5P	0.8699	0.8163	0.9281
hu_nr.siblings	1.013	1.006	1.021
hu_last_born1	1.009	0.9517	1.07

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 24507 should be 1 or 24504

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/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: 68724) 
## 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: 12205) 
##                  Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
## sd(Intercept)        0.27      0.00     0.27     0.28       1269 1.00
## sd(hu_Intercept)     0.63      0.01     0.60     0.66       1008 1.01
## 
## Population-Level Effects: 
##                           Estimate Est.Error l-95% CI u-95% CI Eff.Sample
## Intercept                     2.12      0.03     2.06     2.17        947
## paternalage                   0.01      0.01    -0.01     0.03       1308
## birth_cohort1675M1680         0.02      0.02    -0.01     0.05       1249
## birth_cohort1680M1685         0.05      0.02     0.02     0.08       1156
## birth_cohort1685M1690         0.05      0.02     0.02     0.08       1008
## birth_cohort1690M1695         0.05      0.02     0.02     0.08        966
## birth_cohort1695M1700         0.04      0.02     0.00     0.07        838
## birth_cohort1700M1705         0.02      0.02    -0.01     0.05        733
## birth_cohort1705M1710         0.00      0.02    -0.03     0.03        709
## birth_cohort1710M1715        -0.01      0.02    -0.04     0.02        768
## birth_cohort1715M1720        -0.04      0.02    -0.07    -0.01        737
## birth_cohort1720M1725        -0.04      0.02    -0.07    -0.01        720
## birth_cohort1725M1730        -0.07      0.02    -0.10    -0.04        687
## birth_cohort1730M1735        -0.05      0.02    -0.08    -0.02        660
## birth_cohort1735M1740        -0.06      0.02    -0.09    -0.03        648
## male1                         0.11      0.00     0.10     0.12       3000
## maternalage.factor1420       -0.01      0.01    -0.03     0.01       3000
## maternalage.factor3550        0.00      0.01    -0.01     0.02       3000
## paternalage.mean             -0.03      0.01    -0.05    -0.01       1286
## paternal_loss01              -0.04      0.02    -0.08     0.01       1917
## paternal_loss15              -0.02      0.02    -0.05     0.02       1413
## paternal_loss510             -0.01      0.01    -0.04     0.02       1079
## paternal_loss1015            -0.02      0.01    -0.04     0.01       1099
## paternal_loss1520            -0.02      0.01    -0.05     0.00       1179
## paternal_loss2025            -0.02      0.01    -0.04     0.00       1287
## paternal_loss2530            -0.01      0.01    -0.03     0.01       1237
## paternal_loss3035            -0.02      0.01    -0.04     0.00       1424
## paternal_loss3540            -0.02      0.01    -0.04     0.00       1840
## paternal_loss4045             0.00      0.01    -0.02     0.02       3000
## paternal_lossunclear         -0.05      0.01    -0.08    -0.03       1207
## maternal_loss01              -0.05      0.02    -0.10    -0.01       3000
## maternal_loss15              -0.03      0.02    -0.06     0.00       1833
## maternal_loss510              0.01      0.01    -0.01     0.04       1681
## maternal_loss1015            -0.01      0.01    -0.03     0.02       1772
## maternal_loss1520             0.01      0.01    -0.02     0.03       1868
## maternal_loss2025            -0.02      0.01    -0.04     0.00       1819
## maternal_loss2530            -0.01      0.01    -0.03     0.01       2198
## maternal_loss3035            -0.01      0.01    -0.03     0.01       2407
## maternal_loss3540             0.00      0.01    -0.02     0.02       2320
## maternal_loss4045             0.00      0.01    -0.01     0.02       3000
## maternal_lossunclear         -0.02      0.01    -0.04     0.01       1390
## older_siblings1              -0.02      0.01    -0.04    -0.01       3000
## older_siblings2              -0.03      0.01    -0.04    -0.01       1861
## older_siblings3              -0.03      0.01    -0.05    -0.01       1744
## older_siblings4              -0.02      0.01    -0.04     0.00       1625
## older_siblings5P             -0.04      0.01    -0.06    -0.01       1529
## nr.siblings                   0.01      0.00     0.00     0.01       1849
## last_born1                    0.00      0.01    -0.01     0.02       3000
## hu_Intercept                 -0.80      0.09    -0.97    -0.63        663
## hu_paternalage                0.13      0.04     0.05     0.21       1810
## hu_birth_cohort1675M1680      0.08      0.07    -0.06     0.21        791
## hu_birth_cohort1680M1685      0.23      0.07     0.10     0.37        726
## hu_birth_cohort1685M1690      0.34      0.07     0.20     0.47        622
## hu_birth_cohort1690M1695      0.08      0.07    -0.05     0.22        625
## hu_birth_cohort1695M1700      0.10      0.07    -0.02     0.23        534
## hu_birth_cohort1700M1705      0.25      0.06     0.13     0.38        519
## hu_birth_cohort1705M1710      0.15      0.06     0.03     0.28        461
## hu_birth_cohort1710M1715      0.43      0.06     0.31     0.56        484
## hu_birth_cohort1715M1720      0.29      0.06     0.16     0.40        420
## hu_birth_cohort1720M1725      0.31      0.06     0.19     0.43        459
## hu_birth_cohort1725M1730      0.67      0.06     0.55     0.79        414
## hu_birth_cohort1730M1735      0.72      0.06     0.61     0.85        393
## hu_birth_cohort1735M1740      0.57      0.06     0.45     0.68        428
## hu_male1                      0.44      0.02     0.40     0.47       3000
## hu_maternalage.factor1420     0.05      0.04    -0.03     0.12       3000
## hu_maternalage.factor3550     0.03      0.03    -0.03     0.09       3000
## hu_paternalage.mean          -0.14      0.04    -0.22    -0.06       1773
## hu_paternal_loss01            0.58      0.07     0.43     0.72       3000
## hu_paternal_loss15            0.31      0.05     0.21     0.42       3000
## hu_paternal_loss510           0.29      0.05     0.20     0.39       1811
## hu_paternal_loss1015          0.17      0.04     0.09     0.26       1585
## hu_paternal_loss1520          0.28      0.04     0.20     0.37       1717
## hu_paternal_loss2025          0.18      0.04     0.10     0.26       1604
## hu_paternal_loss2530          0.18      0.04     0.11     0.25       1693
## hu_paternal_loss3035          0.14      0.04     0.07     0.22       1557
## hu_paternal_loss3540          0.11      0.04     0.03     0.18       2191
## hu_paternal_loss4045          0.07      0.04    -0.01     0.15       3000
## hu_paternal_lossunclear       0.37      0.04     0.29     0.45       1571
## hu_maternal_loss01            1.06      0.07     0.92     1.21       3000
## hu_maternal_loss15            0.39      0.05     0.30     0.49       3000
## hu_maternal_loss510           0.26      0.05     0.17     0.35       3000
## hu_maternal_loss1015          0.25      0.04     0.17     0.34       3000
## hu_maternal_loss1520          0.20      0.04     0.12     0.28       3000
## hu_maternal_loss2025          0.17      0.04     0.09     0.25       3000
## hu_maternal_loss2530          0.05      0.04    -0.02     0.13       3000
## hu_maternal_loss3035          0.09      0.04     0.01     0.16       3000
## hu_maternal_loss3540          0.08      0.03     0.02     0.15       3000
## hu_maternal_loss4045         -0.01      0.03    -0.08     0.06       3000
## hu_maternal_lossunclear       0.22      0.04     0.15     0.30       3000
## hu_older_siblings1           -0.06      0.03    -0.12     0.00       3000
## hu_older_siblings2           -0.11      0.03    -0.17    -0.04       2283
## hu_older_siblings3           -0.17      0.04    -0.24    -0.09       1866
## hu_older_siblings4           -0.19      0.04    -0.27    -0.10       1790
## hu_older_siblings5P          -0.27      0.05    -0.37    -0.17       1591
## hu_nr.siblings                0.02      0.00     0.01     0.03       1905
## hu_last_born1                 0.00      0.03    -0.06     0.05       3000
##                           Rhat
## Intercept                 1.00
## paternalage               1.00
## birth_cohort1675M1680     1.00
## birth_cohort1680M1685     1.00
## birth_cohort1685M1690     1.00
## 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.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.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.01
## 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
## hu_Intercept              1.00
## hu_paternalage            1.00
## hu_birth_cohort1675M1680  1.00
## hu_birth_cohort1680M1685  1.00
## hu_birth_cohort1685M1690  1.00
## hu_birth_cohort1690M1695  1.00
## hu_birth_cohort1695M1700  1.00
## hu_birth_cohort1700M1705  1.00
## hu_birth_cohort1705M1710  1.00
## hu_birth_cohort1710M1715  1.00
## hu_birth_cohort1715M1720  1.01
## hu_birth_cohort1720M1725  1.00
## hu_birth_cohort1725M1730  1.01
## hu_birth_cohort1730M1735  1.00
## hu_birth_cohort1735M1740  1.01
## hu_male1                  1.00
## hu_maternalage.factor1420 1.00
## hu_maternalage.factor3550 1.00
## hu_paternalage.mean       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_paternal_lossunclear   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_maternal_lossunclear   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	8.297	7.857	8.771
paternalage	1.011	0.9908	1.031
birth_cohort1675M1680	1.021	0.9899	1.054
birth_cohort1680M1685	1.049	1.016	1.083
birth_cohort1685M1690	1.052	1.018	1.088
birth_cohort1690M1695	1.049	1.016	1.084
birth_cohort1695M1700	1.036	1.004	1.069
birth_cohort1700M1705	1.022	0.9899	1.054
birth_cohort1705M1710	0.9983	0.9673	1.03
birth_cohort1710M1715	0.9881	0.9576	1.019
birth_cohort1715M1720	0.9606	0.9312	0.9914
birth_cohort1720M1725	0.9589	0.93	0.9875
birth_cohort1725M1730	0.9361	0.9068	0.9645
birth_cohort1730M1735	0.9522	0.9235	0.9811
birth_cohort1735M1740	0.9417	0.9128	0.9711
male1	1.117	1.108	1.126
maternalage.factor1420	0.992	0.9738	1.011
maternalage.factor3550	1.002	0.9856	1.017
paternalage.mean	0.9706	0.9489	0.9934
paternal_loss01	0.9625	0.9226	1.005
paternal_loss15	0.9829	0.9523	1.016
paternal_loss510	0.9897	0.9629	1.019
paternal_loss1015	0.9842	0.9602	1.009
paternal_loss1520	0.9761	0.9526	1.001
paternal_loss2025	0.9777	0.9571	0.9997
paternal_loss2530	0.9877	0.9674	1.009
paternal_loss3035	0.9776	0.9588	0.9972
paternal_loss3540	0.9782	0.96	0.9965
paternal_loss4045	0.9969	0.978	1.016
paternal_lossunclear	0.9477	0.9239	0.9721
maternal_loss01	0.9489	0.9062	0.9926
maternal_loss15	0.9725	0.9434	1.003
maternal_loss510	1.011	0.9863	1.037
maternal_loss1015	0.9916	0.9673	1.016
maternal_loss1520	1.006	0.9819	1.031
maternal_loss2025	0.9805	0.9583	1.004
maternal_loss2530	0.987	0.9666	1.007
maternal_loss3035	0.9919	0.9723	1.011
maternal_loss3540	1.001	0.9837	1.018
maternal_loss4045	1.004	0.9869	1.022
maternal_lossunclear	0.9808	0.9572	1.006
older_siblings1	0.9783	0.9629	0.993
older_siblings2	0.9745	0.9576	0.9908
older_siblings3	0.9667	0.9478	0.9856
older_siblings4	0.9781	0.9566	0.9992
older_siblings5P	0.9629	0.9378	0.9882
nr.siblings	1.007	1.005	1.01
last_born1	1.001	0.9863	1.016
hu_Intercept	0.4484	0.3787	0.5329
hu_paternalage	1.139	1.053	1.234
hu_birth_cohort1675M1680	1.081	0.9463	1.236
hu_birth_cohort1680M1685	1.258	1.101	1.443
hu_birth_cohort1685M1690	1.398	1.224	1.601
hu_birth_cohort1690M1695	1.085	0.948	1.243
hu_birth_cohort1695M1700	1.104	0.9757	1.257
hu_birth_cohort1700M1705	1.29	1.141	1.462
hu_birth_cohort1705M1710	1.165	1.032	1.323
hu_birth_cohort1710M1715	1.54	1.37	1.748
hu_birth_cohort1715M1720	1.331	1.179	1.498
hu_birth_cohort1720M1725	1.358	1.208	1.539
hu_birth_cohort1725M1730	1.948	1.732	2.207
hu_birth_cohort1730M1735	2.064	1.843	2.331
hu_birth_cohort1735M1740	1.76	1.571	1.981
hu_male1	1.546	1.495	1.595
hu_maternalage.factor1420	1.049	0.9743	1.129
hu_maternalage.factor3550	1.03	0.9713	1.092
hu_paternalage.mean	0.8659	0.8004	0.9409
hu_paternal_loss01	1.777	1.539	2.048
hu_paternal_loss15	1.369	1.236	1.525
hu_paternal_loss510	1.342	1.225	1.472
hu_paternal_loss1015	1.19	1.091	1.295
hu_paternal_loss1520	1.329	1.225	1.448
hu_paternal_loss2025	1.194	1.101	1.291
hu_paternal_loss2530	1.196	1.112	1.287
hu_paternal_loss3035	1.154	1.074	1.242
hu_paternal_loss3540	1.111	1.032	1.194
hu_paternal_loss4045	1.073	0.9897	1.158
hu_paternal_lossunclear	1.45	1.337	1.567
hu_maternal_loss01	2.882	2.514	3.348
hu_maternal_loss15	1.481	1.345	1.636
hu_maternal_loss510	1.295	1.187	1.412
hu_maternal_loss1015	1.29	1.186	1.402
hu_maternal_loss1520	1.218	1.126	1.33
hu_maternal_loss2025	1.185	1.092	1.286
hu_maternal_loss2530	1.054	0.9768	1.141
hu_maternal_loss3035	1.09	1.013	1.17
hu_maternal_loss3540	1.085	1.017	1.157
hu_maternal_loss4045	0.989	0.9221	1.057
hu_maternal_lossunclear	1.246	1.159	1.344
hu_older_siblings1	0.9402	0.8855	0.9975
hu_older_siblings2	0.898	0.841	0.9603
hu_older_siblings3	0.8478	0.7871	0.9127
hu_older_siblings4	0.831	0.7646	0.9022
hu_older_siblings5P	0.7646	0.6921	0.8461
hu_nr.siblings	1.02	1.011	1.028
hu_last_born1	0.9985	0.9435	1.054

Paternal age effect

pander::pander(paternal_age_10y_effect(model))

effect	median_estimate	ci_95	ci_80
estimate father 25y	5.69	[5.41;5.97]	[5.51;5.87]
estimate father 35y	5.52	[5.20;5.83]	[5.31;5.72]
percentage change	-3.00	[-6.08; 0.24]	[-4.97;-0.90]
OR/IRR	1.01	[0.99;1.03]	[1.00;1.02]
OR hurdle	1.14	[1.05;1.23]	[1.08;1.20]

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 24509 should be 1 or 24505

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/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 67 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: 68724) 
## 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.1      0.09     0.00     0.32        522
## sds(hu_spaternalage_1)      0.8      0.46     0.24     1.97        915
##                        Rhat
## sds(spaternalage_1)    1.01
## sds(hu_spaternalage_1) 1.00
## 
## Group-Level Effects: 
## ~idParents (Number of levels: 12205) 
##                  Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
## sd(Intercept)        0.27      0.00     0.27     0.28        848    1
## sd(hu_Intercept)     0.63      0.01     0.60     0.65       1049    1
## 
## Population-Level Effects: 
##                           Estimate Est.Error l-95% CI u-95% CI Eff.Sample
## Intercept                     2.15      0.05     2.06     2.26        885
## birth_cohort1675M1680         0.02      0.02    -0.01     0.05        734
## birth_cohort1680M1685         0.05      0.02     0.02     0.08        743
## birth_cohort1685M1690         0.05      0.02     0.02     0.08        702
## birth_cohort1690M1695         0.05      0.02     0.02     0.08        558
## birth_cohort1695M1700         0.04      0.02     0.00     0.07        552
## birth_cohort1700M1705         0.02      0.02    -0.01     0.05        475
## birth_cohort1705M1710         0.00      0.02    -0.03     0.03        448
## birth_cohort1710M1715        -0.01      0.02    -0.04     0.02        424
## birth_cohort1715M1720        -0.04      0.02    -0.07    -0.01        425
## birth_cohort1720M1725        -0.04      0.02    -0.07    -0.01        405
## birth_cohort1725M1730        -0.07      0.02    -0.10    -0.04        423
## birth_cohort1730M1735        -0.05      0.02    -0.08    -0.02        399
## birth_cohort1735M1740        -0.06      0.02    -0.09    -0.03        440
## male1                         0.11      0.00     0.10     0.12       3000
## maternalage.factor1420       -0.01      0.01    -0.03     0.01       3000
## maternalage.factor3550        0.00      0.01    -0.01     0.02       2284
## paternalage.mean             -0.03      0.01    -0.05    -0.01       1255
## paternal_loss01              -0.04      0.02    -0.08     0.00       1144
## paternal_loss15              -0.02      0.02    -0.05     0.01        928
## paternal_loss510             -0.01      0.01    -0.04     0.02        768
## paternal_loss1015            -0.02      0.01    -0.04     0.01        845
## paternal_loss1520            -0.02      0.01    -0.05     0.00        745
## paternal_loss2025            -0.02      0.01    -0.04     0.00        795
## paternal_loss2530            -0.01      0.01    -0.03     0.01        760
## paternal_loss3035            -0.02      0.01    -0.04     0.00        958
## paternal_loss3540            -0.02      0.01    -0.04     0.00        951
## paternal_loss4045             0.00      0.01    -0.02     0.02       1532
## paternal_lossunclear         -0.05      0.01    -0.08    -0.03        792
## maternal_loss01              -0.05      0.02    -0.10    -0.01       1719
## maternal_loss15              -0.03      0.02    -0.06     0.00       1444
## maternal_loss510              0.01      0.01    -0.01     0.04       1237
## maternal_loss1015            -0.01      0.01    -0.03     0.02       1200
## maternal_loss1520             0.01      0.01    -0.02     0.03       1282
## maternal_loss2025            -0.02      0.01    -0.04     0.00       1258
## maternal_loss2530            -0.01      0.01    -0.03     0.01       1325
## maternal_loss3035            -0.01      0.01    -0.03     0.01       1375
## maternal_loss3540             0.00      0.01    -0.02     0.02       1762
## maternal_loss4045             0.00      0.01    -0.01     0.02       3000
## maternal_lossunclear         -0.02      0.01    -0.05     0.01        819
## older_siblings1              -0.02      0.01    -0.04    -0.01       1970
## older_siblings2              -0.02      0.01    -0.04    -0.01       1504
## older_siblings3              -0.03      0.01    -0.05    -0.01       1340
## older_siblings4              -0.02      0.01    -0.04     0.00       1275
## older_siblings5P             -0.04      0.01    -0.06    -0.01       1260
## nr.siblings                   0.01      0.00     0.00     0.01        883
## last_born1                    0.00      0.01    -0.01     0.02       3000
## spaternalage_1                0.02      0.02    -0.01     0.09        396
## hu_Intercept                 -0.20      0.18    -0.55     0.17        529
## hu_birth_cohort1675M1680      0.08      0.07    -0.05     0.21        382
## hu_birth_cohort1680M1685      0.23      0.07     0.10     0.36        379
## hu_birth_cohort1685M1690      0.34      0.07     0.21     0.48        357
## hu_birth_cohort1690M1695      0.10      0.07    -0.03     0.23        290
## hu_birth_cohort1695M1700      0.12      0.07    -0.01     0.25        255
## hu_birth_cohort1700M1705      0.27      0.06     0.15     0.39        272
## hu_birth_cohort1705M1710      0.16      0.06     0.04     0.29        266
## hu_birth_cohort1710M1715      0.44      0.06     0.32     0.56        241
## hu_birth_cohort1715M1720      0.30      0.06     0.18     0.42        257
## hu_birth_cohort1720M1725      0.32      0.06     0.21     0.44        251
## hu_birth_cohort1725M1730      0.68      0.06     0.57     0.80        244
## hu_birth_cohort1730M1735      0.74      0.06     0.62     0.85        249
## hu_birth_cohort1735M1740      0.57      0.06     0.46     0.69        245
## hu_male1                      0.44      0.02     0.40     0.47       3000
## hu_maternalage.factor1420     0.05      0.04    -0.02     0.13       3000
## hu_maternalage.factor3550     0.04      0.03    -0.02     0.09       1567
## hu_paternalage.mean          -0.17      0.04    -0.26    -0.09        729
## hu_paternal_loss01            0.58      0.07     0.44     0.72       1406
## hu_paternal_loss15            0.31      0.05     0.21     0.42       1088
## hu_paternal_loss510           0.29      0.05     0.20     0.38        876
## hu_paternal_loss1015          0.17      0.04     0.08     0.26        888
## hu_paternal_loss1520          0.28      0.04     0.19     0.36        670
## hu_paternal_loss2025          0.16      0.04     0.09     0.24        740
## hu_paternal_loss2530          0.16      0.04     0.09     0.24        690
## hu_paternal_loss3035          0.12      0.04     0.05     0.20        874
## hu_paternal_loss3540          0.09      0.04     0.02     0.16        974
## hu_paternal_loss4045          0.06      0.04    -0.02     0.13        876
## hu_paternal_lossunclear       0.36      0.04     0.29     0.44        699
## hu_maternal_loss01            1.06      0.07     0.93     1.20       3000
## hu_maternal_loss15            0.39      0.05     0.30     0.49       1359
## hu_maternal_loss510           0.26      0.04     0.18     0.35       1497
## hu_maternal_loss1015          0.26      0.04     0.17     0.34       1466
## hu_maternal_loss1520          0.20      0.04     0.12     0.29       1664
## hu_maternal_loss2025          0.17      0.04     0.09     0.25       1374
## hu_maternal_loss2530          0.05      0.04    -0.02     0.13       1629
## hu_maternal_loss3035          0.09      0.04     0.01     0.15       1518
## hu_maternal_loss3540          0.08      0.03     0.02     0.14       1569
## hu_maternal_loss4045         -0.01      0.03    -0.07     0.05       3000
## hu_maternal_lossunclear       0.22      0.04     0.14     0.30       1614
## hu_older_siblings1           -0.08      0.03    -0.14    -0.02       1527
## hu_older_siblings2           -0.15      0.03    -0.21    -0.08       1053
## hu_older_siblings3           -0.22      0.04    -0.30    -0.14        830
## hu_older_siblings4           -0.26      0.05    -0.35    -0.17        759
## hu_older_siblings5P          -0.35      0.05    -0.46    -0.25        721
## hu_nr.siblings                0.02      0.00     0.01     0.03       1114
## hu_last_born1                 0.00      0.03    -0.06     0.06       3000
## hu_spaternalage_1             0.11      0.11    -0.12     0.34        922
##                           Rhat
## Intercept                 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.02
## birth_cohort1715M1720     1.02
## birth_cohort1720M1725     1.02
## birth_cohort1725M1730     1.02
## birth_cohort1730M1735     1.02
## birth_cohort1735M1740     1.02
## male1                     1.00
## maternalage.factor1420    1.00
## maternalage.factor3550    1.00
## paternalage.mean          1.00
## paternal_loss01           1.01
## paternal_loss15           1.01
## paternal_loss510          1.01
## paternal_loss1015         1.01
## paternal_loss1520         1.01
## paternal_loss2025         1.01
## paternal_loss2530         1.01
## paternal_loss3035         1.01
## paternal_loss3540         1.01
## paternal_loss4045         1.00
## paternal_lossunclear      1.01
## maternal_loss01           1.00
## maternal_loss15           1.00
## maternal_loss510          1.00
## maternal_loss1015         1.00
## maternal_loss1520         1.00
## maternal_loss2025         1.00
## maternal_loss2530         1.00
## maternal_loss3035         1.00
## maternal_loss3540         1.00
## maternal_loss4045         1.00
## maternal_lossunclear      1.01
## 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
## spaternalage_1            1.02
## hu_Intercept              1.01
## hu_birth_cohort1675M1680  1.02
## hu_birth_cohort1680M1685  1.02
## hu_birth_cohort1685M1690  1.02
## hu_birth_cohort1690M1695  1.02
## hu_birth_cohort1695M1700  1.03
## hu_birth_cohort1700M1705  1.03
## hu_birth_cohort1705M1710  1.03
## hu_birth_cohort1710M1715  1.03
## hu_birth_cohort1715M1720  1.02
## hu_birth_cohort1720M1725  1.03
## hu_birth_cohort1725M1730  1.03
## hu_birth_cohort1730M1735  1.03
## hu_birth_cohort1735M1740  1.03
## hu_male1                  1.00
## hu_maternalage.factor1420 1.00
## hu_maternalage.factor3550 1.00
## hu_paternalage.mean       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_paternal_lossunclear   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_maternal_lossunclear   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.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	8.628	7.824	9.56
birth_cohort1675M1680	1.021	0.9906	1.053
birth_cohort1680M1685	1.049	1.016	1.082
birth_cohort1685M1690	1.052	1.018	1.088
birth_cohort1690M1695	1.049	1.017	1.082
birth_cohort1695M1700	1.036	1.004	1.068
birth_cohort1700M1705	1.022	0.9906	1.054
birth_cohort1705M1710	0.9979	0.9685	1.03
birth_cohort1710M1715	0.988	0.9575	1.021
birth_cohort1715M1720	0.9604	0.9317	0.9899
birth_cohort1720M1725	0.9581	0.929	0.9879
birth_cohort1725M1730	0.9357	0.9073	0.9649
birth_cohort1730M1735	0.9518	0.9222	0.9821
birth_cohort1735M1740	0.9413	0.9131	0.9714
male1	1.118	1.108	1.127
maternalage.factor1420	0.9915	0.9729	1.01
maternalage.factor3550	1.001	0.9863	1.017
paternalage.mean	0.9706	0.9484	0.9925
paternal_loss01	0.9626	0.9243	1.005
paternal_loss15	0.9822	0.9513	1.015
paternal_loss510	0.9894	0.963	1.018
paternal_loss1015	0.9844	0.9593	1.012
paternal_loss1520	0.9765	0.9532	1
paternal_loss2025	0.978	0.9566	0.9999
paternal_loss2530	0.9879	0.9672	1.009
paternal_loss3035	0.9779	0.9585	0.997
paternal_loss3540	0.9782	0.96	0.9964
paternal_loss4045	0.9967	0.9776	1.015
paternal_lossunclear	0.948	0.9239	0.9725
maternal_loss01	0.9496	0.9067	0.9947
maternal_loss15	0.9727	0.9433	1.004
maternal_loss510	1.012	0.9859	1.039
maternal_loss1015	0.9925	0.969	1.018
maternal_loss1520	1.006	0.9838	1.03
maternal_loss2025	0.9813	0.9596	1.004
maternal_loss2530	0.9877	0.9676	1.009
maternal_loss3035	0.9921	0.9733	1.011
maternal_loss3540	1.001	0.9848	1.018
maternal_loss4045	1.004	0.9877	1.02
maternal_lossunclear	0.9806	0.9557	1.006
older_siblings1	0.9793	0.9642	0.995
older_siblings2	0.976	0.9579	0.9936
older_siblings3	0.9683	0.9478	0.9886
older_siblings4	0.9795	0.9561	1.004
older_siblings5P	0.9641	0.9375	0.9919
nr.siblings	1.007	1.005	1.01
last_born1	1.002	0.987	1.017
spaternalage_1	1.023	0.9882	1.095
hu_Intercept	0.8198	0.5781	1.182
hu_birth_cohort1675M1680	1.078	0.951	1.228
hu_birth_cohort1680M1685	1.26	1.102	1.44
hu_birth_cohort1685M1690	1.411	1.235	1.623
hu_birth_cohort1690M1695	1.104	0.9682	1.263
hu_birth_cohort1695M1700	1.123	0.9873	1.285
hu_birth_cohort1700M1705	1.307	1.161	1.481
hu_birth_cohort1705M1710	1.178	1.042	1.338
hu_birth_cohort1710M1715	1.556	1.38	1.757
hu_birth_cohort1715M1720	1.347	1.202	1.521
hu_birth_cohort1720M1725	1.38	1.228	1.554
hu_birth_cohort1725M1730	1.979	1.766	2.229
hu_birth_cohort1730M1735	2.092	1.866	2.351
hu_birth_cohort1735M1740	1.777	1.583	1.998
hu_male1	1.545	1.496	1.597
hu_maternalage.factor1420	1.056	0.9804	1.136
hu_maternalage.factor3550	1.039	0.9807	1.099
hu_paternalage.mean	0.8407	0.7748	0.9095
hu_paternal_loss01	1.781	1.56	2.045
hu_paternal_loss15	1.37	1.234	1.523
hu_paternal_loss510	1.336	1.217	1.466
hu_paternal_loss1015	1.182	1.085	1.293
hu_paternal_loss1520	1.317	1.213	1.434
hu_paternal_loss2025	1.179	1.09	1.273
hu_paternal_loss2530	1.175	1.091	1.266
hu_paternal_loss3035	1.132	1.055	1.219
hu_paternal_loss3540	1.09	1.015	1.17
hu_paternal_loss4045	1.057	0.9773	1.143
hu_paternal_lossunclear	1.436	1.336	1.557
hu_maternal_loss01	2.89	2.536	3.315
hu_maternal_loss15	1.483	1.35	1.637
hu_maternal_loss510	1.297	1.193	1.415
hu_maternal_loss1015	1.295	1.189	1.411
hu_maternal_loss1520	1.224	1.129	1.335
hu_maternal_loss2025	1.188	1.097	1.283
hu_maternal_loss2530	1.055	0.9769	1.136
hu_maternal_loss3035	1.089	1.013	1.167
hu_maternal_loss3540	1.084	1.022	1.154
hu_maternal_loss4045	0.9901	0.9281	1.055
hu_maternal_lossunclear	1.248	1.155	1.349
hu_older_siblings1	0.922	0.8694	0.9765
hu_older_siblings2	0.863	0.8072	0.9257
hu_older_siblings3	0.8011	0.742	0.8655
hu_older_siblings4	0.774	0.7074	0.8449
hu_older_siblings5P	0.703	0.6328	0.7801
hu_nr.siblings	1.023	1.015	1.031
hu_last_born1	1	0.945	1.06
hu_spaternalage_1	1.121	0.8893	1.411

Paternal age effect

pander::pander(paternal_age_10y_effect(model))

effect	median_estimate	ci_95	ci_80
estimate father 25y	5.79	[5.50;6.09]	[5.60;5.99]
estimate father 35y	5.33	[5.00;5.69]	[5.10;5.56]
percentage change	-7.98	[-11.85; -4.12]	[-10.49; -5.23]
OR/IRR	1.02	[0.99;1.10]	[1.00;1.06]
OR hurdle	1.12	[0.89;1.41]	[0.97;1.28]

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 24527 should be 1 or 24523

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/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	299887	992.2
m2	293615	957.5
m3	293619	957.7
m4	293606	957.7
m1 - m2	6272	299.1
m1 - m3	6268	299.3
m1 - m4	6281	299.3
m2 - m3	-3.66	17.39
m2 - m4	9.83	19.87
m3 - m4	13.49	18.8

Watanabe-Akaike information criterion

	WAIC	SE
m1	299887	992.2
m2	292565	950.5
m3	292572	950.7
m4	292556	950.7
m1 - m2	7322	299.6
m1 - m3	7315	299.9
m1 - m4	7331	299.9
m2 - m3	-7.05	13.1
m2 - m4	8.81	16.28
m3 - m4	15.86	14.79