Dear all,
I would like to test for phylogenetic signal in a discrete trait with 3 states
(trait “str”) and another discrete trait with 7 states (trait “hab") using
Pagel’s lambda. I have tried to apply the function “fitDiscrete” in geiger to
my data, but not sure whether this is correct or not:
I first compared which of the 3 models (ER, SYM, ARD) fits to the trait data,
and ER model is accepted for “str”, while ARD for “hab". Then I used the same
function but set transform=“lambda” to get the likelihood based on the original
tree, and compared it with the likelihood based on a lambda 0 tree, which was
transformed to a star-like tree using the function “rescale”. For the “str"
trait, the str_ER_lambda$opt says lambda is 1 (original tree), while the
str_ER_lambda_0$opt says lambda approximates to 0 (lambda 0 tree). It’s clear
that “str” shows phylogenetic signal, as indicated by the likelihood ratio
test. But in the other trait “hab”, although the likelihood ratio test rejects
lambda 0 tree, hab_ARD_lambda$opt$lambda is 0.216 (original tree), and
hab_ARD_lambda_0$opt$lambda is 0.545 (lambda 0 tree). What does the lambda
value mean in the returned list of fitDiscrete? Is the lambda value here 0.216
based on the original tree still indicates the “hab” trait with phylogenetic
signal, because the likelihood ratio test says so? Following are the scripts
for testing “str” trait. It would be very helpful if I can receive some
suggestions from you. Also it would be nice to know whether there is any other
method to test phylogenetic signal for categorical traits. Thanks!
# the phylogenetic tree in use
phy2<-phy.chr.str4
phy2_0 <- rescale(phy.chr.str4, model = "lambda", 0) # lambda transform tree to
star like
# load trait data
trait<-read.csv("traitdata.csv",row.names=1)
# lambda 1 and lambda 0 trees combined with trait data
trait_geiger <- treedata(phy2,trait)
trait_geiger_0 <- treedata(phy2_0,trait)
# Comparison of different models for trait evolution
# stratification
str_ER <- fitDiscrete(trait_geiger$phy, trait_geiger$data[,2], type="discrete",
model = "ER", niter = 1000)
str_SYM <- fitDiscrete(trait_geiger$phy, trait_geiger$data[,2],
type="discrete", model = "SYM", niter = 1000)
str_ARD <- fitDiscrete(trait_geiger$phy, trait_geiger$data[,2],
type="discrete", model = "ARD", niter = 1000)
str_d_ER_vs_SYM <- abs(2*(str_ER$opt$lnL-str_SYM$opt$lnL)) # 1.543
str_p_value_ER_vs_SYM <- pchisq(str_d_ER_vs_SYM, 3-1, lower.tail=FALSE) # 0.462
-- Not rejects ER-model
str_d_ER_vs_ARD <- abs(2*(str_ER$opt$lnL-str_ARD$opt$lnL)) # 8.855
str_p_value_ER_vs_ARD <- pchisq(str_d_ER_vs_ARD, 6-1, lower.tail=FALSE) # 0.114
-- Not rejects ER-model
str_ER # we accept the ER model for "stratification"
# phylogenetic signal test for trait stratification, using Pagel's lambda
str_ER_lambda <- fitDiscrete(trait_geiger$phy,trait_geiger$data[,2],
type="discrete", model = "ER", transform = "lambda", niter = 1000) #AICc =
129.778
str_ER_lambda_0 <- fitDiscrete(trait_geiger_0$phy,trait_geiger_0$data[,2],
type="discrete", model = "ER", transform = "lambda", niter = 1000) #AICc =
167.906
str_d_our_tree_vs_lambda_0 <-
abs(2*(str_ER_lambda_0$opt$lnL-str_ER_lambda$opt$lnL)) # 38.127
str_p_our_tree_vs_lambda_0 <- pchisq(str_d_our_tree_vs_lambda_0, 1,
lower.tail=FALSE) # 6.626e-10 -- Rejects no signal model, stratification has
phylogenetic signal
All the best
Ting-Wen
--
Ting-Wen Chen
J.F. Blumenbach Institute of Zoology and Anthropology
Georg August University Goettingen
Berliner Str. 28
D-37073 Goettingen, Germany
Tel: +49-55139-10943
_______________________________________________
R-sig-phylo mailing list - R-sig-phylo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-phylo
Searchable archive at http://www.mail-archive.com/r-sig-phylo@r-project.org/
