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!