[R-sig-phylo] erroneous results with rTraitCont simulations of the OU process

Emmanuel Paradis Emmanuel.Paradis at ird.fr
Fri Mar 4 03:31:24 CET 2011


Hi Carl,

?rTraitCont says:

           By default the following formula is used:

           x(t'') = x(t') - alpha l (x(t') - theta) + sigma    l epsilon

           where l (= t'' - t') is the branch length, and epsilon ~ N(0,
           1). If alpha > 1, this may lead to chaotic oscillations. Thus
           an alternative parameterisation is used if ‘linear = FALSE’:

           x(t'') = x(t') - (1 - exp(-alpha l)) * (x(t') -
               theta) + sigma l epsilon

I recently found an alternative formulation in the physical literature 
that seems similar to the second one above, but I need to check that and 
shall possibly program it.

Cheers,

Emmanuel

Carl Boettiger wrote on 04/03/2011 05:15:
> Dear list,
> 
> The simulation of OU processes using rTraitCont seems to return
> erroneous results for non-small values of alpha:
> 
> Consider this example:
> 
> require(geiger)
> require(TreeSim)
> tree <- sim.bd.taxa(100,1,1,0)[[1]]
> 
> x <- rTraitCont(tree, model="OU", sigma=1, alpha=10, theta=0, root.value=0)
> mean(x)
> [1] -24.13565
> 
> 
> The mean should near theta=0.  Smaller values of alpha don't have this
> problem.
>  x <- rTraitCont(tree, model="OU", sigma=1, alpha=3, theta=0, root.value=0)
>> mean(x)
> [1] 0.01494487
> 
> 
> I don't seem to receive any warnings about this.  Moreover, it is not
> clear why this should be an issue.  Performing simulations of the same
> model on the same tree in ouch does not give this error.  Perhaps I've
> missed something in my  call to rTraitCont?
> 
> 
> Thanks for your time and advice!
> 
> Carl
> 

-- 
Emmanuel Paradis
IRD, Jakarta, Indonesia
http://ape.mpl.ird.fr/



More information about the R-sig-phylo mailing list