Meta-analytical model comparing relationship between two effect sizes

I have been fitting meta-analytical models that compare the relationship between paired effect sizes (log response ratios).

So for context, I am assessing the relationship between two different ecological response metrics (abundance responses and body condition responses). I’m interested in whether body condition responses can be used to predict abundance responses.

So for each species within each dataset I have the log response ratio and associated variance for abundance (AB_LRR, AB_LRR_se) and the log response ratio and associated variance for body condition (BC_LRR, BC_LRR_se).

The model structure I’ve been using is:

mod1 ← brm(AB_LRR | se(AB_LRR_se) ~ me(BC_LRR, BC_LRR_se) + (1|study_ID) + (1|species)

I’m also interested in investigating whether some covariates may influence this relationship (e.g. habitat type).

mod2 ← brm(AB_LRR|se(AB_LRR_se) ~ me(BC_LRR, BC_LRR_se)*habitat_type + (1|study_ID) + (1|species)

My first question is whether the me() function is adequate to use for the body condition effect size variance. If not, is there a function that is appropriate?

Secondly, I seem to be getting some weird results in my interactions models (e.g. mod2). Say I have 12 habitat types but only 6 of them have large sample sizes. If I run the model on all of the data I’m getting positive slopes between body condition and abundance, but if I subset to only run the model on the 6 habitats with large sample sizes, I’m suddenly getting all negative slopes. So presumably something funky is going on with my model structure.

Any insight would be hugely appreciated!

Operating System: Windows
brms Version: