Power analysis for Ornstein Uhlenbeck model

source("https://tonyjhwueng.info/phyrates/ComparingOU.r")

library(ape)
library(geiger)
#install.packages("mvMORPH")
library(mvMORPH)

# Can look at package mvMORPH
#https://cran.r-project.org/web/packages/mvMORPH/vignettes/tutorial_mvMORPH.pdf

# we use
#install.packages("diversitree")
#https://www.zoology.ubc.ca/~fitzjohn/diversitree.docs/sim.character.html
library(diversitree)

num.spe.array<-c(16,32,64,128)
sigma<-1
sigma1.array<-rep(sigma,5)
sigma2.array<-sigma*c(1.25,1.5,2.0,3.0,4.0)
simstrait<-1000

for(taxaIndex in 1: length(num.spe.array)){
  #taxaIndex<-1
  taxa<-num.spe.array[taxaIndex]
  print(taxa)
  #taxa<-compute.brlen(stree(taxa,"balanced"))#balanced tree
  tree<-rcoal(taxa)
  #plot(tree)
  for(s2Index in 1:length(sigma2.array)){
   # s2Index<-1
    sigma2<-sigma2.array[s2Index]
    
  for(traitIndex in 1:simstrait){
    #print(traitIndex)
    #traitIndex<-1
#?sim.char 
  #x1<-sim.char(phy=tree,par=sigma, model="BM")
  #x2<-sim.char(phy=tree,par=sigma2, model="BM")
    alpha<-0.2
    theta<-0
   x1<- sim.character(tree, pars=c(sigma,alpha,theta), x0=0, model="ou", br=NULL)
   x2<- sim.character(tree, pars=c(sigma2,alpha,theta), x0=0, model="ou", br=NULL)
  x<-cbind(x1,x2)
  rownames(x)<-tree$tip.label
  try(assign(paste("taxa",taxa,"sgmratio",s2Index,"sim",traitIndex,sep="") ,CompareRates.multTrait(phy=tree,x=x,TraitCov=T,ms.err=NULL,ms.cov=NULL)))
  }
 }
}


save.image("simsOU.rda")

#nohup R CMD BATCH simsBM.R >/dev/null &

load(url("https://tonyjhwueng.info/phyrates/simsOUa0.5.rda"))

TypeIIerr<-array(0,c(length(num.spe.array),length(sigma2.array)))
rownames(TypeIIerr)<-num.spe.array
colnames(TypeIIerr)<-sigma2.array
TypeIIerr

for(taxaIndex in 1: length(num.spe.array)){
  #taxaIndex<-1
  taxa<-num.spe.array[taxaIndex]
#  print(taxa)
  #compute.brnlen(rtree(taxa,))
#  tree<-rcoal(taxa)
  #plot(tree)
  for(s2Index in 1:length(sigma2.array)){
    # s2Index<-1
    sigma2<-sigma2.array[s2Index]
    for(traitIndex in 1:simstrait){
      try(if(get(paste("taxa",taxa,"sgmratio",s2Index,"sim",traitIndex,sep=""))$Prob>0.05){
        TypeIIerr[taxaIndex,s2Index]<-TypeIIerr[taxaIndex,s2Index]+1
      })
    }
  }
}   
TypeIIerr/simstrait

power.seq<-1-TypeIIerr/simstrait
power.seq
plot(power.seq[1,],type="b",ylim=c(0,1),xaxt="n",xlab="ratio",ylab="Power",main="OU")
axis(1, at=1:length(sigma2.array),labels=sigma2.array, las=2)
points(1:length(sigma2.array),power.seq[2,],type="b",col="blue",pch=3)
points(1:length(sigma2.array),power.seq[3,],type="b",col="red",pch=4)
points(1:length(sigma2.array),power.seq[4,],type="b",col="purple",pch=5)