Problems fitting multivariate model with option set_rescor

Interfaces brms
1

Hello,
I am struggling to fit a trivariate (i.e. three response variables) mixed model, one of whose random/group effects includes a phylogeny. For reasons that I can explain if needed, I need to obtain the residual sd and correlations of the response variables. The problem is that whenever I set set_rescor(T), the model estimates the sd of the random/group effects to be equal to zero.
When I set set_rescor(F), the model fits very well but I do not obtain the residual correlations of the response variables that I need.

These are the three identical formulae for each response variable: 
formula.LeafN3 = bf(log(SpsLeaf.N.)~Presin.s+pH.H2O.s+ Clay.s+ Corg.s + (1|Morphology)+ (1|Site.Plot)+  (1+Presin.s|Spsabrv.)+(1+Presin.s|p|gr(Spsname1, cov=inv.phylomatrix)))
formula.LeafP3 = bf(log(SpsLeaf.P.)~Presin.s+pH.H2O.s+ Clay.s+ Corg.s + (1|Morphology)+ (1|Site.Plot)+ (1+Presin.s|Spsabrv.)+(1+Presin.s|p|gr(Spsname1, cov=inv.phylomatrix)))
formula.LeafN.P3 = bf(log(SpsLeafN.P)~Presin.s+pH.H2O.s+ Clay.s+ Corg.s + (1|Morphology)+ (1|Site.Plot)+ (1+Presin.s|Spsabrv.)+(1+Presin.s|p|gr(Spsname1, cov=inv.phylomatrix)))

And the model is:
multiv5= brm(formula.LeafN3+formula.LeafP3+formula.LeafN.P3  + set_rescor(T), data=DF,
             family=gaussian(), data2=list(inv.phylomatrix=inv.phylomatrix),
             chains=3,  warmup = 1e3, thin=20, iter=1e5, cores = 3, prior=prior.m1, 
             control = list(adapt_delta = 0.9,max_treedepth=15),
             threads = threading(2, grainsize = 100))
The priors I set were:
prior.m1=prior(normal(0,1), class = "b")+
         prior(normal(0,2), class = "Intercept")+
         prior(lkj(0.5), class="cor")

All explanatory variables were standardized to mean=0 and sd=1 prior to the analysis. The dataset has 207 rows.
Before posting this message, I have tried many simplifications of the above model by deleting many fixed/population effects, by reducing the group/random effects, by increasing the warmup, the number of iterations (up to 10⁶), and by changing a few other parameters involved in the HMC algorithm such as the adapt_delta and max_treedepth). I am also trying to fit this model on a cluster of the Faculty of Agronomy and, after 5 days, the results are seemingly identical to those I am describing. I’d then guess that it is not a matter of “numerical brute force”.
In all cases, I obtained the same frustrating result oof not being able to estimate either the sd of the group/random effects (when I set set_rescor(T)), or the residual correlations (when I set_rescor(F)).
I would appreciate any suggestions that you might have.
Cheers

Pablo

Please also provide the following information in addition to your question:

  • Operating System: Linux ubuntu 22.04
  • brms Version:2.22_0
  • R version 4.5.0
    My machine is a Pentium I7 with 32GB RAM and the processor is AMD® Ryzen 5 2600 six-core processor × 12

Looking forward to your topic!

2

Do you have multiple observations per species?

You write cov=inv.phylomatrix, but brms expects the phylogenetic VCV, not its inverse.

If you are not getting warnings about the fitting it is probably an issue with the model specification rather than the estimation itself.