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)
The rpqa.1
dataset contains only those participants where paternal age is known and the birthdate is between 1630 and 1750.
All of the following models are the same as our main model m3, except for the noted changes to test robustness.
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 = 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).
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 |
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] |
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)
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")
These plots were made to diagnose misfit and nonconvergence.
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")
Did the 6 chains converge?
stanplot(model, pars = "^b_[^I]", type = 'rhat')
stanplot(model, pars = "^b", type = 'neff_hist')
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)
}
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 = "")
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 = 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).
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 |
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] |
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)
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")
These plots were made to diagnose misfit and nonconvergence.
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")
Did the 6 chains converge?
stanplot(model, pars = "^b_[^I]", type = 'rhat')
stanplot(model, pars = "^b", type = 'neff_hist')
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)
}
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 = "")