Power analysis for Early Burst model

#source("~/Dropbox/MingHanHsu/rcodeMHH/ComparingEB.r")
source("https://tonyjhwueng.info/phyrates/ComparingEB.r")

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

# Can look at package mvMORPH
#https://rdrr.io/cran/mvMORPH/man/mvSIM.html
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")
    
    
  x1<- mvSIM(tree, nsim = 1, error = NULL, model = c("EB"),param = list(sigma = sigma, theta=0,beta = 0.1))
  x2<- mvSIM(tree, nsim = 1, error = NULL, model = c("EB"),param = list(sigma = sigma2,theta=0, beta = 0.1))
  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("simsEB.rda")

load(url("https://tonyjhwueng.info/phyrates/simsEB.rda"))
TypeIIerr<-array(0,c(length(num.spe.array),length(sigma2.array)))
rownames(TypeIIerr)<-num.spe.array
colnames(TypeIIerr)<-sigma2.array
TypeIIerr

##     1.25 1.5 2 3 4
## 16     0   0 0 0 0
## 32     0   0 0 0 0
## 64     0   0 0 0 0
## 128    0   0 0 0 0

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

##      1.25   1.5     2     3     4
## 16  0.727 0.706 0.611 0.359 0.242
## 32  0.674 0.573 0.395 0.183 0.086
## 64  0.618 0.465 0.223 0.095 0.057
## 128 0.552 0.294 0.089 0.057 0.032

power.seq<-1-TypeIIerr/simstrait
power.seq

##      1.25   1.5     2     3     4
## 16  0.273 0.294 0.389 0.641 0.758
## 32  0.326 0.427 0.605 0.817 0.914
## 64  0.382 0.535 0.777 0.905 0.943
## 128 0.448 0.706 0.911 0.943 0.968

plot(power.seq[1,],type="b",ylim=c(0,1),xaxt="n",xlab="ratio",ylab="Power",main="EB")
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)