# Controlling for phylogenetic signal - what is statistically appropriate?

I am currently collaborating with a fellow PhD student. We are both in the same Biology department, but my collaborator is more of a natural historian, so I am handling the statistical side of things.

He has published a few papers using only Brownian models and Pagel's $$\lambda$$ (with a different collaborator). However there seems to be a plurality of different methods available - Grafen, Blomberg, and Martin, to name a few. OU models seem inappropriate for our data since we have relatively small phylogenies (Cooper 2016). The literature seems to indicate that Pagel's $$\lambda$$ is more robust than Blomberg's $$\kappa$$, and in general an okay method to check for phylogenetic signal. I've also found that some people throw everything and the kitchen sink at their data, and then compare log-likelihoods, AIC, BIC, with log-likelihood ratios for any nested comparisons.

My first question is then, should you have any a priori assumptions of which method will be appropriate for your data?

I am still new to PCM so I used old R scripts and an AmNat paper (from 2019) as my reference. I am also using the same phylogenies as that AmNat paper. From that paper, and from those scripts, only two models were used, a Brownian model (which is essentially $$\lambda$$ = 1, anyways) and an estimate of $$\lambda$$. They compared the two models, chose the more appropriate model by log-likelihood ratio and that was it.

My second question is, shouldn't you always compare your models to a model with fixed $$\lambda$$=0?

For example, I have the following output in R,

#Brownian Model

pglsModel_BM <- gls(sum_dep ~ ContGroup, correlation = corPagel(1, phy = UltTree, fixed = TRUE), data = temp, method = "ML")

#Estimated Lambda Model

pglsModel_E <- gls(sum_dep ~ ContGroup, correlation = corPagel(0.50, phy = UltTree, fixed = FALSE), data = temp, method = "ML")

###Output of comparison is...

Model df      AIC      BIC    logLik   Test  L.Ratio p-value
pglsModel_BM     1  3 528.4584 534.5344 -261.2292
pglsModel_E      2  4 528.0454 536.1468 -260.0227 1 vs 2 2.412992  0.1203

###ANOVA output of preferred model

Denom. DF: 54
numDF  F-value p-value
(Intercept)     1 7.039370  0.0104
ContGroup       1 6.480427  0.0138



And the Brownian model got the go ahead. It seems that Brownian is treated as a null model, but I can't wrap my head around why $$\lambda =0$$ isn't also a null model or the null model. Moving forward with my assumption,

###Adding in a lambda=0 model

pglsModel_0<- gls(sum_dep ~ ContGroup, correlation = corPagel(0, phy = UltTree, fixed = TRUE), data = temp, method = "ML")

###Using anova for model comparison

Model df      AIC      BIC    logLik
pglsModel_0      1  3 524.0962 530.1723 -259.0481
pglsModel_BM     2  3 528.4584 534.5344 -261.2292

###checking out the 0 model...

Denom. DF: 54
numDF   F-value p-value
(Intercept)     1 289.67228  <.0001
ContGroup       1   0.21234  0.6468


I realize the differences are rather small, but all methods point to $$\lambda =0$$ as the 'more likely' model.

My third question is -- although, possibly answered by the second question -- do we assume that there must be some amount of phylogenetic signal due to shared history, such that it is suffiecient to only test for $$\hat{\lambda}$$ and compare to $$\lambda = 1$$?

Note: Just to be sure, I compared all Brownian models with corBrownian to their $$\lambda = 1$$equivalents, and got the exact same outputs.

I have a tentative answer to my question.

In the publication for the R package phylosignal (Keck, 2016), they state:

To test the presence of phylogenetic signal, the null hypothesis is that trait values are randomly distributed in the phylogeny. Another null hypothesis might be that trait values follow a Brownian motion model but it is less often used and implemented.

So it would seem there are indeed two null hypotheses, however it seems disingenuous to not test for both, especially since it is relatively simple to do. And of the two, should you for some reason be limited to one, testing for $$\lambda = 0$$ (or the log-likelihood of a general linear model) should be your first choice.

Another issue was with corPagel from the package ape. It requires an initial value to estimate Pagel's $$\lambda$$ (Unless you fix the value, of course). Convergence is not guaranteed, and sometimes requires some fine-tuning of the initial value. This set off a bit of a red-flag, so I collected the median lambda from 5000 subsets of my data (I am using a different family/phylogeny/dataset, where the $$\hat{\lambda}$$ model is preferred over a Brownian, but not $$\lambda = 0$$).

I think this is stronger argumentation for a phylogenetic signal, albeit a weak signal.

#output of comparing lambda=0, lambda=1, and median lambda

Model df      AIC      BIC    logLik
pglsModel_0          1  3 636.9876 643.6899 -315.4938
pglsModel_1          2  3 681.9784 688.6807 -337.9892
pglsModel_MEDIAN     3  3 636.4812 643.1836 -315.2406

#output of comparing median lambda to the estimated lambda from the full dataset

Model df      AIC      BIC    logLik   Test   L.Ratio p-value
pglsModel_MEDIAN     1  3 636.4812 643.1836 -315.2406
pglsModel_FULL       2  4 637.7024 646.6388 -314.8512 1 vs 2 0.7788265  0.3775


By the evidence of the distribution of lambdas, and the median model's AIC, BIC, and log-likelihood, we can argue in favor of the median model. Residuals also look normal and random. Confidence intervals should also be obtained. But, in reality, here, the difference here between a weak phylogenetic signal and none at all are vanishingly small.