Hi Stacey,
you can try this function, which is analogous to Blomberg's K for
discrete traits and follows Maddison & Slatkin 1991. An example is
available at the bottom of the syntax.
Cheers,
Enrico
################################## Function phylo.signal.disc
########################################
# This function tests for phylogenetic signal with categorical traits.
# It works similarly to phylo.signal, by randomizing the tip data and
comparing the minimum number of character state changes with a null model.
# Minimum character state change is obtained with parsimony, and the
syntax allows for different evolutionary models.
# To build a matrix of costs of character state transition, see Maddison
& Maddison 2000. MacClade 4 Manual pp.69-72 (unordered parsimony is the
default, cost=NULL).
# Note that this function corresponds to the "Fixed Tree, Character
Radomly Reshuffled" model proposed in Maddison & Slatkin (1991)
Evolution 45:1184.
library(ape)
library(vegan)
library(picante)
library(phangorn)
library(PHYLOGR)
'phylo.signal.disc' <-
function(trait,phy,rep = 999,cost=NULL)
{
lev <- attributes(factor(trait))$levels
if (length(lev) == length(trait))
stop("Are you sure this variable is categorical?")
if(is.null(cost)){
cost1 <- 1-diag(length(lev))
}
else {
if (length(lev) != dim(cost)[1])
stop("Dimensions of the character state transition matrix do
not agree with the number of levels")
cost1<- t(cost)
}
dimnames(cost1) <- list(lev,lev)
trait <- as.numeric(trait)
attributes(trait)$names <- phy$tip
NULL.MODEL <- matrix(NA,rep,1)
obs <- t(data.frame(trait))
obs <- phyDat(t(obs),type="USER",levels=attributes(factor(obs))$levels)
OBS <- parsimony(phy,obs,cost=cost1)
for (i in 1:rep){
null <- sample(as.numeric(trait))
attributes(null)$names <- attributes(trait)$names
null <- t(data.frame(null))
null <-
phyDat(t(null),type="USER",levels=attributes(factor(null))$levels)
NULL.MODEL[i,]<-parsimony(phy,null,cost=cost1)
P.value <- sum(OBS >= NULL.MODEL)/(rep + 1)
}
par(mfrow=c(1,2))
hist(NULL.MODEL,xlab="Transitions.in.Randomizations",xlim=c(min(c(min(NULL.MODEL,OBS-1))),max(NULL.MODEL)+1))
arrows(OBS,rep/10,OBS,0,angle=20,col="red",lwd=4)
phy$tip.label <- rep(".",length(trait))
plot(phy,tip.col=trait+10,cex=250/length(trait),font=1)
title("Character states")
par(mfrow=c(1,1))
OUTPUT1 <- t(data.frame(Number.of.Levels =
length(attributes(factor(trait))$levels),
Evolutionary.Transitions.Observed=OBS,Evolutionary.Transitions.Randomization.Median=median(NULL.MODEL),Evolutionary.Transitions.Randomization.Min=min(NULL.MODEL),Evolutionary.Transitions.Randomization.Max=max(NULL.MODEL),P.value))
if(is.null(cost)){
list(.Randomization.Results=OUTPUT1,.Levels=
lev,.Costs.of.character.state.transition.UNORDERED.PARSIMONY = t(cost1))
}
else {
list(.Randomization.Results=OUTPUT1,.Levels=
lev,.Costs.of.character.state.transition.FROM.ROW.TO.COL = t(cost1))
}
}
###EXAMPLE:
tree <- rcoal(10)
trait <- round(rnorm(10,3,1),0)
phylo.signal.disc(trait,tree)
El 2/3/11 9:09 p.m., sd...@duke.edu escribió:
Hi all,
My lab has been reading literature related to the measurement of
phylogenetic
signal and we are wondering about the application of Blomberg's K
statistic to
discrete data and in particular discrete, unordered characters. We
have seen
it applied to the latter in the literature, but it seems like this would
violate basic model assumptions under Brownian motion.
We'd be curious to hear your opinions on using K to measure signal in
discrete
traits and to be pointed to any releveant literature.
Thanks!
Stacey
_______________________________________________
R-sig-phylo mailing list
R-sig-phylo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-phylo
--
************************************************************************
Enrico L. Rezende
Departament de Genètica i de Microbiologia
Facultat de Biociències, Edifici Cn
Universitat Autònoma de Barcelona
08193 Bellaterra (Barcelona)
SPAIN
Telephone: +34 93 581 4705
Fax: +34 93 581 2387
E-mail: enrico.reze...@uab.cat
_______________________________________________
R-sig-phylo mailing list
R-sig-phylo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-phylo
