[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