[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