Sweden (1947-2000), main models

Loading details

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

make_path = function(file) {
    get_coefficient_path(file, "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.1 dataset contains all 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.

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

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.

m1: No sibling comparison

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

Model summary

Full summary

model_summary = summary(model, use_cache = FALSE, priors = TRUE)
print(model_summary)
##  Family: poisson (log) 
## Formula: children ~ paternalage + birth_cohort + male + maternalage.factor + paternal_loss + maternal_loss + older_siblings + nr.siblings + last_born 
##    Data: model_data (Number of observations: 1408177) 
## Samples: 6 chains, each with iter = 800; warmup = 300; thin = 1; 
##          total post-warmup samples = 3000
##    WAIC: Not computed
##  
## Priors: 
## b ~ normal(0,5)
## 
## Population-Level Effects: 
##                        Estimate Est.Error l-95% CI u-95% CI Eff.Sample
## Intercept                  0.73      0.00     0.72     0.74       3000
## paternalage               -0.05      0.00    -0.05    -0.04       3000
## birth_cohort1950M1955      0.00      0.00     0.00     0.00       3000
## birth_cohort1955M1960      0.00      0.00     0.00     0.01       3000
## male1                     -0.06      0.00    -0.07    -0.06       3000
## maternalage.factor1420     0.06      0.00     0.05     0.06       3000
## maternalage.factor3561     0.00      0.00    -0.01     0.00       3000
## paternal_loss01            0.11      0.05     0.00     0.22       1327
## paternal_loss15            0.03      0.02     0.00     0.07       1893
## paternal_loss510          -0.02      0.01    -0.03     0.00       2062
## paternal_loss1015          0.00      0.01    -0.01     0.01       3000
## paternal_loss1520          0.00      0.00     0.00     0.01       3000
## paternal_loss2025          0.00      0.00    -0.01     0.00       3000
## paternal_loss2530          0.00      0.00     0.00     0.01       3000
## paternal_loss3035          0.00      0.00     0.00     0.00       3000
## paternal_loss3540          0.00      0.00    -0.01     0.00       3000
## paternal_loss4045         -0.01      0.00    -0.01     0.00       3000
## paternal_lossunclear      -0.06      0.00    -0.06    -0.05       3000
## maternal_loss01           -0.20      0.07    -0.34    -0.06       1665
## maternal_loss15           -0.06      0.03    -0.11     0.00       1497
## maternal_loss510          -0.01      0.01    -0.04     0.02       1759
## maternal_loss1015         -0.01      0.01    -0.02     0.01       2587
## maternal_loss1520          0.00      0.01    -0.02     0.01       2838
## maternal_loss2025          0.01      0.00     0.00     0.02       3000
## maternal_loss2530          0.01      0.00     0.00     0.01       3000
## maternal_loss3035          0.01      0.00     0.00     0.01       3000
## maternal_loss3540          0.00      0.00    -0.01     0.00       3000
## maternal_loss4045          0.00      0.00    -0.01     0.01       3000
## maternal_lossunclear      -0.02      0.00    -0.02    -0.02       3000
## older_siblings1            0.02      0.00     0.01     0.02       3000
## older_siblings2            0.02      0.00     0.01     0.02       3000
## older_siblings3            0.01      0.00     0.00     0.02       3000
## older_siblings4           -0.01      0.00    -0.02     0.00       3000
## older_siblings5P          -0.06      0.01    -0.07    -0.05       2350
## nr.siblings                0.04      0.00     0.04     0.04       2470
## last_born1                 0.01      0.00     0.01     0.01       3000
##                        Rhat
## Intercept                 1
## paternalage               1
## birth_cohort1950M1955     1
## birth_cohort1955M1960     1
## male1                     1
## maternalage.factor1420    1
## maternalage.factor3561    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

fixed_eff = data.frame(model_summary$fixed, check.names = F)
fixed_eff$Est.Error = fixed_eff$Eff.Sample = fixed_eff$Rhat = NULL
fixed_eff$OR = exp(fixed_eff$Estimate)
fixed_eff$OR_low = exp(fixed_eff$`l-95% CI`)
fixed_eff$OR_high = exp(fixed_eff$`u-95% CI`)
pander::pander(fixed_eff)
  Estimate l-95% CI u-95% CI OR OR_low OR_high
Intercept 0.7324 0.7238 0.741 2.08 2.062 2.098
paternalage -0.04706 -0.04951 -0.04454 0.954 0.9517 0.9564
birth_cohort1950M1955 0.0007506 -0.002472 0.003975 1.001 0.9975 1.004
birth_cohort1955M1960 0.001589 -0.001701 0.005026 1.002 0.9983 1.005
male1 -0.06361 -0.06598 -0.06126 0.9384 0.9361 0.9406
maternalage.factor1420 0.05691 0.0522 0.06142 1.059 1.054 1.063
maternalage.factor3561 -0.00131 -0.005571 0.002988 0.9987 0.9944 1.003
paternal_loss01 0.114 0.002498 0.2187 1.121 1.003 1.244
paternal_loss15 0.0314 -0.004265 0.06798 1.032 0.9957 1.07
paternal_loss510 -0.01675 -0.03403 0.0008674 0.9834 0.9665 1.001
paternal_loss1015 -0.003914 -0.0146 0.006414 0.9961 0.9855 1.006
paternal_loss1520 0.004768 -0.003012 0.01275 1.005 0.997 1.013
paternal_loss2025 -0.001746 -0.00844 0.004828 0.9983 0.9916 1.005
paternal_loss2530 0.003258 -0.001987 0.008665 1.003 0.998 1.009
paternal_loss3035 0.0002228 -0.004601 0.004953 1 0.9954 1.005
paternal_loss3540 -0.003974 -0.008541 0.0003649 0.996 0.9915 1
paternal_loss4045 -0.008584 -0.01299 -0.004322 0.9915 0.9871 0.9957
paternal_lossunclear -0.05531 -0.05867 -0.05186 0.9462 0.943 0.9495
maternal_loss01 -0.1959 -0.3429 -0.05812 0.8221 0.7097 0.9435
maternal_loss15 -0.05531 -0.1132 -0.0001841 0.9462 0.893 0.9998
maternal_loss510 -0.009533 -0.03653 0.01708 0.9905 0.9641 1.017
maternal_loss1015 -0.006085 -0.02159 0.009571 0.9939 0.9786 1.01
maternal_loss1520 -0.00468 -0.01746 0.006948 0.9953 0.9827 1.007
maternal_loss2025 0.008751 -0.001207 0.01849 1.009 0.9988 1.019
maternal_loss2530 0.00568 -0.002706 0.01384 1.006 0.9973 1.014
maternal_loss3035 0.005608 -0.001537 0.01246 1.006 0.9985 1.013
maternal_loss3540 -0.003965 -0.009709 0.001744 0.996 0.9903 1.002
maternal_loss4045 -0.0006131 -0.006377 0.005184 0.9994 0.9936 1.005
maternal_lossunclear -0.0204 -0.02348 -0.01729 0.9798 0.9768 0.9829
older_siblings1 0.01541 0.01239 0.01836 1.016 1.012 1.019
older_siblings2 0.01884 0.01459 0.02317 1.019 1.015 1.023
older_siblings3 0.01066 0.004143 0.01705 1.011 1.004 1.017
older_siblings4 -0.005764 -0.01521 0.003415 0.9943 0.9849 1.003
older_siblings5P -0.06164 -0.07247 -0.05062 0.9402 0.9301 0.9506
nr.siblings 0.03829 0.03718 0.03936 1.039 1.038 1.04
last_born1 0.01138 0.008642 0.01416 1.011 1.009 1.014

Paternal age effect

pander::pander(paternal_age_10y_effect(model))
effect median_estimate ci_95 ci_80
estimate father 25y 1.92 [1.91;1.93] [1.92;1.93]
estimate father 35y 1.83 [1.83;1.84] [1.83;1.84]
percentage change -4.6 [-4.83;-4.36] [-4.75;-4.44]
OR/IRR 0.95 [0.95;0.96] [0.95;0.96]

Marginal effect plots

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

plot of chunk unnamed-chunk-11plot of chunk unnamed-chunk-11plot of chunk unnamed-chunk-11plot of chunk unnamed-chunk-11plot of chunk unnamed-chunk-11plot of chunk unnamed-chunk-11plot of chunk unnamed-chunk-11plot of chunk unnamed-chunk-11plot of chunk unnamed-chunk-11

Coefficient plot

Coefficient estimates (95% and 80% credibility).

mcmc_intervals(as.matrix(model$fit), regex_pars = "b_[^I]")
plot of chunk unnamed-chunk-12

plot of chunk unnamed-chunk-12

mcmc_areas(as.matrix(model$fit), regex_pars = "b")
plot of chunk unnamed-chunk-12

plot of chunk unnamed-chunk-12

Diagnostics

These plots were made to diagnose misfit and nonconvergence.

Posterior predictive checks

brms::pp_check(model, re_formula = NA, type = "dens_overlay")
plot of chunk unnamed-chunk-13

plot of chunk unnamed-chunk-13

brms::pp_check(model, re_formula = NA, type = "hist")
plot of chunk unnamed-chunk-13

plot of chunk unnamed-chunk-13

Rhat

Did the 6 chains converge?

stanplot(model, pars = "^b_[^I]", type = 'rhat')
plot of chunk unnamed-chunk-14

plot of chunk unnamed-chunk-14

Effective sample size by average sample size

stanplot(model, pars = "^b", type = 'ess')
plot of chunk unnamed-chunk-15

plot of chunk unnamed-chunk-15

Monte Carlo SE

stanplot(model, pars = "^b", type = 'mcse')
plot of chunk unnamed-chunk-16

plot of chunk unnamed-chunk-16

Trace plots

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

Further plots

stanplot(model, pars = "^b", type = 'diag')
plot of chunk unnamed-chunk-18

plot of chunk unnamed-chunk-18

File name

coefs/swed/m1_children_linear_noranef.rds

Cluster script

opts_chunk$set(echo = FALSE)
clusterscript = str_replace(basename(model_filename), "\\.rds",".html")
cat("[Cluster script](" , clusterscript, ")", sep = "")

Cluster script

m2: Sibling comparison, no paternal age effect

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

Model summary

Full summary

model_summary = summary(model, use_cache = FALSE, priors = TRUE)
print(model_summary)
##  Family: poisson (log) 
## Formula: children ~ 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: 1408177) 
## Samples: 6 chains, each with iter = 1500; warmup = 1000; 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: 884975) 
##               Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
## sd(Intercept)        0         0        0     0.01         67 1.07
## 
## Population-Level Effects: 
##                        Estimate Est.Error l-95% CI u-95% CI Eff.Sample
## Intercept                  0.73      0.00     0.72     0.74       3000
## birth_cohort1950M1955      0.00      0.00     0.00     0.00       3000
## birth_cohort1955M1960      0.00      0.00     0.00     0.00       3000
## male1                     -0.06      0.00    -0.07    -0.06       3000
## maternalage.factor1420     0.06      0.00     0.05     0.06       3000
## maternalage.factor3561    -0.01      0.00    -0.01    -0.01       3000
## paternalage.mean          -0.04      0.00    -0.05    -0.04       3000
## paternal_loss01            0.11      0.05     0.02     0.21       3000
## paternal_loss15            0.03      0.02    -0.01     0.06       3000
## paternal_loss510          -0.02      0.01    -0.04     0.00       3000
## paternal_loss1015         -0.01      0.01    -0.02     0.00       3000
## paternal_loss1520          0.00      0.00    -0.01     0.01       3000
## paternal_loss2025          0.00      0.00    -0.01     0.00       3000
## paternal_loss2530          0.00      0.00     0.00     0.01       3000
## paternal_loss3035          0.00      0.00    -0.01     0.00       3000
## paternal_loss3540         -0.01      0.00    -0.01     0.00       3000
## paternal_loss4045         -0.01      0.00    -0.01     0.00       3000
## paternal_lossunclear      -0.05      0.00    -0.06    -0.05       3000
## maternal_loss01           -0.20      0.07    -0.34    -0.06       3000
## maternal_loss15           -0.06      0.03    -0.11     0.00       3000
## maternal_loss510          -0.01      0.01    -0.04     0.01       3000
## maternal_loss1015         -0.01      0.01    -0.02     0.01       3000
## maternal_loss1520          0.00      0.01    -0.02     0.01       3000
## maternal_loss2025          0.01      0.01     0.00     0.02       3000
## maternal_loss2530          0.01      0.00     0.00     0.01       3000
## maternal_loss3035          0.01      0.00     0.00     0.01       3000
## maternal_loss3540          0.00      0.00    -0.01     0.00       3000
## maternal_loss4045          0.00      0.00    -0.01     0.00       3000
## maternal_lossunclear      -0.02      0.00    -0.02    -0.02       3000
## older_siblings1            0.00      0.00     0.00     0.00       3000
## older_siblings2           -0.01      0.00    -0.02    -0.01       3000
## older_siblings3           -0.03      0.00    -0.04    -0.03       3000
## older_siblings4           -0.06      0.00    -0.07    -0.05       3000
## older_siblings5P          -0.14      0.01    -0.15    -0.12       3000
## nr.siblings                0.05      0.00     0.04     0.05       2555
## last_born1                 0.01      0.00     0.01     0.02       2871
##                        Rhat
## Intercept                 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

fixed_eff = data.frame(model_summary$fixed, check.names = F)
fixed_eff$Est.Error = fixed_eff$Eff.Sample = fixed_eff$Rhat = NULL
fixed_eff$OR = exp(fixed_eff$Estimate)
fixed_eff$OR_low = exp(fixed_eff$`l-95% CI`)
fixed_eff$OR_high = exp(fixed_eff$`u-95% CI`)
pander::pander(fixed_eff)
  Estimate l-95% CI u-95% CI OR OR_low OR_high
Intercept 0.7262 0.7174 0.735 2.067 2.049 2.086
birth_cohort1950M1955 -6.501e-05 -0.003164 0.00277 0.9999 0.9968 1.003
birth_cohort1955M1960 -6.963e-05 -0.003467 0.003152 0.9999 0.9965 1.003
male1 -0.0636 -0.06611 -0.06116 0.9384 0.936 0.9407
maternalage.factor1420 0.05861 0.05428 0.06309 1.06 1.056 1.065
maternalage.factor3561 -0.009328 -0.01344 -0.005119 0.9907 0.9867 0.9949
paternalage.mean -0.0445 -0.04695 -0.04204 0.9565 0.9541 0.9588
paternal_loss01 0.115 0.01517 0.215 1.122 1.015 1.24
paternal_loss15 0.02863 -0.008693 0.06492 1.029 0.9913 1.067
paternal_loss510 -0.01971 -0.0375 -0.002419 0.9805 0.9632 0.9976
paternal_loss1015 -0.006282 -0.01692 0.003985 0.9937 0.9832 1.004
paternal_loss1520 0.002715 -0.005341 0.01071 1.003 0.9947 1.011
paternal_loss2025 -0.003603 -0.01005 0.002664 0.9964 0.99 1.003
paternal_loss2530 0.001496 -0.004131 0.007162 1.001 0.9959 1.007
paternal_loss3035 -0.001137 -0.005895 0.003822 0.9989 0.9941 1.004
paternal_loss3540 -0.005074 -0.009499 -0.0006149 0.9949 0.9905 0.9994
paternal_loss4045 -0.009404 -0.01382 -0.004898 0.9906 0.9863 0.9951
paternal_lossunclear -0.05453 -0.05795 -0.05132 0.9469 0.9437 0.95
maternal_loss01 -0.1958 -0.3364 -0.05657 0.8222 0.7143 0.945
maternal_loss15 -0.05708 -0.1138 -0.002132 0.9445 0.8925 0.9979
maternal_loss510 -0.01076 -0.03655 0.01449 0.9893 0.9641 1.015
maternal_loss1015 -0.006463 -0.02209 0.00835 0.9936 0.9782 1.008
maternal_loss1520 -0.004837 -0.01664 0.007031 0.9952 0.9835 1.007
maternal_loss2025 0.008192 -0.002057 0.01836 1.008 0.9979 1.019
maternal_loss2530 0.005365 -0.003211 0.01374 1.005 0.9968 1.014
maternal_loss3035 0.005176 -0.001906 0.01239 1.005 0.9981 1.012
maternal_loss3540 -0.004266 -0.01004 0.001536 0.9957 0.99 1.002
maternal_loss4045 -0.001028 -0.006462 0.00426 0.999 0.9936 1.004
maternal_lossunclear -0.01928 -0.02222 -0.01636 0.9809 0.978 0.9838
older_siblings1 -0.0001088 -0.003114 0.002863 0.9999 0.9969 1.003
older_siblings2 -0.01167 -0.01612 -0.007456 0.9884 0.984 0.9926
older_siblings3 -0.03246 -0.03889 -0.02612 0.9681 0.9619 0.9742
older_siblings4 -0.06004 -0.0693 -0.05055 0.9417 0.9331 0.9507
older_siblings5P -0.135 -0.1459 -0.124 0.8737 0.8643 0.8833
nr.siblings 0.04507 0.04397 0.04613 1.046 1.045 1.047
last_born1 0.01446 0.0117 0.01714 1.015 1.012 1.017

Paternal age effect

pander::pander(paternal_age_10y_effect(model))

Marginal effect plots

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

plot of chunk unnamed-chunk-24plot of chunk unnamed-chunk-24plot of chunk unnamed-chunk-24plot of chunk unnamed-chunk-24plot of chunk unnamed-chunk-24plot of chunk unnamed-chunk-24plot of chunk unnamed-chunk-24plot of chunk unnamed-chunk-24plot of chunk unnamed-chunk-24

Coefficient plot

Coefficient estimates (95% and 80% credibility).

mcmc_intervals(as.matrix(model$fit), regex_pars = "b_[^I]")
plot of chunk unnamed-chunk-25

plot of chunk unnamed-chunk-25

mcmc_areas(as.matrix(model$fit), regex_pars = "b")
plot of chunk unnamed-chunk-25

plot of chunk unnamed-chunk-25

Diagnostics

These plots were made to diagnose misfit and nonconvergence.

Posterior predictive checks

brms::pp_check(model, re_formula = NA, type = "dens_overlay")
plot of chunk unnamed-chunk-26

plot of chunk unnamed-chunk-26

brms::pp_check(model, re_formula = NA, type = "hist")
plot of chunk unnamed-chunk-26

plot of chunk unnamed-chunk-26

Rhat

Did the 6 chains converge?

stanplot(model, pars = "^b_[^I]", type = 'rhat')
plot of chunk unnamed-chunk-27

plot of chunk unnamed-chunk-27

Effective sample size by average sample size

stanplot(model, pars = "^b", type = 'ess')
plot of chunk unnamed-chunk-28

plot of chunk unnamed-chunk-28

Monte Carlo SE

stanplot(model, pars = "^b", type = 'mcse')
plot of chunk unnamed-chunk-29

plot of chunk unnamed-chunk-29

Trace plots

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

Further plots

stanplot(model, pars = "^b", type = 'diag')
plot of chunk unnamed-chunk-31

plot of chunk unnamed-chunk-31

File name

coefs/swed/m2_children_linear_no_diff.rds

Cluster script

opts_chunk$set(echo = FALSE)
clusterscript = str_replace(basename(model_filename), "\\.rds",".html")
cat("[Cluster script](" , clusterscript, ")", sep = "")

Cluster script

m3: Sibling comparison, linear paternal age effect

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

Model summary

Full summary

model_summary = summary(model, use_cache = FALSE, priors = TRUE)
print(model_summary)
##  Family: poisson (log) 
## Formula: children ~ 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: 1408177) 
## 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: 884975) 
##               Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
## sd(Intercept)        0         0        0     0.01        169 1.02
## 
## Population-Level Effects: 
##                        Estimate Est.Error l-95% CI u-95% CI Eff.Sample
## Intercept                  0.73      0.00     0.72     0.74       3000
## paternalage               -0.05      0.00    -0.06    -0.05       3000
## birth_cohort1950M1955      0.00      0.00     0.00     0.00       3000
## birth_cohort1955M1960      0.00      0.00     0.00     0.01       3000
## male1                     -0.06      0.00    -0.07    -0.06       3000
## maternalage.factor1420     0.06      0.00     0.05     0.06       3000
## maternalage.factor3561     0.00      0.00    -0.01     0.00       3000
## paternalage.mean           0.01      0.00     0.00     0.01       3000
## paternal_loss01            0.11      0.05     0.01     0.22       3000
## paternal_loss15            0.03      0.02    -0.01     0.07       3000
## paternal_loss510          -0.02      0.01    -0.03     0.00       3000
## paternal_loss1015          0.00      0.01    -0.01     0.01       3000
## paternal_loss1520          0.00      0.00     0.00     0.01       3000
## paternal_loss2025          0.00      0.00    -0.01     0.00       3000
## paternal_loss2530          0.00      0.00     0.00     0.01       3000
## paternal_loss3035          0.00      0.00     0.00     0.01       3000
## paternal_loss3540          0.00      0.00    -0.01     0.00       3000
## paternal_loss4045         -0.01      0.00    -0.01     0.00       3000
## paternal_lossunclear      -0.06      0.00    -0.06    -0.05       3000
## maternal_loss01           -0.20      0.07    -0.34    -0.06       3000
## maternal_loss15           -0.06      0.03    -0.11     0.00       3000
## maternal_loss510          -0.01      0.01    -0.04     0.02       3000
## maternal_loss1015         -0.01      0.01    -0.02     0.01       3000
## maternal_loss1520          0.00      0.01    -0.02     0.01       3000
## maternal_loss2025          0.01      0.01     0.00     0.02       3000
## maternal_loss2530          0.01      0.00     0.00     0.01       3000
## maternal_loss3035          0.01      0.00     0.00     0.01       3000
## maternal_loss3540          0.00      0.00    -0.01     0.00       3000
## maternal_loss4045          0.00      0.00    -0.01     0.00       3000
## maternal_lossunclear      -0.02      0.00    -0.02    -0.02       3000
## older_siblings1            0.02      0.00     0.01     0.02       3000
## older_siblings2            0.02      0.00     0.02     0.03       3000
## older_siblings3            0.02      0.00     0.01     0.03       3000
## older_siblings4            0.00      0.01    -0.01     0.01       3000
## older_siblings5P          -0.05      0.01    -0.07    -0.04       3000
## nr.siblings                0.04      0.00     0.04     0.04       3000
## last_born1                 0.01      0.00     0.01     0.01       3000
##                        Rhat
## Intercept                 1
## paternalage               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

fixed_eff = data.frame(model_summary$fixed, check.names = F)
fixed_eff$Est.Error = fixed_eff$Eff.Sample = fixed_eff$Rhat = NULL
fixed_eff$OR = exp(fixed_eff$Estimate)
fixed_eff$OR_low = exp(fixed_eff$`l-95% CI`)
fixed_eff$OR_high = exp(fixed_eff$`u-95% CI`)
pander::pander(fixed_eff)
  Estimate l-95% CI u-95% CI OR OR_low OR_high
Intercept 0.7311 0.7225 0.7398 2.077 2.06 2.096
paternalage -0.05282 -0.06076 -0.04518 0.9485 0.941 0.9558
birth_cohort1950M1955 0.0008848 -0.002276 0.004109 1.001 0.9977 1.004
birth_cohort1955M1960 0.001805 -0.001487 0.005156 1.002 0.9985 1.005
male1 -0.0636 -0.06598 -0.06111 0.9384 0.9361 0.9407
maternalage.factor1420 0.05699 0.05241 0.06152 1.059 1.054 1.063
maternalage.factor3561 -0.0005762 -0.00518 0.003765 0.9994 0.9948 1.004
paternalage.mean 0.006093 -0.00192 0.01404 1.006 0.9981 1.014
paternal_loss01 0.115 0.0113 0.2163 1.122 1.011 1.241
paternal_loss15 0.03188 -0.005482 0.06891 1.032 0.9945 1.071
paternal_loss510 -0.0166 -0.03425 0.000544 0.9835 0.9663 1.001
paternal_loss1015 -0.003869 -0.01413 0.006115 0.9961 0.986 1.006
paternal_loss1520 0.004711 -0.00317 0.01238 1.005 0.9968 1.012
paternal_loss2025 -0.001669 -0.008168 0.004476 0.9983 0.9919 1.004
paternal_loss2530 0.003309 -0.002232 0.008929 1.003 0.9978 1.009
paternal_loss3035 0.0002486 -0.00474 0.005001 1 0.9953 1.005
paternal_loss3540 -0.003985 -0.00847 0.0003979 0.996 0.9916 1
paternal_loss4045 -0.00856 -0.01317 -0.00394 0.9915 0.9869 0.9961
paternal_lossunclear -0.05533 -0.05867 -0.05208 0.9462 0.943 0.9493
maternal_loss01 -0.1991 -0.3426 -0.06329 0.8195 0.7099 0.9387
maternal_loss15 -0.056 -0.1125 -0.0002241 0.9455 0.8936 0.9998
maternal_loss510 -0.00949 -0.03583 0.0167 0.9906 0.9648 1.017
maternal_loss1015 -0.005822 -0.0218 0.01011 0.9942 0.9784 1.01
maternal_loss1520 -0.004533 -0.01646 0.00713 0.9955 0.9837 1.007
maternal_loss2025 0.008532 -0.001182 0.01845 1.009 0.9988 1.019
maternal_loss2530 0.005649 -0.002824 0.01408 1.006 0.9972 1.014
maternal_loss3035 0.005573 -0.00139 0.01257 1.006 0.9986 1.013
maternal_loss3540 -0.003918 -0.009787 0.001973 0.9961 0.9903 1.002
maternal_loss4045 -0.0007113 -0.00606 0.004642 0.9993 0.994 1.005
maternal_lossunclear -0.02044 -0.02341 -0.01756 0.9798 0.9769 0.9826
older_siblings1 0.01738 0.01336 0.02154 1.018 1.013 1.022
older_siblings2 0.02266 0.01575 0.02961 1.023 1.016 1.03
older_siblings3 0.01614 0.006488 0.02588 1.016 1.007 1.026
older_siblings4 0.001064 -0.01209 0.01423 1.001 0.988 1.014
older_siblings5P -0.05232 -0.06885 -0.03542 0.949 0.9335 0.9652
nr.siblings 0.03742 0.0358 0.039 1.038 1.036 1.04
last_born1 0.01095 0.008286 0.01367 1.011 1.008 1.014

Paternal age effect

pander::pander(paternal_age_10y_effect(model))
effect median_estimate ci_95 ci_80
estimate father 25y 1.93 [1.92;1.94] [1.92;1.93]
estimate father 35y 1.83 [1.82;1.84] [1.82;1.83]
percentage change -5.14 [-5.9;-4.42] [-5.63;-4.66]
OR/IRR 0.95 [0.94;0.96] [0.94;0.95]

Marginal effect plots

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

plot of chunk unnamed-chunk-37plot of chunk unnamed-chunk-37plot of chunk unnamed-chunk-37plot of chunk unnamed-chunk-37plot of chunk unnamed-chunk-37plot of chunk unnamed-chunk-37plot of chunk unnamed-chunk-37plot of chunk unnamed-chunk-37plot of chunk unnamed-chunk-37plot of chunk unnamed-chunk-37

Coefficient plot

Coefficient estimates (95% and 80% credibility).

mcmc_intervals(as.matrix(model$fit), regex_pars = "b_[^I]")
plot of chunk unnamed-chunk-38

plot of chunk unnamed-chunk-38

mcmc_areas(as.matrix(model$fit), regex_pars = "b")
plot of chunk unnamed-chunk-38

plot of chunk unnamed-chunk-38

Diagnostics

These plots were made to diagnose misfit and nonconvergence.

Posterior predictive checks

brms::pp_check(model, re_formula = NA, type = "dens_overlay")
plot of chunk unnamed-chunk-39

plot of chunk unnamed-chunk-39

brms::pp_check(model, re_formula = NA, type = "hist")
plot of chunk unnamed-chunk-39

plot of chunk unnamed-chunk-39

Rhat

Did the 6 chains converge?

stanplot(model, pars = "^b_[^I]", type = 'rhat')
plot of chunk unnamed-chunk-40

plot of chunk unnamed-chunk-40

Effective sample size by average sample size

stanplot(model, pars = "^b", type = 'ess')
plot of chunk unnamed-chunk-41

plot of chunk unnamed-chunk-41

Monte Carlo SE

stanplot(model, pars = "^b", type = 'mcse')
plot of chunk unnamed-chunk-42

plot of chunk unnamed-chunk-42

Trace plots

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

Further plots

stanplot(model, pars = "^b", type = 'diag')
plot of chunk unnamed-chunk-44

plot of chunk unnamed-chunk-44

File name

coefs/swed/m3_children_linear.rds

Cluster script

opts_chunk$set(echo = FALSE)
clusterscript = str_replace(basename(model_filename), "\\.rds",".html")
cat("[Cluster script](" , clusterscript, ")", sep = "")

Cluster script

m4: Sibling comparison, nonlinear paternal age effect

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

Model summary

Full summary

model_summary = summary(model, use_cache = FALSE, priors = TRUE)
print(model_summary)
##  Family: poisson (log) 
## Formula: children ~ s(paternalage) + birth_cohort + male + maternalage.factor + paternalage.mean + paternal_loss + maternal_loss + older_siblings + nr.siblings + last_born + (1 | idParents) 
##    Data: model_data (Number of observations: 1408177) 
## 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)
## sds ~ student_t(3, 0, 10)
## 
## Spline Effects: 
##                     Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
## sds(spaternalage_1)     0.13      0.07     0.06     0.37         35 1.14
## 
## Group-Level Effects: 
## ~idParents (Number of levels: 884975) 
##               Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
## sd(Intercept)        0         0        0     0.01        152 1.05
## 
## Population-Level Effects: 
##                        Estimate Est.Error l-95% CI u-95% CI Eff.Sample
## Intercept                  0.56      0.01     0.53     0.59       3000
## birth_cohort1950M1955      0.00      0.00     0.00     0.00       3000
## birth_cohort1955M1960      0.00      0.00     0.00     0.01       3000
## male1                     -0.06      0.00    -0.07    -0.06        273
## maternalage.factor1420     0.05      0.00     0.04     0.05       3000
## maternalage.factor3561    -0.01      0.00    -0.01     0.00        118
## paternalage.mean           0.01      0.00     0.00     0.02       3000
## paternal_loss01            0.10      0.05     0.00     0.20       3000
## paternal_loss15            0.03      0.02    -0.01     0.06       3000
## paternal_loss510          -0.02      0.01    -0.04    -0.01       3000
## paternal_loss1015         -0.01      0.01    -0.02     0.00       3000
## paternal_loss1520          0.00      0.00    -0.01     0.01       3000
## paternal_loss2025         -0.01      0.00    -0.01     0.00        122
## paternal_loss2530          0.00      0.00    -0.01     0.01       2499
## paternal_loss3035          0.00      0.00    -0.01     0.00       3000
## paternal_loss3540          0.00      0.00    -0.01     0.00       3000
## paternal_loss4045         -0.01      0.00    -0.01     0.00       3000
## paternal_lossunclear      -0.06      0.00    -0.06    -0.06       3000
## maternal_loss01           -0.20      0.07    -0.34    -0.06       3000
## maternal_loss15           -0.06      0.03    -0.12     0.00       3000
## maternal_loss510          -0.01      0.01    -0.04     0.02       3000
## maternal_loss1015         -0.01      0.01    -0.02     0.01       3000
## maternal_loss1520         -0.01      0.01    -0.02     0.01       3000
## maternal_loss2025          0.01      0.00     0.00     0.02       3000
## maternal_loss2530          0.00      0.00     0.00     0.01       3000
## maternal_loss3035          0.00      0.00     0.00     0.01       3000
## maternal_loss3540          0.00      0.00    -0.01     0.00        223
## maternal_loss4045          0.00      0.00    -0.01     0.00        477
## maternal_lossunclear      -0.02      0.00    -0.02    -0.02       1602
## older_siblings1            0.02      0.00     0.02     0.03       3000
## older_siblings2            0.03      0.00     0.02     0.04       3000
## older_siblings3            0.02      0.00     0.01     0.03       3000
## older_siblings4            0.01      0.01    -0.01     0.02       3000
## older_siblings5P          -0.05      0.01    -0.07    -0.03       3000
## nr.siblings                0.04      0.00     0.04     0.04       3000
## last_born1                 0.01      0.00     0.01     0.01        593
## spaternalage_1            -0.04      0.02    -0.08     0.00        148
##                        Rhat
## Intercept              1.01
## birth_cohort1950M1955  1.00
## birth_cohort1955M1960  1.01
## male1                  1.02
## maternalage.factor1420 1.00
## maternalage.factor3561 1.03
## paternalage.mean       1.01
## paternal_loss01        1.01
## paternal_loss15        1.01
## paternal_loss510       1.00
## paternal_loss1015      1.00
## paternal_loss1520      1.01
## paternal_loss2025      1.04
## paternal_loss2530      1.00
## paternal_loss3035      1.00
## paternal_loss3540      1.00
## paternal_loss4045      1.00
## paternal_lossunclear   1.00
## maternal_loss01        1.01
## maternal_loss15        1.00
## maternal_loss510       1.03
## maternal_loss1015      1.01
## maternal_loss1520      1.02
## maternal_loss2025      1.01
## maternal_loss2530      1.02
## maternal_loss3035      1.01
## maternal_loss3540      1.02
## maternal_loss4045      1.01
## maternal_lossunclear   1.01
## older_siblings1        1.01
## older_siblings2        1.01
## older_siblings3        1.00
## older_siblings4        1.00
## older_siblings5P       1.02
## nr.siblings            1.00
## last_born1             1.00
## spaternalage_1         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

fixed_eff = data.frame(model_summary$fixed, check.names = F)
fixed_eff$Est.Error = fixed_eff$Eff.Sample = fixed_eff$Rhat = NULL
fixed_eff$OR = exp(fixed_eff$Estimate)
fixed_eff$OR_low = exp(fixed_eff$`l-95% CI`)
fixed_eff$OR_high = exp(fixed_eff$`u-95% CI`)
pander::pander(fixed_eff)
  Estimate l-95% CI u-95% CI OR OR_low OR_high
Intercept 0.5593 0.5339 0.5866 1.749 1.706 1.798
birth_cohort1950M1955 0.001172 -0.001897 0.004376 1.001 0.9981 1.004
birth_cohort1955M1960 0.002719 -0.0005264 0.006003 1.003 0.9995 1.006
male1 -0.06348 -0.06602 -0.06104 0.9385 0.9361 0.9408
maternalage.factor1420 0.04514 0.04028 0.04978 1.046 1.041 1.051
maternalage.factor3561 -0.007883 -0.01281 -0.002873 0.9921 0.9873 0.9971
paternalage.mean 0.008075 -0.0001993 0.01598 1.008 0.9998 1.016
paternal_loss01 0.1042 -0.002038 0.2026 1.11 0.998 1.225
paternal_loss15 0.02511 -0.01044 0.05866 1.025 0.9896 1.06
paternal_loss510 -0.02381 -0.04104 -0.006164 0.9765 0.9598 0.9939
paternal_loss1015 -0.01018 -0.02101 2.37e-05 0.9899 0.9792 1
paternal_loss1520 -0.0005987 -0.008226 0.007141 0.9994 0.9918 1.007
paternal_loss2025 -0.005697 -0.01248 0.001058 0.9943 0.9876 1.001
paternal_loss2530 0.0002694 -0.005244 0.005709 1 0.9948 1.006
paternal_loss3035 -0.001436 -0.006377 0.003384 0.9986 0.9936 1.003
paternal_loss3540 -0.0047 -0.009246 -0.0003081 0.9953 0.9908 0.9997
paternal_loss4045 -0.008531 -0.0129 -0.004146 0.9915 0.9872 0.9959
paternal_lossunclear -0.05946 -0.06285 -0.05618 0.9423 0.9391 0.9454
maternal_loss01 -0.2011 -0.3358 -0.05844 0.8178 0.7148 0.9432
maternal_loss15 -0.05905 -0.1165 -0.0003786 0.9427 0.89 0.9996
maternal_loss510 -0.01214 -0.03698 0.01543 0.9879 0.9637 1.016
maternal_loss1015 -0.007337 -0.02176 0.007747 0.9927 0.9785 1.008
maternal_loss1520 -0.00575 -0.01665 0.006421 0.9943 0.9835 1.006
maternal_loss2025 0.007788 -0.001614 0.01703 1.008 0.9984 1.017
maternal_loss2530 0.004332 -0.003913 0.01282 1.004 0.9961 1.013
maternal_loss3035 0.004257 -0.002523 0.01143 1.004 0.9975 1.011
maternal_loss3540 -0.004594 -0.01082 0.0009524 0.9954 0.9892 1.001
maternal_loss4045 -0.001348 -0.006585 0.00421 0.9987 0.9934 1.004
maternal_lossunclear -0.02138 -0.02436 -0.01845 0.9789 0.9759 0.9817
older_siblings1 0.02158 0.01761 0.02545 1.022 1.018 1.026
older_siblings2 0.02872 0.02197 0.03521 1.029 1.022 1.036
older_siblings3 0.02216 0.01296 0.0317 1.022 1.013 1.032
older_siblings4 0.006698 -0.006324 0.02004 1.007 0.9937 1.02
older_siblings5P -0.04837 -0.06563 -0.03255 0.9528 0.9365 0.968
nr.siblings 0.03664 0.03503 0.03819 1.037 1.036 1.039
last_born1 0.01137 0.008597 0.01414 1.011 1.009 1.014
spaternalage_1 -0.0413 -0.07926 -0.004906 0.9595 0.9238 0.9951

Paternal age effect

pander::pander(paternal_age_10y_effect(model))
effect median_estimate ci_95 ci_80
estimate father 25y 1.94 [1.93;1.95] [1.94;1.95]
estimate father 35y 1.81 [1.8;1.82] [1.8;1.82]
percentage change -6.88 [-7.62;-6.03] [-7.4;-6.32]
OR/IRR 0.96 [0.92;1] [0.93;0.98]

Marginal effect plots

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

plot of chunk unnamed-chunk-50plot of chunk unnamed-chunk-50plot of chunk unnamed-chunk-50plot of chunk unnamed-chunk-50plot of chunk unnamed-chunk-50plot of chunk unnamed-chunk-50plot of chunk unnamed-chunk-50plot of chunk unnamed-chunk-50plot of chunk unnamed-chunk-50plot of chunk unnamed-chunk-50

Coefficient plot

Coefficient estimates (95% and 80% credibility).

mcmc_intervals(as.matrix(model$fit), regex_pars = "b_[^I]")
plot of chunk unnamed-chunk-51

plot of chunk unnamed-chunk-51

mcmc_areas(as.matrix(model$fit), regex_pars = "b")
plot of chunk unnamed-chunk-51

plot of chunk unnamed-chunk-51

Diagnostics

These plots were made to diagnose misfit and nonconvergence.

Posterior predictive checks

brms::pp_check(model, re_formula = NA, type = "dens_overlay")
plot of chunk unnamed-chunk-52

plot of chunk unnamed-chunk-52

brms::pp_check(model, re_formula = NA, type = "hist")
plot of chunk unnamed-chunk-52

plot of chunk unnamed-chunk-52

Rhat

Did the 6 chains converge?

stanplot(model, pars = "^b_[^I]", type = 'rhat')
plot of chunk unnamed-chunk-53

plot of chunk unnamed-chunk-53

Effective sample size by average sample size

stanplot(model, pars = "^b", type = 'ess')
plot of chunk unnamed-chunk-54

plot of chunk unnamed-chunk-54

Monte Carlo SE

stanplot(model, pars = "^b", type = 'mcse')
plot of chunk unnamed-chunk-55

plot of chunk unnamed-chunk-55

Trace plots

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

Further plots

stanplot(model, pars = "^b", type = 'diag')
plot of chunk unnamed-chunk-57

plot of chunk unnamed-chunk-57

File name

coefs/swed/m4_children_nonlinear.rds

Cluster script

opts_chunk$set(echo = FALSE)
clusterscript = str_replace(basename(model_filename), "\\.rds",".html")
cat("[Cluster script](" , clusterscript, ")", sep = "")

Cluster script

Model comparison

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

##              LOOIC      SE
## m1      4598045.58 1277.67
## m2      4598206.12 1277.86
## m3      4598045.15 1277.67
## m4      4597745.82 1276.88
## m1 - m2    -160.55   21.28
## m1 - m3       0.43    2.75
## m1 - m4     299.76   33.75
## m2 - m3     160.97   23.86
## m2 - m4     460.31   41.08
## m3 - m4     299.34   33.61

plot of chunk unnamed-chunk-6plot of chunk unnamed-chunk-6plot of chunk unnamed-chunk-6plot of chunk unnamed-chunk-6

##               WAIC      SE
## m1      4598045.58 1277.67
## m2      4598206.12 1277.86
## m3      4598045.15 1277.67
## m4      4597745.82 1276.88
## m1 - m2    -160.54   21.28
## m1 - m3       0.43    2.75
## m1 - m4     299.76   33.75
## m2 - m3     160.97   23.86
## m2 - m4     460.31   41.08
## m3 - m4     299.33   33.61