For reasons I won't get into in detail, I just wrote a simple function to conduct Ornstein-Uhlenbeck simulations by finely discretizing time and then using a simple difference equation version of the OU model to simulate as follows:

```
multiOU<-function(tree,alpha,sig2,theta=NULL,a0=NULL,internal=TRUE,...){
if(hasArg(dt)) dt<-list(...)$dt
else dt<-1/1000*max(nodeHeights(tree))
ss<-sort(unique(c(getStates(tree,"tips"),getStates(tree,"nodes"))))
if(is.null(theta)) theta<-setNames(rep(0,length(ss)),ss)
if(is.null(a0)) a0<-0
tree<-reorder(tree,"cladewise")
root<-Ntip(tree)+1
S<-matrix(NA,nrow(tree$edge),2)
S[which(tree$edge[,1]==root),1]<-a0
for(i in 1:nrow(tree$edge)){
x1<-S[i,1]
for(j in 1:length(tree$maps[[i]])){
t<-tree$maps[[i]][j]
ALPHA<-alpha[names(t)]
SIG2<-sig2[names(t)]
THETA<-theta[names(t)]
t<-c(rep(dt,floor(t/dt)),t%%dt)
for(k in 1:length(t))
x1<-x1+ALPHA*(THETA-x1)*t[k]+
rnorm(n=1,sd=sqrt(SIG2*t[k]))
}
S[i,2]<-x1
if(any(tree$edge[,1]==tree$edge[i,2]))
S[which(tree$edge[,1]==tree$edge[i,2]),1]<-S[i,2]
}
x<-setNames(c(S[1,1],S[,2]),c(tree$edge[1,1],
tree$edge[,2]))[as.character(1:(Ntip(tree)+tree$Nnode))]
names(x)[1:Ntip(tree)]<-tree$tip.label
if(!internal) x<-x[tree$tip.label]
x
}
```

In the following, I'll use it to simulate:

```
library(phytools)
set.seed(116)
cols<-setNames(c("black","blue","red"),letters[1:3])
tree<-pbtree(n=100,scale=2)
Q<-matrix(c(-1,1,0,1,-2,1,0,1,-1),3,3,
dimnames=list(letters[1:3],letters[1:3]))
tree<-sim.history(tree,Q,anc="a")
```

```
## Done simulation(s).
```

```
sig2<-setNames(c(1,2,1),letters[1:3])
alpha<-setNames(c(2,2,2),letters[1:3])
theta<-setNames(c(10,0,-10),letters[1:3])
x<-multiOU(tree,alpha,sig2,theta)
phenogram(tree,x,ftype="off",colors=cols,ylim=c(-11,11),
xlim=c(0,2))
abline(h=c(-10,0,10),lty="dashed",col=sapply(c("red","blue",
"black"),make.transparent,alpha=0.5),lwd=3,lend=2)
add.simmap.legend(colors=cols,prompt=FALSE,x=0,y=-5)
```

I'm going to add (importantly) that the same simulation is done without
the discrete difference equation approximation by `OUwie.sim`

in the powerful OUwie
package by Jeremy Beaulieu & Brian O'Meara.

Finally, let's fit the appropriate 'OUwie' (that is, flexible OU) model to these data:

```
library(OUwie)
data<-data.frame(Genus_species=tree$tip.label,
Reg=getStates(tree,"tips"),X=as.numeric(x[tree$tip.label]))
OUwie(tree,data,model="OUMV",simmap.tree=TRUE)
```

```
## Initializing...
## Finished. Begin thorough search...
## Finished. Summarizing results.
```

```
##
## Fit
## lnL AIC AICc model ntax
## -153.169 320.3379 321.5553 OUMV 100
##
##
## Rates
## a b c
## alpha 1.6410055 1.641006 1.641006
## sigma.sq 0.8331012 12.695801 73.799702
##
## Optima
## a b c
## estimate 9.5391941 -2.5029839 -14.22271
## se 0.1269122 0.8906037 2.51031
##
## Arrived at a reliable solution
```

One of our σ^{2} estimates is way off - but aside from that,
we seem to have more or less captured the generating process. Cool.

## No comments:

## Post a Comment