Extremely slow phylogenetic mixed model in brms

Modeling
, ,
1

Hi folks,

I have the following multiresponse phylogenetic mixed model in brms :

model <- brms::brm(brms::brmsformula(mvbind(t1, t2, t3) ~ state + (1|p|gr(binominal, cov = corrmat)) + (1|q|taxa)) + set_rescor(TRUE), data = fred4, data2 = list(corrmat = corrmat), chains = 4, cores = 4, threads = 6, iter = 10000)

The issue with this is that after 2 days and 15 hours, it hasn’t even progressed past the 1st iteration. Pasting the output below:

SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 1).

SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 2).


SAMPLINGSAMPLING FOR MODEL ' FOR MODEL 'anon_modelanon_model' NOW (CHAIN ' NOW (CHAIN 34).
).
Chain 3:
Chain 3: Gradient evaluation took 1.26292 seconds
Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 12629.2 seconds.
Chain 3: Adjust your expectations accordingly!
Chain 3:
Chain 3:
Chain 1:
Chain 1: Gradient evaluation took 1.61108 seconds
Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 16110.8 seconds.
Chain 1: Adjust your expectations accordingly!
Chain 1:
Chain 1:
Chain 2:
Chain 2: Gradient evaluation took 1.47388 seconds
Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 14738.8 seconds.
Chain 2: Adjust your expectations accordingly!
Chain 2:
Chain 2:
Chain 4:
Chain 4: Gradient evaluation took 1.78846 seconds
Chain 4: 1000 transitions using 10 leapfrog steps per transition would take 17884.6 seconds.
Chain 4: Adjust your expectations accordingly!
Chain 4:
Chain 4:
Chain 3: Iteration:    1 / 10000 [  0%]  (Warmup)
Chain 4: Iteration:    1 / 10000 [  0%]  (Warmup)
Chain 2: Iteration:    1 / 10000 [  0%]  (Warmup)
Chain 1: Iteration:    1 / 10000 [  0%]  (Warmup)

I don’t believe this model is that complex to be taking this much time? Is it?

  • My phylogeny has 1301 plant species
  • t1, t2 and t3 are logged and standardized continuous traits
  • state is a categorical trait with 6 different characters
  • taxa and binominal are identical columns specifying the species names

The following model finished within 43 hours on the same machine:

brms::brm(brms::brmsformula(mvbind(t1, t2, t3) ~ state + (1|gr(binominal, cov = corrmat)) + (1|taxa)) + set_rescor(TRUE), data = fred4, data2 = list(corrmat = corrmat), chains = 4, cores = 4, threads = 6, iter = 10000)

The only difference between the slow one and this one is that the slow one treats the random effects as correlated ((1|p|gr(binominal, cov = corrmat)), ((1|q|taxa)) whereas the other one did not. Is this enough to make the runtime this long? Is this to be expected? Or is there anything I’m doing wrong? Is there a way to speed up the execution (setting priors, using more cores, more chains…?)?

These models were run on an Ubuntu 24.04.4 LTS VM with 32 vCPUs and 128 GBs of RAM. Installed R version is 4.6.1 and brms version is 2.23.0.

I’d appreciate any help with this. Also please feel free to ask if more info is needed. Thanks!

2

My first thought is to question whether you really need 10,000 iterations per chain? Unless you’re getting convergence issues due to poor warmup, it’s often a better idea to run more, shorter chains to get the N_{eff} you’re looking for. I’d first try seeing how well shorter chains work for your non-correlated model; cutting to two chains & increasing the threading may also help during debugging. It’s really hard to work with alternate versions of these models if they take two days to run. If cutting the the iterations still takes a long time, I’d think about testing on a subset of the phylogeny.

It looks like you’re going with the default LKJ(1) prior distribution for your random effect correlations. That is equivalent to a uniform prior across all valid correlation matrices, which is probably a bit stronger than you want to go with. Your uncorrelated model is effectively fixing your prior correlation to an identity matrix; the LKJ(\eta) distribution approaches this as \eta gets larger. I would consider trying a tighter prior here; I’ve found LKJ(2) is a decent place to start for light regularization, though there are other folks here who understand this distribution better than me.