---
title: "phynb2simb0b1box"
output: html_document
date: "2023-08-14"
---
```{r}
rm(list=ls())
load("~/Dropbox/FCU/Teaching/Mentoring/2016Spring-2021Fall_AssociateProfessor/2020Fall/ChiYuWu/rcodeCYW/phynb2sim3.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
dim(tree.coef.array)
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=""))
print(outtable)

library(xtable)
xtable(outtable)

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

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