rm(list=ls())
library(ggpubr)
## 載入需要的套件:ggplot2
library(dplyr)
## 
## 載入套件:'dplyr'
## 下列物件被遮斷自 'package:stats':
## 
##     filter, lag
## 下列物件被遮斷自 'package:base':
## 
##     intersect, setdiff, setequal, union
library(nonpar)

setwd("~/Dropbox/JournalSubmission/0phyregDiversity/rcode/sims/")
load("fulloutput.rda")

#load(url("https://tonyjhwueng.info/pcmreg/fulloutput.rda"))
####################
###full.raw.output##
####################

full.raw.output<-as.data.frame(full.raw.output)
full.raw.output$intercept<-as.numeric(full.raw.output$intercept)
full.raw.output$slope<-as.numeric(full.raw.output$slope)
full.raw.output$aic<-as.numeric(full.raw.output$aic)
full.raw.output$r.squared<-as.numeric(full.raw.output$r.squared)
full.raw.output$rsq<-as.numeric(full.raw.output$r.squared)
full.raw.output$sigma2<-as.numeric(full.raw.output$sigma2)
full.raw.output$mse<-as.numeric(full.raw.output$mse)
full.raw.output$size<-as.factor(full.raw.output$size)
full.raw.output$tree<-as.factor(full.raw.output$tree)
full.raw.output$model<-as.factor(full.raw.output$model)

####################
###full.sim.output##
####################

full.sim.output<-as.data.frame(full.sim.output)
full.sim.output$intercept<-as.numeric(full.sim.output$intercept)
full.sim.output$slope<-as.numeric(full.sim.output$slope)
full.sim.output$sigma2<-as.numeric(full.sim.output$sigma2)
full.sim.output$aic<-as.numeric(full.sim.output$aic)
full.sim.output$r.squared<-as.numeric(full.sim.output$r.squared)
full.sim.output$rsq<-as.numeric(full.sim.output$r.squared)
full.sim.output$mse<-as.numeric(full.sim.output$mse)
full.sim.output$size<-as.factor(full.sim.output$size)
full.sim.output$sims<-as.numeric(full.sim.output$sims)
full.sim.output$collevel<-as.factor(full.sim.output$collevel)
full.sim.output$coltype<-as.factor(full.sim.output$coltype)
full.sim.output$tree<-as.factor(full.sim.output$tree)
full.sim.output$model<-as.factor(full.sim.output$model)
full.sim.output<-full.sim.output%>%
  mutate(colbil=ifelse(collevel %in%c(1:2),1, 0))
full.sim.output$colbil<-as.factor(full.sim.output$colbil)

result0<-round(rbind(
c(median(full.raw.output$intercept,na.rm=T),
median(full.raw.output$slope,na.rm=T),
median(full.raw.output$sigma2,na.rm=T),
median(full.raw.output$aic,na.rm=T),
median(full.raw.output$r.squared,na.rm=T),
median(full.raw.output$mse,na.rm=T)),
c(median(full.sim.output$intercept,na.rm=T),
median(full.sim.output$slope,na.rm=T),
median(full.sim.output$sigma2,na.rm=T),
median(full.sim.output$aic,na.rm=T),
median(full.sim.output$r.squared,na.rm=T),
median(full.sim.output$mse,na.rm=T))
),3)

colnames(result0)<-c(
  "intercept", 
  "slope",    
  "sigma2",    
  "aic",      
  "r.squared", 
  "mse"      
)
rownames(result0)<-c("raw","sims")
library(xtable)
xtable(result0,digits=3)
## % latex table generated in R 4.2.0 by xtable 1.8-4 package
## % Sat Sep 10 20:58:25 2022
## \begin{table}[ht]
## \centering
## \begin{tabular}{rrrrrrr}
##   \hline
##  & intercept & slope & sigma2 & aic & r.squared & mse \\ 
##   \hline
## raw & 4.992 & 2.981 & 25.742 & 1114.414 & 0.324 & 2.317 \\ 
##   sims & 4.997 & 2.999 & 26.638 & 1117.378 & 0.364 & 2.328 \\ 
##    \hline
## \end{tabular}
## \end{table}
load("~/Dropbox/JournalSubmission/0phyregDiversity/rcode/sims/outputpcont.rda")

dim(output.test.p.value.array)
## [1]  4  7 11  2  8  6
dim(output.test.confint.array)
## [1]  4  7 11  2  8  6  2
head(result.intercept)
##      [,1]                               
## [1,] "BalancedTreeBM450clade1"          
## [2,] "BalancedTreeEB150clade3"          
## [3,] "BalancedTreeEB150clade5"          
## [4,] "BalancedTreeEB200clade3"          
## [5,] "BalancedTreeOUfixedRoot500clade3" 
## [6,] "BalancedTreeOUrandomRoot150clade1"
head(result.slope)
##      [,1]                    
## [1,] "BalancedTreeBM20root8" 
## [2,] "BalancedTreeBM50clade3"
## [3,] "BalancedTreeBM50clade4"
## [4,] "BalancedTreeBM50clade5"
## [5,] "BalancedTreeBM50clade6"
## [6,] "BalancedTreeBM50clade7"
head(result.sigma2)
##      [,1]                    
## [1,] "BalancedTreeBM20clade1"
## [2,] "BalancedTreeBM20clade2"
## [3,] "BalancedTreeBM20clade3"
## [4,] "BalancedTreeBM20clade4"
## [5,] "BalancedTreeBM20clade5"
## [6,] "BalancedTreeBM20clade6"
head(result.aic)
##      [,1]                    
## [1,] "BalancedTreeBM20clade1"
## [2,] "BalancedTreeBM20clade2"
## [3,] "BalancedTreeBM20clade3"
## [4,] "BalancedTreeBM20clade4"
## [5,] "BalancedTreeBM20clade5"
## [6,] "BalancedTreeBM20clade6"
head(result.rsq)
##      [,1]                    
## [1,] "BalancedTreeBM20clade1"
## [2,] "BalancedTreeBM20clade2"
## [3,] "BalancedTreeBM20clade3"
## [4,] "BalancedTreeBM20clade4"
## [5,] "BalancedTreeBM20clade5"
## [6,] "BalancedTreeBM20clade6"
head(result.mse)
##      [,1]                    
## [1,] "BalancedTreeBM20clade1"
## [2,] "BalancedTreeBM20clade2"
## [3,] "BalancedTreeBM20clade3"
## [4,] "BalancedTreeBM20clade4"
## [5,] "BalancedTreeBM20clade5"
## [6,] "BalancedTreeBM20clade6"
parameas<-c("intercept","slope","sigma2","aic","rsq","mse")

##################
# remove outlier #
# test   again   #   
##################

removeoutlier<-function(data=data){
  quartiles <- quantile(data, probs=c(.25, .75), na.rm = FALSE)
  IQR <- IQR(data)
  Lower <- quartiles[1] - 1.5*IQR
  Upper <- quartiles[2] + 1.5*IQR 
  data_no_outlier <- subset(data, data > Lower & data < Upper)
  return(list(data_no_outlier=data_no_outlier,  pos= (data > Lower & data < Upper) ))
}

##############
# intercept #
#############

full.raw.output.intercept.nooutlier<-removeoutlier(full.raw.output$intercept)
full.sim.output.intercept.nooutlier<-removeoutlier(full.sim.output$intercept)

df.tree.intercept<-data.frame(
  intercept= c(full.raw.output.intercept.nooutlier$data_no_outlier,full.sim.output.intercept.nooutlier$data_no_outlier),
  treetype = c(full.raw.output$tree[full.raw.output.intercept.nooutlier$pos],full.sim.output$tree[full.sim.output.intercept.nooutlier$pos]),
  datatype = c(rep("raw",length(full.raw.output.intercept.nooutlier$data_no_outlier)),rep("sim",length(full.sim.output.intercept.nooutlier$data_no_outlier)))
)





##############
# slope #
#############

full.raw.output.slope.nooutlier<-removeoutlier(full.raw.output$slope)
full.sim.output.slope.nooutlier<-removeoutlier(full.sim.output$slope)

df.tree.slope<-data.frame(
  slope= c(full.raw.output.slope.nooutlier$data_no_outlier,full.sim.output.slope.nooutlier$data_no_outlier),
  treetype = c(full.raw.output$tree[full.raw.output.slope.nooutlier$pos],full.sim.output$tree[full.sim.output.slope.nooutlier$pos]),
  datatype = c(rep("raw",length(full.raw.output.slope.nooutlier$data_no_outlier)),rep("sim",length(full.sim.output.slope.nooutlier$data_no_outlier)))
)



##############
#  sigma2   #
#############

full.raw.output.sigma2.nooutlier<-removeoutlier(full.raw.output$sigma2)
full.sim.output.sigma2.nooutlier<-removeoutlier(full.sim.output$sigma2)

df.tree.sigma2<-data.frame(
  sigma2= c(full.raw.output.sigma2.nooutlier$data_no_outlier,full.sim.output.sigma2.nooutlier$data_no_outlier),
  treetype = c(full.raw.output$tree[full.raw.output.sigma2.nooutlier$pos],full.sim.output$tree[full.sim.output.sigma2.nooutlier$pos]),
  datatype = c(rep("raw",length(full.raw.output.sigma2.nooutlier$data_no_outlier)),rep("sim",length(full.sim.output.sigma2.nooutlier$data_no_outlier)))
)


##############
#  aic   #
#############

full.raw.output.aic.nooutlier<-removeoutlier(full.raw.output$aic)
full.sim.output.aic.nooutlier<-removeoutlier(full.sim.output$aic)

df.tree.aic<-data.frame(
  aic= c(full.raw.output.aic.nooutlier$data_no_outlier,full.sim.output.aic.nooutlier$data_no_outlier),
  treetype = c(full.raw.output$tree[full.raw.output.aic.nooutlier$pos],full.sim.output$tree[full.sim.output.aic.nooutlier$pos]),
  datatype = c(rep("raw",length(full.raw.output.aic.nooutlier$data_no_outlier)),rep("sim",length(full.sim.output.aic.nooutlier$data_no_outlier)))
)


##############
#  r.squared #
#############

full.raw.output.r.squared.nooutlier<-removeoutlier(full.raw.output$r.squared)
full.sim.output.r.squared.nooutlier<-removeoutlier(full.sim.output$r.squared)

df.tree.r.squared<-data.frame(
  r.squared= c(full.raw.output.r.squared.nooutlier$data_no_outlier,full.sim.output.r.squared.nooutlier$data_no_outlier),
  treetype = c(full.raw.output$tree[full.raw.output.r.squared.nooutlier$pos],full.sim.output$tree[full.sim.output.r.squared.nooutlier$pos]),
  datatype = c(rep("raw",length(full.raw.output.r.squared.nooutlier$data_no_outlier)),rep("sim",length(full.sim.output.r.squared.nooutlier$data_no_outlier)))
)


##############
#   mse      #
#############

full.raw.output.mse.nooutlier<-removeoutlier(full.raw.output$mse)
full.sim.output.mse.nooutlier<-removeoutlier(full.sim.output$mse)

df.tree.mse<-data.frame(
  mse= c(full.raw.output.mse.nooutlier$data_no_outlier,full.sim.output.mse.nooutlier$data_no_outlier),
  treetype = c(full.raw.output$tree[full.raw.output.mse.nooutlier$pos],full.sim.output$tree[full.sim.output.mse.nooutlier$pos]),
  datatype = c(rep("raw",length(full.raw.output.mse.nooutlier$data_no_outlier)),rep("sim",length(full.sim.output.mse.nooutlier$data_no_outlier)))
)



# p.tree.intercep <- ggboxplot(df.tree.intercept, x = "treetype", y = "intercept", palette =c("#00AFBB", "#E7B800","green3","cyan"),add = "jitter", shape = "treetype",color="treetype")
# 
# my_comparisons <- list(c("BalancedTree", "birthdeath"),c("BalancedTree", "pbtree"),c("BalancedTree", "rcoal"),c("rcoal", "birthdeath"),c("rcoal", "pbtree"),c("birthdeath", "pbtree"))
# p.tree.intercep <- p.tree.intercep + stat_compare_means(comparisons = my_comparisons)+ stat_compare_means(label.y = 50)
# p.tree.intercep
#tree and model

colnames(full.sim.output)
##  [1] "intercept" "slope"     "sigma2"    "aic"       "r.squared" "mse"      
##  [7] "tree"      "model"     "size"      "coltype"   "collevel"  "sims"     
## [13] "rsq"       "colbil"
"intercept" 
## [1] "intercept"
"slope"    
## [1] "slope"
"sigma2"    
## [1] "sigma2"
"aic"      
## [1] "aic"
"r.squared" 
## [1] "r.squared"
"mse"      
## [1] "mse"
"tree"      
## [1] "tree"
"model"    
## [1] "model"
"size"      
## [1] "size"
"coltype"  
## [1] "coltype"
"collevel"  
## [1] "collevel"
"sims"  
## [1] "sims"
######################################################
##############  coltype and colbil ###################
######################################################


######################################################
################## intercept coltype and colbil ############
######################################################

df.intercept.raw<-full.raw.output%>%
  group_by(tree,model,size)%>%
  summarise(intercept=mean(intercept))
## `summarise()` has grouped output by 'tree', 'model'. You can override using the
## `.groups` argument.
head(df.intercept.raw)
## # A tibble: 6 × 4
## # Groups:   tree, model [1]
##   tree         model size  intercept
##   <fct>        <fct> <fct>     <dbl>
## 1 BalancedTree BM    100        4.65
## 2 BalancedTree BM    150        5.08
## 3 BalancedTree BM    20         4.80
## 4 BalancedTree BM    200        5.70
## 5 BalancedTree BM    250        6.05
## 6 BalancedTree BM    300        5.42
ls(full.sim.output)
##  [1] "aic"       "colbil"    "collevel"  "coltype"   "intercept" "model"    
##  [7] "mse"       "r.squared" "rsq"       "sigma2"    "sims"      "size"     
## [13] "slope"     "tree"
df.intercept.sim<-full.sim.output %>%
  group_by(tree,model,size,coltype,colbil)%>%
  summarise(intercept=mean(intercept))
## `summarise()` has grouped output by 'tree', 'model', 'size', 'coltype'. You can
## override using the `.groups` argument.
head(df.intercept.sim)
## # A tibble: 6 × 6
## # Groups:   tree, model, size, coltype [3]
##   tree         model size  coltype colbil intercept
##   <fct>        <fct> <fct> <fct>   <fct>      <dbl>
## 1 BalancedTree BM    100   clade   0           4.95
## 2 BalancedTree BM    100   clade   1           4.92
## 3 BalancedTree BM    100   root    0           4.91
## 4 BalancedTree BM    100   root    1           4.89
## 5 BalancedTree BM    150   clade   0           4.90
## 6 BalancedTree BM    150   clade   1           4.94
coltype.intercept<-t.test(df.intercept.sim$intercept[df.intercept.sim$coltype=="clade"],df.intercept.sim$intercept[df.intercept.sim$coltype=="root"],alternative = c("two.sided"))

colbil.intercept<-t.test(df.intercept.sim$intercept[df.intercept.sim$colbil==0],df.intercept.sim$intercept[df.intercept.sim$colbil==1],alternative = c("two.sided"))

coltypebil.intercept<-rbind(
  paste("(",round(coltype.intercept$conf.int[1],3),",",round(coltype.intercept$conf.int[2],3),"),", round(coltype.intercept$p.value,2),sep="")
  ,
  paste("(",round(colbil.intercept$conf.int[1],3),",",round(colbil.intercept$conf.int[2],3),"),", round(colbil.intercept$p.value,2),sep="")
)


######################################################
################## slope coltype and colbil ############
######################################################

df.slope.raw<-full.raw.output%>%
  group_by(tree,model,size)%>%
  summarise(slope=mean(slope))
## `summarise()` has grouped output by 'tree', 'model'. You can override using the
## `.groups` argument.
head(df.slope.raw)
## # A tibble: 6 × 4
## # Groups:   tree, model [1]
##   tree         model size  slope
##   <fct>        <fct> <fct> <dbl>
## 1 BalancedTree BM    100    4.42
## 2 BalancedTree BM    150    3.18
## 3 BalancedTree BM    20     3.08
## 4 BalancedTree BM    200    1.47
## 5 BalancedTree BM    250    1.80
## 6 BalancedTree BM    300    2.66
ls(full.sim.output)
##  [1] "aic"       "colbil"    "collevel"  "coltype"   "intercept" "model"    
##  [7] "mse"       "r.squared" "rsq"       "sigma2"    "sims"      "size"     
## [13] "slope"     "tree"
df.slope.sim<-full.sim.output %>%
  group_by(tree,model,size,coltype,colbil)%>%
  summarise(slope=mean(slope))
## `summarise()` has grouped output by 'tree', 'model', 'size', 'coltype'. You can
## override using the `.groups` argument.
head(df.slope.sim)
## # A tibble: 6 × 6
## # Groups:   tree, model, size, coltype [3]
##   tree         model size  coltype colbil slope
##   <fct>        <fct> <fct> <fct>   <fct>  <dbl>
## 1 BalancedTree BM    100   clade   0       2.82
## 2 BalancedTree BM    100   clade   1       2.83
## 3 BalancedTree BM    100   root    0       2.83
## 4 BalancedTree BM    100   root    1       2.80
## 5 BalancedTree BM    150   clade   0       3.16
## 6 BalancedTree BM    150   clade   1       3.09
coltype.slope<-t.test(df.slope.sim$slope[df.slope.sim$coltype=="clade"],df.slope.sim$slope[df.slope.sim$coltype=="root"],alternative = c("two.sided"))

colbil.slope<-t.test(df.slope.sim$slope[df.slope.sim$colbil==0],df.slope.sim$slope[df.slope.sim$colbil==1],alternative = c("two.sided"))

coltypebil.slope<-rbind(
  paste("(",round(coltype.slope$conf.int[1],3),",",round(coltype.slope$conf.int[2],3),"),", round(coltype.slope$p.value,2),sep="")
  ,
  paste("(",round(colbil.slope$conf.int[1],3),",",round(colbil.slope$conf.int[2],3),"),", round(colbil.slope$p.value,2),sep="")
)

######################################################
################## sigma2 coltype and colbil ##########
######################################################

df.sigma2.raw<-full.raw.output%>%
  group_by(tree,model,size)%>%
  summarise(sigma2=mean(sigma2))
## `summarise()` has grouped output by 'tree', 'model'. You can override using the
## `.groups` argument.
head(df.sigma2.raw)
## # A tibble: 6 × 4
## # Groups:   tree, model [1]
##   tree         model size  sigma2
##   <fct>        <fct> <fct>  <dbl>
## 1 BalancedTree BM    100    117. 
## 2 BalancedTree BM    150    196. 
## 3 BalancedTree BM    20      14.1
## 4 BalancedTree BM    200    243. 
## 5 BalancedTree BM    250    290. 
## 6 BalancedTree BM    300    343.
ls(full.sim.output)
##  [1] "aic"       "colbil"    "collevel"  "coltype"   "intercept" "model"    
##  [7] "mse"       "r.squared" "rsq"       "sigma2"    "sims"      "size"     
## [13] "slope"     "tree"
df.sigma2.sim<-full.sim.output %>%
  group_by(tree,model,size,coltype,colbil)%>%
  summarise(sigma2=mean(sigma2))
## `summarise()` has grouped output by 'tree', 'model', 'size', 'coltype'. You can
## override using the `.groups` argument.
head(df.sigma2.sim)
## # A tibble: 6 × 6
## # Groups:   tree, model, size, coltype [3]
##   tree         model size  coltype colbil sigma2
##   <fct>        <fct> <fct> <fct>   <fct>   <dbl>
## 1 BalancedTree BM    100   clade   0       104. 
## 2 BalancedTree BM    100   clade   1        68.8
## 3 BalancedTree BM    100   root    0        98.1
## 4 BalancedTree BM    100   root    1       108. 
## 5 BalancedTree BM    150   clade   0       172. 
## 6 BalancedTree BM    150   clade   1       114.
coltype.sigma2<-t.test(df.sigma2.sim$sigma2[df.sigma2.sim$coltype=="clade"],df.sigma2.sim$sigma2[df.sigma2.sim$coltype=="root"],alternative = c("two.sided"))

colbil.sigma2<-t.test(df.sigma2.sim$sigma2[df.sigma2.sim$colbil==0],df.sigma2.sim$sigma2[df.sigma2.sim$colbil==1],alternative = c("two.sided"))

coltypebil.sigma2<-rbind(
  paste("(",round(coltype.sigma2$conf.int[1],3),",",round(coltype.sigma2$conf.int[2],3),"),", round(coltype.sigma2$p.value,2),sep="")
  ,
  paste("(",round(colbil.sigma2$conf.int[1],3),",",round(colbil.sigma2$conf.int[2],3),"),", round(colbil.sigma2$p.value,2),sep="")
)


######################################################
################## rsq coltype and colbil ############
######################################################

df.rsq.raw<-full.raw.output%>%
  group_by(tree,model,size)%>%
  summarise(rsq=mean(r.squared))
## `summarise()` has grouped output by 'tree', 'model'. You can override using the
## `.groups` argument.
head(df.rsq.raw)
## # A tibble: 6 × 4
## # Groups:   tree, model [1]
##   tree         model size      rsq
##   <fct>        <fct> <fct>   <dbl>
## 1 BalancedTree BM    100   0.171  
## 2 BalancedTree BM    150   0.0675 
## 3 BalancedTree BM    20    0.527  
## 4 BalancedTree BM    200   0.00956
## 5 BalancedTree BM    250   0.0111 
## 6 BalancedTree BM    300   0.0201
ls(full.sim.output)
##  [1] "aic"       "colbil"    "collevel"  "coltype"   "intercept" "model"    
##  [7] "mse"       "r.squared" "rsq"       "sigma2"    "sims"      "size"     
## [13] "slope"     "tree"
df.rsq.sim<-full.sim.output %>%
  group_by(tree,model,size,coltype,colbil)%>%
  summarise(rsq=mean(r.squared))
## `summarise()` has grouped output by 'tree', 'model', 'size', 'coltype'. You can
## override using the `.groups` argument.
head(df.rsq.sim)
## # A tibble: 6 × 6
## # Groups:   tree, model, size, coltype [3]
##   tree         model size  coltype colbil    rsq
##   <fct>        <fct> <fct> <fct>   <fct>   <dbl>
## 1 BalancedTree BM    100   clade   0      0.0911
## 2 BalancedTree BM    100   clade   1      0.124 
## 3 BalancedTree BM    100   root    0      0.0914
## 4 BalancedTree BM    100   root    1      0.0872
## 5 BalancedTree BM    150   clade   0      0.0701
## 6 BalancedTree BM    150   clade   1      0.0869
coltype.rsq<-t.test(df.rsq.sim$rsq[df.rsq.sim$coltype=="clade"],df.rsq.sim$rsq[df.rsq.sim$coltype=="root"],alternative = c("two.sided"))

colbil.rsq<-t.test(df.rsq.sim$rsq[df.rsq.sim$colbil==0],df.rsq.sim$rsq[df.rsq.sim$colbil==1],alternative = c("two.sided"))

coltypebil.rsq<-rbind(
  paste("(",round(coltype.rsq$conf.int[1],3),",",round(coltype.rsq$conf.int[2],3),"),", round(coltype.rsq$p.value,2),sep="")
  ,
  paste("(",round(colbil.rsq$conf.int[1],3),",",round(colbil.rsq$conf.int[2],3),"),", round(colbil.rsq$p.value,2),sep="")
)



######################################################
################## aic coltype and colbil ############
######################################################

df.aic.raw<-full.raw.output%>%
  group_by(tree,model,size)%>%
  summarise(aic=mean(aic))
## `summarise()` has grouped output by 'tree', 'model'. You can override using the
## `.groups` argument.
head(df.aic.raw)
## # A tibble: 6 × 4
## # Groups:   tree, model [1]
##   tree         model size     aic
##   <fct>        <fct> <fct>  <dbl>
## 1 BalancedTree BM    100    459. 
## 2 BalancedTree BM    150    702. 
## 3 BalancedTree BM    20      85.1
## 4 BalancedTree BM    200    923. 
## 5 BalancedTree BM    250   1130. 
## 6 BalancedTree BM    300   1358.
ls(full.sim.output)
##  [1] "aic"       "colbil"    "collevel"  "coltype"   "intercept" "model"    
##  [7] "mse"       "r.squared" "rsq"       "sigma2"    "sims"      "size"     
## [13] "slope"     "tree"
df.aic.sim<-full.sim.output %>%
  group_by(tree,model,size,coltype,colbil)%>%
  summarise(aic=mean(aic))
## `summarise()` has grouped output by 'tree', 'model', 'size', 'coltype'. You can
## override using the `.groups` argument.
head(df.aic.sim)
## # A tibble: 6 × 6
## # Groups:   tree, model, size, coltype [3]
##   tree         model size  coltype colbil   aic
##   <fct>        <fct> <fct> <fct>   <fct>  <dbl>
## 1 BalancedTree BM    100   clade   0       450.
## 2 BalancedTree BM    100   clade   1       485.
## 3 BalancedTree BM    100   root    0       470.
## 4 BalancedTree BM    100   root    1       453.
## 5 BalancedTree BM    150   clade   0       687.
## 6 BalancedTree BM    150   clade   1       762.
coltype.aic<-t.test(df.aic.sim$aic[df.aic.sim$coltype=="clade"],df.aic.sim$aic[df.aic.sim$coltype=="root"],alternative = c("two.sided"))

colbil.aic<-t.test(df.aic.sim$aic[df.aic.sim$colbil==0],df.aic.sim$aic[df.aic.sim$colbil==1],alternative = c("two.sided"))

coltypebil.aic<-rbind(
  paste("(",round(coltype.aic$conf.int[1],3),",",round(coltype.aic$conf.int[2],3),"),", round(coltype.aic$p.value,2),sep="")
  ,
  paste("(",round(colbil.aic$conf.int[1],3),",",round(colbil.aic$conf.int[2],3),"),", round(colbil.aic$p.value,2),sep="")
)




######################################################
################## mse coltype and colbil ############
######################################################

df.mse.raw<-full.raw.output%>%
  group_by(tree,model,size)%>%
  summarise(mse=mean(mse))
## `summarise()` has grouped output by 'tree', 'model'. You can override using the
## `.groups` argument.
head(df.mse.raw)
## # A tibble: 6 × 4
## # Groups:   tree, model [1]
##   tree         model size    mse
##   <fct>        <fct> <fct> <dbl>
## 1 BalancedTree BM    100    2.61
## 2 BalancedTree BM    150    2.35
## 3 BalancedTree BM    20     1.45
## 4 BalancedTree BM    200    3.28
## 5 BalancedTree BM    250    3.08
## 6 BalancedTree BM    300    2.31
ls(full.sim.output)
##  [1] "aic"       "colbil"    "collevel"  "coltype"   "intercept" "model"    
##  [7] "mse"       "r.squared" "rsq"       "sigma2"    "sims"      "size"     
## [13] "slope"     "tree"
df.mse.sim<-full.sim.output %>%
  group_by(tree,model,size,coltype,colbil)%>%
  summarise(mse=mean(mse))
## `summarise()` has grouped output by 'tree', 'model', 'size', 'coltype'. You can
## override using the `.groups` argument.
head(df.mse.sim)
## # A tibble: 6 × 6
## # Groups:   tree, model, size, coltype [3]
##   tree         model size  coltype colbil   mse
##   <fct>        <fct> <fct> <fct>   <fct>  <dbl>
## 1 BalancedTree BM    100   clade   0       2.92
## 2 BalancedTree BM    100   clade   1       2.87
## 3 BalancedTree BM    100   root    0       3.24
## 4 BalancedTree BM    100   root    1       3.07
## 5 BalancedTree BM    150   clade   0       2.91
## 6 BalancedTree BM    150   clade   1       2.81
t.test(df.mse.sim$mse[df.mse.sim$coltype=="clade"],df.mse.sim$mse[df.mse.sim$coltype=="root"],alternative = c("two.sided"))
## 
##  Welch Two Sample t-test
## 
## data:  df.mse.sim$mse[df.mse.sim$coltype == "clade"] and df.mse.sim$mse[df.mse.sim$coltype == "root"]
## t = 0.13023, df = 1141.5, p-value = 0.8964
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -2.751953  3.143235
## sample estimates:
## mean of x mean of y 
##  5.673452  5.477810
t.test(df.mse.sim$mse[df.mse.sim$colbil==0],df.mse.sim$mse[df.mse.sim$colbil==1],alternative = c("two.sided"))
## 
##  Welch Two Sample t-test
## 
## data:  df.mse.sim$mse[df.mse.sim$colbil == 0] and df.mse.sim$mse[df.mse.sim$colbil == 1]
## t = -0.21882, df = 1128, p-value = 0.8268
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -3.276331  2.618858
## sample estimates:
## mean of x mean of y 
##  5.411263  5.739999
coltype.mse<-t.test(df.mse.sim$mse[df.mse.sim$coltype=="clade"],df.mse.sim$mse[df.mse.sim$coltype=="root"],alternative = c("two.sided"))

colbil.mse<-t.test(df.mse.sim$mse[df.mse.sim$colbil==0],df.mse.sim$mse[df.mse.sim$colbil==1],alternative = c("two.sided"))

coltypebil.mse<-rbind(
  paste("(",round(coltype.mse$conf.int[1],3),",",round(coltype.mse$conf.int[2],3),"),", round(coltype.mse$p.value,2),sep="")
  ,
  paste("(",round(colbil.mse$conf.int[1],3),",",round(colbil.mse$conf.int[2],3),"),", round(colbil.mse$p.value,2),sep="")
)


#################################################




coltypebil.full<-cbind(
  coltypebil.intercept,
  coltypebil.slope,
  coltypebil.sigma2,
  coltypebil.aic,
  coltypebil.rsq,
  coltypebil.mse
)

rownames(coltypebil.full)<-c("type","level")
colnames(coltypebil.full)<-c("intercept","slope","sigma2","aic","rsq","mse")


xtable(t(coltypebil.full))
## % latex table generated in R 4.2.0 by xtable 1.8-4 package
## % Sat Sep 10 20:58:27 2022
## \begin{table}[ht]
## \centering
## \begin{tabular}{rll}
##   \hline
##  & type & level \\ 
##   \hline
## intercept & (-0.014,0.023),0.63 & (-0.023,0.013),0.61 \\ 
##   slope & (-0.006,0.021),0.28 & (-0.019,0.008),0.44 \\ 
##   sigma2 & (-531.447,15.805),0.06 & (-179.965,367.877),0.5 \\ 
##   aic & (-77.071,101.029),0.79 & (-91.487,86.618),0.96 \\ 
##   rsq & (0.011,0.086),0.01 & (-0.065,0.01),0.15 \\ 
##   mse & (-2.752,3.143),0.9 & (-3.276,2.619),0.83 \\ 
##    \hline
## \end{tabular}
## \end{table}