rm(list=ls())
source("~/Dropbox/FCU/Teaching/Mentoring/2022Spring-Current_Professor/2022Fall/YourueiMin/rcode_dcj/empricaldataanalysis/mainPhyarima.R")

ntaxa<-8
stree<-stree(ntaxa)
tree<-compute.brlen(stree(ntaxa,type="balanced"))
#tree<-reorder(tree,"postorder")
fastBM(tree,internal=T,sig2=1)->phen
phen[1:ntaxa]->y
trait<-round(y,2)
trait
edgelength<-tree$edge.length
#names(edgelength)<-edge_labels
#edgelength
betas.rrphylo <- compute_betas_rrphylo(tree, y)$betas.rrphylo
betas.rrphylo

# Rename node labels
node_labels <- paste("nd", (ntaxa+1):(2*ntaxa-1),sep="")
plot(tree,  show.tip.label = FALSE , edge.width = 2, cex = 3)
title(main = "Raw Tree", cex.main = 1.5)
nodelabels(node_labels, cex = 1.5)
# Add tip labels
tiplabels(paste("y",1:ntaxa,sep=""), cex = 1.5)
# Get the number of edges
num_edges <- Nedge(tree)
# Create edge labels
edge_labels <- paste0("eg", seq_len(num_edges),sep="")
edgelabels(edge_labels, cex = 1.5)

df.trait<-data.frame(trait=t(trait))
colnames(df.trait)<-paste("y",1:length(trait),sep="")
df.trait
xtable(df.trait)

df.edgelen<-data.frame(edgelength=round(t(edgelength),3))
colnames(df.edgelen)<-paste("eg",1:length(df.edgelen),sep="")
df.edgelen
xtable(df.edgelen,digits=3)

df.betas.rrphylo<-data.frame(betas.rrphylo=t(round(betas.rrphylo,3)))
colnames(df.betas.rrphylo)<-paste("nd",1:length(df.betas.rrphylo),sep="")
df.betas.rrphylo
xtable(df.betas.rrphylo,digits=3)



#######################################
############ Formal Analysis ##########
#######################################

mfit<-modelfit(trait=y, tree=tree) 
empmodelout<-empmodel(mfit)

empmodelout$parammleset
empmodelout$fitset

xtable(empmodelout$parammleset,digits=3)
xtable(empmodelout$fitset,digits=3)

testrate2subtrees_output <- testrate2subtrees(tree, y) 
tworatedf<- treeARanova(testrate2subtrees_output)
tworatedf
xtable(tworatedf,digits=3)