rm(list=ls())
setwd("~/Dropbox/FCU/Teaching/Mentoring/2022Spring-Current_Professor/2022Fall/YourueiMin/rcode_dcj/simulation")
# Load required libraries
#install.packages("ggplot2")
library(ggplot2)
library(reshape2) # for melting the data
library(xtable)
#install.packages("gridExtra")
library(gridExtra)
# wirte a function to plot the power curve
T1PowerCurve<-function(power_matrix=power_matrix,T1err.all=T1err.all,alt.phi=alt.phi,taxa.all=taxa.all,treetype=treetype){
power <- data.frame(
alt.phi = rep(alt.phi, length(taxa.all)),
power = c(power_matrix),
taxa = rep(taxa.all, each = length(alt.phi))
)
power
type_I_error <- T1err.all
type_I_error
# Plot
p1<-ggplot(power, aes(x = alt.phi, y = power, color = as.factor(taxa), group = taxa)) + geom_line(size = 1) + geom_hline(yintercept = type_I_error, linetype = "dashed", color = "red", size = 0.5) + labs(title = "Power Across Alt.Phi for Different Taxa", x = "Alt.Phi", y = "Power", color = "Taxa") + ggtitle(paste("T1 Power Curve for ",treetype," Tree",sep="")) + scale_color_discrete(name = "Taxa") + theme_minimal()
print(p1)
return(p1)
}
treetype.array<-c("coal","purebirth","balanced","birthdeath")
taxa.all<-c(16,32,64,128)
power_matrix_full<-NULL
T1err.df_full<-NULL
for (treetype in treetype.array){
filetoload<-paste("onerateT1power_",treetype,"Tree.RData",sep="")
print(filetoload)
load(filetoload)
print(T1err.all)
print(power_matrix)
power_matrix_full<-rbind(power_matrix_full,power_matrix)
# p1<-T1PowerCurve(power_matrix=power_matrix,T1err.all=T1err.all, alt.phi=alt_phi_values,taxa.all=taxa.all,treetype=treetype)
# print(p1)
#---given power_matrix of rows names taxa.all and colname phi.alt
# plot the power curve using alt.phi as x axis and power as y axis, plot the line each taxa across alt.phi, each line represent a taxa for each tree type
# add T1err.all as horizontal line,
# make latex table for each tree type
# Reshape data for ggplot
power_melted <- melt(power_matrix, id.vars = "taxa.all", variable.name = "alt.phi", value.name = "power")
power_melted
colnames(power_melted)<-c("taxa.all","alt.phi","power")
T1err.df <- data.frame(taxa.all = names(T1err.all), T1err = unname(T1err.all))
# Convert taxa.all to a factor
power_melted$taxa.all <- as.factor(power_melted$taxa.all)
T1err.df$taxa.all <- as.factor(T1err.df$taxa.all)
T1err.df_full<-rbind(T1err.df_full,T1err.df)
# Define colors (you can change these to any colors you want)
my_colors <- rainbow(length(unique(power_melted$taxa.all)))
# Plot
assign(paste("p.",treetype,sep = ""),
ggplot(power_melted, aes(x = as.numeric(as.character(alt.phi)), y = power, group = taxa.all)) + geom_line(aes(color = taxa.all), size = 1) + geom_hline(data = T1err.df, aes(yintercept = T1err, color = taxa.all), linetype = "dashed", size = 0.5) + scale_color_manual(values = my_colors) + labs(title = paste("Power curve: ", treetype," tree",sep="") , x = expression(phi[a]), y = "Power value", color = "Taxa") + theme_minimal()
)
}
## [1] "onerateT1power_coalTree.RData"
## 16 32 64 128
## 0.066 0.041 0.043 0.055
## 0.05 0.2 0.35 0.5 0.65 0.8 0.95
## 16 0.050 0.211 0.481 0.746 0.924 0.97 0.991
## 32 0.084 0.355 0.792 0.965 0.997 1.00 1.000
## 64 0.081 0.616 0.970 1.000 1.000 1.00 1.000
## 128 0.124 0.903 1.000 1.000 1.000 1.00 1.000
## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
## [1] "onerateT1power_purebirthTree.RData"
## 16 32 64 128
## 0.056 0.057 0.058 0.037
## 0.05 0.2 0.35 0.5 0.65 0.8 0.95
## 16 0.057 0.184 0.443 0.750 0.907 0.971 0.997
## 32 0.071 0.343 0.807 0.971 0.998 1.000 1.000
## 64 0.086 0.642 0.965 0.999 1.000 1.000 1.000
## 128 0.125 0.884 1.000 1.000 1.000 1.000 1.000
## [1] "onerateT1power_balancedTree.RData"
## 16 32 64 128
## 0.058 0.061 0.057 0.053
## 0.05 0.2 0.35 0.5 0.65 0.8 0.95
## 16 0.064 0.197 0.490 0.748 0.924 0.969 0.996
## 32 0.074 0.313 0.766 0.968 0.997 1.000 1.000
## 64 0.085 0.604 0.971 1.000 1.000 1.000 1.000
## 128 0.118 0.861 1.000 1.000 1.000 1.000 1.000
## [1] "onerateT1power_birthdeathTree.RData"
## 16 32 64 128
## 0.054 0.043 0.053 0.066
## 0.05 0.2 0.35 0.5 0.65 0.8 0.95
## 16 0.075 0.210 0.465 0.753 0.910 0.98 0.99
## 32 0.074 0.346 0.763 0.963 0.996 1.00 1.00
## 64 0.085 0.611 0.968 1.000 1.000 1.00 1.00
## 128 0.164 0.893 1.000 1.000 1.000 1.00 1.00
grid.arrange(p.coal,p.purebirth,p.balanced,p.birthdeath,ncol=2)
fulltable<-cbind(T1err.df_full,power_matrix_full)
rownames(fulltable)<-NULL
fulltable
## taxa.all T1err 0.05 0.2 0.35 0.5 0.65 0.8 0.95
## 1 16 0.066 0.050 0.211 0.481 0.746 0.924 0.970 0.991
## 2 32 0.041 0.084 0.355 0.792 0.965 0.997 1.000 1.000
## 3 64 0.043 0.081 0.616 0.970 1.000 1.000 1.000 1.000
## 4 128 0.055 0.124 0.903 1.000 1.000 1.000 1.000 1.000
## 5 16 0.056 0.057 0.184 0.443 0.750 0.907 0.971 0.997
## 6 32 0.057 0.071 0.343 0.807 0.971 0.998 1.000 1.000
## 7 64 0.058 0.086 0.642 0.965 0.999 1.000 1.000 1.000
## 8 128 0.037 0.125 0.884 1.000 1.000 1.000 1.000 1.000
## 9 16 0.058 0.064 0.197 0.490 0.748 0.924 0.969 0.996
## 10 32 0.061 0.074 0.313 0.766 0.968 0.997 1.000 1.000
## 11 64 0.057 0.085 0.604 0.971 1.000 1.000 1.000 1.000
## 12 128 0.053 0.118 0.861 1.000 1.000 1.000 1.000 1.000
## 13 16 0.054 0.075 0.210 0.465 0.753 0.910 0.980 0.990
## 14 32 0.043 0.074 0.346 0.763 0.963 0.996 1.000 1.000
## 15 64 0.053 0.085 0.611 0.968 1.000 1.000 1.000 1.000
## 16 128 0.066 0.164 0.893 1.000 1.000 1.000 1.000 1.000
# make latex table for each tree type
xtable(fulltable,digits=3)
## % latex table generated in R 4.3.2 by xtable 1.8-4 package
## % Sun Jan 7 10:26:55 2024
## \begin{table}[ht]
## \centering
## \begin{tabular}{rlrrrrrrrr}
## \hline
## & taxa.all & T1err & 0.05 & 0.2 & 0.35 & 0.5 & 0.65 & 0.8 & 0.95 \\
## \hline
## 1 & 16 & 0.066 & 0.050 & 0.211 & 0.481 & 0.746 & 0.924 & 0.970 & 0.991 \\
## 2 & 32 & 0.041 & 0.084 & 0.355 & 0.792 & 0.965 & 0.997 & 1.000 & 1.000 \\
## 3 & 64 & 0.043 & 0.081 & 0.616 & 0.970 & 1.000 & 1.000 & 1.000 & 1.000 \\
## 4 & 128 & 0.055 & 0.124 & 0.903 & 1.000 & 1.000 & 1.000 & 1.000 & 1.000 \\
## 5 & 16 & 0.056 & 0.057 & 0.184 & 0.443 & 0.750 & 0.907 & 0.971 & 0.997 \\
## 6 & 32 & 0.057 & 0.071 & 0.343 & 0.807 & 0.971 & 0.998 & 1.000 & 1.000 \\
## 7 & 64 & 0.058 & 0.086 & 0.642 & 0.965 & 0.999 & 1.000 & 1.000 & 1.000 \\
## 8 & 128 & 0.037 & 0.125 & 0.884 & 1.000 & 1.000 & 1.000 & 1.000 & 1.000 \\
## 9 & 16 & 0.058 & 0.064 & 0.197 & 0.490 & 0.748 & 0.924 & 0.969 & 0.996 \\
## 10 & 32 & 0.061 & 0.074 & 0.313 & 0.766 & 0.968 & 0.997 & 1.000 & 1.000 \\
## 11 & 64 & 0.057 & 0.085 & 0.604 & 0.971 & 1.000 & 1.000 & 1.000 & 1.000 \\
## 12 & 128 & 0.053 & 0.118 & 0.861 & 1.000 & 1.000 & 1.000 & 1.000 & 1.000 \\
## 13 & 16 & 0.054 & 0.075 & 0.210 & 0.465 & 0.753 & 0.910 & 0.980 & 0.990 \\
## 14 & 32 & 0.043 & 0.074 & 0.346 & 0.763 & 0.963 & 0.996 & 1.000 & 1.000 \\
## 15 & 64 & 0.053 & 0.085 & 0.611 & 0.968 & 1.000 & 1.000 & 1.000 & 1.000 \\
## 16 & 128 & 0.066 & 0.164 & 0.893 & 1.000 & 1.000 & 1.000 & 1.000 & 1.000 \\
## \hline
## \end{tabular}
## \end{table}