Hi all.
Here is some code that I think will do what Dr. Felsenstein has
suggested for a single character. It returns a vector containing the
simulated states at all internal or terminal nodes, numbered according
to the numbers used in tree$edge. Happy to hear if this works. See below:
# conducts BM simulation by ascending through the tree from the root
node to the tips
# by Liam J. Revell
sim.bm<-function(tree,bm.rate,root.state){
# generate random normal deviates
dev<-rnorm(n=length(tree$edge.length))
# give them appropriate variances
dev<-dev*sqrt(bm.rate*tree$edge.length)
# sort the tree edges, edge.lengths, and devs
ranks<-order(tree$edge[,1])
tree$edge<-tree$edge[ranks,]
tree$edge.length<-tree$edge.length[ranks]
dev<-dev[ranks]
# add them up
state<-matrix(0,length(tree$edge.length),1);
for(j in 1:length(tree$edge.length)){
if(tree$edge[j,1]==length(tree$tip)+1)
state[j]<-root.state+dev[j]
else
state[j]<-state[match(tree$edge[j,1],tree$edge[,2])]+dev[j]
}
# bookkeeping
rownames(state)<-tree$edge[,2]
state<-state[order(as.numeric(rownames(state))),]
sim.bm<-state
}
Sincerely, Liam
Liam J. Revell
NESCent, Duke University
web: http://anolis.oeb.harvard.edu/~liam/
NEW email: lrev...@nescent.org
Joe Felsenstein wrote:
Luke wrote:
The simulations in Geiger do draw from a mvn distribution but - in an
unnecessary step - draw random numbers for every branch and use matrix
multiplication to get the tip values. This is kind of dumb and needlessly
slow but does give you the right answer. Carl is correct that simply drawing
from a mvn distribution is faster and equivalent. If I had time I'd rewrite
the function - maybe I will!
Without knowing exactly what your R code is doing, I just want to add to
these (correct) ideas that there is one simple way of getting tip and
interior node values without much more effort. Work up the tree,
drawing a set of normals (one for each character) to get the net changes in
that branch. Keep going until you reach the tips. (If the characters
are not independent, you also need to have on hand a matrix which
contains the square root of the covariances among characters, and multiply
the values at each node by that.)
This way you don't need a nodes x nodes covariance matrix.
Joe
----
Joe Felsenstein j...@gs.washington.edu
Department of Genome Sciences and Department of Biology,
University of Washington, Box 355065, Seattle, WA 98195-5065 USA
_______________________________________________
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
