This online material is a supplement to our paper “Quantifying macro-evolutionary patterns of trait mean and variance with phylogenetic location-scale models”. You will see how to fit phylogenetic location-scale models (PLSMs) in a Bayesian framework using brms package in R. All the code and related data are available in this repository here.
In this online material, we will first introduce how to fit phylogenetic location-only models (nomral phylgoenetic mixed-effects model) using brms.
Then, we will illustrate how to fit three models described in our paper, subtitled as:
Different Trait Variance in Two Groups (Model 3)
Co-evolution of Mean and Variance (Model 4)
Co-evolution of Two Traits (Model 5)
All results presented in this online material are for illustrative purposes only (more careful thoughts are required to draw biological conclusions).
Additionally, we will show how to obtain phylogenetic heritability and macro-evolvability from the fitted brms object.
Also, importantly, there is a paper, which focuses on location-scale meta-analytic and meta-regression. You can find the sister paper’s tutorial here; also, the preprint of this sister paper is available on DOI:
Nakagawa, S., Mizuno, A., Morrison, K., Ricolfi, L., Williams, C., Drobniak, S. M., … & Yang, Y. (2025). Location‐Scale Meta‐Analysis and Meta‐Regression as a Tool to Capture Large‐Scale Changes in Biological and Methodological Heterogeneity: A Spotlight on Heteroscedasticity. Global Change Biology, 31(5), e70204.
This tutorial has several tips to run location-scale models using brms.
Furthermore, there are two papers, which teaches Bayesian generalised linear modelling.
The first one’s tutorial here; also, the preprint is available on DOI:
Halliwell, B., Holland, B. R., & Yates, L. A. (2025). Multi‐response phylogenetic mixed models: concepts and application. Biological Reviews, 100(3), 1294-1316.
The second one’s tutorial here; also, the preprint is available on EcoEvoRxiv:
Mizuno, A., Drobniak, S.M., Williams, C., Lagisz, M. and Nakagawa, S., 2025. Promoting the use of phylogenetic multinomial generalised mixed-effects model to understand the evolution of discrete traits. EcoEvoRxiv. 2025.
Also, these tutorials teach you how to diagnose Bayesian models and how to interpret the results.
Our tutorial uses R statistical software and existing R packages, which you will first need to download and install.
If the packages are archived in CRAN, use install.packages() to install them. For example, to install the brms , you can execute install.packages("brms") in the console (bottom left pane of R Studio).
Version information of each package is listed at the end of this tutorial.
# attempt to load or install necessary packages
if(!require(pacman)) install.packages("pacman")
pacman::p_load(tidyverse,
ggplot2,
brms,
bayesplot,
tidybayes,
here,
ape,
patchwork,
pander,
stringr,
tibble,
dplyr
)
html_out <- knitr::is_html_output() We also provide some additional helper functions to help visualize the results (mainly extracted from the brms object). If you want to use these custom functions in your own data, you’ll need to change the variable names according to your own data (check out the R code and you’ll see what we mean).
# Function to get variable names dynamically
get_variables_dynamic <- function(model, pattern) {
variables <- get_variables(model)
variables[grep(pattern, variables)]
}
rename_vars <- function(variable) {
# b
variable <- gsub("b_Intercept", "l_int", variable)
variable <- gsub("b_sigma_Intercept", "s_int", variable)
variable <- gsub("b_cmass", "l_cmass", variable)
variable <- gsub("b_sigma_cmass", "s_cmass", variable)
variable <- gsub("b_forest", "l_contrast", variable)
variable <- gsub("b_sigma_forest", "s_contrast", variable)
# mod3
variable <- gsub("b_cbeakwidth_Intercept", "l_width_int", variable)
variable <- gsub("b_sigma_cbeakwidth_Intercept", "s_width_int", variable)
variable <- gsub("b_cbeakdepth_Intercept", "l_depth_int", variable)
variable <- gsub("b_sigma_cbeakdepth_Intercept", "s_depth_int", variable)
variable <- gsub("b_cbeakwidth_cmass", "l_width_cmass", variable)
variable <- gsub("b_sigma_cbeakwidth_cmass", "s_width_cmass", variable)
variable <- gsub("b_cbeakdepth_cmass", "l_depth_cmass", variable)
variable <- gsub("b_sigma_cbeakdepth_cmass", "s_depth_cmass", variable)
# sd
variable <- gsub("sd_Phylo__Intercept", "l_sd", variable)
variable <- gsub("sd_Phylo__sigma_Intercept", "s_sd", variable)
# mod3
variable <- gsub("sd_Phylo__cbeakwidth_Intercept", "l_width_sd", variable)
variable <- gsub("sd_Phylo__sigma_cbeakwidth_Intercept", "s_width_sd", variable)
variable <- gsub("sd_Phylo__cbeakdepth_Intercept", "l_depth_sd", variable)
variable <- gsub("sd_Phylo__sigma_cbeakdepth_Intercept", "s_depth_sd", variable)
# corr
variable <- gsub("cor_Phylo__Intercept__sigma_Intercept", "ls_cor", variable)
# mod3
variable <- gsub("cor_Phylo__cbeakwidth_Intercept__sigma_cbeakwidth_Intercept", "ls_width-width_cor", variable)
variable <- gsub("cor_Phylo__cbeakwidth_Intercept__cbeakdepth_Intercept", "ll_width-depth_cor", variable)
variable <- gsub("cor_Phylo__sigma_cbeakwidth_Intercept__cbeakdepth_Intercept", "sl_width-depth_cor", variable)
variable <- gsub("cor_Phylo__cbeakwidth_Intercept__sigma_cbeakdepth_Intercept", "ls_width-depth_cor", variable)
variable <- gsub("cor_Phylo__sigma_cbeakwidth_Intercept__sigma_cbeakdepth_Intercept", "ss_width-depth_cor", variable)
variable <- gsub("cor_Phylo__cbeakdepth_Intercept__sigma_cbeakdepth_Intercept", "ls_depth-depth_cor", variable)
return(variable)
}
# Function to visualize fixed effects
visualize_fixed_effects <- function(model) {
fixed_effect_vars <- get_variables_dynamic(model, "^b_")
if (length(fixed_effect_vars) == 0) {
message("No fixed effects found")
return(NULL)
}
tryCatch({
fixed_effects_samples <- model %>%
spread_draws(!!!syms(fixed_effect_vars)) %>%
pivot_longer(cols = all_of(fixed_effect_vars), names_to = ".variable", values_to = ".value") %>%
mutate(.variable = rename_vars(.variable))
ggplot(fixed_effects_samples, aes(x = .value, y = .variable)) +
stat_halfeye(
normalize = "xy",
point_interval = "mean_qi",
fill = "lightcyan3",
color = "lightcyan4"
) +
geom_vline(xintercept = 0, linetype = "dashed", color = "#005") +
labs(y = "Fixed effects", x = "Posterior values") +
theme_classic()
}, error = function(e) {
message("Error in visualize_fixed_effects: ", e$message)
return(NULL)
})
}
# Function to visualize random effects
visualize_random_effects <- function(model) {
random_effect_vars <- get_variables_dynamic(model, "^sd_")
random_effect_vars <- random_effect_vars[random_effect_vars != "sd_es_ID__Intercept"]
if (length(random_effect_vars) == 0) {
message("No random effects found")
return(NULL)
}
tryCatch({
random_effects_samples <- model %>%
spread_draws(!!!syms(random_effect_vars)) %>%
pivot_longer(cols = all_of(random_effect_vars), names_to = ".variable", values_to = ".value") %>%
mutate(.variable = rename_vars(.variable)) #%>%
#mutate(.value = .value) # leave SD as it is
ggplot(random_effects_samples, aes(x = .value, y = .variable)) +
stat_halfeye(
normalize = "xy",
point_interval = "mean_qi",
fill = "olivedrab3",
color = "olivedrab4"
) +
geom_vline(xintercept = 0, linetype = "dashed", color = "#005") +
labs(y = "Random effects (SD)", x = "Posterior values") +
theme_classic()
}, error = function(e) {
message("Error in visualize_random_effects: ", e$message)
return(NULL)
})
}
# Function to visualize correlations
visualize_correlations <- function(model) {
correlation_vars <- get_variables_dynamic(model, "^cor_")
if (length(correlation_vars) == 0) {
message("No correlations found")
return(NULL)
}
tryCatch({
correlation_samples <- model %>%
spread_draws(!!!syms(correlation_vars)) %>%
pivot_longer(cols = all_of(correlation_vars), names_to = ".variable", values_to = ".value") %>%
mutate(.variable = rename_vars(.variable))
ggplot(correlation_samples, aes(x = .value, y = .variable)) +
stat_halfeye(
normalize = "xy",
fill = "#FF6347",
color = "#8B3626"
) +
geom_vline(xintercept = 0, linetype = "dashed", color = "#005") +
labs(y = "Correlations", x = "Posterior values") +
theme_classic()
}, error = function(e) {
message("Error in visualize_correlations: ", e$message)
return(NULL)
})
}We use the data from Avonet:
J. A. Tobias, C. Sheard, A. L. Pigot, A. J. Devenish, J. Yang, F. Sayol, M. H. Neate-Clegg, N. Alioravainen, T. L. Weeks, R. A. Barber, et al. Avonet: morphological, ecological and geographical data for all birds. Ecology Letters, 25(3):581–597, 2022
As described in the paper, we only use parrot species (N = 354) and we use phylogenetic trees form Jetz et al. (2012):
W. Jetz, G. H. Thomas, J. B. Joy, K. Hartmann, and A. O. Mooers. The global diversity of birds in space and time. Nature, 491(7424):444–448, 2012
# load data
dat <- read.csv(here("data", "Psittaciformes", "Psittaciformes_354spp.csv"))
#str(dat)
# reading nexus file (use ape package)
tree_all <- read.nexus(here("data", "Psittaciformes", "Psittaciformes_354spp_100.nex"))
tree <- tree_all[[1]]
# Turning phylogenetic tree into phylognetic correlation matrix
A <- vcv.phylo(tree, corr = TRUE)
# center data of interest
dat$cbeak_length <- scale(log(dat$Beak.Length_Culmen), center = TRUE, scale = FALSE)
dat$cbeak_width <- scale(log(dat$Beak.Width), center = TRUE, scale = FALSE)
dat$cbeak_depth <- scale(log(dat$Beak.Depth), center = TRUE, scale = FALSE)
dat$cmass <- scale(log(dat$Mass), center = TRUE, scale = FALSE)
dat$crange_size <- scale(log(dat$Range.Size), center = TRUE, scale = FALSE)
# phynotipic correlation among three beak traits
# cor(dat$cbeak_length, dat$cbeak_width)
# cor(dat$cbeak_length, dat$cbeak_depth)
# cor(dat$cbeak_width, dat$cbeak_depth)
# create a variable which says whether they live in Forest or not
dat$forest <- ifelse(dat$Habitat == "Forest", 1, 0)We use the following variables in our models:
cbeak_length: log-transformed and centred beak length (Beak.Length_Culmen)
cbeak_width: log-transformed and centred beak width (Beak.Width)
cbeak_depth: log-transformed and centred beak depth (Beak.Depth)
cmass: log-transformed and centred body mass (Mass)
crange_size: log-transformed and centred range size (Range.Size)
forest: a binary variable indicating whether the species lives in the forest or other habits (created based on Habitat)
Before adding a scale component it is helpful to see the “classic” phylogenetic generalised linear mixed model (PGLMM) that acts only on the mean.
Here we predict centred log-range size (crange_size) from centred body mass (cmass), while allowing a Brownian-motion random effect (Phylo) on the location part and letting the residual variance stay constant across species.
# creating the formula
form0 <- bf(cbeak_length ~ 1 + cmass + forest +
(1 | gr(Phylo, cov = A)) # phylogenetic random effect
)
# setting priors
prior0 <- default_prior(form0,
data = dat,
data2 = list(A = A),
family = gaussian()
)
# running the model
mod0 <- brm(form0,
data = dat,
data2 = list(A = A),
chains = 2,
cores = 2,
iter = 8000,
warmup = 5000,
prior = prior0,
control = list(adapt_delta = 0.99, max_treedepth = 15)
)
# time taken to run the model
# 179.714 seconds (Total)
saveRDS(mod0, here("Rdata", "mod0.rds"))Check the results:
mod0 <- readRDS(here("Rdata", "mod0.rds"))
summary(mod0)
## Family: gaussian
## Links: mu = identity; sigma = identity
## Formula: cbeak_length ~ 1 + cmass + forest + (1 | gr(Phylo, cov = A))
## Data: dat (Number of observations: 353)
## Draws: 2 chains, each with iter = 8000; warmup = 5000; thin = 1;
## total post-warmup draws = 6000
##
## Multilevel Hyperparameters:
## ~Phylo (Number of levels: 353)
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept) 0.26 0.03 0.20 0.32 1.00 605 1151
##
## Regression Coefficients:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept 0.02 0.14 -0.25 0.30 1.00 2006 2917
## cmass 0.34 0.02 0.31 0.37 1.00 1634 3672
## forest 0.03 0.02 -0.01 0.06 1.00 4329 4685
##
## Further Distributional Parameters:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma 0.09 0.01 0.08 0.11 1.00 633 1294
##
## 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).# Fixed and random effects
plots_mod0 <- list(
visualize_fixed_effects(mod0),
visualize_random_effects(mod0)
)
plots_mod0[[1]] / plots_mod0[[2]] + plot_layout(heights = c(2, 1))
Posterior distributions of different types of parameters; the vertical dashed line indicates zero, aiding the interpretation of effect direction with thick lines showing 66% credible intervals and thin whiskers 95% credible intervals (not spanning zero being interpreted as ‘statistical significant’).
In this first phylogenetic location-scale model (PLSM), we will fit a location-scale model to investigate whether the mean and variance of beak length (cbeak_length) differs between species living in the forest and those living in other habitats.
The brms syntax for Model 3 is as follows:
# creating the formula
form1 <- bf(cbeak_length ~ 1 + cmass + forest +
(1 | gr(Phylo, cov = A)), # how to model phylogenetic correlation
sigma ~ 1 + cmass + forest
)
# setting priors
prior1 <- default_prior(form1,
data = dat,
data2 = list(A = A),
family = gaussian()
)
# running the model
mod1 <- brm(form1,
data = dat,
data2 = list(A = A),
chains = 2,
cores = 2,
iter = 80000,
warmup = 20000,
prior = prior1,
control = list(adapt_delta = 0.99, max_treedepth = 15)
)
# time taken to run the model
# 8673 seconds (Total)
# saving
saveRDS(mod1, here("Rdata", "mod1.rds"))Check the results:
# loading and checking
mod1 <- readRDS(here("Rdata", "mod1.rds"))
summary(mod1)
## Warning: There were 329 divergent transitions after warmup. Increasing
## adapt_delta above 0.99 may help. See
## http://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup
## Family: gaussian
## Links: mu = identity; sigma = log
## Formula: cbeak_length ~ 1 + cmass + forest + (1 | gr(Phylo, cov = A))
## sigma ~ 1 + cmass + forest
## Data: dat (Number of observations: 353)
## Draws: 2 chains, each with iter = 80000; warmup = 20000; thin = 1;
## total post-warmup draws = 120000
##
## Multilevel Hyperparameters:
## ~Phylo (Number of levels: 353)
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept) 0.27 0.03 0.22 0.33 1.00 3122 6870
##
## Regression Coefficients:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept 0.00 0.14 -0.27 0.28 1.00 24037 39809
## sigma_Intercept -3.35 0.47 -4.65 -2.71 1.00 1852 1237
## cmass 0.36 0.01 0.33 0.39 1.00 28036 51686
## forest 0.02 0.01 -0.01 0.05 1.00 21346 54450
## sigma_cmass -0.22 0.08 -0.39 -0.07 1.00 7401 15081
## sigma_forest 1.00 0.44 0.39 2.23 1.00 2049 1263
##
## 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).We can visualize the results (posterior distributions) of parameters (i.e., regression coefficients, variance and standard deviation components, and correlations) using the custom functions we defined earlier.
# getting plots
plots_mod1 <- list(
visualize_fixed_effects(mod1),
visualize_random_effects(mod1)
)
plots_mod1[[1]] / plots_mod1[[2]] + plot_layout(heights = c(2, 1))
Posterior distributions of different types of parameters; the vertical dashed line indicates zero, aiding the interpretation of effect direction with thick lines showing 66% credible intervals and thin whiskers 95% credible intervals (not spanning zero being interpreted as ‘statistical significant’).
Trait variability itself is structured: larger species are less variable beak length.
Forest dwellers are more variable, pointing to environment-specific heteroscedasticity athough the difference in the location model was limited.
In this second model, we will fit a phylogenetic location-scale model to investigate whether the mean and variance of beak length (cbeak_length) are evolving together among parrots, after controlling body mass (cmass).
The brms syntax for Model 4 is as follows:
# creating the formula
# |p| indicates we are going to estimate correlations between these two random effects
form2 <- bf(crange_size ~1 + cmass + (1|p|gr(Phylo, cov = A)),
sigma ~ 1 + cmass + (1|p|gr(Phylo, cov = A))
)
# setting priors
prior2 <- default_prior(form2,
data = dat,
data2 = list(A = A),
family = gaussian()
)
# running the model
mod2 <- brm(form2,
data = dat,
data2 = list(A = A),
chains = 2,
cores = 2,
iter = 6000,
warmup = 3000,
prior = prior2,
control = list(adapt_delta = 0.99, max_treedepth = 15)
)
# time taken to run the model
# 14265.3 seconds (Total)
# saving
saveRDS(mod2, here("Rdata", "mod2.rds"))Check the results:
# loading and checking
mod2 <- readRDS(here("Rdata", "mod2.rds"))
summary(mod2)
## Family: gaussian
## Links: mu = identity; sigma = log
## Formula: crange_size ~ 1 + cmass + (1 | p | gr(Phylo, cov = A))
## sigma ~ 1 + cmass + (1 | p | gr(Phylo, cov = A))
## Data: dat (Number of observations: 351)
## Draws: 2 chains, each with iter = 6000; warmup = 3000; thin = 1;
## total post-warmup draws = 6000
##
## Multilevel Hyperparameters:
## ~Phylo (Number of levels: 351)
## Estimate Est.Error l-95% CI u-95% CI Rhat
## sd(Intercept) 3.11 0.40 2.36 3.92 1.00
## sd(sigma_Intercept) 0.78 0.15 0.51 1.10 1.00
## cor(Intercept,sigma_Intercept) -0.94 0.05 -1.00 -0.81 1.00
## Bulk_ESS Tail_ESS
## sd(Intercept) 1417 2702
## sd(sigma_Intercept) 1612 4149
## cor(Intercept,sigma_Intercept) 2677 3439
##
## Regression Coefficients:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept 0.17 1.38 -2.59 2.86 1.00 3439 3368
## sigma_Intercept 0.63 0.38 -0.12 1.36 1.00 3231 3588
## cmass -0.25 0.23 -0.72 0.20 1.00 3378 3552
## sigma_cmass 0.18 0.08 0.03 0.33 1.00 3569 4099
##
## 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).As above, using the custom functions:
# getting plots
plots_mod2 <- list(
visualize_fixed_effects(mod2),
visualize_random_effects(mod2),
visualize_correlations(mod2)
)
plots_mod2[[1]] / plots_mod2[[2]] /plots_mod2[[3]] + plot_layout(heights = c(2, 1, 1))
Posterior distributions of different types of parameters; the vertical dashed line indicates zero, aiding the interpretation of effect direction with thick lines showing 66% credible intervals and thin whiskers 95% credible intervals (not spanning zero being interpreted as ‘statistical significant’).
The range size’s mean and variance are evolving together, as indicated by the strong negative correlation
That is, larger species have larger range size but smaller variance in range size.
In this third model, we will fit a phylogenetic bivariate location-scale model to investigate whether the mean and variance of beak length (cbeak_length) and beak width (cbeak_width) are evolving together among parrots, after controlling body mass (cmass).
The brms syntax for Model 5 is as follows:
# creating the formula
# |p| indicates we are going to estimate correlations between these 4 random effects
form3A <- bf(cbeak_width ~1 + cmass + (1|p|gr(Phylo, cov = A)),
sigma ~ 1 + cmass + (1|p|gr(Phylo, cov = A))
)
form3B <- bf(cbeak_depth ~1 + cmass + (1|p|gr(Phylo, cov = A)),
sigma ~ 1 + cmass + (1|p|gr(Phylo, cov = A))
)
form3 <- form3A + form3B + set_rescor(TRUE)
# setting prior
prior3 <- default_prior(form3,
data = dat,
data2 = list(A = A),
family = gaussian()
)
# fit model
mod3 <- brm(form3,
data = dat,
data2 = list(A = A),
chains = 2,
cores = 2,
iter = 15000,
warmup = 5000,
prior = prior3,
control = list(adapt_delta = 0.99, max_treedepth = 15)
)
summary(mod3)
# time taken to run the model
# 439397 seconds (Total)
# saving
saveRDS(mod3, here("Rdata", "mod3.rds"))Check the results:
# loading and checking
mod3 <- readRDS(here("Rdata", "mod3.rds"))
summary(mod3)
## Warning: There were 7 divergent transitions after warmup. Increasing
## adapt_delta above 0.99 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: cbeak_width ~ 1 + cmass + (1 | p | gr(Phylo, cov = A))
## sigma ~ 1 + cmass + (1 | p | gr(Phylo, cov = A))
## cbeak_depth ~ 1 + cmass + (1 | p | gr(Phylo, cov = A))
## sigma ~ 1 + cmass + (1 | p | gr(Phylo, cov = A))
## Data: dat (Number of observations: 354)
## Draws: 2 chains, each with iter = 15000; warmup = 5000; thin = 1;
## total post-warmup draws = 20000
##
## Multilevel Hyperparameters:
## ~Phylo (Number of levels: 354)
## Estimate Est.Error
## sd(cbeakwidth_Intercept) 0.23 0.02
## sd(sigma_cbeakwidth_Intercept) 0.76 0.16
## sd(cbeakdepth_Intercept) 0.27 0.02
## sd(sigma_cbeakdepth_Intercept) 0.94 0.17
## cor(cbeakwidth_Intercept,sigma_cbeakwidth_Intercept) 0.23 0.16
## cor(cbeakwidth_Intercept,cbeakdepth_Intercept) 0.89 0.03
## cor(sigma_cbeakwidth_Intercept,cbeakdepth_Intercept) 0.28 0.16
## cor(cbeakwidth_Intercept,sigma_cbeakdepth_Intercept) 0.36 0.16
## cor(sigma_cbeakwidth_Intercept,sigma_cbeakdepth_Intercept) 0.82 0.13
## cor(cbeakdepth_Intercept,sigma_cbeakdepth_Intercept) 0.28 0.15
## l-95% CI u-95% CI
## sd(cbeakwidth_Intercept) 0.19 0.27
## sd(sigma_cbeakwidth_Intercept) 0.46 1.10
## sd(cbeakdepth_Intercept) 0.23 0.32
## sd(sigma_cbeakdepth_Intercept) 0.64 1.32
## cor(cbeakwidth_Intercept,sigma_cbeakwidth_Intercept) -0.10 0.52
## cor(cbeakwidth_Intercept,cbeakdepth_Intercept) 0.82 0.94
## cor(sigma_cbeakwidth_Intercept,cbeakdepth_Intercept) -0.04 0.57
## cor(cbeakwidth_Intercept,sigma_cbeakdepth_Intercept) 0.02 0.64
## cor(sigma_cbeakwidth_Intercept,sigma_cbeakdepth_Intercept) 0.48 0.98
## cor(cbeakdepth_Intercept,sigma_cbeakdepth_Intercept) -0.03 0.56
## Rhat Bulk_ESS
## sd(cbeakwidth_Intercept) 1.00 3184
## sd(sigma_cbeakwidth_Intercept) 1.00 3910
## sd(cbeakdepth_Intercept) 1.00 2349
## sd(sigma_cbeakdepth_Intercept) 1.00 2727
## cor(cbeakwidth_Intercept,sigma_cbeakwidth_Intercept) 1.00 6929
## cor(cbeakwidth_Intercept,cbeakdepth_Intercept) 1.00 3649
## cor(sigma_cbeakwidth_Intercept,cbeakdepth_Intercept) 1.00 3773
## cor(cbeakwidth_Intercept,sigma_cbeakdepth_Intercept) 1.00 5411
## cor(sigma_cbeakwidth_Intercept,sigma_cbeakdepth_Intercept) 1.00 2684
## cor(cbeakdepth_Intercept,sigma_cbeakdepth_Intercept) 1.00 5414
## Tail_ESS
## sd(cbeakwidth_Intercept) 7271
## sd(sigma_cbeakwidth_Intercept) 7196
## sd(cbeakdepth_Intercept) 4902
## sd(sigma_cbeakdepth_Intercept) 4447
## cor(cbeakwidth_Intercept,sigma_cbeakwidth_Intercept) 11541
## cor(cbeakwidth_Intercept,cbeakdepth_Intercept) 8642
## cor(sigma_cbeakwidth_Intercept,cbeakdepth_Intercept) 7289
## cor(cbeakwidth_Intercept,sigma_cbeakdepth_Intercept) 10329
## cor(sigma_cbeakwidth_Intercept,sigma_cbeakdepth_Intercept) 3927
## cor(cbeakdepth_Intercept,sigma_cbeakdepth_Intercept) 9976
##
## Regression Coefficients:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
## cbeakwidth_Intercept -0.05 0.13 -0.30 0.21 1.00 6878
## sigma_cbeakwidth_Intercept -1.99 0.46 -2.90 -1.06 1.00 7701
## cbeakdepth_Intercept -0.05 0.15 -0.34 0.24 1.00 7549
## sigma_cbeakdepth_Intercept -1.85 0.53 -2.88 -0.80 1.00 9697
## cbeakwidth_cmass 0.32 0.01 0.29 0.34 1.00 8547
## sigma_cbeakwidth_cmass -0.10 0.09 -0.26 0.08 1.00 6114
## cbeakdepth_cmass 0.36 0.02 0.33 0.39 1.00 7883
## sigma_cbeakdepth_cmass -0.14 0.10 -0.36 0.06 1.00 3795
## Tail_ESS
## cbeakwidth_Intercept 9798
## sigma_cbeakwidth_Intercept 10102
## cbeakdepth_Intercept 11079
## sigma_cbeakdepth_Intercept 11162
## cbeakwidth_cmass 12524
## sigma_cbeakwidth_cmass 10313
## cbeakdepth_cmass 13214
## sigma_cbeakdepth_cmass 4953
##
## Residual Correlations:
## Estimate Est.Error l-95% CI u-95% CI Rhat
## rescor(cbeakwidth,cbeakdepth) 0.74 0.05 0.65 0.82 1.00
## Bulk_ESS Tail_ESS
## rescor(cbeakwidth,cbeakdepth) 6469 11074
##
## 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).Again, using the custom functions:
# getting plots
plots_mod3 <- list(
visualize_fixed_effects(mod3),
visualize_random_effects(mod3),
visualize_correlations(mod3)
)
plots_mod3[[1]] / plots_mod3[[2]] /plots_mod3[[3]]
Posterior distributions of different types of parameters; the vertical dashed line indicates zero, aiding the interpretation of effect direction with thick lines showing 66% credible intervals and thin whiskers 95% credible intervals (not spanning zero being interpreted as ‘statistical significant’).
This table links each conceptual cell in Figure 2 (the main text) to the corresponding distributional brms formulae and the parameters to report. Use the same group-level ID (e.g., a) in random effects you want correlated across parts/traits.
Implementing the Figure 2 cells in brms (pseudo-code) |
||
|---|---|---|
| Figure 2 cell / biological question | Minimal distributional formulae | |
| Within-trait mean–variance (ρ[a(ls)]) | Do lineages with larger means show greater or lower dispersion in the *same* trait? | bf(
y ~ 1 + X + (1 | a | gr(species, cov = A)),
sigma ~ 1 + X + (1 | a | gr(species, cov = A))
)
family = gaussian()
|
| Across-trait mean–mean (ρ[a(l1 l2)]) | Do two traits coevolve in their means (integration vs. trade-off)? | bf1 <- bf(
y1 ~ 1 + X1 + (1 | a | gr(species, cov = A)),
sigma ~ 1 + X1 + (1 | a | gr(species, cov = A))
)
bf2 <- bf(
y2 ~ 1 + X2 + (1 | a | gr(species, cov = A)),
sigma ~ 1 + X2 + (1 | a | gr(species, cov = A))
)
brm(bf1 + bf2, data2 = list(A = A), family = gaussian(), ...)
|
| Across-trait variance–variance (ρ[a(s1 s2)]) | Do two traits co-diverge (or contra-diverge) in dispersion? | # same multivariate structure as above; ensure both traits have distributional parts
bf1 <- bf(
y1 ~ 1 + X1 + (1 | a | gr(species, cov = A)),
sigma ~ 1 + X1 + (1 | a | gr(species, cov = A))
)
bf2 <- bf(
y2 ~ 1 + X2 + (1 | a | gr(species, cov = A)),
sigma ~ 1 + X2 + (1 | a | gr(species, cov = A))
)
|
| Across-trait mean–variance (ρ[a(l1 s2)], ρ[a(l2 s1)]) | Does a shift in the mean of trait 1 relax or constrain variability in trait 2 (and vice versa)? | # same multivariate distributional model; cross-part/trait correlations learned via shared ID `a`
brm(bf1 + bf2, data2 = list(A = A), family = gaussian(), ...)
|
Below, we show how to calculate phylogenetic heritability and macro-evolvability from the fitted brms object with and without a fixed effect(s). As you will see, phylogenetic heritability is fairly straightforward to obtain but macro-evolvability is a bit more complex as it is not clear what scale to calculate it on (e.g., log-scale or original scale). In fact, the same issue should be relevant to phylogenetic heritability as well but we do not consider this here.
We first create a phylogenetic location-scale (double-hierarchical) model for cmass without any fixed effects.
The brms syntax for this phylogenetic location-scale model is as follows:
# creating the formula
# |p| indicates we are going to estimate correlations between these 2 random effects
form4 <- bf(crange_size ~ 1 + (1 | p | gr(Phylo, cov = A)),
sigma ~ 1 + (1 | p | gr(Phylo, cov = A))
)
# setting prior
prior4 <- default_prior(form4,
data = dat,
data2 = list(A = A),
family = gaussian()
)
# fit model
mod4 <- brm(form4,
data = dat,
data2 = list(A = A),
chains = 2,
cores = 2,
iter = 8000,
warmup = 1000,
prior = prior4,
control = list(adapt_delta = 0.99, max_treedepth = 15)
)
# time taken to run the model
# 8523.23 seconds (Total)
# saving
saveRDS(mod4, here("Rdata", "mod4.rds"))Check the results:
# loading and checking
mod4 <- readRDS(here("Rdata", "mod4.rds"))
summary(mod4)
## Family: gaussian
## Links: mu = identity; sigma = log
## Formula: crange_size ~ 1 + (1 | p | gr(Phylo, cov = A))
## sigma ~ 1 + (1 | p | gr(Phylo, cov = A))
## Data: dat (Number of observations: 351)
## Draws: 2 chains, each with iter = 8000; warmup = 1000; thin = 1;
## total post-warmup draws = 14000
##
## Multilevel Hyperparameters:
## ~Phylo (Number of levels: 351)
## Estimate Est.Error l-95% CI u-95% CI Rhat
## sd(Intercept) 3.14 0.43 2.34 4.02 1.00
## sd(sigma_Intercept) 0.81 0.15 0.53 1.14 1.00
## cor(Intercept,sigma_Intercept) -0.92 0.06 -0.99 -0.77 1.00
## Bulk_ESS Tail_ESS
## sd(Intercept) 2357 4994
## sd(sigma_Intercept) 2820 4818
## cor(Intercept,sigma_Intercept) 3935 4577
##
## Regression Coefficients:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept -0.12 1.40 -2.94 2.62 1.00 3716 6128
## sigma_Intercept 0.72 0.39 -0.07 1.50 1.00 3684 5326
##
## 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).Using Equations 14, 15 & 17 in the paper, we can calculate the phylogenetic heritability of the trait crange_size (the location part):
# sigma2_p (total phenotypic variance)
# posterior draws
post <- as_draws_df(mod4)
# Equation 17: first calculate sigma2_e_bar (average residual variance)
sigma2_e_bar <- exp(2*post$b_sigma_Intercept + 2*post$sd_Phylo__sigma_Intercept^2)
# Equation 15 (no fixed effects)
sigma2_p <- post$sd_Phylo__Intercept^2 + sigma2_e_bar
# Equation 14
H2l <- post$sd_Phylo__Intercept^2 / sigma2_p
#summary(H2)
quantile(H2l, probs = c(0.025, 0.5, 0.975)) %>% pander()| 2.5% | 50% | 97.5% |
|---|---|---|
| 0.07655 | 0.3863 | 0.7532 |
We can see the median phylogenetic heritablity is around 39%.
Using Equations 18-21 in the paper, we can calculate the phylogenetic heritability of the trait crange_size (the scale part):
# Equation 21 - sigma2_sigma2_e (variance of the residual variance)
sigma2_sigma2_e <- (exp(4*post$sd_Phylo__sigma_Intercept^2) - 1) * exp(4*(post$b_sigma_Intercept + post$sd_Phylo__sigma_Intercept^2))
# Equation 19 - sigma2_a_star
sigma2_as_star <- sigma2_sigma2_e
# Equation 20 - sigma2_sigma2_p (variance of the total phenotypic variance)
sigma2_sigma2_p = sigma2_p ^ 2 + 3 * sigma2_sigma2_e
# Equation 18
H2s <- sigma2_as_star / sigma2_sigma2_p
#summary(H2)
quantile(H2s, probs = c(0.025, 0.5, 0.975)) %>% pander()| 2.5% | 50% | 97.5% |
|---|---|---|
| 0.1537 | 0.3081 | 0.3323 |
Here we found a median phylogenetic heritability of around 31%.
Using Equations 6 in the paper, we can calculate the macro-evolvability (CV_Al) of the trait crange_size (log-transformed trait; the location part), but the difficulty is that we need to calculate the phylogenetic variance on the original scale. We can do this by using the delta method. The details of this calculation are provided in this earlier work (see the supplementary material):
O’Dea RE, Noble DW, Nakagawa S. Unifying individual differences in personality, predictability and plasticity: a practical guide. Methods in Ecology and Evolution. 2022 Feb;13(2):278-93.
# exp(b0) is the median of the trait
ln_b0 <- post$b_Intercept + mean(log(dat$Range.Size), na.rm = T)
sigma2_al <- post$sd_Phylo__Intercept^2
mean_b0 <- exp(ln_b0 + 0.5 * (sigma2_al + sigma2_e_bar) )
# conversion formula from log-scale to original scale
sd_al <- sqrt( (exp(sigma2_al) - 1) * exp(2*ln_b0 + sigma2_al) )
CV_Al <- sd_al / mean_b0
quantile(CV_Al, probs = c(0.025, 0.5, 0.975)) %>% pander()| 2.5% | 50% | 97.5% |
|---|---|---|
| 3.659e-24 | 0.06627 | 58.41 |
We can see the median macro-evolvability is around 0.07
# alternatively, we can get CV on ln-scale
CV_Al2 <- sigma2_al / ln_b0
quantile(CV_Al2, probs = c(0.025, 0.5, 0.975)) %>% pander()| 2.5% | 50% | 97.5% |
|---|---|---|
| 0.4634 | 0.8535 | 1.549 |
If we get the CV on the log-scale, we can see the median macro-evolvability is around 0.85 (nearly 10 times bigger); this CV for the location is probably more comparable to CV for the scale below (as it turns, our CV for the scale is not estimable on the original scale).
We now estimate the macro-evolvability of the trait crange_size (the scale part) using Equation 24 (we assume that this CV_As is the same as the value for the original scale):
CV_As <- sqrt(exp(4*post$sd_Phylo__sigma_Intercept^2) - 1)
quantile(CV_As, probs = c(0.025, 0.5, 0.975)) %>% pander()| 2.5% | 50% | 97.5% |
|---|---|---|
| 1.46 | 3.423 | 13.26 |
We get the median macro-evolvability of around 3.4. This is macro-evolvability of variance but we can get macro-evolvability of the standard deviation using Equation 25.
# Alternatively, we can get CV_As on the SD scale
CV_As2 <- sqrt(exp(post$sd_Phylo__sigma_Intercept^2) - 1)
quantile(CV_As2, probs = c(0.025, 0.5, 0.975)) %>% pander()| 2.5% | 50% | 97.5% |
|---|---|---|
| 0.5746 | 0.9425 | 1.627 |
We get the median macro-evolvability of around 0.94, which is not that different from the macro-evolvability of the mean estimated on the log scale.
Also, we can use Equation 23:
CV_As2 <- sqrt(sigma2_as_star) / sigma2_e_bar
quantile(CV_As2, probs = c(0.025, 0.5, 0.975), na.rm = T) %>% pander()| 2.5% | 50% | 97.5% |
|---|---|---|
| 1.46 | 3.423 | 13.26 |
This estimation, as expected, matches with that from Equation 24.
We also try to get the macro-evolvability for the scale part on the original scale
# original scale conversion does not seem to make sense so commented out
sigma2_as_star_orig <- sqrt( (exp(sigma2_as_star) - 1) * exp(2*ln_b0 + sigma2_as_star) )
sigma2_e_bar_orig <- sqrt( (exp(sigma2_e_bar) - 1) * exp(2*ln_b0 + sigma2_e_bar) )
CV_As2_orig <- sqrt(sigma2_as_star_orig) / sigma2_e_bar_orig
quantile(CV_As2_orig, probs = c(0.025, 0.5, 0.975), na.rm = T) %>% pander()| 2.5% | 50% | 97.5% |
|---|---|---|
| 11483244 | Inf | Inf |
We get the median macro-evolvability of a Inf value so this does not seem to work.
Here we use mod2 as an example to see the phylogenetic heritability and macro-evolvability of the trait crange_size after controlling for body mass cmass.
Using Equations 14-17 in the paper, we can calculate the phylogenetic heritability of the trait crange_size (the location part) after controlling for body mass cmass:
# sigma2_p (total phenotypic variance)
# posterior draws
post <- as_draws_df(mod2)
# Equation 17: first calculate sigma2_e_bar (average residual variance)
sigma2_e_bar <- exp(2*post$b_sigma_Intercept + 2*post$sd_Phylo__sigma_Intercept^2)
# Equation 16 - sigma2_f (variance of the fixed effects)
sigma2_fl <- sapply(post$b_cmass, function(b) var(b*dat$cmass))
# Equation 15 (no fixed effects)
sigma2_p <- post$sd_Phylo__Intercept^2 + sigma2_fl + sigma2_e_bar
# Equation 14
H2l <- post$sd_Phylo__Intercept^2 / sigma2_p
#summary(H2)
quantile(H2l, probs = c(0.025, 0.5, 0.975)) %>% pander()| 2.5% | 50% | 97.5% |
|---|---|---|
| 0.1075 | 0.4399 | 0.763 |
We can see the median phylogenetic heritablity is around 44%.
Using Equations 18-21 in the paper, we can calculate the phylogenetic heritability of the trait crange_size (the scale part):
# Equation 22 - sigma2_fs
sigma2_fs <- sapply(post$b_sigma_cmass, function(b) var(b*dat$cmass))
# Equation 21 - sigma2_sigma2_e (variance of the residual variance)
sigma2_sigma2_e <- (exp(4*(post$sd_Phylo__sigma_Intercept^2 + sigma2_fs)) - 1) * exp(4*(post$b_sigma_Intercept + post$sd_Phylo__sigma_Intercept^2 + sigma2_fs))
# Equation 19 - sigma2_a_star
sigma2_as_star <- sigma2_sigma2_e * (post$sd_Phylo__sigma_Intercept^2 /(post$sd_Phylo__sigma_Intercept^2 + sigma2_fs))
# Equation 20 - sigma2_sigma2_p (variance of the total phenotypic variance)
sigma2_sigma2_p = sigma2_p ^ 2 + 3 * sigma2_sigma2_e
# Equation 18
H2s <- sigma2_as_star / sigma2_sigma2_p
#summary(H2)
quantile(H2s, probs = c(0.025, 0.5, 0.975)) %>% pander()| 2.5% | 50% | 97.5% |
|---|---|---|
| 0.1391 | 0.2877 | 0.326 |
Here we found a median phylogenetic heritability of around 29%.
Using Equations 6 in the paper, we can calculate the macro-evolvability (CV_Al) of the trait crange_size (log-transformed trait; the location part), but we first need to back-transform as we did above.
# exp(b0) is the median of the trait
ln_b0 <- post$b_Intercept + mean(log(dat$Range.Size), na.rm = T)
sigma2_al <- post$sd_Phylo__Intercept^2
mean_b0 <- exp(ln_b0 + 0.5 * (sigma2_al + sigma2_e_bar) )
# (exp(sigma_ln^2) - 1) * exp(2 * mu_ln + sigma_ln^2)
sd_al <- sqrt( (exp(sigma2_al) - 1) * exp(2*ln_b0 + sigma2_al) )
CV_Al <- sd_al / mean_b0
quantile(CV_Al, probs = c(0.025, 0.5, 0.975)) %>% pander()| 2.5% | 50% | 97.5% |
|---|---|---|
| 2.556e-16 | 0.3042 | 66.77 |
We can see the median macro-evolvability is around 0.30.
# alternatively, we can get CV on ln-scale
CV_Al2 <- sigma2_al / ln_b0
quantile(CV_Al2, probs = c(0.025, 0.5, 0.975)) %>% pander()| 2.5% | 50% | 97.5% |
|---|---|---|
| 0.4576 | 0.8214 | 1.449 |
If we get the CV on the log-scale, we can see the median macro-evolvability is around 0.82, which is higher than what we got on the original scale (like above).
We now estimate the macro-evolvability of the trait crange_size (the scale part) using Equation 24 (we assume that this CV_As is the same as the value for the original scale):
# Using Equation 24
CV_As <- sqrt(exp(4*post$sd_Phylo__sigma_Intercept^2) - 1)
quantile(CV_As, probs = c(0.025, 0.5, 0.975)) %>% pander()| 2.5% | 50% | 97.5% |
|---|---|---|
| 1.34 | 3.177 | 11.31 |
We get the median macro-evolvability of around 3.2. As above, this is macro-evolvability of variance but we can get macro-evolvability of the standard deviation using Equation 25.
# Alternatively, we can get CV_As on the SD scale
CV_As2 <- sqrt(exp(post$sd_Phylo__sigma_Intercept^2) - 1)
quantile(CV_As2, probs = c(0.025, 0.5, 0.975)) %>% pander()| 2.5% | 50% | 97.5% |
|---|---|---|
| 0.5414 | 0.9083 | 1.539 |
We get the median macro-evolvability of around 0.91, which is not that different from the macro-evolvability of the mean estimated on the log scale.
Additionally, we can use Equation 23:
CV_As <- sqrt(sigma2_as_star) / sigma2_e_bar
quantile(CV_As, probs = c(0.025, 0.5, 0.975), na.rm = T) %>% pander()| 2.5% | 50% | 97.5% |
|---|---|---|
| 1.461 | 3.628 | 14.12 |
This estimation matches the result from Equation 24, unlike the one above. This is because we have now a fixed effect which is taken into account for obtaining the phylogenetic variance on the scale part: sigma2_as_star (we will not try to get the macro-evolvability on the original scale as it does not seem to work).
# putting all plots together
(plots_mod1[[1]] + plots_mod1[[2]]) +
(plots_mod2[[1]] / plots_mod2[[2]] /plots_mod2[[3]] ) +
(plots_mod3[[1]] / plots_mod3[[2]] /plots_mod3[[3]]) + plot_layout(heights = c(1, 3), widths = c(1, 1)) +
# put A, B, C
plot_annotation(tag_levels = 'A')
sessionInfo() %>% pander()R version 4.5.0 (2025-04-11)
Platform: aarch64-apple-darwin20
locale: en_US.UTF-8||en_US.UTF-8||en_US.UTF-8||C||en_US.UTF-8||en_US.UTF-8
attached base packages: stats, graphics, grDevices, utils, datasets, methods and base
other attached packages: gt(v.1.0.0), pander(v.0.6.6), patchwork(v.1.3.0), ape(v.5.8-1), here(v.1.0.1), tidybayes(v.3.0.7), bayesplot(v.1.12.0), brms(v.2.22.0), Rcpp(v.1.0.14), lubridate(v.1.9.4), forcats(v.1.0.0), stringr(v.1.5.1), dplyr(v.1.1.4), purrr(v.1.0.4), readr(v.2.1.5), tidyr(v.1.3.1), tibble(v.3.2.1), ggplot2(v.3.5.2), tidyverse(v.2.0.0) and pacman(v.0.5.1)
loaded via a namespace (and not attached): tidyselect(v.1.2.1), svUnit(v.1.0.6), farver(v.2.1.2), loo(v.2.8.0), fastmap(v.1.2.0), TH.data(v.1.1-3), tensorA(v.0.36.2.1), digest(v.0.6.37), estimability(v.1.5.1), timechange(v.0.3.0), lifecycle(v.1.0.4), StanHeaders(v.2.32.10), survival(v.3.8-3), magrittr(v.2.0.3), posterior(v.1.6.1), compiler(v.4.5.0), sass(v.0.4.10), rlang(v.1.1.6), tools(v.4.5.0), yaml(v.2.3.10), knitr(v.1.50), labeling(v.0.4.3), bridgesampling(v.1.1-2), htmlwidgets(v.1.6.4), curl(v.6.2.2), pkgbuild(v.1.4.7), xml2(v.1.3.8), plyr(v.1.8.9), multcomp(v.1.4-28), abind(v.1.4-8), withr(v.3.0.2), grid(v.4.5.0), stats4(v.4.5.0), xtable(v.1.8-4), colorspace(v.2.1-1), inline(v.0.3.21), emmeans(v.1.11.1), scales(v.1.3.0), MASS(v.7.3-65), cli(v.3.6.4), mvtnorm(v.1.3-3), rmarkdown(v.2.29), generics(v.0.1.4), RcppParallel(v.5.1.10), rstudioapi(v.0.17.1), reshape2(v.1.4.4), tzdb(v.0.5.0), commonmark(v.1.9.5), rstan(v.2.32.7), splines(v.4.5.0), parallel(v.4.5.0), base64enc(v.0.1-3), matrixStats(v.1.5.0), vctrs(v.0.6.5), V8(v.6.0.4), Matrix(v.1.7-3), sandwich(v.3.1-1), jsonlite(v.2.0.0), litedown(v.0.7), hms(v.1.1.3), arrayhelpers(v.1.1-0), ggdist(v.3.3.3), glue(v.1.8.0), codetools(v.0.2-20), distributional(v.0.5.0), stringi(v.1.8.7), gtable(v.0.3.6), QuickJSR(v.1.7.0), munsell(v.0.5.1), pillar(v.1.10.2), htmltools(v.0.5.8.1), Brobdingnag(v.1.2-9), R6(v.2.6.1), rprojroot(v.2.0.4), evaluate(v.1.0.3), lattice(v.0.22-7), markdown(v.2.0), backports(v.1.5.0), rstantools(v.2.4.0), gridExtra(v.2.3), coda(v.0.19-4.1), nlme(v.3.1-168), checkmate(v.2.3.2), xfun(v.0.52), zoo(v.1.8-14) and pkgconfig(v.2.0.3)