Hi Jonathan.
The thing to do in this case is to estimate K with within-species
variability. This is described in Ives et al. (2007; Syst. Biol.) and
implemented in the phytools function phylosig. This will give an
unbiased estimate of K. (K estimated when within-species variability is
ignored is downwardly biased.)
To do this, you need the species means for each species and a vector of
the standard error of the means. If we cannot estimate the standard
error of the mean for some species (because n=1), then we can just use
the mean variance. To get the means and standard errors (assuming your
raw data is in a vector, x, with species names), we can just do the
following:
# get the mean by species
temp<-aggregate(x,by=list(names(x)),mean)
xbar<-temp[,2]; names(xbar)<-temp[,1]
# get the variance by species
temp<-aggregate(x,by=list(names(x)),var)
xvar<-temp[,2]; names(xvar)<-temp[,1]
# replace NA with mean variance
xvar[is.na(xvar)]<-mean(xvar,na.rm=TRUE)
# get the N per species
n<-as.vector(table(names(y)))
# compute the standard errors
se<-sqrt(xvar/n)
# compute K
K<-phylosig(tree,xbar,se=se)
Hopefully that works for you.
All the best, Liam
--
Liam J. Revell
University of Massachusetts Boston
web: http://faculty.umb.edu/liam.revell/
email: liam.rev...@umb.edu
blog: http://phytools.blogspot.com
On 6/6/2012 6:12 PM, Jonathan Benstead wrote:
Dear r-sig-phylo members - Colleagues and I came across an old post about
calculating the Blomberg K statistic in studies where there is more than one
data point for some species, by using a tree that includes zero-branch lengths
for those multiple studies. This seems ideal for our data, which has several
species with up to 14 data points for trait values. We have prepared a tree
with the necessary zero-branch lengths and our tip labels match our row labels.
However, when I run the code I get the following error message:
Error in Initialize.corSymm(X[[1L]], ...) :
Initial values for corSymm must be between -1 and 1
I found two references to this problem among old posts, but nothing that helped
solve the problem, except that it might be something to do with the tree. As
someone else has pointed out, the problem goes away if we change the zero
branch lengths to a positive number, but we don't want to give these branches
an arbitrary length. Is there another way to do this?
Here's our code:
library(ape)
library(picante)
Fish = read.table("Fish.csv", header = TRUE, sep = ",")
AuthP2<-as.matrix(as.matrix(Fish)[,1])
tree1 = read.nexus("consensus_constraint_polytomies.nex")
tree2 = drop.tip(tree1, "Erpetoichthys_calabaricus_AY442348_")
names(AuthP2)<-tree2$tip.label
Kcalc(AuthP2[tree2$tip.label], tree2)
Any help getting this analysis to work would be greatly appreciated.
Jon
Jon Benstead
Associate Professor
Aquatic Biology Program
University of Alabama, Box 870206
1124 Bevill Building, 201 7th Ave.
Tuscaloosa, Alabama 35487-0206
Lab homepage: http://bama.ua.edu/~jbenstead
Iceland blog: http://icelandstreams.blogspot.com/
Tel: 205-348-9034
Fax: 205-348-1403
_______________________________________________
R-sig-phylo mailing list
R-sig-phylo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-phylo
_______________________________________________
R-sig-phylo mailing list
R-sig-phylo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-phylo
