Heterogeneous residual covariance structure OR zero error multivariate model

Interfaces brms
1

I have a (Gaussian linear) multivariate model of tree traits with a single 8-vector of traits measured on each of about 1000 trees over about 90 locations and 5 species. After adjusting for the effect of study from which the data was collected, I’m more or less interested in the covariance structure of traits as attributed to variability across species, to variability across locations (or species “subpopulations”) within species and to variability across trees within locations within species (3 nested levels). Again, there is one 8-vector observation per tree, so, as you see in my code, below, I will get one tree level random effect vector per observation vector. Thus, I set_rescor(FALSE) as the correlation among traits at the tree level is handled by the tree random effects specification via the |id3| being shared among the 8 traits’ tree random effects. But then this estimates a non-correlated error vector, which, I might be distinguishable (identifiable) from the correlated error (via the tree effects), but mixing is poor and takes about 4.5 hours on my Mac Studio with M1 Ultra.

So, now I want a model that omits the uncorrelated error vectors (as in non-gaussian models), i.e., with zero error, along with the tree level effects (which act as the correlated error) (and with the other effects). I can get this “zero (non-correlated) error” behavior if I omit the tree level effects and set_rescor(TRUE) (results look great (not shown)), but I eventually want to model these effects (whether via tree level effects or residuals) as a function of species (a different covariance component per species). I see how to do that via random effects with (1|id3|gr(sp:loc:tree, by=sp)) (I’ve done this with the other random effects (not shown in code below)), but I don’t see how to get this via a heterogeneous correlated residual structure specification. So, I can get what I want via random effects, but seem to have to settle for the uncorrelated residuals, which, as I mentioned, seem to cause a mixing (perhaps weak identifiability) problem.

So, again, (1) can I specify zero errors or (2) can I specify correlated residuals to vary with the level of species (analogous to the random effects specification (1|id3|gr(sp:loc:tree, id=sp))?

I’ve seen discussions related to (2) but could not make much hay from them. Help appreciated.

Best,

Jay

form1<- bf(mvbind(d13C,Cond,SLA,SS_conc,
                  ST_conc,WD,diameter,mean_rw) ~
               study + (1|id1|sp) + (1|id2|sp:loc) + (1|id3|sp:loc:tree))
fit1<- brm(form1 + set_rescor(FALSE), data=trait.dat, chains=4, cores=4, threads=4,
           iter=8000, warmup=4000, seed=24601, backend="cmdstanr")
  • Operating System: Mac OS Sonoma 14.5
  • brms Version: 2.21.0
1 Like
2

Hi,
not sure I understand your problem correctly, but I think you can fudge the model you are aiming at by modelling the correlations with per-observation (i.e. per row) random effects and setting a fixed but small residual standard error - basically you will be adding a fixed diagonal matrix to the covariance from your random effects. There’s a tradeoff - the smaller the fixed residual error is, the closer you are to your target model but you are also more likely to encounter sampling difficulties.

If that sounds good, fixed residuals error can be achieved in two ways:

  1. via the se addition term
  2. by putting a constant prior on the relevant parameter

Does that sound like an acceptable solution?

3

Thanks, Martin.

I’ve tried both the fixed diagonal approach, using the constant function, as you suggest, as well as using set_rescor(TRUE) without the otherwise explicit random effect in a formula. Both work fine (but can’t get different rescor by group…).

Since my post, I’ve moved to a Stan model that gives me more flexibility and, incidentally, allows me to easily implement hierarchical centering. I’m able to identify (more or less) the uncorrelated (“nugget”) variability from the correlated variability in my Stan model.

Still, as a long-time nlme/lme4 user, I like brms a lot.

Best,

Jay

PS Apologies for delayed reply