Specifying a nested phylogenetic model

Interfaces brms
1

Hi there,

I am trying to correctly specify a nested phylogenetic model in brms where the nested random term needs to be constrained by the phylogeny. i.e. within my model I have the term +(1|Region/Species) and cov_ranef=list(Species = A) where A is the phylo structure (I followed the phylo brms vignette).

Where my concern arises is: if I run the phylogenetic signal hypothesis test I have to add in sd_Region:Species__Intercept etc for it to work on my nested structure but as I haven’t specified the nested structure in the cov_ranef argument I’m convinced I must be doing something wrong. I’ve tried inputting it as cov_ranef=list(Region:Species = A) but this doesn’t work.

Any help/advice would be much appreciated as I trying to carefully learn how to use Bayesian GLMM methods correctly.

Thanks in advance,
Liam

Please also provide the following information in addition to your question:

  • Operating System: Windows 10
  • brms Version: 2.3.1
2

The Region:Species needs its own phylo matrix of the right dimension. I am not an expert in these kinds of models, so not sure how to best expand A to be reasonable for Region:Species. To then specifiy that in cov_ranef, go for cov_ranef=list("Region:Species" = B), where B is the matrix you created.

3

Thank you for your reply, Paul, I really appreciate it. I guess the most simplistic way to expand it would be a matter of specifying the rownames of ‘B’ to match the region:species grouping levels (given it is nested it should have the same length? or there is at least no overlap in species). Is there a way to extract the grouping levels as a list to match the species rownames? I hope that makes sense :)

4

The levels of region:species are given by <region>_<species>. If you have already fitted you model, you can see the levels in <model>$data[["region:species"]].

5

Following your advice, I re-labelled the phylogenetic matrix with the region_species labels. However when I attempt to re-fit the model, it returns the following error:
Error: Covariance matrix of grouping factor 'Region:Species' is not symmetric.

This is confusing as the original ‘A’ matrix wasn’t symmetric when fitted with a non-nested Species term so I don’t understand why it would need to be symmetric for Region:Species but not Species. I am definitely out of my depth here…

6

These matrices always need to be symmetric. I am sure, A was symmetric before (after all its a correlation matrix, right?).

7

Apologises, this is my error. If I run isSymmetric(A, check.attributes = FALSE) this returns as TRUE for the original but as I ordered the matrix so Species was alphabetical so I could easily change the Species names to Region_Species names, it must’ve messed with the symmetry. I will try some other methods of changing the names and hopefully will find a solution.

8

Ok so I used the following code to re-gig the original correlation matrix to fit a nested structure

Data frame = bee_mean
Phylogeny = bee.tree
Phylogenetic correlation matrix set up following brms phylogenetic vignette = B

  1. Join region and species into single character string called Rspecies that matches <Region_Species>
    bee_mean$Rspecies=paste(bee_mean$Region,bee_mean$Species,sep="_")

  2. Match the order of Rspecies to the bee.tree tip order (as it is in the same order as B)
    region_species=unique(bee_mean[ order(match(bee_mean$Species, bee.tree$tip.label)), ]$Rspecies)

  3. Make rownames of
    rownames(B)=unique(region_species)

  4. Check matrix symmetry
    isSymmetric(B,check.attributes=FALSE)
    [1] TRUE

I do wonder if there is anything I am overlooking as I am new to Bayesian GLMM methods and any input regarding if it is statistically plausible to specify a phylogenetic model in this way would be much appreciated!