Hi all,
In rTraitCont, the argument model can be a function, eg:
f <- function(x, l) x + rnorm(1, 1, sqrt(l * 0.1))
which is a way to model a trend. Here's an example of what you can do:
tr <- rbdtree(.1, .05)
tr <- makeNodeLabel(tr)
x <- rTraitCont(tr, model = f, ancestor = TRUE)
bt <- branching.times(tr)
plot(-bt, x[names(bt)])
You may add the tip values with:
n <- Ntip(tr)
points(rep(0, n), x[1:n], pch = 2)
You can also draw lines linking the ancestors with their descendants in
this "morpho-time-space" since they can be identified in the edge matrix.
If you don't like this model, just change f. As a side note:
f <- function(x, l) x + rnorm(1, 0, sqrt(l * 0.1))
x <- rTraitCont(tr, model = f)
simulates the same model (ie, BM) than:
x <- rTraitCont(tr)
HTH,
Emmanuel
Hunt, Gene wrote on 17/02/2011 01:25:
Hi,
To Liam's extremely helpful and practical reply, I'd add that one can also use
theory to generate tip data for a trend. For a given tree, all tip values can
be considered a single draw from a multivariate normal distribution, with each
dimension representing a tip taxon. For this MVN, the vector of means is equal
to the directionality term (mu) times the root-to-tip distance for each
terminal taxon, and the variance-covariance is that for BM. So the following
two lines also generate tip data for a trend model for phylogeny (phy):
V<- vcv(phy)
x<- mvrnorm(n=1, diag(V)*mu, V)
Setting n to a higher number will generate more replicates. Liam in his blog
talks about why the above code can be slow, but I think it is useful to point
out theory. BM can be generated the same way with the multivariate mean equal
to the ancestral value repeated as many times as there are tips in the tree.
Best,
Gene
On 2/16/11 1:04 PM, "Liam J. Revell" <liam.rev...@umb.edu> wrote:
Hi Dave.
I believe BM with a trend should just have a non-zero mean for the
random normal changes along individual branches of the tree. It was
easy to add this to a function that I already have for fast Brownian
motion simulation. I will post this update to my website:
http://anolis.oeb.harvard.edu/~liam/R-phylogenetics/ but it is so simple
that I have also pasted the code below.
- Liam
# Simulates BM evolution more quickly.
# A trend can be simulated by mu!=0.
# mu=0 is standard BM; mu<0 downward trend; mu>0 upward trend.
# Written by Liam J. Revell 2011
fastBM.trend<-function(tree,a=0,mu=0,sig2=1,internal=FALSE){
n<-length(tree$tip)
# first simulate changes along each branch
x<-rnorm(n=length(tree$edge.length),mean=mu,sd=sqrt(sig2*tree$edge.length))
# now add them up
y<-matrix(0,nrow(tree$edge),ncol(tree$edge))
for(i in 1:length(x)){
if(tree$edge[i,1]==(n+1))
y[i,1]<-a
else
y[i,1]<-y[match(tree$edge[i,1],tree$edge[,2]),2]
y[i,2]<-y[i,1]+x[i]
}
rm(x); x<-c(y[1,1],y[,2])
names(x)<-c(n+1,tree$edge[,2])
x<-x[as.character(1:(n+tree$Nnode))]
names(x)[1:n]<-tree$tip.label
# return simulated data
if(internal==TRUE)
return(x) # include internal nodes
else
return(x[1:length(tree$tip.label)]) # tip nodes only
}
--
Liam J. Revell
University of Massachusetts Boston
web: http://www.bio.umb.edu/facstaff/faculty_Revell.html
(new) email: liam.rev...@umb.edu
(new) blog: http://phytools.blogspot.com
On 2/16/2011 12:28 PM, David Bapst wrote:
All-
For a class presentation I was working on, I wanted other students to use
model fitting with fitContinuous() on simulated traits, under different
models of trait evolution. I was surprised to discover there seemed to be no
easy way to simulate active directional trend mechanisms, unlike other trait
evolution models.
One idea I had was to kludge it.
tree<-rtree(20)
trait1<-,rTraitCont(tree,sigma=1)+2*diag(vcv.phylo(tree))
An alternative would be an OU with the root far from the theta value. Takes
some playing with the parameters to get something that looks like a trend in
a traitgram.
trait2<-rTraitCont(tree,model="OU",sigma=1,alpha=0.5,root.value=-10)
Both of these were strongly preferred as 'trend' in fitContinuous() compared
to BM and OU1, but I was wondering if there were any alternatives I was
missing.
-Dave
_______________________________________________
R-sig-phylo mailing list
R-sig-phylo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-phylo
--
Gene Hunt
Curator, Department of Paleobiology
National Museum of Natural History
Smithsonian Institution [NHB, MRC 121]
P.O. Box 37012
Washington DC 20013-7012
Phone: 202-633-1331 Fax: 202-786-2832
http://paleobiology.si.edu/staff/individuals/hunt.cfm
_______________________________________________
R-sig-phylo mailing list
R-sig-phylo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-phylo
--
Emmanuel Paradis
IRD, Jakarta, Indonesia
http://ape.mpl.ird.fr/
_______________________________________________
R-sig-phylo mailing list
R-sig-phylo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-phylo
