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

Mon Mar 7 04:47:13 CET 2011

Dear Carl,

Thanks for the reply too. I have designed rTraitCont (and rTraitDisc) as
an alternative to the methods based on the expected var-covar.
It allows you to specify separately parameters (selective regimes, rates
of evolution, and /or optima) on each branch of the tree at your will,
and you can output the ancestral values of each run.

Carl Boettiger wrote on 05/03/2011 02:15:
> Dear Emmanuel, List,
>
> Thank you for the reply, I hadn't realized that rTraitCont simulates
> the trait evolution numerically from the SDE.  I believe it would be
> more efficient and accurate to use the analytic solution,

Generally speaking, I don't see why this should be so. If you simulate a
BM in continuous time on t units, or if you draw random normal variates
with variance proportional to t, it's exactly the same.

> which is
> what sim.char and simulate in geiger and ouch do.   Since both BM and
> OU are linear sdes, the solutions are normals.  Putting them on a tree
> makes the distribution of tips multivariate normal, as in Hansen's
> 1997 paper.  I believe it should be much faster as well.

Maybe for small trees but this will scale as O(n^2), whereas
rTraitCont's algorithm will scale as O(n) because it is proportional to
the number of branches.

> Luke Harmon points out this can be done in sim.char by rescaling the tree first
>
> x <- sim.char(ouTree(tree, alpha=50), model.matrix-matrix(sigma) )
> mean(x)
>
> (though as he says simply rescaling this doesn't let you simulate from
> an arbitrary root state and trait optimum theta).
>
> In general you can simulate the OU model (with arbitrary root state
> X_0 and theta) by drawing random numbers from the multivariate
> distribution satisfying
> E(x_t)=x_0 e^{-\theta t}+\mu(1-e^{-\theta t})
> Cov(x_s,x_t) = \frac{\sigma^2}{2\theta}\left( e^{-\theta(t-s)} -
> e^{-\theta(t+s)} \right).
> for s < t (so that min(s, t) = s),
>
> The rTraitCont approach could be used to simulate/numerically
> intergrate SDEs even when they can't be solved analytically, but is
> unnecessary here -- and I believe can lead to errors even for alpha
> less then one?

Yes, if you don't have the correct formulae. In the case of BM, this is
quite straightforward (see my comment above). In the case of OU, I'm not
so sure of it, that's why I was looking at the literature recently.

> A good discussion of these issues is in
> ﻿Iacus, S.M. Simulation and Inference for Stochastic Differential
> Equations With R Examples. (Springer: New York, 2008).

Thanks for the reference. The one I mentioned in my previous mail is:

Daniel T. Gillespie (1996) Exact numerical simulation of the
Ornstein-Uhlenbeck process and its integral. Physical Review E 54(2):
2084-2091.

Cheers,

Emmanuel

> Thanks for the reply, hope this is helpful!
>
> -Carl
>
>
> On Thu, Mar 3, 2011 at 6:31 PM, Emmanuel Paradis
> <Emmanuel.Paradis at ird.fr> wrote:
>> 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
>>>
>> --
>> IRD, Jakarta, Indonesia
>> http://ape.mpl.ird.fr/
>>
>
>
>

--