# Multilevel multiresponse phylogenetic mixed models with brms

**URL:** <https://discourse.mc-stan.org/t/multilevel-multiresponse-phylogenetic-mixed-models-with-brms/41411>\
**Category:** Modeling\
**Tags:** loo, brms, specification, techniques\
**Created:** [July 1, 2026, 12:18am UTC](https://discourse.mc-stan.org/t/multilevel-multiresponse-phylogenetic-mixed-models-with-brms/41411 "2026-07-01T00:18:27Z")\
**Posts on this page:** 5\
**Page:** 1

<div class="post-metadata">

**Author:** ![Anoban](https://avatars.discourse-cdn.com/v4/letter/a/df705f/32.png) [@Anoban](https://discourse.mc-stan.org/u/Anoban)\
**Post date:** [July 1, 2026, 12:18am UTC](https://discourse.mc-stan.org/t/multilevel-multiresponse-phylogenetic-mixed-models-with-brms/41411/1 "2026-07-01T00:18:27Z")

</div>

Hi folks!

I’m a 2nd year PhD student studying root ecological strategies and mycorrhizal outsourcing. I’m interested in understanding how the evolution of three response variables (three root traits) is influenced by a categorical trait (mycorrhizal states). I’m relatively new to this space (both phylogenetics and MR-PMM) and would appreciate your guidance on the following aspects:

After reading through some [vignettes](https://paulbuerkner.com/brms/articles/brms_phylogenetics.html)and papers, I’ve put together the following formula:

```r
brms::brmsformula(mvbind(rt1, rt2, rt3) ~ state + (1|gr(species, cov = corrmat)) + (1|p|taxa)) + set_rescor(TRUE)

```

The `rt` variables are the root traits and state is the categorical trait (with 6 states). The first random intercept is the phylogenetically structured trait correlation between species and the second random intercept is the phylogenetically unstructured trait correlations within and between species (`taxa` is identical to `species`). My data has records for 1,300 unique species (5,200 total records - many species have repeated measurements).

First question: `(1|p|taxa)` fits a random intercept over species and accounts for correlations between the three root traits and the random intercepts based on species, as mentioned [here](https://paulbuerkner.com/brms/articles/brms_multivariate.html). Is that correct?

Second question: The variance covariance matrix used here is generated by `ape::vcv.phylo` which uses a BM process while the rest of my analysis is based on OU models. This [discussion](https://discourse.mc-stan.org/t/rescaling-phylogenetic-tree-for-non-brownian-phylogenetic-models-in-brms/29637) is basically the problem I have at the moment. Is there a way to convert my phylogeny into a vcv matrix based on an OU process without knowing how to accurately parametrize `alpha` and `theta`? I have three different continuous traits and how do I decide these two parameters which could vary depending on the traits?

Third question: If I were to also include higher taxonomic ranks in addition to the species names (say genus names and family names) so I could see the (potentially) hierarchical nature of random effects, would the following formula be appropriate? I intend to compare these using the `LOO` method eventually to see how their fits compare to the first one above?

1. 

```r
brms::brmsformula(mvbind(rt1, rt2, rt3) ~ state + (1|gr(species, cov = corrmat)) + (1|p|taxa)) * (1|q|genus) + set_rescor(TRUE)

```

2. 

```r
brms::brmsformula(mvbind(rt1, rt2, rt3) ~ state + (1|gr(species, cov = corrmat)) + (1|p|taxa)) * (1|q|genus) * (1|r|family) + set_rescor(TRUE)

```

Allowing interaction between these new random intercepts with `taxa` since these are taxonomically hierarchically strucured. Any help will be greatly appreciated!. Thanks!

---

<div class="post-metadata">

**Author:** ![Ax3man](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/ax3man/32/11894_2.png) [@Ax3man](https://discourse.mc-stan.org/u/Ax3man)\
**Post date:** [July 6, 2026, 6:54pm UTC](https://discourse.mc-stan.org/t/multilevel-multiresponse-phylogenetic-mixed-models-with-brms/41411/2 "2026-07-06T18:54:32Z")

</div>

> [@Anoban](#):
>
> First question: `(1|p|taxa)` fits a random intercept over species and accounts for correlations between the three root traits and the random intercepts based on species, as mentioned [here](https://paulbuerkner.com/brms/articles/brms_multivariate.html). Is that correct?

Yes. It estimates standard-deviations, reflecting the between-species variation of each root trait, and correlations, reflecting the similarity of root traits across species (but only the non-phylogenetic part).

In most cases, I would recommend also fitting the phylogenetic correlations, like so:

`(1 | p | gr(species, cov = corrmat)) + (1 | q | taxa)`

> [@Anoban](#):
>
> Second question: The variance covariance matrix used here is generated by `ape::vcv.phylo` which uses a BM process while the rest of my analysis is based on OU models. This [discussion](https://discourse.mc-stan.org/t/rescaling-phylogenetic-tree-for-non-brownian-phylogenetic-models-in-brms/29637) is basically the problem I have at the moment. Is there a way to convert my phylogeny into a vcv matrix based on an OU process without knowing how to accurately parametrize `alpha` and `theta`? I have three different continuous traits and how do I decide these two parameters which could vary depending on the traits?

I don’t think there has been any development in supporting OU models, so the discussion you link to is still relevant. `brms` does currently supports known covariance matrices, not covariance structures that depend on additional parameters to be estimated, like OU.

> [@Anoban](#):
>
> Third question: If I were to also include higher taxonomic ranks in addition to the species names (say genus names and family names) so I could see the (potentially) hierarchical nature of random effects, would the following formula be appropriate? I intend to compare these using the `LOO` method eventually to see how their fits compare to the first one above?

I don’t think the `*` notation in your formulas is supported syntax, to get nested effects you’d write something like `(1|p|taxa) + (1|q|taxa:genus) + (1|r|taxa:genus:family)`.

But regardless, I am not sure if this is a useful model to fit? Categorizations of what comprises a genus or family is arbitrary, and the actual expected similarities are captured by the branch lengths of phylogeny. The higher order random effects would start “fighting” with the phylogenetically structured random effect. How would you interpret a better fit for one of these models? I guess it would mean that the evolutionary model for the VCV is misspecified?

---

<div class="post-metadata">

**Author:** ![Anoban](https://avatars.discourse-cdn.com/v4/letter/a/df705f/32.png) [@Anoban](https://discourse.mc-stan.org/u/Anoban)\
**Post date:** [July 7, 2026, 10:38pm UTC](https://discourse.mc-stan.org/t/multilevel-multiresponse-phylogenetic-mixed-models-with-brms/41411/3 "2026-07-07T22:38:50Z")

</div>

Hey, thank you so much for the reply. I’m going to change the phylogenetic random effect to `(1|p|gr(species, cov = corrmat))`, as you advised. What would happen if I estimated the OU parameters `alpha` and `sigma.sq` using other means e.g. `geiger::fitContinuous`, used those values to transform the phylogeny using `geiger:::rescale.phylo` for each continuous trait independently and fitted the same model using a verbose version:

`fm1 <- brms::brmsformula(rt1 ~ state + (1|p|gr(species, cov = corr_rt1)) + (1|q|taxa))`  
`fm2 <- brms::brmsformula(rt2 ~ state + (1|p|gr(species, cov = corr_rt2)) + (1|q|taxa))`  
`fm3 <- brms::brmsformula(rt3 ~ state + (1|p|gr(species, cov = corr_rt3)) + (1|q|taxa))`

where each correlation matrix is branch length transformed based on an OU model, using estimates of `alpha` and `sigma.sq` for each continuous trait, independently, and finally combined these three elements like below:

`model <- brms::brm(fm1 + fm2 + fm3 + set_rescor(TRUE))`

Would the above be a sensible approach?

And thanks for the clarification and you are right about the third question. My supervisor also had similar concerns. No longer planning to include those as the hierarchical relationships I anticipated to result from fitting nested random effects will already be accounted by the phylogeny. Your help is very much appreciated!

---

<div class="post-metadata">

**Author:** ![Ax3man](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/ax3man/32/11894_2.png) [@Ax3man](https://discourse.mc-stan.org/u/Ax3man)\
**Post date:** [July 8, 2026, 7:19pm UTC](https://discourse.mc-stan.org/t/multilevel-multiresponse-phylogenetic-mixed-models-with-brms/41411/4 "2026-07-08T19:19:56Z")

</div>

I would have a few possible concerns with that approach:

- I am not sure if you can model the correlations when the covariance matrices are not the same. I am both not sure if `brms` allows it this way, and I am not sure what the correlations mean when the branch length rescaling is not the same between traits.
- The OU parameters estimated by geiger are not conditioned on an effect of `state`, like they would be in a single complete model. If state has a large effect, then you’d want the model the OU process after accounting for state, which this approach does not do. You could also be interested in whether state affects the OU parameter directly. For this, ou may want to look at e.g. the `OUwie` package, where you can vary the optimum and pull strength of the OU process by different values of `state`. It’s not Bayesian and univariate, but could be relevant for some of your questions?
- Your solution would discard any uncertainty around the OU parameters (with your data size this may be a trivial concern, but OU can be hard to identify).

Perhaps one could try feeding the brms generated stancode for the BM model into an LMM, and ask it to change into an OU model with estimated parameters. You would have to do a bunch of test cases though, e.g. validate the converted Stan model against `OUwie`/`geiger` on simulated data with known alpha before trusting it on the real dataset.

---

<div class="post-metadata">

**Author:** ![Anoban](https://avatars.discourse-cdn.com/v4/letter/a/df705f/32.png) [@Anoban](https://discourse.mc-stan.org/u/Anoban)\
**Post date:** [July 10, 2026, 2:48am UTC](https://discourse.mc-stan.org/t/multilevel-multiresponse-phylogenetic-mixed-models-with-brms/41411/5 "2026-07-10T02:48:26Z")

</div>

Hey, once again thanks for your time and help. I realize that using different `vcv` matrices won’t work, as you mentioned, the correlation structures won’t be identical across all the trait pairs. `OUwie`, particularly `hOUwie` models are already a significant part of my work (working together with one of it’s authors), however, they are univariate (like you said) - which only allow one continuous trait and one categorical trait to be tested for evolutionary correlations. That’s why I started experimenting with MR-PMM models to look into multi trait coordinations. I’m going to just use the BM based `ape::vcv.phylo` and going to explicitly declare that it assumes a different model of evolution, compared to the rest of the analysys. Thanks again!