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 and papers, I’ve put together the following formula:
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. 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 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?
Allowing interaction between these new random intercepts with taxa since these are taxonomically hierarchically strucured. Any help will be greatly appreciated!. Thanks!
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:
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.
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?
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:
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!
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.
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!