First load brms:
Loading required package: RcppLoading 'brms' package (version 2.21.0). Useful instructions
can be found by typing help('brms'). A more detailed introduction
to the package is available through vignette('brms_overview').
Attaching package: 'brms'The following object is masked from 'package:stats':
arFitting a multivariate model in brms involves several new consideration above those for fitting univariate models. First, we need to create two models/objects with the function bf fitting the desired univariate model structure for each response variable (here bwt and tarsus). It is the equivalent of writing mvbf(bwt, tarsus), but the advantage to create two distinct model is to specific different model structure (fixed or random effect) for each response variable.
Then, the two objects/models are added into a third model to quantify all the estimates in addition to their covariance. Contrary to MCMCglmm or asreml-R, brms directly estimate the covariance and the correlation in its outputs. Our most basic model can be specified as:
bf_bwt <- bf(bwt ~ 1 + (1 | a | gr(animal, cov = Amat)))
bf_tarsus <- bf(tarsus ~ 1 + (1 | a | gr(animal, cov = Amat)))
brms_m2.1 <- brm(
bf_bwt + bf_tarsus + set_rescor(TRUE),
data = gryphon,
data2 = list(Amat = Amat),
chains = 2, cores = 2, iter = 1000
)
save(brms_m2.1, file = "r-obj/brms_m2_1.rda")Again we have provided the data from one such run. It can be accessed using the code:
Warning: Parts of the model have not converged (some Rhats are > 1.05). Be
careful when analysing the results! We recommend running more iterations and/or
setting stronger priors. Family: MV(gaussian, gaussian)
Links: mu = identity; sigma = identity
mu = identity; sigma = identity
Formula: bwt ~ 1 + (1 | p | gr(animal, cov = Amat))
tarsus ~ 1 + (1 | p | gr(animal, cov = Amat))
Data: gryphon (Number of observations: 683)
Draws: 2 chains, each with iter = 1000; warmup = 500; thin = 1;
total post-warmup draws = 1000
Multilevel Hyperparameters:
~animal (Number of levels: 683)
Estimate Est.Error l-95% CI u-95% CI Rhat
sd(bwt_Intercept) 1.81 0.21 1.41 2.20 1.06
sd(tarsus_Intercept) 3.44 0.43 2.49 4.25 1.05
cor(bwt_Intercept,tarsus_Intercept) 0.38 0.14 0.08 0.62 1.02
Bulk_ESS Tail_ESS
sd(bwt_Intercept) 31 192
sd(tarsus_Intercept) 61 173
cor(bwt_Intercept,tarsus_Intercept) 101 232
Regression Coefficients:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
bwt_Intercept 7.49 0.16 7.20 7.79 1.00 608 839
tarsus_Intercept 20.47 0.30 19.92 21.03 1.00 868 803
Further Distributional Parameters:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sigma_bwt 1.97 0.16 1.66 2.28 1.06 27 172
sigma_tarsus 4.24 0.30 3.63 4.82 1.04 72 162
Residual Correlations:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
rescor(bwt,tarsus) 0.39 0.09 0.21 0.55 1.02 95 179
Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).plot(brms_m2.1, ask = FALSE)
VarCorr(brms_m2.1)$animal
$animal$sd
Estimate Est.Error Q2.5 Q97.5
bwt_Intercept 1.808171 0.2050233 1.412824 2.204805
tarsus_Intercept 3.438368 0.4283612 2.491218 4.245264
$animal$cor
, , bwt_Intercept
Estimate Est.Error Q2.5 Q97.5
bwt_Intercept 1.0000000 0.0000000 1.00000000 1.0000000
tarsus_Intercept 0.3814062 0.1380014 0.07581464 0.6209038
, , tarsus_Intercept
Estimate Est.Error Q2.5 Q97.5
bwt_Intercept 0.3814062 0.1380014 0.07581464 0.6209038
tarsus_Intercept 1.0000000 0.0000000 1.00000000 1.0000000
$animal$cov
, , bwt_Intercept
Estimate Est.Error Q2.5 Q97.5
bwt_Intercept 3.311473 0.7430185 1.9960721 4.861167
tarsus_Intercept 2.440166 1.0901689 0.3870783 4.668720
, , tarsus_Intercept
Estimate Est.Error Q2.5 Q97.5
bwt_Intercept 2.440166 1.090169 0.3870783 4.66872
tarsus_Intercept 12.005688 2.918741 6.2061701 18.02226
$residual__
$residual__$sd
Estimate Est.Error Q2.5 Q97.5
bwt 1.970532 0.1597581 1.658782 2.276074
tarsus 4.244704 0.2984518 3.632824 4.820109
$residual__$cor
, , bwt
Estimate Est.Error Q2.5 Q97.5
bwt 1.0000000 0.00000000 1.0000000 1.0000000
tarsus 0.3888754 0.08510488 0.2127907 0.5526631
, , tarsus
Estimate Est.Error Q2.5 Q97.5
bwt 0.3888754 0.08510488 0.2127907 0.5526631
tarsus 1.0000000 0.00000000 1.0000000 1.0000000
$residual__$cov
, , bwt
Estimate Est.Error Q2.5 Q97.5
bwt 3.908493 0.6282892 2.751557 5.180511
tarsus 3.289995 0.9305960 1.572647 5.147133
, , tarsus
Estimate Est.Error Q2.5 Q97.5
bwt 3.289995 0.930596 1.572647 5.147133
tarsus 18.106495 2.530138 13.197409 23.233452It is also possible to calculate the heritability for each trait using the function ‘as.mcmc’
v_animal <- (VarCorr(brms_m2.1, summary = FALSE)$animal$sd)^2
v_r <- (VarCorr(brms_m2.1, summary = FALSE)$residual$sd)^2
h.bwt.2 <- as.mcmc(v_animal[, 1] / (v_animal[, 1] + v_r[, 1]))
h.tarsus.2 <- as.mcmc(v_animal[, 2] / (v_animal[, 2] + v_r[, 2]))
summary(h.bwt.2)
Iterations = 1:1000
Thinning interval = 1
Number of chains = 1
Sample size per chain = 1000
1. Empirical mean and standard deviation for each variable,
plus standard error of the mean:
Mean SD Naive SE Time-series SE
0.457051 0.090878 0.002874 0.011675
2. Quantiles for each variable:
2.5% 25% 50% 75% 97.5%
0.2878 0.3926 0.4596 0.5254 0.6297 summary(h.tarsus.2)
Iterations = 1:1000
Thinning interval = 1
Number of chains = 1
Sample size per chain = 1000
1. Empirical mean and standard deviation for each variable,
plus standard error of the mean:
Mean SD Naive SE Time-series SE
0.397237 0.087350 0.002762 0.009971
2. Quantiles for each variable:
2.5% 25% 50% 75% 97.5%
0.2174 0.3390 0.3982 0.4553 0.5682 plot(h.bwt.2)
plot(h.tarsus.2)
It is also possible to extract the correlation. Just to remember it is an example, the correlation distribution is skewed to 1 due to a weak prior and model parameters. Note, since
cor_g <- as.mcmc((VarCorr(brms_m2.1, summary = FALSE)$animal$cor[, 1, 2]))
cor_res <- as.mcmc((VarCorr(brms_m2.1, summary = FALSE)$residual$cor[, 1, 2]))
summary(cor_g)
Iterations = 1:1000
Thinning interval = 1
Number of chains = 1
Sample size per chain = 1000
1. Empirical mean and standard deviation for each variable,
plus standard error of the mean:
Mean SD Naive SE Time-series SE
0.381406 0.138001 0.004364 0.014946
2. Quantiles for each variable:
2.5% 25% 50% 75% 97.5%
0.07581 0.30354 0.39041 0.47497 0.62090 summary(cor_res)
Iterations = 1:1000
Thinning interval = 1
Number of chains = 1
Sample size per chain = 1000
1. Empirical mean and standard deviation for each variable,
plus standard error of the mean:
Mean SD Naive SE Time-series SE
0.388875 0.085105 0.002691 0.009215
2. Quantiles for each variable:
2.5% 25% 50% 75% 97.5%
0.2128 0.3319 0.3913 0.4484 0.5527 plot(cor_g)
plot(cor_res)
Here we can plot the genetic correlation by extraction the breeding values or BLUP.
Here, some simple code to plot the genetic correlation.
plot(tarsus_Intercept ~ bwt_Intercept, bl_m2.1,
xlab = "", ylab = "",
xlim = c(min(bl_m2.1$bwt_Intercept_lo), max(bl_m2.1$bwt_Intercept_up)),
ylim = c(min(bl_m2.1$tarsus_Intercept_lo), max(bl_m2.1$tarsus_Intercept_up)),
las = 1.2, type = "n"
)
with(
bl_m2.1,
segments(
x0 = bwt_Intercept, y0 = tarsus_Intercept_lo,
x1 = bwt_Intercept, y1 = tarsus_Intercept_up,
col = "black"
)
)
with(bl_m2.1, segments(
x0 = bwt_Intercept_lo, y0 = tarsus_Intercept,
x1 = bwt_Intercept_up, y1 = tarsus_Intercept,
col = "black"
))
points(tarsus_Intercept ~ bwt_Intercept, bl_m2.1, pch = 16, col = "red", cex = 1.5)
points(tarsus_Intercept ~ bwt_Intercept, bl_m2.1, pch = 1, col = rgb(0, 0, 0, 0.3), cex = c(1.5))
mtext("btw (BV±CI)", side = 1, line = 2.4)
mtext("tarsus (BV±CI)", side = 2, line = 2, las = 3)
Fixed and random effects can be added just as for the univariate case. Given that our full model of bwt from tutorial 1 had sex as a fixed effect as well as random effects of byear and mother, we could specify a bivariate formulation of this using the following code (including a line to save the output):
bf_bwt_2 <- bf(bwt ~ 1 + sex + (1 | a | gr(animal, cov = Amat)) + (1 | b | byear) + (1 | c | mother))
bf_tarsus_2 <- bf(tarsus ~ 1 + sex + (1 | a | gr(animal, cov = Amat)) + (1 | b | byear) + (1 | c | mother))
brms_m2.2 <- brm(
bf_bwt_2 + bf_tarsus_2 + set_rescor(TRUE),
data = gryphon,
data2 = list(Amat = Amat),
chains = 2, cores = 2, iter = 1000
)
save(brms_m2.2, file = "r-obj/brms_m2_2.rda")Again we have provided the data from one such run. It can be accessed using the code:
Warning: Parts of the model have not converged (some Rhats are > 1.05). Be
careful when analysing the results! We recommend running more iterations and/or
setting stronger priors.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 Family: MV(gaussian, gaussian)
Links: mu = identity; sigma = identity
mu = identity; sigma = identity
Formula: bwt ~ 1 + sex + (1 | a | gr(animal, cov = Amat)) + (1 | b | byear) + (1 | c | mother)
tarsus ~ 1 + sex + (1 | a | gr(animal, cov = Amat)) + (1 | b | byear) + (1 | c | mother)
Data: gryphon (Number of observations: 683)
Draws: 2 chains, each with iter = 1000; warmup = 500; thin = 1;
total post-warmup draws = 1000
Multilevel Hyperparameters:
~animal (Number of levels: 683)
Estimate Est.Error l-95% CI u-95% CI Rhat
sd(bwt_Intercept) 1.31 0.21 0.86 1.69 1.06
sd(tarsus_Intercept) 2.88 0.47 1.88 3.70 1.01
cor(bwt_Intercept,tarsus_Intercept) 0.60 0.17 0.18 0.89 1.09
Bulk_ESS Tail_ESS
sd(bwt_Intercept) 54 58
sd(tarsus_Intercept) 52 162
cor(bwt_Intercept,tarsus_Intercept) 25 30
~byear (Number of levels: 34)
Estimate Est.Error l-95% CI u-95% CI Rhat
sd(bwt_Intercept) 0.99 0.17 0.71 1.39 1.00
sd(tarsus_Intercept) 2.02 0.34 1.45 2.78 1.00
cor(bwt_Intercept,tarsus_Intercept) 0.01 0.22 -0.44 0.44 1.01
Bulk_ESS Tail_ESS
sd(bwt_Intercept) 424 509
sd(tarsus_Intercept) 525 699
cor(bwt_Intercept,tarsus_Intercept) 478 492
~mother (Number of levels: 352)
Estimate Est.Error l-95% CI u-95% CI Rhat
sd(bwt_Intercept) 1.14 0.12 0.90 1.36 1.01
sd(tarsus_Intercept) 2.09 0.29 1.54 2.67 1.01
cor(bwt_Intercept,tarsus_Intercept) -0.64 0.20 -0.97 -0.24 1.02
Bulk_ESS Tail_ESS
sd(bwt_Intercept) 370 764
sd(tarsus_Intercept) 134 420
cor(bwt_Intercept,tarsus_Intercept) 84 147
Regression Coefficients:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
bwt_Intercept 6.28 0.24 5.80 6.73 1.00 453 703
tarsus_Intercept 20.39 0.52 19.45 21.39 1.00 755 760
bwt_sex2 2.05 0.17 1.71 2.37 1.00 1097 715
tarsus_sex2 0.11 0.42 -0.67 0.90 1.00 780 578
Further Distributional Parameters:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sigma_bwt 1.40 0.16 1.05 1.68 1.04 59 60
sigma_tarsus 3.73 0.32 3.14 4.32 1.00 55 176
Residual Correlations:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
rescor(bwt,tarsus) 0.89 0.07 0.72 0.98 1.50 4 24
Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).plot(brms_m2.2, ask = FALSE)
VarCorr(brms_m2.2)$animal
$animal$sd
Estimate Est.Error Q2.5 Q97.5
bwt_Intercept 1.306921 0.2109780 0.8555702 1.687254
tarsus_Intercept 2.876731 0.4718524 1.8827128 3.699873
$animal$cor
, , bwt_Intercept
Estimate Est.Error Q2.5 Q97.5
bwt_Intercept 1.0000000 0.0000000 1.0000000 1.0000000
tarsus_Intercept 0.5961514 0.1727943 0.1849725 0.8929607
, , tarsus_Intercept
Estimate Est.Error Q2.5 Q97.5
bwt_Intercept 0.5961514 0.1727943 0.1849725 0.8929607
tarsus_Intercept 1.0000000 0.0000000 1.0000000 1.0000000
$animal$cov
, , bwt_Intercept
Estimate Est.Error Q2.5 Q97.5
bwt_Intercept 1.752510 0.5433004 0.7320013 2.846827
tarsus_Intercept 2.357852 1.0462644 0.3960956 4.526731
, , tarsus_Intercept
Estimate Est.Error Q2.5 Q97.5
bwt_Intercept 2.357852 1.046264 0.3960956 4.526731
tarsus_Intercept 8.498003 2.672465 3.5446115 13.689062
$byear
$byear$sd
Estimate Est.Error Q2.5 Q97.5
bwt_Intercept 0.9904858 0.1696752 0.7102391 1.390260
tarsus_Intercept 2.0166804 0.3366954 1.4463862 2.777005
$byear$cor
, , bwt_Intercept
Estimate Est.Error Q2.5 Q97.5
bwt_Intercept 1.00000000 0.0000000 1.0000000 1.0000000
tarsus_Intercept 0.01351405 0.2206186 -0.4367809 0.4412319
, , tarsus_Intercept
Estimate Est.Error Q2.5 Q97.5
bwt_Intercept 0.01351405 0.2206186 -0.4367809 0.4412319
tarsus_Intercept 1.00000000 0.0000000 1.0000000 1.0000000
$byear$cov
, , bwt_Intercept
Estimate Est.Error Q2.5 Q97.5
bwt_Intercept 1.00982294 0.3583697 0.5044397 1.932823
tarsus_Intercept 0.06412895 0.4880116 -0.8559486 1.092317
, , tarsus_Intercept
Estimate Est.Error Q2.5 Q97.5
bwt_Intercept 0.06412895 0.4880116 -0.8559486 1.092317
tarsus_Intercept 4.18025049 1.4249595 2.0920330 7.711755
$mother
$mother$sd
Estimate Est.Error Q2.5 Q97.5
bwt_Intercept 1.137249 0.1175495 0.8968349 1.358078
tarsus_Intercept 2.088602 0.2865291 1.5425683 2.673791
$mother$cor
, , bwt_Intercept
Estimate Est.Error Q2.5 Q97.5
bwt_Intercept 1.0000000 0.0000000 1.0000000 1.0000000
tarsus_Intercept -0.6413329 0.1968102 -0.9745925 -0.2367393
, , tarsus_Intercept
Estimate Est.Error Q2.5 Q97.5
bwt_Intercept -0.6413329 0.1968102 -0.9745925 -0.2367393
tarsus_Intercept 1.0000000 0.0000000 1.0000000 1.0000000
$mother$cov
, , bwt_Intercept
Estimate Est.Error Q2.5 Q97.5
bwt_Intercept 1.307139 0.2669748 0.8043132 1.8443762
tarsus_Intercept -1.467511 0.3650113 -2.1491979 -0.6718737
, , tarsus_Intercept
Estimate Est.Error Q2.5 Q97.5
bwt_Intercept -1.467511 0.3650113 -2.149198 -0.6718737
tarsus_Intercept 4.444274 1.2018090 2.379517 7.1491575
$residual__
$residual__$sd
Estimate Est.Error Q2.5 Q97.5
bwt 1.396815 0.1597393 1.052042 1.684810
tarsus 3.733614 0.3170549 3.140204 4.319446
$residual__$cor
, , bwt
Estimate Est.Error Q2.5 Q97.5
bwt 1.0000000 0.0000000 1.0000000 1.0000000
tarsus 0.8925119 0.0672286 0.7167862 0.9774545
, , tarsus
Estimate Est.Error Q2.5 Q97.5
bwt 0.8925119 0.0672286 0.7167862 0.9774545
tarsus 1.0000000 0.0000000 1.0000000 1.0000000
$residual__$cov
, , bwt
Estimate Est.Error Q2.5 Q97.5
bwt 1.976584 0.4347085 1.106793 2.838585
tarsus 4.683451 0.9103808 2.846560 6.440550
, , tarsus
Estimate Est.Error Q2.5 Q97.5
bwt 4.683451 0.9103808 2.846560 6.44055
tarsus 14.040300 2.3674563 9.860883 18.65761Evaluation of the statistical support for these genetic and maternal correlations is straightforward. Because we imposed no constraint on their estimation, we can evaluate the extent to which the posterior distributions overlap zero:
cor_g <- as.mcmc((VarCorr(brms_m2.2, summary = FALSE)$animal$cor[, 1, 2]))
cor_res <- as.mcmc((VarCorr(brms_m2.2, summary = FALSE)$residual$cor[, 1, 2]))
cor_mother <- as.mcmc((VarCorr(brms_m2.2, summary = FALSE)$mother$cor[, 1, 2]))
cor_byear <- as.mcmc((VarCorr(brms_m2.2, summary = FALSE)$byear$cor[, 1, 2]))
summary(cor_g)
Iterations = 1:1000
Thinning interval = 1
Number of chains = 1
Sample size per chain = 1000
1. Empirical mean and standard deviation for each variable,
plus standard error of the mean:
Mean SD Naive SE Time-series SE
0.596151 0.172794 0.005464 0.028178
2. Quantiles for each variable:
2.5% 25% 50% 75% 97.5%
0.1850 0.5065 0.6117 0.7064 0.8930 summary(cor_mother)
Iterations = 1:1000
Thinning interval = 1
Number of chains = 1
Sample size per chain = 1000
1. Empirical mean and standard deviation for each variable,
plus standard error of the mean:
Mean SD Naive SE Time-series SE
-0.641333 0.196810 0.006224 0.021355
2. Quantiles for each variable:
2.5% 25% 50% 75% 97.5%
-0.9746 -0.7933 -0.6507 -0.5027 -0.2367 summary(cor_byear)
Iterations = 1:1000
Thinning interval = 1
Number of chains = 1
Sample size per chain = 1000
1. Empirical mean and standard deviation for each variable,
plus standard error of the mean:
Mean SD Naive SE Time-series SE
0.013514 0.220619 0.006977 0.009633
2. Quantiles for each variable:
2.5% 25% 50% 75% 97.5%
-0.43678 -0.13505 0.02408 0.16947 0.44123 summary(cor_res)
Iterations = 1:1000
Thinning interval = 1
Number of chains = 1
Sample size per chain = 1000
1. Empirical mean and standard deviation for each variable,
plus standard error of the mean:
Mean SD Naive SE Time-series SE
0.892512 0.067229 0.002126 0.026496
2. Quantiles for each variable:
2.5% 25% 50% 75% 97.5%
0.7168 0.8620 0.9069 0.9425 0.9775 plot(cor_g)
plot(cor_res)
plot(cor_mother)
plot(cor_byear)
Neither or these posterior distributions overlaps zero, so we can consider them both statistically supported.
cor.est <- rbind(
cbind(summary(cor_g)$statistics[1], summary(cor_g)$quantiles[1], summary(cor_g)$quantiles[5]),
cbind(summary(cor_mother)$statistics[1], summary(cor_mother)$quantiles[1], summary(cor_mother)$quantiles[5]),
cbind(summary(cor_byear)$statistics[1], summary(cor_byear)$quantiles[1], summary(cor_byear)$quantiles[5]),
cbind(summary(cor_res)$statistics[1], summary(cor_res)$quantiles[1], summary(cor_res)$quantiles[5])
)
plot(c(1, 2, 3, 4) ~ cor.est[, 1], xlim = c(-1, 1), ylim = c(0, 5), xlab = "", ylab = "", cex = 2, yaxt = "n")
segments(y0 = 1, x0 = cor.est[1, 1] - cor.est[1, 2], y1 = 1, x1 = cor.est[1, 1] + cor.est[1, 2], lwd = 2)
segments(y0 = 2, x0 = cor.est[2, 1] - cor.est[2, 2], y1 = 2, x1 = cor.est[2, 1] + cor.est[2, 2], lwd = 2)
segments(y0 = 3, x0 = cor.est[3, 1] - cor.est[3, 2], y1 = 3, x1 = cor.est[3, 1] + cor.est[3, 2], lwd = 2)
segments(y0 = 4, x0 = cor.est[4, 1] - cor.est[4, 2], y1 = 4, x1 = cor.est[4, 1] + cor.est[4, 2], lwd = 2)
mtext("Correlation (±CI)", side = 1, las = 1, adj = 0.4, line = 3, cex = 1.6)
axis(2, at = 1, labels = c("genetic"), las = 2, cex.axis = 1)
axis(2, at = 2, labels = c("mother"), las = 2, cex.axis = 1)
axis(2, at = 3, labels = c("year"), las = 2, cex.axis = 1)
axis(2, at = 4, labels = c("residual"), las = 2, cex.axis = 1)
Note, brms estimates the correlation and also the covariance. We can also recalculate the correlation directly from the covariance.To facilitate the extraction of the different parameter, we can the fucntion as_draws_df
cov_g <- (VarCorr(brms_m2.2, summary = FALSE)$animal$cov)[, 1, 2]
cov_res <- (VarCorr(brms_m2.2, summary = FALSE)$residual$cov)[, 1, 2]
cov_mother <- (VarCorr(brms_m2.2, summary = FALSE)$mother$cov)[, 1, 2]
cov_byear <- (VarCorr(brms_m2.2, summary = FALSE)$byear$cov)[, 1, 2]
var.est <- as_draws_df(brms_m2.2, variable = c("sd", "sigma"), regex = TRUE)
var.est <- var.est^2
cor_g_2 <- as.mcmc(cov_g / sqrt(var.est[1] * var.est[2]))
cor_byear_2 <- as.mcmc(cov_byear / sqrt(var.est[3] * var.est[4]))
cor_mother_2 <- as.mcmc(cov_g / sqrt(var.est[5] * var.est[6]))
cor_res_2 <- as.mcmc(cov_res / sqrt(var.est[7] * var.est[8]))
summary(cor_g_2)
Iterations = 1:1000
Thinning interval = 1
Number of chains = 1
Sample size per chain = 1000
1. Empirical mean and standard deviation for each variable,
plus standard error of the mean:
Mean SD Naive SE Time-series SE
0.596151 0.172794 0.005464 0.028178
2. Quantiles for each variable:
2.5% 25% 50% 75% 97.5%
0.1850 0.5065 0.6117 0.7064 0.8930 summary(cor_byear_2)
Iterations = 1:1000
Thinning interval = 1
Number of chains = 1
Sample size per chain = 1000
1. Empirical mean and standard deviation for each variable,
plus standard error of the mean:
Mean SD Naive SE Time-series SE
0.013514 0.220619 0.006977 0.009633
2. Quantiles for each variable:
2.5% 25% 50% 75% 97.5%
-0.43678 -0.13505 0.02408 0.16947 0.44123 summary(cor_mother_2)
Iterations = 1:1000
Thinning interval = 1
Number of chains = 1
Sample size per chain = 1000
1. Empirical mean and standard deviation for each variable,
plus standard error of the mean:
Mean SD Naive SE Time-series SE
1.01862 0.46337 0.01465 0.06126
2. Quantiles for each variable:
2.5% 25% 50% 75% 97.5%
0.1514 0.6969 1.0310 1.3241 2.0032 summary(cor_res_2)
Iterations = 1:1000
Thinning interval = 1
Number of chains = 1
Sample size per chain = 1000
1. Empirical mean and standard deviation for each variable,
plus standard error of the mean:
Mean SD Naive SE Time-series SE
0.892512 0.067229 0.002126 0.026496
2. Quantiles for each variable:
2.5% 25% 50% 75% 97.5%
0.7168 0.8620 0.9069 0.9425 0.9775 plot(cor_g_2)
plot(cor_byear_2)
plot(cor_mother_2)
plot(cor_res_2)
As in the tutorial 1, it is possible to partition the variance-covariance matrix between groups (here sex)
bf_bwt_3 <- bf(bwt ~ 1 + sex + ((1 | a | gr(animal, cov = Amat, by = sex))) + (1 | b | byear) + (1 | c | mother))
bf_tarsus_3 <- bf(tarsus ~ 1 + sex + (1 | a | gr(animal, cov = Amat, by = sex)) + (1 | b | byear) + (1 | c | mother))
brms_m2.3 <- brm(
bf_bwt_3 + bf_tarsus_3 + set_rescor(TRUE),
data = gryphon,
data2 = list(Amat = Amat),
chains = 2, cores = 2, iter = 1000
)
save(brms_m2.3, file = "r-obj/brms_m2_3.rda")Again we have provided the data from one such run. It can be accessed using the code:
Warning: Parts of the model have not converged (some Rhats are > 1.05). Be
careful when analysing the results! We recommend running more iterations and/or
setting stronger priors.Warning: There were 6 divergent transitions after warmup. Increasing
adapt_delta above 0.8 may help. See
http://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup Family: MV(gaussian, gaussian)
Links: mu = identity; sigma = identity
mu = identity; sigma = identity
Formula: bwt ~ 1 + sex + ((1 | a | gr(animal, cov = Amat, by = sex))) + (1 | b | byear) + (1 | c | mother)
tarsus ~ 1 + sex + (1 | a | gr(animal, cov = Amat, by = sex)) + (1 | b | byear) + (1 | c | mother)
Data: gryphon (Number of observations: 683)
Draws: 2 chains, each with iter = 1000; warmup = 500; thin = 1;
total post-warmup draws = 1000
Multilevel Hyperparameters:
~animal (Number of levels: 683)
Estimate Est.Error l-95% CI
sd(bwt_Intercept:sex1) 1.05 0.28 0.48
sd(tarsus_Intercept:sex1) 1.57 0.78 0.14
sd(bwt_Intercept:sex2) 1.23 0.24 0.71
sd(tarsus_Intercept:sex2) 3.24 0.48 2.17
cor(bwt_Intercept:sex1,tarsus_Intercept:sex1) 0.32 0.43 -0.78
cor(bwt_Intercept:sex2,tarsus_Intercept:sex2) 0.70 0.12 0.45
u-95% CI Rhat Bulk_ESS Tail_ESS
sd(bwt_Intercept:sex1) 1.58 1.27 6 39
sd(tarsus_Intercept:sex1) 3.03 1.08 24 117
sd(bwt_Intercept:sex2) 1.65 1.24 7 46
sd(tarsus_Intercept:sex2) 4.13 1.05 37 79
cor(bwt_Intercept:sex1,tarsus_Intercept:sex1) 0.90 1.07 25 202
cor(bwt_Intercept:sex2,tarsus_Intercept:sex2) 0.90 1.09 28 44
~byear (Number of levels: 34)
Estimate Est.Error l-95% CI u-95% CI Rhat
sd(bwt_Intercept) 0.97 0.14 0.70 1.26 1.01
sd(tarsus_Intercept) 2.03 0.34 1.47 2.80 1.00
cor(bwt_Intercept,tarsus_Intercept) 0.01 0.21 -0.41 0.41 1.01
Bulk_ESS Tail_ESS
sd(bwt_Intercept) 282 292
sd(tarsus_Intercept) 324 361
cor(bwt_Intercept,tarsus_Intercept) 183 393
~mother (Number of levels: 352)
Estimate Est.Error l-95% CI u-95% CI Rhat
sd(bwt_Intercept) 1.18 0.11 0.96 1.39 1.02
sd(tarsus_Intercept) 2.04 0.34 1.33 2.65 1.09
cor(bwt_Intercept,tarsus_Intercept) -0.63 0.21 -0.97 -0.23 1.08
Bulk_ESS Tail_ESS
sd(bwt_Intercept) 170 352
sd(tarsus_Intercept) 24 60
cor(bwt_Intercept,tarsus_Intercept) 22 100
Regression Coefficients:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
bwt_Intercept 6.27 0.22 5.82 6.68 1.00 195 123
tarsus_Intercept 20.35 0.48 19.35 21.21 1.00 494 578
bwt_sex2 2.04 0.17 1.71 2.36 1.01 265 384
tarsus_sex2 0.14 0.39 -0.60 0.88 1.00 691 621
Further Distributional Parameters:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sigma_bwt 1.49 0.15 1.16 1.74 1.33 5 81
sigma_tarsus 3.88 0.27 3.27 4.35 1.04 42 85
Residual Correlations:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
rescor(bwt,tarsus) 0.87 0.05 0.78 0.94 1.04 22 94
Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).plot(brms_m2.3, ask = FALSE)
VarCorr(brms_m2.3)$animal
$animal$sd
Estimate Est.Error Q2.5 Q97.5
bwt_Intercept:sex1 1.051791 0.2757021 0.4836926 1.575837
tarsus_Intercept:sex1 1.573820 0.7777373 0.1389929 3.025426
bwt_Intercept:sex2 1.232722 0.2436415 0.7075357 1.650628
tarsus_Intercept:sex2 3.237363 0.4773725 2.1667846 4.126560
$animal$cor
, , bwt_Intercept:sex1
Estimate Est.Error Q2.5 Q97.5
bwt_Intercept:sex1 1.0000000 0.0000000 1.0000000 1.0000000
tarsus_Intercept:sex1 0.3232502 0.4324034 -0.7811983 0.9007145
bwt_Intercept:sex2 0.0000000 0.0000000 0.0000000 0.0000000
tarsus_Intercept:sex2 0.0000000 0.0000000 0.0000000 0.0000000
, , tarsus_Intercept:sex1
Estimate Est.Error Q2.5 Q97.5
bwt_Intercept:sex1 0.3232502 0.4324034 -0.7811983 0.9007145
tarsus_Intercept:sex1 1.0000000 0.0000000 1.0000000 1.0000000
bwt_Intercept:sex2 0.0000000 0.0000000 0.0000000 0.0000000
tarsus_Intercept:sex2 0.0000000 0.0000000 0.0000000 0.0000000
, , bwt_Intercept:sex2
Estimate Est.Error Q2.5 Q97.5
bwt_Intercept:sex1 0.0000000 0.0000000 0.0000000 0.0000000
tarsus_Intercept:sex1 0.0000000 0.0000000 0.0000000 0.0000000
bwt_Intercept:sex2 1.0000000 0.0000000 1.0000000 1.0000000
tarsus_Intercept:sex2 0.6975925 0.1213692 0.4474181 0.9047623
, , tarsus_Intercept:sex2
Estimate Est.Error Q2.5 Q97.5
bwt_Intercept:sex1 0.0000000 0.0000000 0.0000000 0.0000000
tarsus_Intercept:sex1 0.0000000 0.0000000 0.0000000 0.0000000
bwt_Intercept:sex2 0.6975925 0.1213692 0.4474181 0.9047623
tarsus_Intercept:sex2 1.0000000 0.0000000 1.0000000 1.0000000
$animal$cov
, , bwt_Intercept:sex1
Estimate Est.Error Q2.5 Q97.5
bwt_Intercept:sex1 1.1821990 0.5768155 0.2339592 2.483263
tarsus_Intercept:sex1 0.7957577 0.9482719 -0.4879668 3.035395
bwt_Intercept:sex2 0.0000000 0.0000000 0.0000000 0.000000
tarsus_Intercept:sex2 0.0000000 0.0000000 0.0000000 0.000000
, , tarsus_Intercept:sex1
Estimate Est.Error Q2.5 Q97.5
bwt_Intercept:sex1 0.7957577 0.9482719 -0.48796685 3.035395
tarsus_Intercept:sex1 3.0811806 2.5183190 0.01931907 9.153208
bwt_Intercept:sex2 0.0000000 0.0000000 0.00000000 0.000000
tarsus_Intercept:sex2 0.0000000 0.0000000 0.00000000 0.000000
, , bwt_Intercept:sex2
Estimate Est.Error Q2.5 Q97.5
bwt_Intercept:sex1 0.000000 0.0000000 0.0000000 0.000000
tarsus_Intercept:sex1 0.000000 0.0000000 0.0000000 0.000000
bwt_Intercept:sex2 1.578907 0.5865683 0.5006069 2.724572
tarsus_Intercept:sex2 2.842372 1.0153859 1.1593760 5.086448
, , tarsus_Intercept:sex2
Estimate Est.Error Q2.5 Q97.5
bwt_Intercept:sex1 0.000000 0.000000 0.000000 0.000000
tarsus_Intercept:sex1 0.000000 0.000000 0.000000 0.000000
bwt_Intercept:sex2 2.842372 1.015386 1.159376 5.086448
tarsus_Intercept:sex2 10.708178 3.017245 4.694964 17.028497
$byear
$byear$sd
Estimate Est.Error Q2.5 Q97.5
bwt_Intercept 0.9676965 0.1434955 0.7018779 1.258518
tarsus_Intercept 2.0290144 0.3382466 1.4650518 2.801932
$byear$cor
, , bwt_Intercept
Estimate Est.Error Q2.5 Q97.5
bwt_Intercept 1.000000000 0.0000000 1.0000000 1.000000
tarsus_Intercept 0.009103073 0.2077977 -0.4096021 0.410902
, , tarsus_Intercept
Estimate Est.Error Q2.5 Q97.5
bwt_Intercept 0.009103073 0.2077977 -0.4096021 0.410902
tarsus_Intercept 1.000000000 0.0000000 1.0000000 1.000000
$byear$cov
, , bwt_Intercept
Estimate Est.Error Q2.5 Q97.5
bwt_Intercept 0.95700691 0.2875804 0.4926327 1.583869
tarsus_Intercept 0.04233908 0.4457863 -0.8475401 1.014920
, , tarsus_Intercept
Estimate Est.Error Q2.5 Q97.5
bwt_Intercept 0.04233908 0.4457863 -0.8475401 1.014920
tarsus_Intercept 4.23119588 1.4453420 2.1463767 7.850826
$mother
$mother$sd
Estimate Est.Error Q2.5 Q97.5
bwt_Intercept 1.175168 0.1083817 0.958963 1.388926
tarsus_Intercept 2.038471 0.3445355 1.327549 2.648997
$mother$cor
, , bwt_Intercept
Estimate Est.Error Q2.5 Q97.5
bwt_Intercept 1.0000000 0.0000000 1.0000000 1.0000000
tarsus_Intercept -0.6301934 0.2113836 -0.9652044 -0.2283133
, , tarsus_Intercept
Estimate Est.Error Q2.5 Q97.5
bwt_Intercept -0.6301934 0.2113836 -0.9652044 -0.2283133
tarsus_Intercept 1.0000000 0.0000000 1.0000000 1.0000000
$mother$cov
, , bwt_Intercept
Estimate Est.Error Q2.5 Q97.5
bwt_Intercept 1.392755 0.2562907 0.9196105 1.929117
tarsus_Intercept -1.434727 0.3664843 -2.0769067 -0.674943
, , tarsus_Intercept
Estimate Est.Error Q2.5 Q97.5
bwt_Intercept -1.434727 0.3664843 -2.076907 -0.674943
tarsus_Intercept 4.273951 1.3914013 1.762387 7.017188
$residual__
$residual__$sd
Estimate Est.Error Q2.5 Q97.5
bwt 1.489342 0.1489641 1.162192 1.741559
tarsus 3.875633 0.2650604 3.268802 4.345400
$residual__$cor
, , bwt
Estimate Est.Error Q2.5 Q97.5
bwt 1.0000000 0.00000000 1.0000000 1.0000000
tarsus 0.8685184 0.04534008 0.7772864 0.9414488
, , tarsus
Estimate Est.Error Q2.5 Q97.5
bwt 0.8685184 0.04534008 0.7772864 0.9414488
tarsus 1.0000000 0.00000000 1.0000000 1.0000000
$residual__$cov
, , bwt
Estimate Est.Error Q2.5 Q97.5
bwt 2.240307 0.4341334 1.350691 3.033029
tarsus 5.034739 0.7852137 3.272088 6.400539
, , tarsus
Estimate Est.Error Q2.5 Q97.5
bwt 5.034739 0.7852137 3.272088 6.400539
tarsus 15.090721 2.0309529 10.685067 18.882505However, this model is lacking an important and essentiel group-specific partitionning (we do with the asreml-R and MCMCglmm). We need to partition the residual variance (or sigma) as well. Doing so, we will use the argument `sigma`` to parititon the model by sex. To avoid an estimation of the difference between sexes, we need to remove the estimate of the intercept at the sigma level.
bf_bwt_4 <- bf(bwt ~ 1 + sex + ((1 | a | gr(animal, cov = Amat, by = sex))) + (1 | b | byear) + (1 | c | mother), sigma ~ sex - 1)
bf_tarsus_4 <- bf(tarsus ~ 1 + sex + (1 | a | gr(animal, cov = Amat, by = sex)) + (1 | b | byear) + (1 | c | mother), sigma ~ sex - 1)
brms_m2.4 <- brm(
bf_bwt_4 + bf_tarsus_4 + set_rescor(TRUE),
data = gryphon,
data2 = list(Amat = Amat),
chains = 2, cores = 2, iter = 1000
)
save(brms_m2.4, file = "r-obj/brms_m2_4.rda")Again we have provided the data from one such run. It can be accessed using the code:
Warning: Parts of the model have not converged (some Rhats are > 1.05). Be
careful when analysing the results! We recommend running more iterations and/or
setting stronger priors.Warning: There were 6 divergent transitions after warmup. Increasing
adapt_delta above 0.8 may help. See
http://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup Family: MV(gaussian, gaussian)
Links: mu = identity; sigma = log
mu = identity; sigma = log
Formula: bwt ~ 1 + sex + ((1 | a | gr(animal, cov = Amat, by = sex))) + (1 | b | byear) + (1 | c | mother)
sigma ~ sex - 1
tarsus ~ 1 + sex + (1 | a | gr(animal, cov = Amat, by = sex)) + (1 | b | byear) + (1 | c | mother)
sigma ~ sex - 1
Data: gryphon (Number of observations: 683)
Draws: 2 chains, each with iter = 1000; warmup = 500; thin = 1;
total post-warmup draws = 1000
Multilevel Hyperparameters:
~animal (Number of levels: 683)
Estimate Est.Error l-95% CI
sd(bwt_Intercept:sex1) 0.86 0.35 0.14
sd(tarsus_Intercept:sex1) 1.40 0.82 0.10
sd(bwt_Intercept:sex2) 1.40 0.23 0.94
sd(tarsus_Intercept:sex2) 3.49 0.72 1.96
cor(bwt_Intercept:sex1,tarsus_Intercept:sex1) 0.28 0.51 -0.85
cor(bwt_Intercept:sex2,tarsus_Intercept:sex2) 0.78 0.10 0.56
u-95% CI Rhat Bulk_ESS Tail_ESS
sd(bwt_Intercept:sex1) 1.48 1.00 59 101
sd(tarsus_Intercept:sex1) 3.00 1.07 45 111
sd(bwt_Intercept:sex2) 1.83 1.04 42 40
sd(tarsus_Intercept:sex2) 4.62 1.22 8 50
cor(bwt_Intercept:sex1,tarsus_Intercept:sex1) 0.94 1.05 60 287
cor(bwt_Intercept:sex2,tarsus_Intercept:sex2) 0.97 1.09 30 94
~byear (Number of levels: 34)
Estimate Est.Error l-95% CI u-95% CI Rhat
sd(bwt_Intercept) 0.97 0.15 0.73 1.30 1.01
sd(tarsus_Intercept) 2.01 0.33 1.41 2.70 1.00
cor(bwt_Intercept,tarsus_Intercept) -0.00 0.22 -0.42 0.45 1.02
Bulk_ESS Tail_ESS
sd(bwt_Intercept) 283 412
sd(tarsus_Intercept) 349 554
cor(bwt_Intercept,tarsus_Intercept) 225 256
~mother (Number of levels: 352)
Estimate Est.Error l-95% CI u-95% CI Rhat
sd(bwt_Intercept) 1.19 0.11 0.98 1.42 1.00
sd(tarsus_Intercept) 2.14 0.30 1.57 2.71 1.05
cor(bwt_Intercept,tarsus_Intercept) -0.55 0.21 -0.96 -0.16 1.04
Bulk_ESS Tail_ESS
sd(bwt_Intercept) 279 434
sd(tarsus_Intercept) 46 227
cor(bwt_Intercept,tarsus_Intercept) 46 122
Regression Coefficients:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
bwt_Intercept 6.27 0.23 5.82 6.73 1.01 336 556
tarsus_Intercept 20.39 0.49 19.42 21.37 1.00 384 509
bwt_sex2 2.04 0.17 1.70 2.38 1.00 483 507
sigma_bwt_sex1 0.45 0.12 0.18 0.63 1.00 68 128
sigma_bwt_sex2 0.26 0.17 -0.17 0.52 1.06 38 33
tarsus_sex2 0.10 0.41 -0.69 0.92 1.00 658 659
sigma_tarsus_sex1 1.37 0.08 1.19 1.50 1.04 65 215
sigma_tarsus_sex2 1.24 0.16 0.91 1.50 1.22 7 67
Residual Correlations:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
rescor(bwt,tarsus) 0.84 0.05 0.73 0.93 1.15 11 29
Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).plot(brms_m2.4, ask = FALSE)
VarCorr(brms_m2.4)$animal
$animal$sd
Estimate Est.Error Q2.5 Q97.5
bwt_Intercept:sex1 0.8628389 0.3531696 0.14360929 1.481121
tarsus_Intercept:sex1 1.3985239 0.8163731 0.09977045 3.002609
bwt_Intercept:sex2 1.4023202 0.2325863 0.93896619 1.833688
tarsus_Intercept:sex2 3.4858243 0.7167230 1.96177030 4.620105
$animal$cor
, , bwt_Intercept:sex1
Estimate Est.Error Q2.5 Q97.5
bwt_Intercept:sex1 1.000000 0.0000000 1.0000000 1.0000000
tarsus_Intercept:sex1 0.277338 0.5119501 -0.8479996 0.9398158
bwt_Intercept:sex2 0.000000 0.0000000 0.0000000 0.0000000
tarsus_Intercept:sex2 0.000000 0.0000000 0.0000000 0.0000000
, , tarsus_Intercept:sex1
Estimate Est.Error Q2.5 Q97.5
bwt_Intercept:sex1 0.277338 0.5119501 -0.8479996 0.9398158
tarsus_Intercept:sex1 1.000000 0.0000000 1.0000000 1.0000000
bwt_Intercept:sex2 0.000000 0.0000000 0.0000000 0.0000000
tarsus_Intercept:sex2 0.000000 0.0000000 0.0000000 0.0000000
, , bwt_Intercept:sex2
Estimate Est.Error Q2.5 Q97.5
bwt_Intercept:sex1 0.0000000 0.00000000 0.000000 0.0000000
tarsus_Intercept:sex1 0.0000000 0.00000000 0.000000 0.0000000
bwt_Intercept:sex2 1.0000000 0.00000000 1.000000 1.0000000
tarsus_Intercept:sex2 0.7781659 0.09850178 0.564036 0.9699765
, , tarsus_Intercept:sex2
Estimate Est.Error Q2.5 Q97.5
bwt_Intercept:sex1 0.0000000 0.00000000 0.000000 0.0000000
tarsus_Intercept:sex1 0.0000000 0.00000000 0.000000 0.0000000
bwt_Intercept:sex2 0.7781659 0.09850178 0.564036 0.9699765
tarsus_Intercept:sex2 1.0000000 0.00000000 1.000000 1.0000000
$animal$cov
, , bwt_Intercept:sex1
Estimate Est.Error Q2.5 Q97.5
bwt_Intercept:sex1 0.8690949 0.6002446 0.0206250 2.193721
tarsus_Intercept:sex1 0.6483841 0.9143700 -0.4481203 2.989048
bwt_Intercept:sex2 0.0000000 0.0000000 0.0000000 0.000000
tarsus_Intercept:sex2 0.0000000 0.0000000 0.0000000 0.000000
, , tarsus_Intercept:sex1
Estimate Est.Error Q2.5 Q97.5
bwt_Intercept:sex1 0.6483841 0.914370 -0.448120294 2.989048
tarsus_Intercept:sex1 2.6216677 2.500989 0.009954235 9.015666
bwt_Intercept:sex2 0.0000000 0.000000 0.000000000 0.000000
tarsus_Intercept:sex2 0.0000000 0.000000 0.000000000 0.000000
, , bwt_Intercept:sex2
Estimate Est.Error Q2.5 Q97.5
bwt_Intercept:sex1 0.000000 0.000000 0.0000000 0.000000
tarsus_Intercept:sex1 0.000000 0.000000 0.0000000 0.000000
bwt_Intercept:sex2 2.020544 0.639550 0.8816577 3.362416
tarsus_Intercept:sex2 3.875299 1.298562 1.4247624 6.415724
, , tarsus_Intercept:sex2
Estimate Est.Error Q2.5 Q97.5
bwt_Intercept:sex1 0.000000 0.000000 0.000000 0.000000
tarsus_Intercept:sex1 0.000000 0.000000 0.000000 0.000000
bwt_Intercept:sex2 3.875299 1.298562 1.424762 6.415724
tarsus_Intercept:sex2 12.664149 4.814957 3.848544 21.345370
$byear
$byear$sd
Estimate Est.Error Q2.5 Q97.5
bwt_Intercept 0.9707654 0.1478361 0.7256973 1.298668
tarsus_Intercept 2.0073203 0.3290043 1.4102530 2.699770
$byear$cor
, , bwt_Intercept
Estimate Est.Error Q2.5 Q97.5
bwt_Intercept 1.000000000 0.0000000 1.0000000 1.0000000
tarsus_Intercept -0.001551923 0.2236193 -0.4237849 0.4526618
, , tarsus_Intercept
Estimate Est.Error Q2.5 Q97.5
bwt_Intercept -0.001551923 0.2236193 -0.4237849 0.4526618
tarsus_Intercept 1.000000000 0.0000000 1.0000000 1.0000000
$byear$cov
, , bwt_Intercept
Estimate Est.Error Q2.5 Q97.5
bwt_Intercept 0.9642191 0.3021042 0.5266366 1.686538
tarsus_Intercept 0.0252866 0.4713288 -0.8080314 1.069567
, , tarsus_Intercept
Estimate Est.Error Q2.5 Q97.5
bwt_Intercept 0.0252866 0.4713288 -0.8080314 1.069567
tarsus_Intercept 4.1374703 1.3676889 1.9888135 7.288761
$mother
$mother$sd
Estimate Est.Error Q2.5 Q97.5
bwt_Intercept 1.189902 0.1113270 0.9791374 1.423333
tarsus_Intercept 2.139290 0.2985875 1.5714814 2.708199
$mother$cor
, , bwt_Intercept
Estimate Est.Error Q2.5 Q97.5
bwt_Intercept 1.0000000 0.0000000 1.0000000 1.0000000
tarsus_Intercept -0.5501934 0.2066985 -0.9589938 -0.1591737
, , tarsus_Intercept
Estimate Est.Error Q2.5 Q97.5
bwt_Intercept -0.5501934 0.2066985 -0.9589938 -0.1591737
tarsus_Intercept 1.0000000 0.0000000 1.0000000 1.0000000
$mother$cov
, , bwt_Intercept
Estimate Est.Error Q2.5 Q97.5
bwt_Intercept 1.428247 0.2672833 0.958710 2.0258758
tarsus_Intercept -1.335450 0.3932707 -2.051866 -0.5146954
, , tarsus_Intercept
Estimate Est.Error Q2.5 Q97.5
bwt_Intercept -1.335450 0.3932707 -2.051866 -0.5146954
tarsus_Intercept 4.665626 1.2688941 2.469554 7.3343400Evaluation of the statistical support for these sex-specific correlations is straightforward. Because we imposed no constraint on their estimation, we can evaluate the extent to which the posterior distributions overlap zero or overlap each other:
cor_g_F <- as.mcmc((VarCorr(brms_m2.4, summary = FALSE)$animal$cor[, 1, 2]))
cor_g_M <- as.mcmc((VarCorr(brms_m2.4, summary = FALSE)$animal$cor[, 3, 4]))
summary(cor_g_F)
Iterations = 1:1000
Thinning interval = 1
Number of chains = 1
Sample size per chain = 1000
1. Empirical mean and standard deviation for each variable,
plus standard error of the mean:
Mean SD Naive SE Time-series SE
0.27734 0.51195 0.01619 0.05981
2. Quantiles for each variable:
2.5% 25% 50% 75% 97.5%
-0.84800 -0.02476 0.43226 0.67124 0.93982 summary(cor_g_M)
Iterations = 1:1000
Thinning interval = 1
Number of chains = 1
Sample size per chain = 1000
1. Empirical mean and standard deviation for each variable,
plus standard error of the mean:
Mean SD Naive SE Time-series SE
0.778166 0.098502 0.003115 0.013827
2. Quantiles for each variable:
2.5% 25% 50% 75% 97.5%
0.5640 0.7240 0.7805 0.8424 0.9700 plot(cor_g_F)
plot(cor_g_M)
Here a plot to visualize the overlaps of covariances.
cor.est <- rbind(
cbind(summary(cor_g_F)$statistics[1], summary(cor_g_F)$quantiles[1], summary(cor_g_F)$quantiles[5]),
cbind(summary(cor_g_M)$statistics[1], summary(cor_g_M)$quantiles[1], summary(cor_g_M)$quantiles[5])
)
plot(c(1, 2) ~ cor.est[, 1], xlim = c(0, 1.5), ylim = c(0, 2.5), xlab = "", ylab = "", col = c("red", "blue"), pch = c(16, 17), cex = 2, yaxt = "n")
segments(y0 = 1, x0 = cor.est[1, 2], y1 = 1, x1 = cor.est[1, 3], col = c("red"), lwd = 2)
segments(y0 = 2, x0 = cor.est[2, 2], y1 = 2, x1 = cor.est[2, 3], col = c("blue"), lwd = 2)
mtext("Covariance (±CI)", side = 1, las = 1, adj = 0.4, line = 3, cex = 1.6)
axis(2, at = 1, labels = c("female"), las = 3, cex.axis = 1.6)
axis(2, at = 2, labels = c("male"), las = 3, cex.axis = 1.6)
Here a simple plot of the sex-specific genetic correlation using the BLUPs
bls_m2.4 <- ranef(brms_m2.4)$animal
bl_m2.4 <- as.data.frame(abind::abind(lapply(1:dim(bls_m2.4)[3], function(x) bls_m2.4[, c(1, 3, 4), x])))
colnames(bl_m2.4) <- paste0(rep(dimnames(bls_m2.4)[[3]], each = 3), c("", "_lo", "_up"))
bl_m2.4$id <- rownames(bl_m2.4)
bl_m2.4$sex <- attr(dimnames(bls_m2.4)[[1]], "by")
FEM <- subset(bl_m2.4, sex == "1")
MAL <- subset(bl_m2.4, sex == "2")#
par(mfrow = c(1, 2))
plot(tarsus_Intercept ~ bwt_Intercept, FEM,
xlab = "", ylab = "",
xlim = c(min(FEM$bwt_Intercept_lo), max(FEM$bwt_Intercept_up)),
ylim = c(min(FEM$tarsus_Intercept_lo), max(FEM$tarsus_Intercept_up)),
las = 1.2, type = "n"
)
segments(
x0 = FEM$bwt_Intercept, y0 = FEM$tarsus_Intercept_lo,
x1 = FEM$bwt_Intercept, y1 = FEM$tarsus_Intercept_up,
col = "black"
)
segments(
x0 = FEM$bwt_Intercept_lo, y0 = FEM$tarsus_Intercept,
x1 = FEM$bwt_Intercept_up, y1 = FEM$tarsus_Intercept,
col = "black"
)
points(tarsus_Intercept ~ bwt_Intercept, FEM, pch = 16, col = "red", cex = 1.5)
points(tarsus_Intercept ~ bwt_Intercept, FEM, pch = 1, col = rgb(0, 0, 0, 0.3), cex = c(1.5))
mtext("btw (BV±CI)", side = 1, line = 2.4)
mtext("tarsus (BV±CI)", side = 2, line = 2, las = 3)
#
plot(tarsus_Intercept ~ bwt_Intercept, MAL,
xlab = "", ylab = "",
xlim = c(min(MAL$bwt_Intercept_lo), max(MAL$bwt_Intercept_up)),
ylim = c(min(MAL$tarsus_Intercept_lo), max(MAL$tarsus_Intercept_up)),
las = 1.2, type = "n"
)
segments(
x0 = MAL$bwt_Intercept, y0 = MAL$tarsus_Intercept_lo,
x1 = MAL$bwt_Intercept, y1 = MAL$tarsus_Intercept_up, col = "black"
)
segments(
x0 = MAL$bwt_Intercept_lo, y0 = MAL$tarsus_Intercept,
x1 = MAL$bwt_Intercept_up, y1 = MAL$tarsus_Intercept, col = "black"
)
points(tarsus_Intercept ~ bwt_Intercept, MAL, pch = 17, col = "blue", cex = 1.5)
points(tarsus_Intercept ~ bwt_Intercept, MAL, pch = 1, col = rgb(0, 0, 0, 0.3), cex = c(1.5))
mtext("btw (BV±CI)", side = 1, line = 2.4)
mtext("tarsus (BV±CI)", side = 2, line = 2, las = 3)
Animal models are amazing model. With different group within a population, it is also possible to estimate how much the different groups shared the same genetic via the cross-group genetic covariance. This covariance is essential to understand ontogenic or sexual conflict, which can constraint or enhanced response to evolution. As an example, we estimate the cross-sex genetic correlation r_{fm}
It is important to keep in mind the covariance matrix at the residual level is zero and it is important to avoid estimating the cross-sex residual covariance because no individual switched sex during the experiment.
Note: the way of partitionning variance per sex is a bit different then the previous code “,by=sex”.
This code is faster and also easier to understand. Note, it is possible to play with the | or || to estimate or not covariance between sexes.
bf_bwt_5 <- bf(
bwt ~ 1 + sex + (0 + sex | a | gr(animal, cov = Amat)) + (0 + sex | b | mother) + (0 + sex | c | byear),
sigma ~ sex - 1
)
bf_tarsus_5 <- bf(
tarsus ~ 1 + sex + (0 + sex | a | gr(animal, cov = Amat)) + (0 + sex | b | mother) + (0 + sex | c | byear),
sigma ~ sex - 1
)
brms_m2.5 <- brm(
bf_bwt_5 + bf_tarsus_5 + set_rescor(TRUE),
data = gryphon,
data2 = list(Amat = Amat),
chains = 2, cores = 2, iter = 1000
)
save(brms_m2.5, file = "r-obj/brms_m2_5.rda")Again we have provided the data from one such run. It can be accessed using the code:
Warning: Parts of the model have not converged (some Rhats are > 1.05). Be
careful when analysing the results! We recommend running more iterations and/or
setting stronger priors.Warning: There were 45 divergent transitions after warmup. Increasing
adapt_delta above 0.8 may help. See
http://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup Family: MV(gaussian, gaussian)
Links: mu = identity; sigma = log
mu = identity; sigma = log
Formula: bwt ~ 1 + sex + (0 + sex | a | gr(animal, cov = Amat)) + (0 + sex | b | mother) + (0 + sex | c | byear)
sigma ~ sex - 1
tarsus ~ 1 + sex + (0 + sex | a | gr(animal, cov = Amat)) + (0 + sex | b | mother) + (0 + sex | c | byear)
sigma ~ sex - 1
Data: gryphon (Number of observations: 683)
Draws: 2 chains, each with iter = 1000; warmup = 500; thin = 1;
total post-warmup draws = 1000
Multilevel Hyperparameters:
~animal (Number of levels: 683)
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
sd(bwt_sex1) 1.26 0.30 0.63 1.73 1.06 21
sd(bwt_sex2) 1.08 0.42 0.20 1.77 1.08 18
sd(tarsus_sex1) 2.26 0.72 0.61 3.57 1.04 40
sd(tarsus_sex2) 2.74 1.05 0.61 4.47 1.13 12
cor(bwt_sex1,bwt_sex2) 0.48 0.29 -0.24 0.87 1.02 84
cor(bwt_sex1,tarsus_sex1) 0.57 0.25 -0.07 0.89 1.14 10
cor(bwt_sex2,tarsus_sex1) 0.38 0.38 -0.53 0.91 1.25 7
cor(bwt_sex1,tarsus_sex2) 0.17 0.31 -0.49 0.75 1.04 60
cor(bwt_sex2,tarsus_sex2) 0.52 0.33 -0.37 0.87 1.20 8
cor(tarsus_sex1,tarsus_sex2) 0.44 0.29 -0.30 0.87 1.03 47
Tail_ESS
sd(bwt_sex1) 104
sd(bwt_sex2) 25
sd(tarsus_sex1) 99
sd(tarsus_sex2) 42
cor(bwt_sex1,bwt_sex2) 112
cor(bwt_sex1,tarsus_sex1) 145
cor(bwt_sex2,tarsus_sex1) 67
cor(bwt_sex1,tarsus_sex2) 94
cor(bwt_sex2,tarsus_sex2) 50
cor(tarsus_sex1,tarsus_sex2) 44
~byear (Number of levels: 34)
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
sd(bwt_sex1) 0.80 0.16 0.53 1.16 1.00 394
sd(bwt_sex2) 1.14 0.19 0.81 1.55 1.01 358
sd(tarsus_sex1) 2.23 0.46 1.50 3.18 1.01 297
sd(tarsus_sex2) 2.34 0.49 1.56 3.41 1.01 229
cor(bwt_sex1,bwt_sex2) 0.74 0.15 0.35 0.96 1.01 266
cor(bwt_sex1,tarsus_sex1) -0.11 0.24 -0.55 0.35 1.01 190
cor(bwt_sex2,tarsus_sex1) -0.39 0.20 -0.73 0.00 1.01 410
cor(bwt_sex1,tarsus_sex2) 0.29 0.23 -0.17 0.71 1.00 256
cor(bwt_sex2,tarsus_sex2) 0.29 0.21 -0.16 0.66 1.01 327
cor(tarsus_sex1,tarsus_sex2) 0.52 0.19 0.12 0.84 1.00 285
Tail_ESS
sd(bwt_sex1) 619
sd(bwt_sex2) 653
sd(tarsus_sex1) 559
sd(tarsus_sex2) 239
cor(bwt_sex1,bwt_sex2) 433
cor(bwt_sex1,tarsus_sex1) 603
cor(bwt_sex2,tarsus_sex1) 656
cor(bwt_sex1,tarsus_sex2) 319
cor(bwt_sex2,tarsus_sex2) 474
cor(tarsus_sex1,tarsus_sex2) 600
~mother (Number of levels: 352)
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
sd(bwt_sex1) 1.08 0.15 0.79 1.39 1.01 281
sd(bwt_sex2) 1.33 0.15 1.03 1.62 1.00 233
sd(tarsus_sex1) 2.21 0.40 1.36 2.95 1.01 72
sd(tarsus_sex2) 2.31 0.49 1.38 3.34 1.05 50
cor(bwt_sex1,bwt_sex2) 0.83 0.11 0.57 0.98 1.01 68
cor(bwt_sex1,tarsus_sex1) -0.50 0.24 -0.91 -0.07 1.01 57
cor(bwt_sex2,tarsus_sex1) -0.64 0.17 -0.93 -0.28 1.06 50
cor(bwt_sex1,tarsus_sex2) -0.51 0.22 -0.88 -0.08 1.05 67
cor(bwt_sex2,tarsus_sex2) -0.36 0.26 -0.84 0.11 1.05 54
cor(tarsus_sex1,tarsus_sex2) 0.72 0.16 0.37 0.95 1.01 249
Tail_ESS
sd(bwt_sex1) 518
sd(bwt_sex2) 585
sd(tarsus_sex1) 224
sd(tarsus_sex2) 191
cor(bwt_sex1,bwt_sex2) 296
cor(bwt_sex1,tarsus_sex1) 156
cor(bwt_sex2,tarsus_sex1) 268
cor(bwt_sex1,tarsus_sex2) 124
cor(bwt_sex2,tarsus_sex2) 297
cor(tarsus_sex1,tarsus_sex2) 596
Regression Coefficients:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
bwt_Intercept 6.18 0.22 5.81 6.66 1.01 421 413
tarsus_Intercept 20.24 0.55 19.12 21.25 1.01 497 527
bwt_sex2 2.11 0.24 1.61 2.57 1.01 541 569
sigma_bwt_sex1 0.27 0.20 -0.16 0.56 1.08 18 66
sigma_bwt_sex2 0.31 0.23 -0.35 0.59 1.11 18 22
tarsus_sex2 0.29 0.63 -0.96 1.56 1.00 490 535
sigma_tarsus_sex1 1.26 0.12 1.01 1.47 1.02 50 109
sigma_tarsus_sex2 1.28 0.20 0.80 1.53 1.13 12 42
Residual Correlations:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
rescor(bwt,tarsus) 0.88 0.05 0.71 0.95 1.32 5 33
Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).plot(brms_m2.5, ask = FALSE)
VarCorr(brms_m2.5)$animal
$animal$sd
Estimate Est.Error Q2.5 Q97.5
bwt_sex1 1.256797 0.3027123 0.6293789 1.733715
bwt_sex2 1.077842 0.4168427 0.1958640 1.767136
tarsus_sex1 2.259727 0.7239967 0.6135135 3.567762
tarsus_sex2 2.744785 1.0537686 0.6051135 4.468542
$animal$cor
, , bwt_sex1
Estimate Est.Error Q2.5 Q97.5
bwt_sex1 1.0000000 0.0000000 1.00000000 1.0000000
bwt_sex2 0.4823461 0.2872049 -0.23539931 0.8699222
tarsus_sex1 0.5711397 0.2473353 -0.06616842 0.8936581
tarsus_sex2 0.1685774 0.3102320 -0.48914110 0.7539686
, , bwt_sex2
Estimate Est.Error Q2.5 Q97.5
bwt_sex1 0.4823461 0.2872049 -0.2353993 0.8699222
bwt_sex2 1.0000000 0.0000000 1.0000000 1.0000000
tarsus_sex1 0.3772862 0.3765302 -0.5312556 0.9085651
tarsus_sex2 0.5246223 0.3336679 -0.3678213 0.8692250
, , tarsus_sex1
Estimate Est.Error Q2.5 Q97.5
bwt_sex1 0.5711397 0.2473353 -0.06616842 0.8936581
bwt_sex2 0.3772862 0.3765302 -0.53125561 0.9085651
tarsus_sex1 1.0000000 0.0000000 1.00000000 1.0000000
tarsus_sex2 0.4401433 0.2929178 -0.29616405 0.8720453
, , tarsus_sex2
Estimate Est.Error Q2.5 Q97.5
bwt_sex1 0.1685774 0.3102320 -0.4891411 0.7539686
bwt_sex2 0.5246223 0.3336679 -0.3678213 0.8692250
tarsus_sex1 0.4401433 0.2929178 -0.2961641 0.8720453
tarsus_sex2 1.0000000 0.0000000 1.0000000 1.0000000
$animal$cov
, , bwt_sex1
Estimate Est.Error Q2.5 Q97.5
bwt_sex1 1.6710810 0.7317405 0.39611874 3.005768
bwt_sex2 0.7428820 0.5579922 -0.10600955 2.049110
tarsus_sex1 1.8733421 1.2159805 -0.05815667 4.340959
tarsus_sex2 0.6471034 1.0646022 -1.25297349 2.930019
, , bwt_sex2
Estimate Est.Error Q2.5 Q97.5
bwt_sex1 0.742882 0.5579922 -0.10600955 2.049110
bwt_sex2 1.335327 0.8557561 0.03836336 3.122771
tarsus_sex1 1.105102 1.0991534 -0.68104784 3.388813
tarsus_sex2 2.171388 1.8198507 -0.18718074 5.946047
, , tarsus_sex1
Estimate Est.Error Q2.5 Q97.5
bwt_sex1 1.873342 1.215980 -0.05815667 4.340959
bwt_sex2 1.105102 1.099153 -0.68104784 3.388813
tarsus_sex1 5.630014 3.143741 0.37639882 12.728924
tarsus_sex2 3.150235 2.476892 -0.67170548 8.755673
, , tarsus_sex2
Estimate Est.Error Q2.5 Q97.5
bwt_sex1 0.6471034 1.064602 -1.2529735 2.930019
bwt_sex2 2.1713876 1.819851 -0.1871807 5.946047
tarsus_sex1 3.1502347 2.476892 -0.6717055 8.755673
tarsus_sex2 8.6431609 5.649764 0.3661935 19.967865
$byear
$byear$sd
Estimate Est.Error Q2.5 Q97.5
bwt_sex1 0.7989572 0.1620838 0.5318383 1.156162
bwt_sex2 1.1420876 0.1912205 0.8090083 1.549360
tarsus_sex1 2.2286834 0.4609182 1.4995560 3.183107
tarsus_sex2 2.3428101 0.4941316 1.5596304 3.405106
$byear$cor
, , bwt_sex1
Estimate Est.Error Q2.5 Q97.5
bwt_sex1 1.0000000 0.0000000 1.0000000 1.0000000
bwt_sex2 0.7404024 0.1542690 0.3534945 0.9558481
tarsus_sex1 -0.1137836 0.2357386 -0.5534464 0.3543949
tarsus_sex2 0.2922708 0.2345550 -0.1749420 0.7113421
, , bwt_sex2
Estimate Est.Error Q2.5 Q97.5
bwt_sex1 0.7404024 0.1542690 0.3534945 0.955848060
bwt_sex2 1.0000000 0.0000000 1.0000000 1.000000000
tarsus_sex1 -0.3874250 0.1987378 -0.7323699 0.004258363
tarsus_sex2 0.2937354 0.2059346 -0.1577533 0.659531338
, , tarsus_sex1
Estimate Est.Error Q2.5 Q97.5
bwt_sex1 -0.1137836 0.2357386 -0.5534464 0.354394929
bwt_sex2 -0.3874250 0.1987378 -0.7323699 0.004258363
tarsus_sex1 1.0000000 0.0000000 1.0000000 1.000000000
tarsus_sex2 0.5226217 0.1897518 0.1183176 0.839911077
, , tarsus_sex2
Estimate Est.Error Q2.5 Q97.5
bwt_sex1 0.2922708 0.2345550 -0.1749420 0.7113421
bwt_sex2 0.2937354 0.2059346 -0.1577533 0.6595313
tarsus_sex1 0.5226217 0.1897518 0.1183176 0.8399111
tarsus_sex2 1.0000000 0.0000000 1.0000000 1.0000000
$byear$cov
, , bwt_sex1
Estimate Est.Error Q2.5 Q97.5
bwt_sex1 0.6645776 0.2761781 0.2828520 1.3367115
bwt_sex2 0.6843681 0.2597687 0.2749460 1.2955568
tarsus_sex1 -0.1581409 0.4589758 -1.0281983 0.8339776
tarsus_sex2 0.5456796 0.5013040 -0.3608892 1.5951028
, , bwt_sex2
Estimate Est.Error Q2.5 Q97.5
bwt_sex1 0.6843681 0.2597687 0.2749460 1.295556807
bwt_sex2 1.3408929 0.4593417 0.6544951 2.400516888
tarsus_sex1 -1.0167438 0.6693667 -2.4775184 0.009371017
tarsus_sex2 0.8646682 0.7242669 -0.3384814 2.623722863
, , tarsus_sex1
Estimate Est.Error Q2.5 Q97.5
bwt_sex1 -0.1581409 0.4589758 -1.0281983 0.833977585
bwt_sex2 -1.0167438 0.6693667 -2.4775184 0.009371017
tarsus_sex1 5.1792626 2.3047474 2.2486683 10.132170288
tarsus_sex2 2.7818157 1.5318128 0.5297591 5.970346660
, , tarsus_sex2
Estimate Est.Error Q2.5 Q97.5
bwt_sex1 0.5456796 0.5013040 -0.3608892 1.595103
bwt_sex2 0.8646682 0.7242669 -0.3384814 2.623723
tarsus_sex1 2.7818157 1.5318128 0.5297591 5.970347
tarsus_sex2 5.7326811 2.5639312 2.4324504 11.594758
$mother
$mother$sd
Estimate Est.Error Q2.5 Q97.5
bwt_sex1 1.076447 0.1526711 0.7866154 1.390273
bwt_sex2 1.325206 0.1540539 1.0256350 1.621175
tarsus_sex1 2.214033 0.3976449 1.3647309 2.946910
tarsus_sex2 2.310902 0.4940802 1.3795156 3.339175
$mother$cor
, , bwt_sex1
Estimate Est.Error Q2.5 Q97.5
bwt_sex1 1.0000000 0.0000000 1.0000000 1.00000000
bwt_sex2 0.8260360 0.1110784 0.5676601 0.97826511
tarsus_sex1 -0.5024557 0.2385053 -0.9137423 -0.06798800
tarsus_sex2 -0.5073494 0.2162931 -0.8806394 -0.07677131
, , bwt_sex2
Estimate Est.Error Q2.5 Q97.5
bwt_sex1 0.8260360 0.1110784 0.5676601 0.9782651
bwt_sex2 1.0000000 0.0000000 1.0000000 1.0000000
tarsus_sex1 -0.6373889 0.1744852 -0.9272366 -0.2775698
tarsus_sex2 -0.3593672 0.2577845 -0.8364419 0.1135921
, , tarsus_sex1
Estimate Est.Error Q2.5 Q97.5
bwt_sex1 -0.5024557 0.2385053 -0.9137423 -0.0679880
bwt_sex2 -0.6373889 0.1744852 -0.9272366 -0.2775698
tarsus_sex1 1.0000000 0.0000000 1.0000000 1.0000000
tarsus_sex2 0.7206209 0.1561300 0.3688228 0.9518425
, , tarsus_sex2
Estimate Est.Error Q2.5 Q97.5
bwt_sex1 -0.5073494 0.2162931 -0.8806394 -0.07677131
bwt_sex2 -0.3593672 0.2577845 -0.8364419 0.11359208
tarsus_sex1 0.7206209 0.1561300 0.3688228 0.95184250
tarsus_sex2 1.0000000 0.0000000 1.0000000 1.00000000
$mother$cov
, , bwt_sex1
Estimate Est.Error Q2.5 Q97.5
bwt_sex1 1.182023 0.3330334 0.6187659 1.9328578
bwt_sex2 1.179468 0.2744898 0.6821359 1.7479435
tarsus_sex1 -1.110068 0.4557029 -1.9379754 -0.2103052
tarsus_sex2 -1.223308 0.5582766 -2.3558658 -0.1900476
, , bwt_sex2
Estimate Est.Error Q2.5 Q97.5
bwt_sex1 1.1794683 0.2744898 0.6821359 1.7479435
bwt_sex2 1.7798788 0.4088241 1.0519272 2.6282093
tarsus_sex1 -1.8437371 0.5678929 -2.9129058 -0.7723013
tarsus_sex2 -0.9438083 0.6393327 -2.0044262 0.4581811
, , tarsus_sex1
Estimate Est.Error Q2.5 Q97.5
bwt_sex1 -1.110068 0.4557029 -1.937975 -0.2103052
bwt_sex2 -1.843737 0.5678929 -2.912906 -0.7723013
tarsus_sex1 5.059904 1.7495725 1.862491 8.6842757
tarsus_sex2 3.692391 1.3309441 1.432767 6.5202798
, , tarsus_sex2
Estimate Est.Error Q2.5 Q97.5
bwt_sex1 -1.2233079 0.5582766 -2.355866 -0.1900476
bwt_sex2 -0.9438083 0.6393327 -2.004426 0.4581811
tarsus_sex1 3.6923914 1.3309441 1.432767 6.5202798
tarsus_sex2 5.5841373 2.3245515 1.903063 11.1500916The cross-sex genetic correlation can estimate form the output of the model. For tarsus length at fledging, sexes shared a lot of genetic variance which is commun for a trait with low sexual dimorphism. If the selection is antagonistic between males and females, sexes can not evolve freely from the other sexes and a intralocus sexual conflict can appeared.
cross_sex.cor.btw <- as.mcmc((VarCorr(brms_m2.5, summary = FALSE)$animal$cor[, 1, 2]))
cross_sex.cor.tarsus <- as.mcmc((VarCorr(brms_m2.5, summary = FALSE)$animal$cor[, 3, 4]))
summary(cross_sex.cor.btw)
Iterations = 1:1000
Thinning interval = 1
Number of chains = 1
Sample size per chain = 1000
1. Empirical mean and standard deviation for each variable,
plus standard error of the mean:
Mean SD Naive SE Time-series SE
0.482346 0.287205 0.009082 0.032430
2. Quantiles for each variable:
2.5% 25% 50% 75% 97.5%
-0.2354 0.3433 0.5365 0.6861 0.8699 summary(cross_sex.cor.tarsus)
Iterations = 1:1000
Thinning interval = 1
Number of chains = 1
Sample size per chain = 1000
1. Empirical mean and standard deviation for each variable,
plus standard error of the mean:
Mean SD Naive SE Time-series SE
0.440143 0.292918 0.009263 0.048863
2. Quantiles for each variable:
2.5% 25% 50% 75% 97.5%
-0.2962 0.2950 0.4846 0.6419 0.8720 plot(cross_sex.cor.btw)
plot(cross_sex.cor.tarsus)
Here, some simple code to extract the BLUP.
Here, some simple code to plot the cross-sex genetic correlation.
par(mfrow = c(1, 2))
plot(bwt_sex2 ~ bwt_sex1, bl_m2.5,
xlab = "", ylab = "", las = 1.2, type = "n",
xlim = c(min(bl_m2.5$bwt_sex1_lo), max(bl_m2.5$bwt_sex1_up)),
ylim = c(min(bl_m2.5$bwt_sex2_lo), max(bl_m2.5$bwt_sex2_up))
)
with(bl_m2.5, segments(x0 = bwt_sex1, y0 = bwt_sex2_lo, x1 = bwt_sex1, y1 = bwt_sex2_up, col = "black"))
with(bl_m2.5, segments(x0 = bwt_sex1_lo, y0 = bwt_sex2, x1 = bwt_sex1_up, y1 = bwt_sex2, col = "black"))
points(bwt_sex2 ~ bwt_sex1, bl_m2.5, pch = 16, col = "red", cex = 1.5)
points(bwt_sex2 ~ bwt_sex1, bl_m2.5, pch = 1, col = rgb(0, 0, 0, 0.3), cex = c(1.5))
mtext("bwt male (BV±CI)", side = 1, line = 2.4)
mtext("bwt female (BV±CI)", side = 2, line = 2, las = 3)
plot(tarsus_sex2 ~ tarsus_sex1, bl_m2.5,
xlab = "", ylab = "", las = 1.2, type = "n",
xlim = c(min(bl_m2.5$tarsus_sex1_lo), max(bl_m2.5$tarsus_sex1_up)),
ylim = c(min(bl_m2.5$tarsus_sex2_lo), max(bl_m2.5$tarsus_sex2_up))
)
with(bl_m2.5, segments(x0 = tarsus_sex1, y0 = tarsus_sex2_lo, x1 = tarsus_sex1, y1 = tarsus_sex2_up, col = "black"))
with(bl_m2.5, segments(x0 = tarsus_sex1_lo, y0 = tarsus_sex2, x1 = tarsus_sex1_up, y1 = tarsus_sex2, col = "black"))
points(tarsus_sex2 ~ tarsus_sex1, bl_m2.5, pch = 16, col = "red", cex = 1.5)
points(tarsus_sex2 ~ tarsus_sex1, bl_m2.5, pch = 1, col = rgb(0, 0, 0, 0.3), cex = c(1.5))
mtext("tarsus male (BV±CI)", side = 1, line = 2.4)
mtext("tarsus female (BV±CI)", side = 2, line = 2, las = 3)
Within this model, we also have acces to the rest of the B-matrix. Note, the cross-sex genetic correlation is just the diagonal of the B matrix. For now on, you can explore this matrix and estimate the cross-sex-cross-trait genetic correlation.