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.

3

Hey, thanks for the helpful suggestions. I reran the non-correlated model with just 1000 iterations and it finished in 2.25 hours. Then I reran the correlated model and It took a little more than 12 hours to get to the 100 / 1000 [ 10%] (Warmup) stage. It’s still a big difference, isn’t it?. One thing I forgot to mention in the original post is that my dataset has 5218 rows (many species have multiple measurements). If reducing the iterations is the right approach, what would you reccomend, in terms of chains and iterations?

Regarding you suggestion about priors, could you please expand on that? My background is in botany, so not very knowledgeable in hardcore stats. Sorry. In other words, how do I set this LKJ(2) prioir in the model?

Thanks again for your time!

4

In my experience, models with random effects that have thousands of levels are often slow to fit.

One option is to simplify the model by first summarizing the data to the species level (so take the mean per species for the root traits), then fitting:

brms::brmsformula(
  mvbind(t1, t2, t3) ~ state + (1|p|gr(binominal, cov = corrmat))
) + set_rescor(TRUE)

This will of course no longer take into account the uncertainty associated with the species mean, and weigh all species the same. This matters less if the differences between species are (much) larger than the differences within species.

5

Hey, thanks for the suggestion. (This is actually the correlated model you advised in another thread). The issue with that is the other multivaruiate OU models in this work (e.g. Rphylopars) use the total dataset, not the averages. This became important as, so far the hOUwie results show no support for correlated evolution. One reason I can think of is high intra-specific variance and the other is low phylogenetic signal in the traits. Hence I opted to use the whole dataset for the multivariate models.

And after around 5 days of running all the 4 chains have progressed to the 1000 / 10000 [10%] (Warmup) mark. I anticipate it’s going to take around 50 days in total for the model to finish, given that the sampling iterations are as slow as the warmup iterations. Since this is on a cloud VM, I’m thinking of just waiting it out. Thanks for your help anyways!

6

Yes, sometimes it makes sense to just be patient.

Make sure you are using backend = "cmdstanr" with the latest version installed, it can make a substantial difference compared to Rstan. And consider running much fewer than 10k iterations with more chains (as others have suggested).