rm(list=ls())
load("~/Dropbox/FCU/Teaching/Mentoring/2016Spring-2021Fall_AssociateProfessor/2020Fall/ChiYuWu/rcodeCYW/phynb2sim2_1.rda")
# Define filter ranges for both dimensions
filter.range <- list(c(1, 5), c(3, 7))
# Apply filter to each of the last dimension of tree.coef.array
for (i in 1:2) {
# Apply filter on the 4th dimension
tree.coef.array[,,,i] <- apply(tree.coef.array[,,,i], c(1,2,3), function(x) {
# If the value is in the desired range, keep it, otherwise set it to NA
ifelse(x >= filter.range[[i]][1] & x <= filter.range[[i]][2], x, NA)
})
}
#tree.coef.array[2,1,,]
#
# outtable<-array(NA,c(length(taxa.array),length(treetype.array)*2))
# rownames(outtable)<-taxa.array
# outtable
# dim(tree.coef.array)
# for(treetypeIndex in 1:length(treetype.array)){
# for(taxaIndex in 1:length(taxa.array)){
# outtable[taxaIndex,(2*treetypeIndex-1):(2*treetypeIndex)] <- c(paste(round(mean(tree.coef.array[treetypeIndex,taxaIndex,,1],na.rm=T),3),"(",round(sd(tree.coef.array[treetypeIndex,taxaIndex,,1],na.rm=T),3),")",sep=""),paste(round(mean(tree.coef.array[treetypeIndex,taxaIndex,,2],na.rm=T),3),"(",round(sd(tree.coef.array[treetypeIndex,taxaIndex,,2],na.rm=T),3),")",sep="")
# )
# }
# }
#
# print(paste("r=",r,sep=""))
# print(outtable)
outtable<-array(NA,c(length(taxa.array),length(treetype.array)*2))
rownames(outtable)<-taxa.array
outtable
## [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8]
## 16 NA NA NA NA NA NA NA NA
## 32 NA NA NA NA NA NA NA NA
## 64 NA NA NA NA NA NA NA NA
## 128 NA NA NA NA NA NA NA NA
dim(tree.coef.array)
## [1] 4 4 1000 2
for(treetypeIndex in 1:length(treetype.array)){
for(taxaIndex in 1:length(taxa.array)){
b0_values <- tree.coef.array[treetypeIndex,taxaIndex,,1][is.finite(tree.coef.array[treetypeIndex,taxaIndex,,1])]
b1_values <- tree.coef.array[treetypeIndex,taxaIndex,,2][is.finite(tree.coef.array[treetypeIndex,taxaIndex,,2])]
outtable[taxaIndex,(2*treetypeIndex-1):(2*treetypeIndex)] <- c(paste(round(mean(b0_values, na.rm=T),3),"(",round(sd(b0_values, na.rm=T),3),")",sep=""),
paste(round(mean(b1_values, na.rm=T),3),"(",round(sd(b1_values, na.rm=T),3),")",sep="")
)
}
}
print(paste("r=",r,sep=""))
## [1] "r=10.698"
print(outtable)
## [,1] [,2] [,3] [,4] [,5]
## 16 "2.113(0.864)" "4.693(0.917)" "1.651(0.673)" "5.029(1.343)" "1.878(0.913)"
## 32 "2.14(0.826)" "4.97(0.856)" "2.029(0.823)" "4.613(0.924)" "1.543(0.548)"
## 64 "2.031(0.779)" "4.997(1.008)" "2.312(0.831)" "4.82(0.832)" "1.71(0.686)"
## 128 "2.866(0.802)" "5.174(0.741)" "2.856(0.91)" "5.856(0.677)" "1.79(0.777)"
## [,6] [,7] [,8]
## 16 "4.425(0.929)" "2.714(0.951)" "5.226(1.078)"
## 32 "4.917(1.035)" "2.705(0.891)" "4.699(1.117)"
## 64 "5.255(1.087)" "2.991(0.899)" "4.556(0.918)"
## 128 "4.808(1.052)" "3.143(0.889)" "4.21(1.02)"
library(xtable)
xtable(outtable)
## % latex table generated in R 4.2.3 by xtable 1.8-4 package
## % Mon Aug 14 11:20:18 2023
## \begin{table}[ht]
## \centering
## \begin{tabular}{rllllllll}
## \hline
## & 1 & 2 & 3 & 4 & 5 & 6 & 7 & 8 \\
## \hline
## 16 & 2.113(0.864) & 4.693(0.917) & 1.651(0.673) & 5.029(1.343) & 1.878(0.913) & 4.425(0.929) & 2.714(0.951) & 5.226(1.078) \\
## 32 & 2.14(0.826) & 4.97(0.856) & 2.029(0.823) & 4.613(0.924) & 1.543(0.548) & 4.917(1.035) & 2.705(0.891) & 4.699(1.117) \\
## 64 & 2.031(0.779) & 4.997(1.008) & 2.312(0.831) & 4.82(0.832) & 1.71(0.686) & 5.255(1.087) & 2.991(0.899) & 4.556(0.918) \\
## 128 & 2.866(0.802) & 5.174(0.741) & 2.856(0.91) & 5.856(0.677) & 1.79(0.777) & 4.808(1.052) & 3.143(0.889) & 4.21(1.02) \\
## \hline
## \end{tabular}
## \end{table}
# Prepare data for b0
df_b0 <- data.frame()
for(treetypeIndex in 1:length(treetype.array)){
for(taxaIndex in 1:length(taxa.array)){
b0_values <- tree.coef.array[treetypeIndex,taxaIndex,,1][is.finite(tree.coef.array[treetypeIndex,taxaIndex,,1])]
temp_df_b0 <- data.frame(
taxa = taxa.array[taxaIndex],
treetype = treetype.array[treetypeIndex],
value = b0_values
)
df_b0 <- rbind(df_b0, temp_df_b0)
}
}
head(df_b0)
## taxa treetype value
## 1 16 coal 2.333157
## 2 16 coal 2.333157
## 3 16 coal 2.333157
## 4 16 coal 3.549299
## 5 16 coal 1.110057
## 6 16 coal 1.649923
# Load library
library(ggplot2)
df_b0$taxa<-as.factor(df_b0$taxa)
# Create boxplot for b0 across all tree types
p_b0 <- ggplot(df_b0, aes(x=treetype, y=value, fill=taxa)) +
geom_boxplot() +
labs(x = "Tree Type",
y = expression(beta[0] ~ "Estimates"),
fill = "Taxa Size",
title = "Phylogenetic NB2 Regression:" ~ beta[0] ~ "Estimates") +
theme_minimal() + ylim(3,5)+
theme(plot.title = element_text(hjust = 0.5)) # This line centers the title
# Print the plot
print(p_b0)
## Warning: Removed 4524 rows containing non-finite values (stat_boxplot).
# png(file="phynb2simb0v2.png",
# width=600, height=350)
# print(p_b0)
# dev.off()
# Prepare data for b1
df_b1 <- data.frame()
for(treetypeIndex in 1:length(treetype.array)){
for(taxaIndex in 1:length(taxa.array)){
b1_values <- tree.coef.array[treetypeIndex,taxaIndex,,2][is.finite(tree.coef.array[treetypeIndex,taxaIndex,,2])]
temp_df_b1 <- data.frame(
taxa = taxa.array[taxaIndex],
treetype = treetype.array[treetypeIndex],
value = b1_values
)
df_b1 <- rbind(df_b1, temp_df_b1)
}
}
# Load library
library(ggplot2)
df_b1$taxa<-as.factor(df_b1$taxa)
# Create boxplot for b1 across all tree types
p_b1 <- ggplot(df_b1, aes(x=treetype, y=value, fill=taxa)) +
geom_boxplot() +
labs(x = "Tree Type",
y = expression(beta[1] ~ "Estimates"),
fill = "Taxa Size",
title = "Phylogenetic NB2 Regression:" ~ beta[1] ~ "Estimates") +
theme_minimal() +
theme(plot.title = element_text(hjust = 0.5)) # This line centers the title
# Print the plot
print(p_b1)
library(gridExtra)
grid.arrange(p_b0,p_b1,nrow=1)
## Warning: Removed 4524 rows containing non-finite values (stat_boxplot).
