[R-sig-phylo] Phylogenetic signal
Simon Blomberg
s.blomberg1 at uq.edu.au
Thu Apr 17 05:33:38 CEST 2008
Here's a function I wrote several months ago to calculate K. It's
basically just a direct R translation of the mathematics in the paper,
using GLS. Dan Ackerly has prettied it up quite a bit. I think it was
going into a package at some stage, but Dan et al. were working on that.
I had planned to write a package for the PHYSIG analyses, but well, life
is short...
Cheers,
Simon.
Kcalc <- function(x,phy) {
## Kcalc calculates Blomberg et al's K statistic for continuous
## traits evolving along a phylogeny.
## See: Blomberg, S. P., T. Garland, and A. R. Ives. 2003. Testing for
phylogenetic signal in comparative data: Behavioral traits are more
labile. Evolution 57:717-745.
## For more information, please contact Simon Blomberg:
## S.Blomberg1_at_uq.edu.au
##
## Kcalc wrapper by: David Ackerly, dackerly at berkeley.edu
## Further hacking by SB 7 August 2007 to do some error-checking and to
## handle nonultrametric trees.
require(ape)
if (length(x) != length(phy$tip.label)) stop(
"Data vector and tree contain different numbers of taxa.")
if (!setequal(names(x), phy$tip.label)) warning(
"Taxon names in data vector and names of tip labels do not
match.")
mat <- vcv.phylo(phy)
if (!is.ultrametric(phy)) {
vars <- diag(mat)
weights <- varFixed(~vars)} else weights <- NULL
x <- x[order(phy$tip.label)]
matc <- vcv.phylo(phy, cor=TRUE) # correlation matrix
ntax <- length(phy$tip.label)
# calculate "phylogenetic" mean via gls
fit <- gls(x ~ 1,
correlation=corSymm(matc[lower.tri(matc)],
fixed=TRUE), weights=weights)
ahat <- coef(fit)
#observed
MSE <- fit$sigma^2
MSE0 <- t(x-ahat) %*% (x - ahat)/(ntax-1)
#expected
MSE0.MSE <- 1/(ntax-1) * (sum(diag(mat))-ntax/sum(solve(mat)))
K <- MSE0/MSE / MSE0.MSE
K
}
On Wed, 2008-04-16 at 21:15 +0200, Albert Barberan Torrents wrote:
> Hello,
> I am new in this field and it would be great if you could help me with a
> question. How can I do to calculate phylogenetic signal (K) and its
> significance as described in Blomberg et al. Testing for phylogenetic
> signal in comparative data: Behavioral traits are more labile? Is there
> any working R package or better does phylocom calculate it?
> I found your page really useful for learning R for phylogenetics. Thank you,
>
--
Simon Blomberg, BSc (Hons), PhD, MAppStat.
Lecturer and Consultant Statistician
Faculty of Biological and Chemical Sciences
The University of Queensland
St. Lucia Queensland 4072
Australia
Room 320 Goddard Building (8)
T: +61 7 3365 2506
http://www.uq.edu.au/~uqsblomb
email: S.Blomberg1_at_uq.edu.au
Policies:
1. I will NOT analyse your data for you.
2. Your deadline is your problem.
The combination of some data and an aching desire for
an answer does not ensure that a reasonable answer can
be extracted from a given body of data. - John Tukey.
More information about the R-sig-phylo
mailing list