## Wednesday, February 16, 2011

### BM simulation with a trend

Dave Bapst just posted a query to the R-sig-phylo about simulating Brownian evolution with a trend. I take this to mean that evolutionary changes along branches in the tree are drawn from a random normal distribution with variances sig2*tree\$edge.length and means mu*tree\$edge.length. [Initially, in my first response to the query, I did not multiply mu by tree\$edge.length; but that doesn't make any sense.] I have now added this capability to v0.2 of my function fastBM() which can be downloaded from my R phylogenetics page here.

> x<-fastBM(tree,a=0,mu=0.5,sig2=2,internal=FALSE)

at the command prompt. This will return a n x 1 vector of character values (in x for the n species in tree.

To test out this function, let's try simulating data and then fitting a "trend" model using {geiger}'s fitContinuous() function. [Obviously, since these are stochastic simulations, individual results will vary.]

> tree<-rtree(500) # first simulate stochastic tree
> x<-fastBM(tree,mu=0.5,sig2=2) # now data on the tree
> fitContinuous(phy=tree,data=x,model="trend") # now fit model
Fitting trend model:
\$Trait1
\$Trait1\$lnl
[1] -900.9902
\$Trait1\$mean
[1] 1.087359
\$Trait1\$beta
[1] 1.87744
\$Trait1\$mu
[1] 0.4968745
\$Trait1\$aic
[1] 1807.980
\$Trait1\$aicc
[1] 1808.029
\$Trait1\$k
[1] 3

Seems to work!

Note that Gene Hunt also pointed out (correctly, of course) that this can be done using the {MASS} function mvrnorm(). To see the discussion thread on this topic, check out the archived R-sig-phylo emails here.