Power plot

rm(list=ls())
load("~/Dropbox/FCU/Teaching/Mentoring/2016Spring-2021Fall_AssociateProfessor/2020Fall/MingHanHsu/rcodeMHH/simsBM.rda")

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

bm.power.seq<-power.seq


load("~/Dropbox/FCU/Teaching/Mentoring/2016Spring-2021Fall_AssociateProfessor/2020Fall/MingHanHsu/rcodeMHH/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

ou.power.seq<-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)



load("~/Dropbox/FCU/Teaching/Mentoring/2016Spring-2021Fall_AssociateProfessor/2020Fall/MingHanHsu/rcodeMHH/simsEB.rda")
TypeIIerr<-array(0,c(length(num.spe.array),length(sigma2.array)))
rownames(TypeIIerr)<-num.spe.array
colnames(TypeIIerr)<-sigma2.array


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

eb.power.seq<-power.seq



# 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)




load("~/Dropbox/FCU/Teaching/Mentoring/2016Spring-2021Fall_AssociateProfessor/2020Fall/MingHanHsu/rcodeMHH/simsPMM.rda")
TypeIIerr<-array(0,c(length(num.spe.array),length(sigma2.array)))
rownames(TypeIIerr)<-num.spe.array
colnames(TypeIIerr)<-sigma2.array

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
pmm.power.seq<-power.seq
 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)

save(sigma2.array,
  bm.power.seq,
  ou.power.seq,
  pmm.power.seq,
  eb.power.seq
  ,file="~/Dropbox/FCU/Teaching/Mentoring/2016Spring-2021Fall_AssociateProfessor/2020Fall/MingHanHsu/rcodeMHH/powerplot.RData")

load("~/Dropbox/FCU/Teaching/Mentoring/2016Spring-2021Fall_AssociateProfessor/2020Fall/MingHanHsu/rcodeMHH/powerplot.RData")


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


plot(ou.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),ou.power.seq[2,],type="b",col="blue",pch=3)
points(1:length(sigma2.array),ou.power.seq[3,],type="b",col="red",pch=4)
points(1:length(sigma2.array),ou.power.seq[4,],type="b",col="purple",pch=5)

plot(eb.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),eb.power.seq[2,],type="b",col="blue",pch=3)
points(1:length(sigma2.array),eb.power.seq[3,],type="b",col="red",pch=4)
points(1:length(sigma2.array),eb.power.seq[4,],type="b",col="purple",pch=5)

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

# Load required libraries
library(ggplot2)
library(dplyr)

## 
## Attaching package: 'dplyr'

## The following objects are masked from 'package:stats':
## 
##     filter, lag

## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union

library(tidyr)
library(gridExtra)

## 
## Attaching package: 'gridExtra'

## The following object is masked from 'package:dplyr':
## 
##     combine

# Convert power matrices to data frames and add necessary columns
bm.df <- as.data.frame(t(bm.power.seq)) %>% mutate(sigma = sigma2.array, Model = "BM")
ou.df <- as.data.frame(t(ou.power.seq)) %>% mutate(sigma = sigma2.array, Model = "OU")
eb.df <- as.data.frame(t(eb.power.seq)) %>% mutate(sigma = sigma2.array, Model = "EB")
pmm.df <- as.data.frame(t(pmm.power.seq)) %>% mutate(sigma = sigma2.array, Model = "PMM")

# Bind the data frames together
df <- bind_rows(bm.df, ou.df, eb.df, pmm.df)

# Convert data to long format
df.long <- df %>% pivot_longer(cols = -c(sigma, Model), names_to = "PowerSeq", values_to = "Power")



# convert PowerSeq to a factor with specified levels
df.long$PowerSeq <- factor(df.long$PowerSeq, levels = c("16", "32", "64", "128"))

# Draw a single faceted plot with ordered legend
combined_plot <- df.long %>%
  ggplot(aes(x = sigma, y = Power, color = PowerSeq, shape = PowerSeq)) +
  geom_point() +
  geom_line() +
  facet_wrap(~ Model, scales = "free", ncol = 2) +
  scale_color_manual(values = c("black", "blue", "red", "purple"), labels = c("16", "32", "64", "128")) +
  scale_shape_manual(values = c(19, 3, 4, 5), labels = c("16", "32", "64", "128")) +
  labs(x = expression(sigma[2]/sigma[1]), y = "Power", color = "Sample size", shape = "Sample size") +
  ylim(0, 1) +
  theme_minimal() +
  theme(legend.position = "bottom") # move the legend to the bottom
print(combined_plot)