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}