Québec sensitivity analyses

Loading details

source("0__helpers.R")
opts_chunk$set(warning=TRUE, cache=F,cache.lazy=F,tidy=FALSE,autodep=TRUE,dev=c('png','pdf'),fig.width=20,fig.height=12.5,out.width='1440px',out.height='900px')

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 are the same as our main model m3, except for the noted changes to test robustness.

s1: Mediation via age

Here, we tested whether the effect on reproductive success is mediated by age (mortality). We entered an unknown age for people who did not have a death date on their records and most likely outlived the observation period of the church records.

Model summary

Full summary

model_summary = summary(model, use_cache = FALSE, priors = TRUE)
print(model_summary)
##  Family: hurdle_poisson(log) 
## Formula: children ~ paternalage + age + birth_cohort + male + maternalage.factor + paternalage.mean + paternal_loss + maternal_loss + older_siblings + nr.siblings + last_born + (1 | idParents) 
##          hu ~ paternalage + age + 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 = 1000; warmup = 500; thin = 1; 
##          total post-warmup samples = 3000
##    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.23      0.00     0.23     0.24       1098 1.00
## sd(hu_Intercept)     0.69      0.02     0.65     0.74        899 1.01
## 
## Population-Level Effects: 
##                           Estimate Est.Error l-95% CI u-95% CI Eff.Sample
## Intercept                     1.94      0.03     1.89     1.99       1475
## paternalage                   0.01      0.01    -0.01     0.03       1899
## age50112                      0.35      0.01     0.34     0.37       3000
## age025                       -1.52      0.03    -1.59    -1.46       3000
## ageunknown                    0.03      0.01     0.02     0.05       3000
## birth_cohort1675M1680         0.02      0.02    -0.01     0.05       1298
## birth_cohort1680M1685         0.03      0.02     0.00     0.06       1262
## birth_cohort1685M1690         0.01      0.02    -0.02     0.04       1108
## birth_cohort1690M1695         0.01      0.02    -0.02     0.04        955
## birth_cohort1695M1700         0.01      0.02    -0.02     0.04        932
## birth_cohort1700M1705         0.00      0.02    -0.03     0.03        905
## birth_cohort1705M1710        -0.02      0.02    -0.05     0.01        850
## birth_cohort1710M1715        -0.02      0.02    -0.05     0.01        828
## birth_cohort1715M1720        -0.05      0.02    -0.08    -0.02        835
## birth_cohort1720M1725        -0.05      0.01    -0.08    -0.02        826
## birth_cohort1725M1730        -0.09      0.01    -0.12    -0.06        841
## birth_cohort1730M1735        -0.07      0.01    -0.10    -0.04        820
## birth_cohort1735M1740        -0.08      0.01    -0.11    -0.05        815
## male1                         0.06      0.00     0.05     0.07       3000
## maternalage.factor1420       -0.01      0.01    -0.03     0.01       3000
## maternalage.factor3550        0.00      0.01    -0.02     0.02       3000
## paternalage.mean             -0.03      0.01    -0.05    -0.01       1970
## paternal_loss01              -0.02      0.02    -0.06     0.02       3000
## paternal_loss15              -0.01      0.02    -0.04     0.02       1797
## paternal_loss510              0.01      0.01    -0.01     0.04       1506
## paternal_loss1015             0.00      0.01    -0.02     0.02       1714
## paternal_loss1520             0.00      0.01    -0.02     0.02       1557
## paternal_loss2025             0.00      0.01    -0.02     0.02       1710
## paternal_loss2530             0.00      0.01    -0.02     0.02       1623
## paternal_loss3035            -0.01      0.01    -0.03     0.00       1860
## paternal_loss3540            -0.02      0.01    -0.04     0.00       2257
## paternal_loss4045             0.00      0.01    -0.02     0.02       3000
## paternal_lossunclear         -0.04      0.01    -0.06    -0.02       1495
## maternal_loss01              -0.02      0.02    -0.07     0.02       3000
## maternal_loss15               0.00      0.01    -0.03     0.03       2350
## maternal_loss510              0.01      0.01    -0.01     0.04       2484
## maternal_loss1015             0.00      0.01    -0.02     0.03       2289
## maternal_loss1520             0.02      0.01     0.00     0.04       2218
## maternal_loss2025             0.00      0.01    -0.02     0.02       2668
## maternal_loss2530             0.00      0.01    -0.02     0.02       2274
## maternal_loss3035             0.00      0.01    -0.01     0.02       2563
## maternal_loss3540             0.00      0.01    -0.01     0.02       3000
## maternal_loss4045             0.01      0.01    -0.01     0.03       3000
## maternal_lossunclear         -0.01      0.01    -0.03     0.01       1796
## older_siblings1              -0.03      0.01    -0.04    -0.01       3000
## older_siblings2              -0.03      0.01    -0.05    -0.01       2650
## older_siblings3              -0.04      0.01    -0.06    -0.02       2369
## older_siblings4              -0.03      0.01    -0.05     0.00       2268
## older_siblings5P             -0.04      0.01    -0.07    -0.02       2053
## nr.siblings                   0.01      0.00     0.00     0.01       2189
## last_born1                    0.01      0.01    -0.01     0.02       3000
## hu_Intercept                 -2.04      0.12    -2.27    -1.80       1230
## hu_paternalage                0.04      0.06    -0.08     0.15       1404
## hu_age50112                  -0.73      0.03    -0.80    -0.67       3000
## hu_age025                     4.84      0.05     4.74     4.94       3000
## hu_ageunknown                 2.12      0.04     2.05     2.19       3000
## hu_birth_cohort1675M1680     -0.01      0.09    -0.19     0.18       1173
## hu_birth_cohort1680M1685      0.07      0.10    -0.12     0.25       1036
## hu_birth_cohort1685M1690      0.15      0.10    -0.05     0.33       1181
## hu_birth_cohort1690M1695      0.08      0.09    -0.10     0.25        974
## hu_birth_cohort1695M1700      0.00      0.09    -0.18     0.17        939
## hu_birth_cohort1700M1705      0.03      0.09    -0.14     0.20        931
## hu_birth_cohort1705M1710     -0.01      0.09    -0.18     0.15        876
## hu_birth_cohort1710M1715      0.08      0.09    -0.09     0.24        946
## hu_birth_cohort1715M1720      0.15      0.09    -0.01     0.32        870
## hu_birth_cohort1720M1725      0.12      0.08    -0.03     0.28        862
## hu_birth_cohort1725M1730      0.27      0.08     0.11     0.43        815
## hu_birth_cohort1730M1735      0.33      0.08     0.17     0.48        825
## hu_birth_cohort1735M1740      0.27      0.08     0.11     0.42        808
## hu_male1                      0.50      0.02     0.45     0.55       3000
## hu_maternalage.factor1420    -0.06      0.06    -0.18     0.04       3000
## hu_maternalage.factor3550     0.00      0.04    -0.09     0.08       3000
## hu_paternalage.mean           0.01      0.06    -0.11     0.13       1762
## hu_paternal_loss01            0.20      0.11     0.00     0.41       3000
## hu_paternal_loss15            0.14      0.08    -0.02     0.30       3000
## hu_paternal_loss510           0.07      0.07    -0.08     0.20       2175
## hu_paternal_loss1015          0.01      0.07    -0.12     0.14       1919
## hu_paternal_loss1520          0.11      0.06     0.00     0.23       1749
## hu_paternal_loss2025          0.05      0.06    -0.06     0.16       1545
## hu_paternal_loss2530          0.11      0.06     0.00     0.22       1600
## hu_paternal_loss3035          0.12      0.05     0.01     0.22       1679
## hu_paternal_loss3540          0.09      0.05    -0.02     0.20       1903
## hu_paternal_loss4045          0.08      0.06    -0.03     0.20       3000
## hu_paternal_lossunclear       0.24      0.06     0.14     0.35       1565
## hu_maternal_loss01            0.31      0.11     0.10     0.52       3000
## hu_maternal_loss15            0.00      0.07    -0.13     0.14       3000
## hu_maternal_loss510           0.12      0.06     0.00     0.25       3000
## hu_maternal_loss1015          0.11      0.06     0.00     0.23       3000
## hu_maternal_loss1520         -0.02      0.06    -0.14     0.11       3000
## hu_maternal_loss2025          0.06      0.06    -0.06     0.16       3000
## hu_maternal_loss2530         -0.06      0.06    -0.17     0.04       3000
## hu_maternal_loss3035          0.00      0.05    -0.09     0.10       3000
## hu_maternal_loss3540          0.06      0.05    -0.03     0.15       3000
## hu_maternal_loss4045         -0.04      0.05    -0.13     0.06       3000
## hu_maternal_lossunclear       0.10      0.05     0.00     0.20       3000
## hu_older_siblings1            0.10      0.05     0.00     0.19       3000
## hu_older_siblings2            0.07      0.05    -0.03     0.17       2201
## hu_older_siblings3            0.07      0.06    -0.05     0.18       1919
## hu_older_siblings4            0.07      0.06    -0.06     0.19       1759
## hu_older_siblings5P           0.04      0.08    -0.11     0.20       1541
## hu_nr.siblings                0.00      0.01    -0.01     0.01       2426
## hu_last_born1                -0.03      0.04    -0.11     0.06       3000
##                           Rhat
## Intercept                    1
## paternalage                  1
## age50112                     1
## age025                       1
## ageunknown                   1
## birth_cohort1675M1680        1
## birth_cohort1680M1685        1
## birth_cohort1685M1690        1
## birth_cohort1690M1695        1
## birth_cohort1695M1700        1
## birth_cohort1700M1705        1
## birth_cohort1705M1710        1
## birth_cohort1710M1715        1
## birth_cohort1715M1720        1
## birth_cohort1720M1725        1
## birth_cohort1725M1730        1
## birth_cohort1730M1735        1
## birth_cohort1735M1740        1
## male1                        1
## maternalage.factor1420       1
## maternalage.factor3550       1
## paternalage.mean             1
## paternal_loss01              1
## paternal_loss15              1
## paternal_loss510             1
## paternal_loss1015            1
## paternal_loss1520            1
## paternal_loss2025            1
## paternal_loss2530            1
## paternal_loss3035            1
## paternal_loss3540            1
## paternal_loss4045            1
## paternal_lossunclear         1
## maternal_loss01              1
## maternal_loss15              1
## maternal_loss510             1
## maternal_loss1015            1
## maternal_loss1520            1
## maternal_loss2025            1
## maternal_loss2530            1
## maternal_loss3035            1
## maternal_loss3540            1
## maternal_loss4045            1
## maternal_lossunclear         1
## older_siblings1              1
## older_siblings2              1
## older_siblings3              1
## older_siblings4              1
## older_siblings5P             1
## nr.siblings                  1
## last_born1                   1
## hu_Intercept                 1
## hu_paternalage               1
## hu_age50112                  1
## hu_age025                    1
## hu_ageunknown                1
## hu_birth_cohort1675M1680     1
## hu_birth_cohort1680M1685     1
## hu_birth_cohort1685M1690     1
## hu_birth_cohort1690M1695     1
## hu_birth_cohort1695M1700     1
## hu_birth_cohort1700M1705     1
## hu_birth_cohort1705M1710     1
## hu_birth_cohort1710M1715     1
## hu_birth_cohort1715M1720     1
## hu_birth_cohort1720M1725     1
## hu_birth_cohort1725M1730     1
## hu_birth_cohort1730M1735     1
## hu_birth_cohort1735M1740     1
## hu_male1                     1
## hu_maternalage.factor1420    1
## hu_maternalage.factor3550    1
## hu_paternalage.mean          1
## hu_paternal_loss01           1
## hu_paternal_loss15           1
## hu_paternal_loss510          1
## hu_paternal_loss1015         1
## hu_paternal_loss1520         1
## hu_paternal_loss2025         1
## hu_paternal_loss2530         1
## hu_paternal_loss3035         1
## hu_paternal_loss3540         1
## hu_paternal_loss4045         1
## hu_paternal_lossunclear      1
## hu_maternal_loss01           1
## hu_maternal_loss15           1
## hu_maternal_loss510          1
## hu_maternal_loss1015         1
## hu_maternal_loss1520         1
## hu_maternal_loss2025         1
## hu_maternal_loss2530         1
## hu_maternal_loss3035         1
## hu_maternal_loss3540         1
## hu_maternal_loss4045         1
## hu_maternal_lossunclear      1
## hu_older_siblings1           1
## hu_older_siblings2           1
## hu_older_siblings3           1
## hu_older_siblings4           1
## hu_older_siblings5P          1
## hu_nr.siblings               1
## hu_last_born1                1
## 
## Samples were drawn using sampling(NUTS). For each parameter, Eff.Sample 
## is a crude measure of effective sample size, and Rhat is the potential 
## scale reduction factor on split chains (at convergence, Rhat = 1).

Table of fixed effects

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 6.971 6.635 7.32
paternalage 1.008 0.9875 1.028
age50112 1.425 1.41 1.441
age025 0.2187 0.2044 0.2334
ageunknown 1.033 1.016 1.05
birth_cohort1675M1680 1.017 0.9866 1.049
birth_cohort1680M1685 1.027 0.9956 1.059
birth_cohort1685M1690 1.012 0.9788 1.045
birth_cohort1690M1695 1.012 0.981 1.044
birth_cohort1695M1700 1.006 0.9756 1.038
birth_cohort1700M1705 0.9998 0.97 1.031
birth_cohort1705M1710 0.9827 0.9536 1.013
birth_cohort1710M1715 0.9798 0.9507 1.01
birth_cohort1715M1720 0.9529 0.925 0.9821
birth_cohort1720M1725 0.9478 0.9204 0.9758
birth_cohort1725M1730 0.9173 0.8909 0.9444
birth_cohort1730M1735 0.9357 0.9091 0.9628
birth_cohort1735M1740 0.9235 0.8978 0.951
male1 1.064 1.055 1.073
maternalage.factor1420 0.9881 0.9703 1.007
maternalage.factor3550 0.9993 0.9836 1.016
paternalage.mean 0.9703 0.9492 0.9928
paternal_loss01 0.9787 0.9396 1.019
paternal_loss15 0.9934 0.9647 1.023
paternal_loss510 1.012 0.9856 1.039
paternal_loss1015 0.9991 0.9755 1.022
paternal_loss1520 1.001 0.9777 1.024
paternal_loss2025 0.9971 0.9755 1.018
paternal_loss2530 1.001 0.9812 1.021
paternal_loss3035 0.9851 0.9667 1.004
paternal_loss3540 0.9797 0.9621 0.9977
paternal_loss4045 1.002 0.9833 1.021
paternal_lossunclear 0.9605 0.9383 0.9831
maternal_loss01 0.9777 0.9366 1.023
maternal_loss15 0.9986 0.971 1.028
maternal_loss510 1.014 0.9891 1.041
maternal_loss1015 1.004 0.981 1.027
maternal_loss1520 1.02 0.9982 1.043
maternal_loss2025 1 0.979 1.022
maternal_loss2530 0.998 0.9773 1.018
maternal_loss3035 1.003 0.9852 1.021
maternal_loss3540 1.002 0.9855 1.018
maternal_loss4045 1.009 0.9924 1.026
maternal_lossunclear 0.9883 0.9661 1.012
older_siblings1 0.9714 0.9566 0.986
older_siblings2 0.9686 0.9525 0.9858
older_siblings3 0.9633 0.945 0.9826
older_siblings4 0.9748 0.954 0.9967
older_siblings5P 0.9592 0.9334 0.9843
nr.siblings 1.007 1.004 1.009
last_born1 1.007 0.9927 1.022
hu_Intercept 0.1306 0.1032 0.1655
hu_paternalage 1.039 0.9247 1.162
hu_age50112 0.4802 0.45 0.5134
hu_age025 126.8 114.6 139.6
hu_ageunknown 8.31 7.746 8.914
hu_birth_cohort1675M1680 0.9938 0.8303 1.194
hu_birth_cohort1680M1685 1.068 0.8852 1.287
hu_birth_cohort1685M1690 1.158 0.951 1.395
hu_birth_cohort1690M1695 1.079 0.9028 1.287
hu_birth_cohort1695M1700 0.9953 0.8349 1.187
hu_birth_cohort1700M1705 1.033 0.8712 1.225
hu_birth_cohort1705M1710 0.986 0.8328 1.167
hu_birth_cohort1710M1715 1.079 0.9119 1.272
hu_birth_cohort1715M1720 1.164 0.9885 1.377
hu_birth_cohort1720M1725 1.131 0.9669 1.327
hu_birth_cohort1725M1730 1.308 1.111 1.536
hu_birth_cohort1730M1735 1.39 1.185 1.62
hu_birth_cohort1735M1740 1.307 1.117 1.524
hu_male1 1.648 1.569 1.729
hu_maternalage.factor1420 0.9383 0.8393 1.042
hu_maternalage.factor3550 0.9962 0.9178 1.085
hu_paternalage.mean 1.009 0.8956 1.138
hu_paternal_loss01 1.226 0.9986 1.502
hu_paternal_loss15 1.146 0.9834 1.347
hu_paternal_loss510 1.068 0.9252 1.22
hu_paternal_loss1015 1.012 0.8886 1.15
hu_paternal_loss1520 1.119 0.997 1.258
hu_paternal_loss2025 1.049 0.9394 1.174
hu_paternal_loss2530 1.121 1.003 1.249
hu_paternal_loss3035 1.124 1.014 1.246
hu_paternal_loss3540 1.095 0.9847 1.219
hu_paternal_loss4045 1.085 0.966 1.221
hu_paternal_lossunclear 1.276 1.145 1.424
hu_maternal_loss01 1.367 1.11 1.685
hu_maternal_loss15 1.005 0.8773 1.145
hu_maternal_loss510 1.13 1.004 1.278
hu_maternal_loss1015 1.121 0.9959 1.262
hu_maternal_loss1520 0.9837 0.871 1.111
hu_maternal_loss2025 1.057 0.9437 1.179
hu_maternal_loss2530 0.9408 0.8439 1.045
hu_maternal_loss3035 1.003 0.9134 1.106
hu_maternal_loss3540 1.057 0.9669 1.165
hu_maternal_loss4045 0.9654 0.876 1.062
hu_maternal_lossunclear 1.105 1 1.217
hu_older_siblings1 1.101 1.005 1.209
hu_older_siblings2 1.074 0.9736 1.186
hu_older_siblings3 1.069 0.9547 1.198
hu_older_siblings4 1.067 0.9447 1.209
hu_older_siblings5P 1.043 0.8969 1.218
hu_nr.siblings 0.9979 0.9867 1.009
hu_last_born1 0.9747 0.8938 1.064

Paternal age effect

pander::pander(paternal_age_10y_effect(model))
effect median_estimate ci_95 ci_80
estimate father 25y 5.9 [5.66;6.13] [5.75;6.05]
estimate father 35y 5.92 [5.66;6.17] [5.74;6.08]
percentage change 0.22 [-2.26;2.74] [-1.43;1.93]
OR/IRR 1.01 [0.99;1.03] [0.99;1.02]
OR hurdle 1.04 [0.92;1.16] [0.97;1.12]

Marginal effect plots

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

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

Coefficient plot

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

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

Diagnostics

These plots were made to diagnose misfit and nonconvergence.

Posterior predictive checks

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

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

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

Rhat

Did the 6 chains converge?

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

Effective sample size over average sample size

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

Trace plots

Trace plots are only shown in the case of nonconvergence.

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

File/cluster script name

This model was stored in the file: coefs/rpqa/s1_control_age.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

s2: Mediation via reproductive timing

Here, we tested whether the paternal age effect on reproductive succes is mediated by reproductive timing (as indexed by anchors’ ages at first and last birth). Because age at first and last birth are by definition only available for anchors who had at least one child, this analysis has to be restricted to such anchors. Hence, paternal age effects on mortality until age 1 and 15 cannot, in principle, be mediated by reproductive timing of the anchors.

Model summary

Full summary

model_summary = summary(model, use_cache = FALSE, priors = TRUE)
## Warning: There were 4 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 ~ paternalage + age_at_1st_child + age_at_last_child + birth_cohort + male + maternalage.factor + paternalage.mean + paternal_loss + maternal_loss + older_siblings + nr.siblings + last_born + (1 | idParents) 
##          hu ~ paternalage + age_at_1st_child + age_at_last_child + 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: 32562) 
## Samples: 6 chains, each with iter = 800; warmup = 300; thin = 1; 
##          total post-warmup samples = 3000
##    WAIC: Not computed
##  
## Priors: 
## b ~ normal(0,5)
## sd ~ student_t(3, 0, 5)
## b_hu ~ normal(0,5)
## sd_hu ~ student_t(3, 0, 10)
## 
## Group-Level Effects: 
## ~idParents (Number of levels: 10222) 
##                  Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
## sd(Intercept)        0.04      0.01     0.03     0.05        421 1.03
## sd(hu_Intercept)     4.92      3.94     0.23    14.84       2530 1.00
## 
## Population-Level Effects: 
##                           Estimate Est.Error l-95% CI u-95% CI Eff.Sample
## Intercept                     1.28      0.02     1.24     1.32       1139
## paternalage                   0.02      0.01     0.00     0.04       1933
## age_at_1st_child             -0.60      0.00    -0.61    -0.59       3000
## age_at_last_child             0.56      0.00     0.55     0.56       2474
## birth_cohort1675M1680         0.02      0.01     0.00     0.05        587
## birth_cohort1680M1685         0.06      0.01     0.03     0.09        564
## birth_cohort1685M1690         0.08      0.01     0.05     0.11        561
## birth_cohort1690M1695         0.07      0.01     0.04     0.10        537
## birth_cohort1695M1700         0.10      0.01     0.07     0.12        512
## birth_cohort1700M1705         0.08      0.01     0.05     0.10        497
## birth_cohort1705M1710         0.08      0.01     0.05     0.10        465
## birth_cohort1710M1715         0.07      0.01     0.05     0.10        471
## birth_cohort1715M1720         0.08      0.01     0.05     0.10        472
## birth_cohort1720M1725         0.06      0.01     0.04     0.09        442
## birth_cohort1725M1730         0.04      0.01     0.01     0.06        445
## birth_cohort1730M1735         0.06      0.01     0.03     0.08        440
## birth_cohort1735M1740         0.07      0.01     0.04     0.09        434
## male1                        -0.03      0.00    -0.04    -0.02       3000
## maternalage.factor1420       -0.01      0.01    -0.03     0.00       3000
## maternalage.factor3550       -0.01      0.01    -0.02     0.01       3000
## paternalage.mean             -0.03      0.01    -0.05    -0.01       2054
## paternal_loss01               0.01      0.02    -0.03     0.04       3000
## paternal_loss15              -0.01      0.01    -0.03     0.01       3000
## paternal_loss510             -0.02      0.01    -0.04     0.00       2579
## paternal_loss1015            -0.01      0.01    -0.03     0.01       2416
## paternal_loss1520            -0.02      0.01    -0.03     0.00       2368
## paternal_loss2025            -0.02      0.01    -0.04     0.00       2258
## paternal_loss2530             0.00      0.01    -0.02     0.01       2187
## paternal_loss3035             0.00      0.01    -0.02     0.01       2337
## paternal_loss3540            -0.01      0.01    -0.02     0.01       2300
## paternal_loss4045             0.00      0.01    -0.01     0.02       3000
## paternal_lossunclear         -0.03      0.01    -0.05    -0.02       2084
## maternal_loss01               0.03      0.02    -0.01     0.07       3000
## maternal_loss15               0.02      0.01    -0.01     0.04       3000
## maternal_loss510              0.00      0.01    -0.02     0.02       3000
## maternal_loss1015            -0.01      0.01    -0.02     0.01       3000
## maternal_loss1520             0.00      0.01    -0.02     0.02       3000
## maternal_loss2025             0.00      0.01    -0.02     0.02       3000
## maternal_loss2530             0.00      0.01    -0.02     0.02       3000
## maternal_loss3035            -0.01      0.01    -0.02     0.01       3000
## maternal_loss3540             0.00      0.01    -0.02     0.01       3000
## maternal_loss4045             0.01      0.01    -0.01     0.02       3000
## maternal_lossunclear         -0.02      0.01    -0.04    -0.01       3000
## older_siblings1              -0.02      0.01    -0.03    -0.01       3000
## older_siblings2              -0.02      0.01    -0.03     0.00       2438
## older_siblings3              -0.02      0.01    -0.04     0.00       2103
## older_siblings4              -0.03      0.01    -0.05    -0.01       1978
## older_siblings5P             -0.04      0.01    -0.06    -0.01       1819
## nr.siblings                   0.01      0.00     0.01     0.01       2288
## last_born1                    0.00      0.01    -0.01     0.02       3000
## hu_Intercept                 -1.51      4.88   -10.95     8.13       3000
## hu_paternalage               -4.27      4.44   -13.09     4.24       3000
## hu_age_at_1st_child          -3.34      4.58   -12.38     5.38       3000
## hu_age_at_last_child         -4.58      4.12   -13.28     2.90       3000
## hu_birth_cohort1675M1680      0.03      5.02    -9.95     9.90       3000
## hu_birth_cohort1680M1685     -0.21      4.92    -9.89     9.21       3000
## hu_birth_cohort1685M1690     -0.12      4.89    -9.77     9.56       3000
## hu_birth_cohort1690M1695     -0.13      4.95    -9.70     9.51       3000
## hu_birth_cohort1695M1700     -0.09      4.93    -9.77     9.56       3000
## hu_birth_cohort1700M1705     -0.11      5.00    -9.86     9.50       3000
## hu_birth_cohort1705M1710     -0.17      4.93    -9.61     9.75       3000
## hu_birth_cohort1710M1715     -0.21      5.14    -9.79     9.53       3000
## hu_birth_cohort1715M1720     -0.10      4.77    -9.37     9.11       3000
## hu_birth_cohort1720M1725     -0.14      4.97   -10.07     9.55       3000
## hu_birth_cohort1725M1730     -0.18      4.87    -9.79     9.13       3000
## hu_birth_cohort1730M1735     -0.28      4.80    -9.67     9.04       3000
## hu_birth_cohort1735M1740     -0.13      5.04    -9.94     9.82       3000
## hu_male1                     -0.45      4.82    -9.99     8.85       3000
## hu_maternalage.factor1420    -0.31      4.68    -9.54     8.86       3000
## hu_maternalage.factor3550    -0.20      4.95    -9.85     9.44       3000
## hu_paternalage.mean          -4.29      4.40   -13.23     4.43       3000
## hu_paternal_loss01           -0.28      4.86    -9.81     9.06       3000
## hu_paternal_loss15           -0.11      4.91    -9.77     9.19       3000
## hu_paternal_loss510          -0.06      4.88    -9.66     9.34       3000
## hu_paternal_loss1015         -0.02      4.91    -9.59     9.20       3000
## hu_paternal_loss1520         -0.14      4.83    -9.24     9.60       3000
## hu_paternal_loss2025         -0.03      5.06   -10.14     9.94       3000
## hu_paternal_loss2530          0.06      4.83    -9.32     9.45       3000
## hu_paternal_loss3035         -0.30      5.04   -10.12     9.48       3000
## hu_paternal_loss3540         -0.32      4.89    -9.99     9.08       3000
## hu_paternal_loss4045         -0.05      4.92    -9.92     9.55       3000
## hu_paternal_lossunclear      -0.38      4.86    -9.85     8.93       3000
## hu_maternal_loss01           -0.08      4.79    -9.76     9.06       3000
## hu_maternal_loss15           -0.24      5.04    -9.91     9.77       3000
## hu_maternal_loss510           0.02      5.12    -9.76    10.02       3000
## hu_maternal_loss1015         -0.05      4.92    -9.42     9.52       3000
## hu_maternal_loss1520         -0.07      4.90    -9.51     9.64       3000
## hu_maternal_loss2025          0.01      4.85    -9.61     9.13       3000
## hu_maternal_loss2530         -0.20      4.89    -9.94     9.55       3000
## hu_maternal_loss3035         -0.13      4.98    -9.84     9.45       3000
## hu_maternal_loss3540          0.07      5.06   -10.21    10.07       3000
## hu_maternal_loss4045         -0.06      4.92    -9.47     9.72       3000
## hu_maternal_lossunclear      -0.15      4.82    -9.79     9.03       3000
## hu_older_siblings1           -0.41      5.04   -10.26     9.22       3000
## hu_older_siblings2           -0.11      4.89    -9.84     9.44       3000
## hu_older_siblings3           -0.10      5.10   -10.05     9.89       3000
## hu_older_siblings4           -0.10      4.87    -9.89     9.21       3000
## hu_older_siblings5P           0.01      5.01   -10.09     9.84       3000
## hu_nr.siblings               -3.92      3.30   -11.21     1.23       2316
## hu_last_born1                -0.29      4.88    -9.91     9.24       3000
##                           Rhat
## Intercept                 1.00
## paternalage               1.00
## age_at_1st_child          1.00
## age_at_last_child         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.02
## 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.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.00
## hu_paternalage            1.00
## hu_age_at_1st_child       1.00
## hu_age_at_last_child      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

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.589 3.439 3.746
paternalage 1.022 1.004 1.041
age_at_1st_child 0.5506 0.5459 0.5554
age_at_last_child 1.743 1.733 1.752
birth_cohort1675M1680 1.024 0.9961 1.055
birth_cohort1680M1685 1.062 1.033 1.093
birth_cohort1685M1690 1.079 1.048 1.111
birth_cohort1690M1695 1.071 1.042 1.1
birth_cohort1695M1700 1.101 1.073 1.131
birth_cohort1700M1705 1.081 1.054 1.11
birth_cohort1705M1710 1.079 1.053 1.107
birth_cohort1710M1715 1.077 1.051 1.104
birth_cohort1715M1720 1.079 1.052 1.105
birth_cohort1720M1725 1.066 1.041 1.093
birth_cohort1725M1730 1.038 1.013 1.063
birth_cohort1730M1735 1.058 1.032 1.084
birth_cohort1735M1740 1.071 1.045 1.098
male1 0.974 0.9651 0.9825
maternalage.factor1420 0.9878 0.9724 1.005
maternalage.factor3550 0.994 0.9805 1.008
paternalage.mean 0.9714 0.9537 0.9903
paternal_loss01 1.005 0.9725 1.039
paternal_loss15 0.9914 0.9686 1.015
paternal_loss510 0.9816 0.9624 1.001
paternal_loss1015 0.993 0.9752 1.012
paternal_loss1520 0.9835 0.966 1
paternal_loss2025 0.9817 0.9655 0.9979
paternal_loss2530 0.9974 0.9809 1.014
paternal_loss3035 0.9961 0.9814 1.012
paternal_loss3540 0.9921 0.9771 1.007
paternal_loss4045 1.003 0.9865 1.021
paternal_lossunclear 0.9657 0.9504 0.9805
maternal_loss01 1.03 0.9906 1.07
maternal_loss15 1.016 0.9934 1.039
maternal_loss510 0.9987 0.9795 1.017
maternal_loss1015 0.9938 0.9757 1.012
maternal_loss1520 1.002 0.9842 1.02
maternal_loss2025 0.9974 0.9795 1.015
maternal_loss2530 1 0.9846 1.016
maternal_loss3035 0.9946 0.9799 1.009
maternal_loss3540 0.9967 0.9821 1.011
maternal_loss4045 1.007 0.9934 1.022
maternal_lossunclear 0.9783 0.9642 0.9929
older_siblings1 0.9805 0.9668 0.9942
older_siblings2 0.982 0.9668 0.9978
older_siblings3 0.9804 0.963 0.9982
older_siblings4 0.9734 0.9547 0.9932
older_siblings5P 0.9627 0.9403 0.9864
nr.siblings 1.007 1.005 1.009
last_born1 1.003 0.9892 1.017
hu_Intercept 0.2205 1.764e-05 3402
hu_paternalage 0.01399 2.066e-06 69.74
hu_age_at_1st_child 0.03538 4.189e-06 216
hu_age_at_last_child 0.01021 1.714e-06 18.16
hu_birth_cohort1675M1680 1.031 4.776e-05 19992
hu_birth_cohort1680M1685 0.8085 5.077e-05 9952
hu_birth_cohort1685M1690 0.8898 5.699e-05 14117
hu_birth_cohort1690M1695 0.8816 6.133e-05 13505
hu_birth_cohort1695M1700 0.9179 5.733e-05 14216
hu_birth_cohort1700M1705 0.8966 5.236e-05 13416
hu_birth_cohort1705M1710 0.8423 6.722e-05 17154
hu_birth_cohort1710M1715 0.8116 5.623e-05 13802
hu_birth_cohort1715M1720 0.9007 8.497e-05 9064
hu_birth_cohort1720M1725 0.8701 4.218e-05 14053
hu_birth_cohort1725M1730 0.837 5.584e-05 9233
hu_birth_cohort1730M1735 0.7573 6.315e-05 8441
hu_birth_cohort1735M1740 0.8739 4.798e-05 18375
hu_male1 0.6361 4.6e-05 6986
hu_maternalage.factor1420 0.7349 7.166e-05 7025
hu_maternalage.factor3550 0.8148 5.287e-05 12520
hu_paternalage.mean 0.01364 1.802e-06 84.28
hu_paternal_loss01 0.7569 5.464e-05 8593
hu_paternal_loss15 0.8979 5.713e-05 9829
hu_paternal_loss510 0.9432 6.396e-05 11349
hu_paternal_loss1015 0.9822 6.847e-05 9900
hu_paternal_loss1520 0.8722 9.669e-05 14693
hu_paternal_loss2025 0.9665 3.944e-05 20751
hu_paternal_loss2530 1.059 8.97e-05 12757
hu_paternal_loss3035 0.7389 4.036e-05 13076
hu_paternal_loss3540 0.7235 4.579e-05 8789
hu_paternal_loss4045 0.9556 4.937e-05 13978
hu_paternal_lossunclear 0.6823 5.267e-05 7587
hu_maternal_loss01 0.9234 5.798e-05 8647
hu_maternal_loss15 0.7891 4.959e-05 17436
hu_maternal_loss510 1.022 5.799e-05 22397
hu_maternal_loss1015 0.9509 8.103e-05 13600
hu_maternal_loss1520 0.9292 7.374e-05 15303
hu_maternal_loss2025 1.01 6.726e-05 9198
hu_maternal_loss2530 0.8228 4.837e-05 14048
hu_maternal_loss3035 0.879 5.324e-05 12732
hu_maternal_loss3540 1.073 3.691e-05 23669
hu_maternal_loss4045 0.9374 7.718e-05 16589
hu_maternal_lossunclear 0.8631 5.615e-05 8390
hu_older_siblings1 0.6668 3.5e-05 10115
hu_older_siblings2 0.8959 5.304e-05 12519
hu_older_siblings3 0.9077 4.317e-05 19704
hu_older_siblings4 0.9022 5.054e-05 9995
hu_older_siblings5P 1.007 4.165e-05 18702
hu_nr.siblings 0.01988 1.352e-05 3.405
hu_last_born1 0.7471 4.964e-05 10311

Paternal age effect

pander::pander(paternal_age_10y_effect(model))
effect median_estimate ci_95 ci_80
estimate father 25y 8.06 [7.85;8.27] [7.92;8.19]
estimate father 35y 8.23 [8;8.47] [8.09;8.39]
percentage change 2.24 [0.44;4.06] [1.01;3.42]
OR/IRR 1.02 [1;1.04] [1.01;1.03]
OR hurdle 0.02 [0;69.74] [0;4.2]

Marginal effect plots

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

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

Coefficient plot

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

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

Diagnostics

These plots were made to diagnose misfit and nonconvergence.

Posterior predictive checks

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

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

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

Rhat

Did the 6 chains converge?

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

Effective sample size over average sample size

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

Trace plots

Trace plots are only shown in the case of nonconvergence.

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

File/cluster script name

This model was stored in the file: coefs/rpqa/s2_reproductive_timing.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