Example 3: Phylogenetic Signal
We may test for phylogenetic signal in the data by finding the tree transformation parameter lambda
that maximizes the likelihood of the comparative model (Pagel 1997,1999; Revell 2010).
Again we use a simulated phylogeny and dataset with missing observations from the example in Example 1: Getting Started):
library(Rphylopars)
library(phytools) # for simulating pure-birth phylogenies
set.seed(21) # Set the seed for reproducible results
trait_cov <- matrix(c(4,2,2.2,2,3,1.5,2.2,1.5,2.5),nrow = 3,ncol = 3)
tree <- pbtree(n = 20)
To estimate Pagel's lambda, simply include the argument model="lambda"
:
p_lambda <- phylopars(trait_data = sim_data$trait_data,
tree = sim_data$tree,model = "lambda")
p_lambda # Estimated trait covariance and Pagel's lambda
Now we fit a model assuming a star phylogeny, which is the null hypothesis when testing for significant phylogenetic signal under Pagel's lambda (i.e., lambda=0
).
p_star <- phylopars(trait_data = sim_data$trait_data,tree = sim_data$tree,model = "star")
p_star # Estimated trait covariance, fixed Pagel's lambda = 0
Now we use the likelihood ratio to test for phylogenetic signal:
chi_square <- as.double(2*(logLik(p_lambda) - # 2*(logLik_alt - logLik_null)
logLik(p_star)))
degrees_freedom <- p_lambda$npars - p_star$npars # df = difference in # model parameters
pchisq(q = chi_square,df = degrees_freedom,lower.tail = FALSE) # p-value
The test is significant (p<0.05
), suggesting significant phylogenetic signal in the phylogenetic residuals (Revell 2010). NOTE: this method jointly fits a single lambda
parameter for all traits (and their covariances), as opposed to fitting a different lambda
for each trait (Freckleton 2012). When applied to a single trait, this method is identical to the phylosig
function in phytools
, except that 1) phytools
fits models by maximum likelihood whereas Rphylopars
uses restricted maximum likelihood (which can be overridden by setting REML=FALSE
in the phylopars
function), and 2) lambda
is constrained between 0 and 1 in Rphylopars
whereas phytools
does not set the upper limit of Pagel's lambda to 1.