Phylogenetic signal in brms

Interfaces brms
1

Hello, I have a question of clarification. When constructing a model using the phylogenetic method described in the brms vignette (https://cran.r-project.org/web/packages/brms/vignettes/brms_phylogenetics.html), does the variance covariance matrix assume a phylogenetic signal that perfectly mirrors Brownian motion (the equivalent of lambda = 1), or does it weight the VCV matrix by the estimate of the phylogenetic signal that you can extract using the hypothesis method? In other words, if I wanted to simulate a PGLS, would I need to multiply the VCV by an estimate of lambda, or does the covariance among the residuals that is used to fit the model reflects an estimate of the phylogenetic signal?

Thank you,

Ethan

2 Likes
Heritability from betabinomial GLMM variance estimates using pedigree data
2

Maybe @Ax3man knows the answer.

3

As far as I know brms is estimating the “phylogenetic mixed model”. A good reference is this Am Nat paper by Housworth, Martins & Lynch: https://doi.org/10.1086/380570. From that paper:

We note that the univariate phylogenetic heritability estimator is mathematically equivalent to the phylogenetic correlation estimator lambda examined recently in Freckleton et al. (2002) despite those authors’ claim to the contrary. The Freckleton et al. article makes a very nice companion to our current work because those authors determined the value of lambda & h2 for numerous real phylogenies and data sets from the biology literature.

The h2 (heritability) they talk about there is:

In quantitative genetics, the heritability of a trait is the proportion of the variance in the trait explained by the relationship between individuals as given in the pedigree. Similarly, we define the phylogenetic heritability of a character as the proportion of variance in the character explained by the relationship among taxa as given by the phylogeny. In both cases, the mathematical formula is h2 = sigma_a ^ 2 / (sigma_a ^ 2 + sigma_e ^ 2 ).

See also here: Phylogenetic signal for binomial and lognormal families.

So the model does not assume perfect Brownian motion (lambda = 1) but scales the effect of the phylogeny the same way as Pagel’s lambda model does. The phylogenetic signal explained in the vignette is the phylogenetic heritability that Housworth et al. are talking about, and the same as lambda in a PGLS.

edit:

For completeness, here is how the referenced Freckleton et al. paper (https://doi.org/10.1086/343873) argues that lambda and phylogenetic h2 are not the same:

For traits evolving along the branches of a phylogeny entirely by Brownian motion, the expected values of both lambda and H2 are equal to 1.0. Conversely, both measures are expected to be equal to 0 for traits that have no phylogenetic component. However, the two measures are not identical since lambda is not interpretable as the percentage of variance attributable to phylogenetic effects. This difference between the two parameters arises because one optimizes its estimates on the basis of partitioning the data into the phylogenetic and nonphylogenetic components, whereas the other optimizes its estimates on the total or unpartitioned variance in the trait. Thus, l is defined as the transformation of the phylogeny that makes the unpartitioned trait data best fit the Brownian motion model. This means that lambda can also adopt values > 1.0. A value of lambda > 1.0 can arise if, for instance, traits are more similar than that predicted by Brownian motion, given the phylogeny. As an example, the value of H2 is 1.0 for Lynch’s (1991) data on body size, but lamba = 1.24 . In practice, however, the range of values > 1.0 that lambda can adopt is restricted since the off-diagonal elements in V cannot be larger than the diagonal elements.

So lambda is not constrained to be between 0 and 1, but h2 is. (In my personal and limited experience, lambda > 1 doesn’t come up that often, and you’re more likely to see lambda < 0. Lamda < 0 could be interpreted as closely related species being more different than unrelated species, perhaps because of displacement or competition. Also note that some implementations of PGLS will estimate lambda but restrict it to [0, 1] without telling you.)

5 Likes
4

Thank you! That was very helpful.

1 Like
5

Apologies for resurrecting this post. I just looked through the stancode generated from the brms model given by resp ~ 1 + expl + (1 + expl | gr(sp, cov=phy)) and the way it uses Lcov (the Cholesky decomposition of the phylogenetic variance-covariance matrix) makes me think that there is no Pagel’s lambda or any equivalent—it just uses the covariance matrix assuming \lambda=1.

For instance, for the species-level intercept it calculates vector[N_sp] r_sp_intercept = sd_sp_intercept * (Lcov * z_sp_slope) (I have changed the notation slightly to make this more readable, originally the variable name was r_1_1). If there were a Pagel’s lambda in the model already, there it would be a line above defining matrix[N_sp, N_sp] Lcov = cholesky_decompose(lambda * cov_nodiag + cov_diag)) instead of having it provided in the data block. Of course, if \lambda=1 then the Lcov would be cholesky_decompose(cov_nodiag + cov_diag), which is exactly the Lcov in the brms-generated model.

Am I reading this right? It uses Lcov directly from the data block without any transformation, which makes me think there’s no estimation of phylogenetic signal.

(I’ve just implemented this Pagel’s lambda version of the model, and one chain finished nearly immediately, but the others look like they might continue for days and days. Does anyone know of a more efficient way to implement this?)

6

I’m not super up to speed on the technicalities of lambda values outside of [0,1] or the difference between lambda and h2. Ignoring these details, the lambda model can be thought of as partitioning the variance across species into a phylogenetically structured component defined by the covariance matrix, and a phylogenetically unstructured component.

For Gaussian models with one observation per species and where the phylogenetic effect acts on the intercept, the unstructured component is simply the model residual. Thus, a formula like resp ~ 1 + (1 | gr(sp, cov=phy)) is the Pagel’s lambda model when there is one observation per species.

With multiple observations per species, there are actually two different things that we could hypothetically mean by the Pagel’s lambda model.

  1. We could mean that we have a phylogeny resolved all the way down to the individual level, such that we have a variance component controlled by relatedness including at the intra-specific level, and then a totally unstructured component. Thus, the Pagel’s lambda variance partition is still between a phylogenetically structured part and an observation-level residual. Furthermore, note if the intra-species branches are all very short compared to the inter-species branches, we might approximate this tree by one where the intra-species branches are of length zero, in which case we could use a covariance matrix consisting of the species-level phylogeny.
  2. Alternatively, we could mean that we wish to apply the Pagel’s lambda model to the species means, and then to allow for an additional within-species between-individual variance component. To achieve this, we need to add a species level residual to the model, in addition to the individual level residual. This species-level residual is simply (1 | species), which is pretty neat! Pagel’s lamba here is how the variance is partitioned between the phylogenetically structured random effect of species and the phylogneticaly unstructured random effect of species.

You’re interested in a model that includes phylogentically structured slopes. In this case, the model doesn’t automatically have a variance component at the individual level (akin to the residual) that yields the Pagel’s lambda model automatically. Instead, we need to put it in explicitly with a (expl | sp) random slope term. Now Pagel’s lambda will be the variance partition between the phylogenetically structured random slope term (expl | gr(sp, cov=phy)) and the phylogentically unstructured slope term (expl | sp).

For more, see https://cran.r-project.org/web/packages/brms/vignettes/brms_phylogenetics.html, and particularly the lines dealing with “phylogenetic signal” in that vignette.

2 Likes
7

Thank you @jsocolar, this was really helpful. I remember you left a similar comment on another post about including both a phylogenetically controlled intercept/slope (1|gr(sp,cov=phy)) and a free intercept/slope (1|sp) but I never realized this was to give the model flexibility to fit the intraspecific variation. Presumably anyone trying to do a phylogenetic analysis with multiple representatives for each tip should use this (1|sp) term to avoid overconstraining the model at the intrapsecific level? Or is this decision the equivalent of using a Brownian covariance vs a Pagel’s lambda on Brownian covariance vs using quadratic kernel/OU process vs etc. etc. to fix phylogenetic signal where you choose what model you want to use and explain your choices in the methods? Thanks!

8

There isn’t a strong “should” or “shouldn’t” here. The Pagel’s Lambda model is motivated by the idea that evolution might proceed by brownian motion, but that there might be some component of the variance that isn’t phylogenetic (thus we talk about computing “phylogenetic signal”). This model is phenomenological and may or may not be a good approximation of reality. Nobody has ever claimed that it’s a principled model for trait evolution, but sometimes it’s useful for extracting insights from data.

The random effect of species is a strategy to give some hierarchical structure to the unstructured variance in the intercept, or to inject some unstructured variance into the slopes. When you compute a Pagel’s Lambda style variance partition in the intercept, it’s up to you whether it is more informative to study the variance partition between the structured variance and all unstructured variance (including the intrapsecific variance), or to study the variance partition between the structured variance and the species-level unstructured variance only. In the former case, it’s also up to you whether to imbue the unstructured variance with some taxonomic structure (individuals nested inside of unstructred species-level variance), versus assuming that all of the species-level variance is phylogenetically structured.

The only case where it’s NOT up to you is when there is one observation per species and a Gaussian likelihood. In this setting, the data contain absolutely no information about how to partition variance between a species-level component and an individual-level component.

2 Likes
9

Hi @jsocolar, thank you! this aligns with my understanding of model choice around phylogenetic structure. I appreciate all of your informative answers on this subject!