[R-sig-phylo] Phylogenetic signal

Steve Kembel skembel at berkeley.edu
Thu Apr 17 06:44:28 CEST 2008


Dear Albert,

To complement Simon Blomberg and David Ackerly's Kcalc function, I  
wrote some R code for the test of phylogenetic signal as the variance  
of standardized contrasts vs. a null model of shuffling trait values  
across the tips of the phylogeny (Blomberg et al. 2002). The function  
is called phylosignal, it's included in the R package picante which  
will be released publicly in the next few weeks. In the meantime here  
are the functions and you can also try installing a prerelease version  
of the package using the following commands from within R.

Once you've installed it, the function "phylosignal" will report the K  
statistic and P-value based on variance of contrast values vs. a tip  
shuffling randomization. The documentation for the package is in flux  
so check back soon for more details. In particular be careful since  
the code currently assumes the trait data are sorted in the same order  
as the phylogeny$tip.labels.

Best,
Steve

Link for picante package:
http://picante.r-forge.r-project.org/

Install code:
install.packages("picante",repos="http://R-Forge.R- 
project.org",dependencies=TRUE)
library(picante)
help(phylosignal)
randtree <- rcoal(20)
randtraits <- evolve.brownian(randtree)
phylosignal(randtraits[randtree$tip.label],randtree)

Raw phylosignal code:
Kcalc <- function(x,phy) {
	mat <- vcv.phylo(phy, cor=TRUE) # correlation matrix
	ntax = length(phy$tip.label)
	ntax1 = ntax-1

	dat = data.frame(x)
	names(dat) = 'x'
	# calculate "phylogenetic" mean via gls
	fit <- gls(x ~ 1, data = dat, 	
		correlation=corSymm(mat[lower.tri(mat)], fixed=TRUE))
	ahat <- coef(fit)
	
	#observed
	MSE <- fit$sigma^2
	MSE0 <- t(dat$x - ahat) %*% (dat$x - ahat)/ ntax1

	#expected
	MSE0.MSE <- 1/ ntax1 * 	
		(sum(diag(mat))- ntax/sum(solve(mat)))

	K <- MSE0/MSE / MSE0.MSE
	return(K)
}

pic.variance <- function(x,phy,scaled=TRUE) {
	pics <- pic(x,phy,scaled)
	N <- length(pics)
	sum(pics^2) / (N -1)
}

phylosignal <- function(x,phy,reps=999,...) {

     K <- Kcalc(x,phy)

     if (!is.vector(x)) {
         x.orig <- x
         x <- as.vector(x)
         names(x) <- row.names(x.orig)
     }
	
	obs.var.pic = pic.variance(x,phy,...)
	
	rnd.var.pic <- numeric(reps)
	
	#significance based on tip shuffle
	for (i in 1:reps) {
		tipsh.x <- sample(x)
	    names(tipsh.x) <- names(x)
		rnd.var.pic[i] <- pic.variance(tipsh.x,phy,...)
	}

	var.pics = c(obs.var.pic,rnd.var.pic)
	var.pics.p = rank(var.pics)[1] / (reps + 1)
	var.pics.z = (obs.var.pic - mean(var.pics))/sd(var.pics)
      
data 
.frame 
(K 
,PIC 
.variance 
.obs 
= 
obs 
.var 
.pic,PIC.variance.rnd.mean=mean(var.pics),PIC.variance.P=var.pics.p,  
PIC.variance.Z=var.pics.z)

}



More information about the R-sig-phylo mailing list