Hi all,

I’m having trouble with the running speed of my phylogenetic models and I wondered if people had some advice. The models I’m running are for 9835 species, so they are expected to take a very long time, but unfortunately I’m only allowed 72 hours on my university cluster in one go.

For anyone that doesn’t really use phylo-models, the main speed decrease comes from using a covariance matrix for the random species level effect. For a few species (~1000) it only takes about 30 mins, but often with the full 10k half the models don’t finish in the allotted 72 hours. I have to run a sample of 50 trees, and repeat for sensitivity analysis as well.

Even with 10,000 iterations I’m still getting Rhat values of 1.01, normally for the group level effect, or else I’d just run for less iterations.

So my question: is there anything I can do to speed up this model?

- Already run in parallel with 2 chains, 25 threads per chain so 50 cores total (benchmarking showed even more cores would’ve been quicker but resources are limited).
- Already randomised data so that each chunk should be similar processing speed (I hope).
- Using decomp = “QR”, which does it speed it up a bit.
- Init = 0 seems to work better than random.
- All variables are centered and scaled to the same units already.
- Priors are weakly informative. I didn’t want to set N(0,1) because I get an effect size of ~3 for one trait. However I get the feeling that with this much data the priors probably won’t affect the posterior. Plus with this much data I doubt they would speed it up, but I did put a stronger prior on the group level effect.
- pp_check is really good fit, so seems like models are parametrized correctly.
- I know using “threshold = equidistant” slows down my model, but it fits better if I do.

This is perhaps a rogue suggestion, but could I try and decrease the max_treedepth? I know increasing it slows down sampling speed, but my models are converging. Or should I try and increase tree_depth, but then have less iterations? Same applies for adapt_delta. Benchmarking speed of a sample shows they speed up models, but also looked like convergence was slightly worse.

Should I increase warmup? I guess I don’t need all those iters? Thinning is 20 for file size afterwards.

Anything else you can think of? Using Max Likelihood didn’t speed it up so i stopped trying with that.

```
###############################################################################
# Ordinal brms models on cluster #
###############################################################################
# Packages to load.
library(magrittr)
library(caper)
library(dplyr)
library(effectsize)
library(phytools)
library(brms)
library(graph4lg)
# Clear the workspace.
rm(list=ls())
###############################################################################
#### Read in the data #####
# Functions.
source("Code/functions.R")
# Read in the tree.
model_tree <- read.tree("Data/Trees/prum_trees.tre")[[tree_number]]
# Read in the life history traits.
model_data <- read.csv("Data/eco_traits.csv")
model_data$tree_tip <- gsub(" ", "_", model_data$birdtree_name)
# Drop tips on the tree.
model_tree <- drop.tip(model_tree, setdiff(model_tree$tip.label, model_data$tree_tip))
# Make a covariance matrix, and order data the same.
model_covar <- ape::vcv.phylo(model_tree)
# Reorder the matrix so it's random, to maximise parallel processing speed.
mat_order <- sample(1:nrow(model_covar), size = nrow(model_covar), replace = FALSE)
model_covar <- reorder_mat(model_covar, rownames(model_covar)[mat_order])
row.names(model_data) <- model_data$tree_tip
model_data <- model_data[row.names(model_covar),]
###############################################################################
#### Prepare predictor variables ######
# Set as factor, then re-level for appropriate reference group.
model_data %<>% mutate(
trait_1 = relevel(as.factor(trait_1), ref = "No"),
trait_2 = relevel(as.factor(trait_2), ref = "Weak"),
trait_3 = relevel(as.factor(trait_3), ref = "Primary"),
)
# Center categorical predictors.
model_data %<>% mutate(
trait_1 _c = center_categorical(trait_1),
trait_2 _c = center_categorical(trait_2),
trait_3 _c = center_categorical(trait_3),
)
# Scale continuous predictors to two SD.
model_data %<>% mutate(
trait_4_z = standardize(trait_4, two_sd = TRUE),
)
###############################################################################
#### Set model formula ######
# Model forumla.
# (3 categorical dummy variables, 1 continuous, 2 interaction terms, group level effect).
model_formula <- "response_score ~
trait_1_c + trait_2_c +
trait_3_c + trait_4_z
trait_1_c:trait_3_c +
trait_1_c:trait_4_z+
(1|gr(tree_tip, cov=A))"
# brms formula.
brms_formula <- brmsformula(model_formula,
family = cumulative(threshold = "equidistant"),
decomp = "QR")
# Add weak priors.
normal_priors <- c(prior(normal(0,5), class="Intercept"),
prior(normal(0,5), class="b"),
prior(gamma(2,1), "sd"))
# Model pathway.
model_pathway <- paste0("Results/Models/, example_model, ".rds")
# Run models.
brms_model <- brm(
brms_formula,
data = model_data,
data2 = list(A=model_covar),
prior = normal_priors,
iter = 10000,
warmup = 5000,
chains = 2,
thin = 20,
cores = 50,
init = 0,
file = model_pathway,
normalize = FALSE,
backend = "cmdstanr",
control = list(max_treedepth = 5),
threads = threading(25),
)
### OS + BRMS INFO ###
> Sys.info()
sysname = "Linux"
release = "4.18.0-348.20.1.el8_5.x86_64"
version = "#1 SMP Tue Mar 8 12:56:54 EST 2022"
machine = "x86_64"
> sessionInfo()
R version 4.0.3 (2020-10-10)
Platform: x86_64-conda-linux-gnu (64-bit)
Running under: Red Hat Enterprise Linux 8.5 (Ootpa)
Matrix products: default
BLAS/LAPACK: /rds/general/user/home/anaconda3/envs/env1/lib/libopenblasp-r0.3.12.so
locale:
[1] LC_CTYPE=C.UTF-8 LC_NUMERIC=C LC_TIME=C.UTF-8
[4] LC_COLLATE=C.UTF-8 LC_MONETARY=C.UTF-8 LC_MESSAGES=C.UTF-8
[7] LC_PAPER=C.UTF-8 LC_NAME=C LC_ADDRESS=C
[10] LC_TELEPHONE=C LC_MEASUREMENT=C.UTF-8 LC_IDENTIFICATION=C
attached base packages:
[1] stats graphics grDevices utils datasets methods base
other attached packages:
[1] brms_2.18.0 Rcpp_1.0.8.3
```