First of all, I realized that I omitted a 5th and more relevant factor, the gypsum affinity of the plant, which indeed is the factor explaining more variance in the elemental accumulation in the leaves. So we are interested in including in the model:
-Affinity to gypsum (categoric)
-Soil as PC1 of a PCA with soil features
-MAT (Mean Annual Temperature)
-MAP (Mean Annual Precipitation)
-phylogeny (we have a tree)
And answering your questions:
1.-Yes, we have read about brms, but never used Stan before
2.-We go for linear models for the moment
We need Gamma distribution, since our data are chemical concentrations with variance of data increasing with means. The distribution of our data is hence skewed, with most data concentrated around zero values but a long tail towards increased scores.
4.- We discarded GLS and PhyloANOVA since they do not allow including intraspecific variability (different reps per species). We have used GLMMs (TMB) with phylogenetic effects included as random terms where “Species” is nested within “Family”. Family explains circa 40% of variation and affinity to gypsum around 10%. The other variables are less explicative.
However, we have received the recommendation to include the whole phylogeny in our model using a covariance matrix. For that reason we got interested in brms as a feasible approach. However, we are fearing our models may be too complicated, since they include 5 factors.
In fact that has been our starting point, but we feel our data are more complex than the situation explained in the vignette. Concretely they may fix to the very last point, and not developed:
Phylogenetic models with multiple group-level effects
That is why we are searching for help. We wonder whether our analysis is really feasible with brms or not.
I believe the part about “Phylogenetic models with multiple group-level effects” refers to models that don’t just model the phylogeny as a single effect but also allow other parameters to vary with phylogeny, so (1 + x | gr(phylo, cov = A)) type models instead of (1 | gr(phylo, cov = A)). As far as I know, that’s not very commonly done in “run-of-the-mill” comparative analysis.
It sounds like you want something along the lines of this:
m <- brm(
accumulation ~ 1 + soil + MAT + MAP + gypsum_affinity +
(1 | gr(phylo, cov = A)) + (1 | species) ,
data = your_data, family = 'gamma'
)
Where A is obtained with ape::vcv.phylo(phylo). This satisfies the recommendation you got to include the whole phylogeny.
If this is not what you want, or if it doesn’t work, try to be specific about what your question is.