Power analysis for Phylognetic Mixed Model model

source("https://tonyjhwueng.info/phyrates/ComparingPMM.r")
library(ape)
library(geiger)
#install.packages("mvMORPH")
library(mvMORPH)
#install.packages("phylolm")
library(phylolm)
# Can look at package mvMORPH
#https://cran.r-project.org/web/packages/mvMORPH/vignettes/tutorial_mvMORPH.pdf

# we use
#install.packages("phylocurve")
#library(phylocurve)

#https://rdrr.io/cran/phylocurve/man/sim.traits.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")
    #https://rdrr.io/cran/phylolm/man/rTrait.html
  #  tree<-rcoal(4)
   x1<- rTrait(n=1, phy=tree, model="lambda",parameters=list(lambda = 0.5))
   x2<- rTrait(n=1, phy=tree, model="lambda",parameters=list(lambda = 0.5))
  #x1<-  sim.traits(ntaxa =taxa, ntraits = 1, nreps = 1,nmissing = 0, tree,   model = "lambda", parameters=list(lambda = 0.5), nsim = 1, return.type = "matrix")
  #x2<-  sim.traits(ntaxa =taxa, ntraits = 1, nreps = 1,nmissing = 0, tree,   model = "lambda", parameters=list(lambda = 0.5), nsim = 1, return.type = "matrix")
  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("simsPMMh0.5.rda")

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

load(url("https://tonyjhwueng.info/phyrates/simsPMM.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.264 0.279 0.266 0.255 0.298
## 32  0.214 0.210 0.212 0.211 0.224
## 64  0.081 0.097 0.096 0.102 0.104
## 128 0.070 0.068 0.085 0.075 0.078

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

##      1.25   1.5     2     3     4
## 16  0.736 0.721 0.734 0.745 0.702
## 32  0.786 0.790 0.788 0.789 0.776
## 64  0.919 0.903 0.904 0.898 0.896
## 128 0.930 0.932 0.915 0.925 0.922

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