Sweden (1947-2000), selective episodes
Loading details
source("0__helpers.R")
opts_chunk$set(warning=TRUE, cache = FALSE,tidy=FALSE,dev=c('png'),fig.width=20,fig.height=12.5)
make_path = function(file) {
get_coefficient_path(file, "swed")
}
# options for each chunk calling knit_child
opts_chunk$set(warning=FALSE, message = FALSE, echo = FALSE)Analysis description
Data subset
The swed_subset_survival.1 dataset contains only those participants where paternal age is known and birth years are from 1969 to 2000. We use it for the analyses pertaining to survival. We randomly sampled 100,000 families from the whole dataset.
The swed_subset_children.1 dataset contains only those participants where paternal age is known and birth years are from 1947 to 1959. We use it for the analyses pertaining to marriage, divorce and reproductive success. We randomly sampled 100,000 families from the whole dataset.
Model description
All of the following models have the following in common:
Estimation
We fit all models using brms v. 1.2.0, a Bayesian regression analysis statistical package. brms uses Stan, a probabilistic programming langugage to fit models using Hamiltonian Monte Carlo.
Covariates
As in our main models we adjust for average paternal age in the family, birth cohort (birth years in five equally large bins), for male sex, for age at paternal and maternal loss (0-1, 2-5, 6-10, …, 41-45, 45+, unknown), for maternal age (bins of 14-20, 20-35 and 35-50), for the number of siblings, for the number of older siblings (0-5, 5+) and for being last born.
Model stratification
We added random intercepts for each family (father-mother dyad). We then controlled for the average paternal age in the family. Hence, the paternal age effects in the plot are split into those between families and those within families or between siblings. We are interested in the effect of paternal age between siblings, as this effect cannot be explained by e.g. genetic propensities of the father to reproduce late.
e1: Survival to first year
Here, we predict the probability that the anchor survives the first year of life. All children born to this father are compared, if their death date is known or their survival can be inferred (from later marriage or children).
Model summary
Full summary
model_summary = summary(model, use_cache = FALSE, priors = TRUE)
print(model_summary)## Family: bernoulli(cauchit)
## Formula: survive1y ~ paternalage + birth_cohort + male + maternalage.factor + paternalage.mean + paternal_loss + maternal_loss + older_siblings + nr.siblings + last_born + (1 | idParents)
## Data: model_data (Number of observations: 363744)
## Samples: 6 chains, each with iter = 1000; warmup = 500; thin = 1;
## total post-warmup samples = 3000
## WAIC: Not computed
##
## Priors:
## b ~ normal(0,5)
## sd ~ student_t(3, 0, 5)
##
## Group-Level Effects:
## ~idParents (Number of levels: 200000)
## Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
## sd(Intercept) 21.41 0.8 19.76 22.96 205 1.03
##
## Population-Level Effects:
## Estimate Est.Error l-95% CI u-95% CI Eff.Sample
## Intercept 96.06 4.72 87.17 105.59 463
## paternalage -10.72 1.94 -14.50 -6.88 488
## birth_cohort1970M1975 -5.41 1.94 -9.34 -1.77 251
## birth_cohort1975M1980 -3.12 2.00 -7.06 0.75 317
## birth_cohort1980M1985 4.66 2.20 0.36 8.89 405
## birth_cohort1985M1990 3.90 2.14 -0.41 8.23 327
## birth_cohort1990M1995 5.97 2.20 1.72 10.42 361
## birth_cohort1995M2000 14.75 3.05 8.92 20.72 368
## male1 -4.69 1.00 -6.73 -2.73 337
## maternalage.factor1420 0.42 2.18 -3.58 4.93 578
## maternalage.factor3561 -4.91 1.75 -8.19 -1.47 753
## paternalage.mean 10.23 2.04 6.20 14.17 665
## paternal_loss01 0.33 4.99 -9.73 10.01 3000
## paternal_loss15 -5.30 4.50 -13.41 4.04 1098
## paternal_loss510 3.95 4.33 -4.09 12.77 1661
## paternal_loss1015 3.54 3.74 -3.41 11.17 1277
## paternal_loss1520 -3.51 2.94 -8.89 2.43 708
## paternal_loss2025 -2.70 2.99 -8.26 3.67 464
## paternal_loss2530 0.87 2.68 -4.04 6.38 907
## paternal_loss3035 2.02 2.93 -3.34 7.97 514
## paternal_loss3540 -4.42 2.66 -9.17 1.01 819
## paternal_loss4045 -1.32 4.98 -10.96 8.96 3000
## maternal_loss01 -5.21 5.62 -16.28 5.58 765
## maternal_loss15 1.90 4.72 -7.17 11.19 3000
## maternal_loss510 1.54 4.54 -6.87 10.56 3000
## maternal_loss1015 -0.35 4.18 -8.03 8.20 1203
## maternal_loss1520 -6.18 3.25 -12.29 0.52 941
## maternal_loss2025 -6.09 3.05 -11.88 0.33 951
## maternal_loss2530 -4.50 3.09 -10.11 1.93 978
## maternal_loss3035 -2.91 3.27 -9.07 3.99 792
## maternal_loss3540 -1.60 3.96 -8.63 6.56 848
## maternal_loss4045 0.97 4.70 -8.54 10.14 3000
## older_siblings1 -6.34 1.37 -9.06 -3.75 249
## older_siblings2 0.01 2.01 -3.87 3.98 147
## older_siblings3 14.49 2.50 9.72 19.44 326
## older_siblings4 16.49 2.71 11.00 21.83 602
## older_siblings5P 25.63 3.53 18.62 32.67 625
## nr.siblings -7.96 0.40 -8.76 -7.19 420
## last_born1 -12.79 1.02 -14.73 -10.83 245
## Rhat
## Intercept 1.01
## paternalage 1.01
## birth_cohort1970M1975 1.01
## birth_cohort1975M1980 1.01
## birth_cohort1980M1985 1.00
## birth_cohort1985M1990 1.01
## birth_cohort1990M1995 1.01
## birth_cohort1995M2000 1.01
## male1 1.02
## maternalage.factor1420 1.01
## maternalage.factor3561 1.01
## paternalage.mean 1.01
## paternal_loss01 1.00
## paternal_loss15 1.00
## paternal_loss510 1.00
## paternal_loss1015 1.00
## paternal_loss1520 1.01
## paternal_loss2025 1.01
## paternal_loss2530 1.01
## paternal_loss3035 1.01
## paternal_loss3540 1.01
## paternal_loss4045 1.00
## maternal_loss01 1.00
## maternal_loss15 1.00
## maternal_loss510 1.00
## maternal_loss1015 1.00
## maternal_loss1520 1.00
## maternal_loss2025 1.01
## maternal_loss2530 1.00
## maternal_loss3035 1.00
## maternal_loss3540 1.01
## maternal_loss4045 1.00
## older_siblings1 1.02
## older_siblings2 1.06
## older_siblings3 1.02
## older_siblings4 1.01
## older_siblings5P 1.01
## nr.siblings 1.01
## last_born1 1.03
##
## 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 | 5.236e+41 | 7.169e+37 | 7.219e+45 |
| paternalage | 2.215e-05 | 5.064e-07 | 0.001033 |
| birth_cohort1970M1975 | 0.004475 | 8.75e-05 | 0.1705 |
| birth_cohort1975M1980 | 0.04427 | 0.0008587 | 2.109 |
| birth_cohort1980M1985 | 105.5 | 1.433 | 7230 |
| birth_cohort1985M1990 | 49.61 | 0.6604 | 3745 |
| birth_cohort1990M1995 | 390.9 | 5.559 | 33363 |
| birth_cohort1995M2000 | 2556476 | 7504 | 996516331 |
| male1 | 0.009161 | 0.001189 | 0.06493 |
| maternalage.factor1420 | 1.523 | 0.0278 | 138.5 |
| maternalage.factor3561 | 0.007358 | 0.0002762 | 0.2295 |
| paternalage.mean | 27741 | 491.5 | 1420899 |
| paternal_loss01 | 1.386 | 5.937e-05 | 22327 |
| paternal_loss15 | 0.004974 | 1.502e-06 | 57.06 |
| paternal_loss510 | 51.97 | 0.01677 | 352889 |
| paternal_loss1015 | 34.55 | 0.03319 | 70655 |
| paternal_loss1520 | 0.02978 | 0.0001379 | 11.39 |
| paternal_loss2025 | 0.06721 | 0.0002594 | 39.27 |
| paternal_loss2530 | 2.394 | 0.01751 | 588.8 |
| paternal_loss3035 | 7.52 | 0.03542 | 2895 |
| paternal_loss3540 | 0.01198 | 0.0001041 | 2.735 |
| paternal_loss4045 | 0.2672 | 1.741e-05 | 7818 |
| maternal_loss01 | 0.005475 | 8.498e-08 | 263.8 |
| maternal_loss15 | 6.672 | 0.0007666 | 72578 |
| maternal_loss510 | 4.676 | 0.001035 | 38393 |
| maternal_loss1015 | 0.7025 | 0.0003243 | 3641 |
| maternal_loss1520 | 0.002073 | 4.595e-06 | 1.688 |
| maternal_loss2025 | 0.002258 | 6.936e-06 | 1.387 |
| maternal_loss2530 | 0.01115 | 4.053e-05 | 6.861 |
| maternal_loss3035 | 0.05461 | 0.0001147 | 54.01 |
| maternal_loss3540 | 0.2023 | 0.0001782 | 708.5 |
| maternal_loss4045 | 2.634 | 0.0001961 | 25350 |
| older_siblings1 | 0.001758 | 0.0001159 | 0.02347 |
| older_siblings2 | 1.012 | 0.02092 | 53.47 |
| older_siblings3 | 1966735 | 16620 | 276670549 |
| older_siblings4 | 14519530 | 59968 | 3.018e+09 |
| older_siblings5P | 1.348e+11 | 121468561 | 1.54e+14 |
| nr.siblings | 0.0003476 | 0.0001564 | 0.0007578 |
| last_born1 | 2.8e-06 | 3.996e-07 | 1.978e-05 |
Paternal age effect
pander::pander(paternal_age_10y_effect(model))| effect | median_estimate | ci_95 | ci_80 |
|---|---|---|---|
| estimate father 25y | 1 | [1;1] | [1;1] |
| estimate father 35y | 1 | [1;1] | [1;1] |
| percentage change | -0.05 | [-0.06;-0.03] | [-0.06;-0.03] |
| OR/IRR | 0 | [0;0] | [0;0] |
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/swed/e1_survive1y.rds.
Click the following link to see the script used to generate this model:
opts_chunk$set(echo = FALSE)
clusterscript = str_replace(basename(model_filename), "\\.rds",".html")
cat("[Cluster script](" , clusterscript, ")", sep = "")
e2: Probability of surviving the first 15 years of life
Here, we predict the probability that the anchor survives the first fifteen of life. All children born to this father who lived at least one year are compared, if their death date is known or their survival can be inferred (from later marriage or children).
Model summary
Full summary
model_summary = summary(model, use_cache = FALSE, priors = TRUE)
print(model_summary)## Family: bernoulli(cauchit)
## Formula: surviveR ~ paternalage + birth_cohort + male + maternalage.factor + paternalage.mean + paternal_loss + maternal_loss + older_siblings + nr.siblings + last_born + (1 | idParents)
## Data: model_data (Number of observations: 242223)
## Samples: 6 chains, each with iter = 1500; warmup = 1000; thin = 1;
## total post-warmup samples = 3000
## WAIC: Not computed
##
## Priors:
## b ~ normal(0,70)
## sd ~ student_t(3, 0, 50)
##
## Group-Level Effects:
## ~idParents (Number of levels: 142552)
## Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
## sd(Intercept) 47.62 3.59 41.08 55.12 79 1.05
##
## Population-Level Effects:
## Estimate Est.Error l-95% CI u-95% CI Eff.Sample
## Intercept 191.88 16.82 160.72 227.46 262
## paternalage 26.88 12.18 3.10 51.20 157
## birth_cohort1970M1975 1.05 9.26 -20.20 16.31 114
## birth_cohort1975M1980 18.78 10.56 -3.31 37.26 111
## birth_cohort1980M1985 13.23 10.28 -10.82 30.57 98
## birth_cohort1985M1990 22.18 10.76 -1.21 40.90 136
## male1 -17.14 4.30 -25.90 -8.99 269
## maternalage.factor1420 18.05 8.51 2.10 36.77 199
## maternalage.factor3561 5.65 8.66 -10.18 23.63 289
## paternalage.mean -32.48 12.09 -56.97 -9.09 163
## paternal_loss01 43.19 50.58 -37.96 155.56 3000
## paternal_loss15 -43.37 11.16 -63.46 -19.46 495
## paternal_loss510 -39.25 8.92 -56.47 -21.70 534
## paternal_loss1015 -15.51 13.94 -37.54 17.04 386
## paternal_loss1520 -18.70 12.77 -37.27 15.44 300
## paternal_loss2025 -12.34 11.89 -32.10 15.07 251
## paternal_loss2530 -0.11 15.76 -23.49 37.71 185
## paternal_loss3035 -12.04 11.28 -31.26 13.57 266
## paternal_loss3540 4.80 21.31 -23.80 60.13 250
## paternal_loss4045 42.35 48.40 -37.55 145.72 2223
## maternal_loss01 24.01 54.56 -65.32 144.66 3000
## maternal_loss15 -68.39 11.68 -89.65 -43.70 597
## maternal_loss510 -58.70 9.42 -76.72 -39.30 709
## maternal_loss1015 -43.64 13.69 -65.81 -11.26 166
## maternal_loss1520 10.31 34.31 -41.86 89.87 257
## maternal_loss2025 -9.40 21.18 -39.47 42.96 168
## maternal_loss2530 -24.50 14.92 -45.60 14.77 184
## maternal_loss3035 -17.48 13.60 -39.44 13.27 349
## maternal_loss3540 -14.87 16.56 -40.06 22.70 195
## maternal_loss4045 -26.47 45.12 -84.22 94.82 348
## older_siblings1 9.27 5.53 -1.35 20.14 202
## older_siblings2 8.27 9.53 -10.93 26.37 126
## older_siblings3 10.09 16.53 -19.72 44.73 168
## older_siblings4 13.43 27.50 -28.11 80.78 302
## older_siblings5P 22.35 27.97 -26.56 80.14 200
## nr.siblings -6.86 1.95 -10.63 -3.06 183
## last_born1 17.85 5.26 7.75 28.14 147
## Rhat
## Intercept 1.02
## paternalage 1.04
## birth_cohort1970M1975 1.04
## birth_cohort1975M1980 1.05
## birth_cohort1980M1985 1.05
## birth_cohort1985M1990 1.03
## male1 1.01
## maternalage.factor1420 1.03
## maternalage.factor3561 1.02
## paternalage.mean 1.04
## paternal_loss01 1.00
## paternal_loss15 1.00
## paternal_loss510 1.01
## paternal_loss1015 1.01
## paternal_loss1520 1.02
## paternal_loss2025 1.02
## paternal_loss2530 1.03
## paternal_loss3035 1.01
## paternal_loss3540 1.02
## paternal_loss4045 1.00
## maternal_loss01 1.00
## maternal_loss15 1.01
## maternal_loss510 1.01
## maternal_loss1015 1.04
## maternal_loss1520 1.01
## maternal_loss2025 1.04
## maternal_loss2530 1.04
## maternal_loss3035 1.01
## maternal_loss3540 1.02
## maternal_loss4045 1.01
## older_siblings1 1.03
## older_siblings2 1.04
## older_siblings3 1.04
## older_siblings4 1.03
## older_siblings5P 1.04
## nr.siblings 1.03
## last_born1 1.03
##
## 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 | 2.157e+83 | 6.307e+69 | 6.086e+98 |
| paternalage | 4.718e+11 | 22.29 | 1.729e+22 |
| birth_cohort1970M1975 | 2.864 | 1.691e-09 | 12151048 |
| birth_cohort1975M1980 | 143513462 | 0.03641 | 1.515e+16 |
| birth_cohort1980M1985 | 554314 | 2.003e-05 | 1.891e+13 |
| birth_cohort1985M1990 | 4.296e+09 | 0.2977 | 5.776e+17 |
| male1 | 3.602e-08 | 5.65e-12 | 0.000125 |
| maternalage.factor1420 | 68913143 | 8.176 | 9.34e+15 |
| maternalage.factor3561 | 285.1 | 3.803e-05 | 1.827e+10 |
| paternalage.mean | 7.819e-15 | 1.821e-25 | 0.0001126 |
| paternal_loss01 | 5.696e+18 | 3.265e-17 | 3.639e+67 |
| paternal_loss15 | 1.454e-19 | 2.755e-28 | 3.544e-09 |
| paternal_loss510 | 8.995e-18 | 2.983e-25 | 3.748e-10 |
| paternal_loss1015 | 1.838e-07 | 4.953e-17 | 25182535 |
| paternal_loss1520 | 7.525e-09 | 6.49e-17 | 5076020 |
| paternal_loss2025 | 4.37e-06 | 1.142e-14 | 3496308 |
| paternal_loss2530 | 0.8916 | 6.314e-11 | 2.394e+16 |
| paternal_loss3035 | 5.882e-06 | 2.647e-14 | 781850 |
| paternal_loss3540 | 121 | 4.592e-11 | 1.295e+26 |
| paternal_loss4045 | 2.462e+18 | 4.915e-17 | 1.928e+63 |
| maternal_loss01 | 2.671e+10 | 4.271e-29 | 6.696e+62 |
| maternal_loss15 | 1.991e-30 | 1.165e-39 | 1.047e-19 |
| maternal_loss510 | 3.203e-26 | 4.806e-34 | 8.595e-18 |
| maternal_loss1015 | 1.112e-19 | 2.633e-29 | 1.284e-05 |
| maternal_loss1520 | 30102 | 6.594e-19 | 1.071e+39 |
| maternal_loss2025 | 8.267e-05 | 7.223e-18 | 4.534e+18 |
| maternal_loss2530 | 2.292e-11 | 1.565e-20 | 2588648 |
| maternal_loss3035 | 2.56e-08 | 7.43e-18 | 579005 |
| maternal_loss3540 | 3.487e-07 | 3.99e-18 | 7.187e+09 |
| maternal_loss4045 | 3.204e-12 | 2.663e-37 | 1.52e+41 |
| older_siblings1 | 10624 | 0.259 | 558935740 |
| older_siblings2 | 3894 | 1.79e-05 | 2.845e+11 |
| older_siblings3 | 24211 | 2.74e-09 | 2.678e+19 |
| older_siblings4 | 677823 | 6.196e-13 | 1.204e+35 |
| older_siblings5P | 5.104e+09 | 2.92e-12 | 6.377e+34 |
| nr.siblings | 0.001046 | 2.427e-05 | 0.04667 |
| last_born1 | 56262650 | 2323 | 1.66e+12 |
Paternal age effect
pander::pander(paternal_age_10y_effect(model))| effect | median_estimate | ci_95 | ci_80 |
|---|---|---|---|
| estimate father 25y | 1 | [1;1] | [1;1] |
| estimate father 35y | 1 | [1;1] | [1;1] |
| percentage change | 0.03 | [0;0.06] | [0.01;0.05] |
| OR/IRR | 440625466756.42 | [22.29;1.72900464486011e+22] | [59649.6;2832585833496194560] |
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/swed/e2_surviveR.rds.
Click the following link to see the script used to generate this model:
opts_chunk$set(echo = FALSE)
clusterscript = str_replace(basename(model_filename), "\\.rds",".html")
cat("[Cluster script](" , clusterscript, ")", sep = "")
e3: Probability of ever marrying
Here, we predict the probability that the anchor ever marries. All anchors who reached reproductive age (15) are included.
Model summary
Full summary
model_summary = summary(model, use_cache = FALSE, priors = TRUE)
print(model_summary)## Family: bernoulli(cauchit)
## Formula: ever_married ~ paternalage + birth_cohort + male + maternalage.factor + paternalage.mean + paternal_loss + maternal_loss + older_siblings + nr.siblings + last_born + (1 | idParents)
## Data: model_data (Number of observations: 1392421)
## Samples: 6 chains, each with iter = 800; warmup = 300; thin = 1;
## total post-warmup samples = 3000
## WAIC: Not computed
##
## Priors:
## b ~ normal(0,5)
## sd ~ student_t(3, 0, 5)
##
## Group-Level Effects:
## ~idParents (Number of levels: 874603)
## Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
## sd(Intercept) 0.93 0.01 0.91 0.94 572 1.01
##
## Population-Level Effects:
## Estimate Est.Error l-95% CI u-95% CI Eff.Sample
## Intercept 2.35 0.02 2.31 2.39 3000
## paternalage 0.14 0.02 0.11 0.18 940
## birth_cohort1950M1955 -0.46 0.01 -0.48 -0.45 3000
## birth_cohort1955M1960 -0.89 0.01 -0.90 -0.87 3000
## male1 -0.47 0.01 -0.48 -0.46 3000
## maternalage.factor1420 -0.05 0.01 -0.07 -0.03 3000
## maternalage.factor3561 -0.01 0.01 -0.03 0.01 3000
## paternalage.mean -0.22 0.02 -0.26 -0.19 991
## paternal_loss01 -0.18 0.22 -0.59 0.27 3000
## paternal_loss15 -0.15 0.07 -0.28 -0.01 3000
## paternal_loss510 -0.19 0.03 -0.25 -0.12 3000
## paternal_loss1015 -0.20 0.02 -0.24 -0.15 3000
## paternal_loss1520 -0.17 0.02 -0.20 -0.14 3000
## paternal_loss2025 -0.13 0.01 -0.15 -0.10 3000
## paternal_loss2530 -0.10 0.01 -0.12 -0.08 3000
## paternal_loss3035 -0.10 0.01 -0.12 -0.08 3000
## paternal_loss3540 -0.10 0.01 -0.11 -0.08 3000
## paternal_loss4045 -0.08 0.01 -0.10 -0.06 3000
## paternal_lossunclear 0.00 0.01 -0.02 0.01 3000
## maternal_loss01 -0.51 0.22 -0.93 -0.05 3000
## maternal_loss15 -0.32 0.10 -0.50 -0.12 3000
## maternal_loss510 -0.22 0.05 -0.31 -0.12 3000
## maternal_loss1015 -0.20 0.03 -0.26 -0.14 3000
## maternal_loss1520 -0.18 0.02 -0.23 -0.14 3000
## maternal_loss2025 -0.07 0.02 -0.11 -0.03 3000
## maternal_loss2530 -0.09 0.02 -0.13 -0.06 3000
## maternal_loss3035 -0.07 0.01 -0.10 -0.04 3000
## maternal_loss3540 -0.08 0.01 -0.10 -0.06 3000
## maternal_loss4045 -0.06 0.01 -0.08 -0.03 3000
## maternal_lossunclear 0.04 0.01 0.03 0.06 3000
## older_siblings1 -0.01 0.01 -0.03 0.00 1268
## older_siblings2 -0.03 0.01 -0.06 0.00 949
## older_siblings3 -0.07 0.02 -0.11 -0.03 1082
## older_siblings4 -0.15 0.03 -0.21 -0.10 1066
## older_siblings5P -0.24 0.04 -0.31 -0.17 880
## nr.siblings 0.02 0.00 0.01 0.03 999
## last_born1 0.05 0.01 0.04 0.06 3000
## Rhat
## Intercept 1.00
## paternalage 1.01
## birth_cohort1950M1955 1.00
## birth_cohort1955M1960 1.00
## male1 1.00
## maternalage.factor1420 1.00
## maternalage.factor3561 1.00
## paternalage.mean 1.01
## paternal_loss01 1.00
## paternal_loss15 1.00
## paternal_loss510 1.00
## paternal_loss1015 1.00
## paternal_loss1520 1.00
## paternal_loss2025 1.00
## paternal_loss2530 1.00
## paternal_loss3035 1.00
## paternal_loss3540 1.00
## paternal_loss4045 1.00
## paternal_lossunclear 1.00
## maternal_loss01 1.00
## maternal_loss15 1.00
## maternal_loss510 1.00
## maternal_loss1015 1.00
## maternal_loss1520 1.00
## maternal_loss2025 1.00
## maternal_loss2530 1.00
## maternal_loss3035 1.00
## maternal_loss3540 1.00
## maternal_loss4045 1.00
## maternal_lossunclear 1.00
## older_siblings1 1.00
## older_siblings2 1.01
## older_siblings3 1.01
## older_siblings4 1.01
## older_siblings5P 1.01
## nr.siblings 1.01
## 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 | 10.52 | 10.1 | 10.95 |
| paternalage | 1.152 | 1.115 | 1.193 |
| birth_cohort1950M1955 | 0.6305 | 0.6209 | 0.6401 |
| birth_cohort1955M1960 | 0.4118 | 0.4055 | 0.4186 |
| male1 | 0.6232 | 0.6165 | 0.6297 |
| maternalage.factor1420 | 0.9553 | 0.9363 | 0.9753 |
| maternalage.factor3561 | 0.9896 | 0.9717 | 1.008 |
| paternalage.mean | 0.7989 | 0.7718 | 0.826 |
| paternal_loss01 | 0.8365 | 0.5541 | 1.305 |
| paternal_loss15 | 0.8641 | 0.7526 | 0.989 |
| paternal_loss510 | 0.8309 | 0.7792 | 0.8865 |
| paternal_loss1015 | 0.8224 | 0.7897 | 0.8571 |
| paternal_loss1520 | 0.8452 | 0.8185 | 0.8733 |
| paternal_loss2025 | 0.8817 | 0.8592 | 0.9056 |
| paternal_loss2530 | 0.9054 | 0.8837 | 0.9265 |
| paternal_loss3035 | 0.9025 | 0.8836 | 0.9218 |
| paternal_loss3540 | 0.9089 | 0.8921 | 0.9259 |
| paternal_loss4045 | 0.9245 | 0.9059 | 0.9431 |
| paternal_lossunclear | 0.998 | 0.9836 | 1.013 |
| maternal_loss01 | 0.5983 | 0.3931 | 0.9509 |
| maternal_loss15 | 0.7294 | 0.606 | 0.8865 |
| maternal_loss510 | 0.8035 | 0.7305 | 0.8848 |
| maternal_loss1015 | 0.8208 | 0.7698 | 0.8731 |
| maternal_loss1520 | 0.832 | 0.7949 | 0.8729 |
| maternal_loss2025 | 0.9349 | 0.8985 | 0.9722 |
| maternal_loss2530 | 0.9097 | 0.8802 | 0.9424 |
| maternal_loss3035 | 0.9296 | 0.9025 | 0.957 |
| maternal_loss3540 | 0.9218 | 0.9006 | 0.9435 |
| maternal_loss4045 | 0.9461 | 0.9249 | 0.9676 |
| maternal_lossunclear | 1.046 | 1.032 | 1.06 |
| older_siblings1 | 0.9864 | 0.9699 | 1.003 |
| older_siblings2 | 0.968 | 0.9399 | 0.9954 |
| older_siblings3 | 0.9319 | 0.8938 | 0.9716 |
| older_siblings4 | 0.8588 | 0.8097 | 0.9074 |
| older_siblings5P | 0.7865 | 0.7307 | 0.8459 |
| nr.siblings | 1.018 | 1.011 | 1.026 |
| last_born1 | 1.05 | 1.038 | 1.062 |
Paternal age effect
pander::pander(paternal_age_10y_effect(model))| effect | median_estimate | ci_95 | ci_80 |
|---|---|---|---|
| estimate father 25y | 0.86 | [0.85;0.86] | [0.85;0.86] |
| estimate father 35y | 0.86 | [0.86;0.87] | [0.86;0.86] |
| percentage change | 0.82 | [0.63;1.01] | [0.69;0.95] |
| OR/IRR | 1.15 | [1.11;1.19] | [1.13;1.18] |
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/swed/e3_ever_married.rds.
Click the following link to see the script used to generate this model:
opts_chunk$set(echo = FALSE)
clusterscript = str_replace(basename(model_filename), "\\.rds",".html")
cat("[Cluster script](" , clusterscript, ")", sep = "")
e4: Number of children
Here, we predict the number of children that the anchor had. To separate this effect from previous selective episodes, we include only ever-married anchors and control for their number of spouses (interacted with sex, because men tend to have more additional children from further spouses).
Model summary
Full summary
model_summary = summary(model, use_cache = FALSE, priors = TRUE)
print(model_summary)## Family: poisson(log)
## Formula: children ~ paternalage + birth_cohort + spouses * male + maternalage.factor + paternalage.mean + paternal_loss + maternal_loss + older_siblings + nr.siblings + last_born + (1 | idParents)
## Data: model_data (Number of observations: 1032368)
## Samples: 6 chains, each with iter = 1500; warmup = 500; thin = 1;
## total post-warmup samples = 6000
## WAIC: Not computed
##
## Priors:
## b ~ normal(0,5)
## sd ~ student_t(3, 0, 5)
##
## Group-Level Effects:
## ~idParents (Number of levels: 719113)
## Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
## sd(Intercept) 0 0 0 0 2848 1
##
## Population-Level Effects:
## Estimate Est.Error l-95% CI u-95% CI Eff.Sample
## Intercept 0.57 0.01 0.56 0.58 6000
## paternalage -0.04 0.00 -0.05 -0.03 2828
## birth_cohort1950M1955 0.06 0.00 0.05 0.06 6000
## birth_cohort1955M1960 0.08 0.00 0.08 0.08 6000
## spouses 0.10 0.00 0.09 0.10 6000
## male1 -0.04 0.00 -0.05 -0.03 6000
## maternalage.factor1420 0.03 0.00 0.02 0.03 6000
## maternalage.factor3561 0.01 0.00 0.00 0.01 6000
## paternalage.mean 0.03 0.00 0.02 0.04 2939
## paternal_loss01 0.10 0.06 -0.02 0.22 6000
## paternal_loss15 0.02 0.02 -0.02 0.06 6000
## paternal_loss510 0.00 0.01 -0.02 0.01 6000
## paternal_loss1015 0.01 0.01 -0.01 0.02 6000
## paternal_loss1520 0.00 0.00 0.00 0.01 6000
## paternal_loss2025 0.00 0.00 -0.01 0.01 6000
## paternal_loss2530 0.00 0.00 -0.01 0.00 6000
## paternal_loss3035 0.00 0.00 -0.01 0.00 6000
## paternal_loss3540 0.00 0.00 -0.01 0.00 6000
## paternal_loss4045 0.00 0.00 -0.01 0.00 6000
## paternal_lossunclear -0.01 0.00 -0.02 -0.01 6000
## maternal_loss01 -0.07 0.08 -0.23 0.09 6000
## maternal_loss15 -0.01 0.03 -0.07 0.05 6000
## maternal_loss510 0.02 0.01 -0.01 0.05 6000
## maternal_loss1015 0.00 0.01 -0.02 0.02 6000
## maternal_loss1520 0.00 0.01 -0.01 0.02 6000
## maternal_loss2025 0.01 0.01 0.00 0.02 6000
## maternal_loss2530 0.01 0.00 0.00 0.01 6000
## maternal_loss3035 0.00 0.00 0.00 0.01 6000
## maternal_loss3540 0.00 0.00 -0.01 0.01 6000
## maternal_loss4045 0.00 0.00 -0.01 0.00 6000
## maternal_lossunclear 0.00 0.00 0.00 0.00 6000
## older_siblings1 0.01 0.00 0.00 0.01 3530
## older_siblings2 0.01 0.00 0.00 0.02 2955
## older_siblings3 0.01 0.01 0.00 0.02 2952
## older_siblings4 0.00 0.01 -0.01 0.01 3123
## older_siblings5P -0.04 0.01 -0.06 -0.02 2899
## nr.siblings 0.03 0.00 0.03 0.03 2892
## last_born1 0.00 0.00 0.00 0.00 6000
## spouses:male1 0.05 0.00 0.04 0.05 6000
## Rhat
## Intercept 1
## paternalage 1
## birth_cohort1950M1955 1
## birth_cohort1955M1960 1
## spouses 1
## male1 1
## maternalage.factor1420 1
## maternalage.factor3561 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
## spouses:male1 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 | 1.77 | 1.751 | 1.79 |
| paternalage | 0.9617 | 0.9536 | 0.9702 |
| birth_cohort1950M1955 | 1.058 | 1.054 | 1.061 |
| birth_cohort1955M1960 | 1.084 | 1.08 | 1.088 |
| spouses | 1.1 | 1.096 | 1.104 |
| male1 | 0.9617 | 0.955 | 0.9685 |
| maternalage.factor1420 | 1.03 | 1.024 | 1.035 |
| maternalage.factor3561 | 1.007 | 1.002 | 1.012 |
| paternalage.mean | 1.029 | 1.02 | 1.038 |
| paternal_loss01 | 1.105 | 0.9832 | 1.242 |
| paternal_loss15 | 1.023 | 0.9818 | 1.065 |
| paternal_loss510 | 0.9961 | 0.9773 | 1.015 |
| paternal_loss1015 | 1.006 | 0.9947 | 1.018 |
| paternal_loss1520 | 1.004 | 0.9956 | 1.013 |
| paternal_loss2025 | 0.9983 | 0.9911 | 1.005 |
| paternal_loss2530 | 0.9983 | 0.9925 | 1.004 |
| paternal_loss3035 | 0.9996 | 0.9942 | 1.005 |
| paternal_loss3540 | 0.9982 | 0.9935 | 1.003 |
| paternal_loss4045 | 0.9955 | 0.9905 | 1 |
| paternal_lossunclear | 0.9869 | 0.9833 | 0.9907 |
| maternal_loss01 | 0.9341 | 0.795 | 1.093 |
| maternal_loss15 | 0.9902 | 0.9281 | 1.054 |
| maternal_loss510 | 1.019 | 0.9898 | 1.049 |
| maternal_loss1015 | 0.9996 | 0.9827 | 1.017 |
| maternal_loss1520 | 1.003 | 0.9898 | 1.016 |
| maternal_loss2025 | 1.007 | 0.996 | 1.018 |
| maternal_loss2530 | 1.005 | 0.9964 | 1.015 |
| maternal_loss3035 | 1.005 | 0.9972 | 1.012 |
| maternal_loss3540 | 0.9991 | 0.9928 | 1.005 |
| maternal_loss4045 | 0.9988 | 0.9927 | 1.005 |
| maternal_lossunclear | 0.9985 | 0.9952 | 1.002 |
| older_siblings1 | 1.009 | 1.005 | 1.013 |
| older_siblings2 | 1.012 | 1.005 | 1.02 |
| older_siblings3 | 1.008 | 0.9976 | 1.02 |
| older_siblings4 | 1 | 0.986 | 1.015 |
| older_siblings5P | 0.9589 | 0.9416 | 0.9767 |
| nr.siblings | 1.031 | 1.029 | 1.033 |
| last_born1 | 1.001 | 0.9981 | 1.004 |
| spouses:male1 | 1.048 | 1.043 | 1.054 |
Paternal age effect
pander::pander(paternal_age_10y_effect(model))| effect | median_estimate | ci_95 | ci_80 |
|---|---|---|---|
| estimate father 25y | 2.05 | [2.04;2.06] | [2.05;2.06] |
| estimate father 35y | 1.97 | [1.96;1.99] | [1.96;1.98] |
| percentage change | -3.82 | [-4.64;-2.98] | [-4.38;-3.28] |
| OR/IRR | 0.96 | [0.95;0.97] | [0.96;0.97] |
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/swed/e4_children.rds.
Click the following link to see the script used to generate this model:
opts_chunk$set(echo = FALSE)
clusterscript = str_replace(basename(model_filename), "\\.rds",".html")
cat("[Cluster script](" , clusterscript, ")", sep = "")
e5: Probability of ever divorcing
Here, we predict the probability of ever divorcing. To separate this effect from previous selective episodes, we include only ever-married anchors. Divorce data was only analysed in modern Sweden.
Model summary
Full summary
model_summary = summary(model, use_cache = FALSE, priors = TRUE)
print(model_summary)## Family: bernoulli(cauchit)
## Formula: ever_divorced ~ paternalage + spouses + 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: 1032368)
## Samples: 6 chains, each with iter = 800; warmup = 300; thin = 1;
## total post-warmup samples = 3000
## WAIC: Not computed
##
## Priors:
## b ~ normal(0,5)
## sd ~ student_t(3, 0, 5)
##
## Group-Level Effects:
## ~idParents (Number of levels: 719113)
## Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
## sd(Intercept) 0.81 0.01 0.79 0.83 427 1.01
##
## Population-Level Effects:
## Estimate Est.Error l-95% CI u-95% CI Eff.Sample
## Intercept -0.79 0.03 -0.84 -0.74 3000
## paternalage -0.03 0.02 -0.07 0.02 1573
## spouses 14.53 0.20 14.15 14.91 3000
## birth_cohort1950M1955 0.04 0.01 0.02 0.05 3000
## birth_cohort1955M1960 0.14 0.01 0.12 0.16 3000
## male1 -0.03 0.01 -0.05 -0.02 3000
## maternalage.factor1420 0.16 0.01 0.13 0.18 3000
## maternalage.factor3561 -0.11 0.01 -0.14 -0.09 3000
## paternalage.mean -0.17 0.02 -0.21 -0.12 1566
## paternal_loss01 -0.21 0.34 -0.94 0.37 3000
## paternal_loss15 0.35 0.09 0.17 0.52 3000
## paternal_loss510 0.37 0.04 0.29 0.45 3000
## paternal_loss1015 0.28 0.03 0.23 0.33 3000
## paternal_loss1520 0.24 0.02 0.20 0.28 3000
## paternal_loss2025 0.18 0.02 0.15 0.22 3000
## paternal_loss2530 0.15 0.02 0.12 0.18 3000
## paternal_loss3035 0.13 0.01 0.10 0.16 3000
## paternal_loss3540 0.08 0.01 0.06 0.10 3000
## paternal_loss4045 0.06 0.01 0.04 0.09 3000
## paternal_lossunclear 0.08 0.01 0.06 0.09 3000
## maternal_loss01 0.55 0.27 -0.02 1.07 3000
## maternal_loss15 0.06 0.14 -0.23 0.32 3000
## maternal_loss510 0.22 0.07 0.08 0.34 3000
## maternal_loss1015 0.16 0.04 0.08 0.25 3000
## maternal_loss1520 0.19 0.03 0.13 0.25 3000
## maternal_loss2025 0.15 0.03 0.09 0.20 3000
## maternal_loss2530 0.08 0.02 0.03 0.12 3000
## maternal_loss3035 0.13 0.02 0.09 0.16 3000
## maternal_loss3540 0.07 0.02 0.04 0.10 3000
## maternal_loss4045 0.04 0.02 0.02 0.07 3000
## maternal_lossunclear 0.01 0.01 0.00 0.03 3000
## older_siblings1 -0.03 0.01 -0.05 -0.01 1818
## older_siblings2 -0.03 0.02 -0.06 0.01 1650
## older_siblings3 0.02 0.03 -0.04 0.07 1571
## older_siblings4 0.04 0.04 -0.04 0.11 1608
## older_siblings5P 0.06 0.05 -0.03 0.15 1696
## nr.siblings 0.00 0.00 0.00 0.01 1599
## last_born1 -0.02 0.01 -0.04 -0.01 3000
## Rhat
## Intercept 1
## paternalage 1
## spouses 1
## birth_cohort1950M1955 1
## birth_cohort1955M1960 1
## male1 1
## maternalage.factor1420 1
## maternalage.factor3561 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
##
## 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 | 0.4529 | 0.4302 | 0.4759 |
| paternalage | 0.9711 | 0.9299 | 1.016 |
| spouses | 2043904 | 1402661 | 2995616 |
| birth_cohort1950M1955 | 1.036 | 1.017 | 1.056 |
| birth_cohort1955M1960 | 1.147 | 1.126 | 1.169 |
| male1 | 0.9688 | 0.9559 | 0.9813 |
| maternalage.factor1420 | 1.171 | 1.144 | 1.198 |
| maternalage.factor3561 | 0.893 | 0.8711 | 0.9159 |
| paternalage.mean | 0.8453 | 0.8083 | 0.8833 |
| paternal_loss01 | 0.808 | 0.3913 | 1.451 |
| paternal_loss15 | 1.415 | 1.184 | 1.68 |
| paternal_loss510 | 1.453 | 1.34 | 1.572 |
| paternal_loss1015 | 1.322 | 1.254 | 1.394 |
| paternal_loss1520 | 1.266 | 1.218 | 1.319 |
| paternal_loss2025 | 1.202 | 1.16 | 1.244 |
| paternal_loss2530 | 1.158 | 1.123 | 1.194 |
| paternal_loss3035 | 1.142 | 1.11 | 1.173 |
| paternal_loss3540 | 1.083 | 1.058 | 1.11 |
| paternal_loss4045 | 1.067 | 1.041 | 1.094 |
| paternal_lossunclear | 1.078 | 1.059 | 1.098 |
| maternal_loss01 | 1.73 | 0.9785 | 2.923 |
| maternal_loss15 | 1.058 | 0.7947 | 1.38 |
| maternal_loss510 | 1.24 | 1.087 | 1.408 |
| maternal_loss1015 | 1.179 | 1.085 | 1.278 |
| maternal_loss1520 | 1.214 | 1.14 | 1.29 |
| maternal_loss2025 | 1.157 | 1.099 | 1.219 |
| maternal_loss2530 | 1.081 | 1.035 | 1.129 |
| maternal_loss3035 | 1.134 | 1.093 | 1.175 |
| maternal_loss3540 | 1.068 | 1.036 | 1.102 |
| maternal_loss4045 | 1.046 | 1.015 | 1.078 |
| maternal_lossunclear | 1.012 | 0.9957 | 1.028 |
| older_siblings1 | 0.9732 | 0.9526 | 0.9943 |
| older_siblings2 | 0.9738 | 0.9376 | 1.01 |
| older_siblings3 | 1.016 | 0.9629 | 1.071 |
| older_siblings4 | 1.038 | 0.9649 | 1.115 |
| older_siblings5P | 1.058 | 0.9669 | 1.16 |
| nr.siblings | 1.005 | 0.9956 | 1.014 |
| last_born1 | 0.9755 | 0.9621 | 0.9891 |
Paternal age effect
pander::pander(paternal_age_10y_effect(model))| effect | median_estimate | ci_95 | ci_80 |
|---|---|---|---|
| estimate father 25y | 0.2 | [0.2;0.2] | [0.2;0.2] |
| estimate father 35y | 0.2 | [0.19;0.2] | [0.19;0.2] |
| percentage change | -0.32 | [-0.78;0.17] | [-0.62;0] |
| OR/IRR | 0.97 | [0.93;1.02] | [0.94;1] |
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/swed/e5_divorce.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 = "")All episodes
Here we show the effect of paternal age for each episode.
