rm(list=ls())
#install.packages("rstatix")
library(ggpubr)
## Loading required package: ggplot2
library(dplyr)
##
## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':
##
## filter, lag
## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, union
library(nonpar)
library(ggplot2)
library(gridExtra)
##
## Attaching package: 'gridExtra'
## The following object is masked from 'package:dplyr':
##
## combine
setwd("~/Dropbox/JournalSubmission/0.0SuccessProjects/2022phyregDiversity/rcode/sims/")
load("fulloutput.rda")
load("~/Dropbox/JournalSubmission/0.0SuccessProjects/2022phyregDiversity/rcode/sims/outputpcont.rda")
####################
###full.raw.output##
####################
full.raw.output<-as.data.frame(full.raw.output)
plot(full.raw.output$intercept)
plot(full.raw.output$slope)
dim(full.raw.output)
## [1] 308 9
colnames(full.raw.output)
## [1] "intercept" "slope" "sigma2" "aic" "r.squared" "mse"
## [7] "tree" "model" "size"
head(full.raw.output)
## intercept slope sigma2
## temp.raw.output 4.80309973352035 3.07518810263625 14.0891381228433
## temp.raw.output.1 4.41289062354188 3.73978241135478 38.6038994255435
## temp.raw.output.2 4.64827179569421 4.42202434896682 116.596741729584
## temp.raw.output.3 5.07580018141365 3.17879430256683 196.231932640169
## temp.raw.output.4 5.70068359866706 1.47468886836927 243.107832595785
## temp.raw.output.5 6.04881233760031 1.79841933017121 289.801278304944
## aic r.squared mse
## temp.raw.output 85.0571351335261 0.527048806731055 1.44659670038882
## temp.raw.output.1 211.118052892473 0.319684303757031 2.02040642737371
## temp.raw.output.2 459.216557939856 0.171107843763553 2.61152272832044
## temp.raw.output.3 701.62767476371 0.0674588476749469 2.35230546260522
## temp.raw.output.4 922.525717527539 0.00955801251463874 3.28112727766279
## temp.raw.output.5 1129.70911672311 0.0111203524902557 3.07892268033534
## tree model size
## temp.raw.output BalancedTree BM 20
## temp.raw.output.1 BalancedTree BM 50
## temp.raw.output.2 BalancedTree BM 100
## temp.raw.output.3 BalancedTree BM 150
## temp.raw.output.4 BalancedTree BM 200
## temp.raw.output.5 BalancedTree BM 250
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$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)
str(full.raw.output)
## 'data.frame': 308 obs. of 9 variables:
## $ intercept: num 4.8 4.41 4.65 5.08 5.7 ...
## $ slope : num 3.08 3.74 4.42 3.18 1.47 ...
## $ sigma2 : num 14.1 38.6 116.6 196.2 243.1 ...
## $ aic : num 85.1 211.1 459.2 701.6 922.5 ...
## $ r.squared: num 0.52705 0.31968 0.17111 0.06746 0.00956 ...
## $ mse : num 1.45 2.02 2.61 2.35 3.28 ...
## $ tree : Factor w/ 4 levels "BalancedTree",..: 1 1 1 1 1 1 1 1 1 1 ...
## $ model : Factor w/ 7 levels "BM","delta","EB",..: 1 1 1 1 1 1 1 1 1 1 ...
## $ size : Factor w/ 11 levels "100","150","20",..: 3 10 1 2 4 5 6 7 8 9 ...
ls(full.raw.output)
## [1] "aic" "intercept" "model" "mse" "r.squared" "sigma2"
## [7] "size" "slope" "tree"
####################
###full.sim.output##
####################
head(full.sim.output)
## intercept slope sigma2
## temp.sim.clade.output "4.52970746363137" "4.0402783799179" "11.281103254218"
## temp.sim.root.output "4.61755031685766" "3.73183075438454" "24.2752871835689"
## temp.sim.clade.output "4.03380673311132" "3.67543134070236" "16.6221528686918"
## temp.sim.root.output "4.48323972916464" "3.26835306446401" "23.4297256659566"
## temp.sim.clade.output "4.97649883813443" "3.04383334651632" "8.74081904505687"
## temp.sim.root.output "4.83391799459114" "3.210843499219" "8.91847058127604"
## aic r.squared mse
## temp.sim.clade.output "90.937856921515" "0.596650814845572" "2.38209113851361"
## temp.sim.root.output "97.5337287023863" "0.382736192466198" "2.1248245133016"
## temp.sim.clade.output "98.690002288419" "0.430769786524581" "2.75999961502552"
## temp.sim.root.output "96.8246626874404" "0.286404316904917" "2.82911020382223"
## temp.sim.clade.output "85.835353921787" "0.635626204758499" "2.42642099643274"
## temp.sim.root.output "77.5070412595935" "0.600091773475409" "2.24448910492674"
## tree model size coltype collevel sims
## temp.sim.clade.output "BalancedTree" "BM" "20" "clade" "1" "1"
## temp.sim.root.output "BalancedTree" "BM" "20" "root" "1" "1"
## temp.sim.clade.output "BalancedTree" "BM" "20" "clade" "1" "2"
## temp.sim.root.output "BalancedTree" "BM" "20" "root" "1" "2"
## temp.sim.clade.output "BalancedTree" "BM" "20" "clade" "1" "3"
## temp.sim.root.output "BalancedTree" "BM" "20" "root" "1" "3"
full.sim.output<-as.data.frame(full.sim.output)
dim(full.sim.output)
## [1] 492800 12
colnames(full.sim.output)
## [1] "intercept" "slope" "sigma2" "aic" "r.squared" "mse"
## [7] "tree" "model" "size" "coltype" "collevel" "sims"
head(full.sim.output)
## intercept slope sigma2
## temp.sim.clade.output 4.52970746363137 4.0402783799179 11.281103254218
## temp.sim.root.output 4.61755031685766 3.73183075438454 24.2752871835689
## temp.sim.clade.output.1 4.03380673311132 3.67543134070236 16.6221528686918
## temp.sim.root.output.1 4.48323972916464 3.26835306446401 23.4297256659566
## temp.sim.clade.output.2 4.97649883813443 3.04383334651632 8.74081904505687
## temp.sim.root.output.2 4.83391799459114 3.210843499219 8.91847058127604
## aic r.squared mse
## temp.sim.clade.output 90.937856921515 0.596650814845572 2.38209113851361
## temp.sim.root.output 97.5337287023863 0.382736192466198 2.1248245133016
## temp.sim.clade.output.1 98.690002288419 0.430769786524581 2.75999961502552
## temp.sim.root.output.1 96.8246626874404 0.286404316904917 2.82911020382223
## temp.sim.clade.output.2 85.835353921787 0.635626204758499 2.42642099643274
## temp.sim.root.output.2 77.5070412595935 0.600091773475409 2.24448910492674
## tree model size coltype collevel sims
## temp.sim.clade.output BalancedTree BM 20 clade 1 1
## temp.sim.root.output BalancedTree BM 20 root 1 1
## temp.sim.clade.output.1 BalancedTree BM 20 clade 1 2
## temp.sim.root.output.1 BalancedTree BM 20 root 1 2
## temp.sim.clade.output.2 BalancedTree BM 20 clade 1 3
## temp.sim.root.output.2 BalancedTree BM 20 root 1 3
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$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)
#################
### ggboxplot ###
#################
full.raw.output$tree <- factor(full.raw.output$tree, c("birthdeath", "pbtree", "BalancedTree", "rcoal"))
full.sim.output$tree <- factor(full.sim.output$tree, c("birthdeath", "pbtree", "BalancedTree", "rcoal"))
legend_title <- "Polytomy Level"
###############################
#### This one for a, b, sigma2#
###############################
ylim1 = boxplot.stats(full.sim.output$intercept)$stats[c(1, 5)]
p1.intercept<-ggboxplot(full.sim.output,y="intercept",x="colbil",fill="coltype")+scale_y_continuous(trans='log10') +xlab("polytomy levels")+ylab("values") +ggtitle("Intercept a") + theme(plot.title = element_text(hjust = 0.5))+scale_fill_discrete(labels=c('clade','root'))+ scale_x_discrete(breaks=c("0","1"),labels=c("High", "Low"))+scale_fill_manual(legend_title,values=c("orange","blue"))+ coord_cartesian(ylim = ylim1)
## Scale for 'fill' is already present. Adding another scale for 'fill', which
## will replace the existing scale.
ylim1 = boxplot.stats(full.sim.output$slope)$stats[c(1, 5)]
p1.slope<-ggboxplot(full.sim.output,y="slope",x="colbil",fill="coltype")+scale_y_continuous(trans='log10') +xlab("polytomy levels")+ylab("values") +ggtitle("Slope b") + theme(plot.title = element_text(hjust = 0.5))+scale_fill_discrete(labels=c('clade','root'))+ scale_x_discrete(breaks=c("0","1"),labels=c("High", "Low"))+scale_fill_manual(legend_title,values=c("orange","blue"))+ coord_cartesian(ylim = ylim1)
## Scale for 'fill' is already present. Adding another scale for 'fill', which
## will replace the existing scale.
ylim1 = boxplot.stats(full.sim.output$sigma2)$stats[c(1, 5)]
p1.sigma2<-ggboxplot(full.sim.output,y="sigma2",x="colbil",fill="coltype")+xlab("polytomy levels")+ylab("values") +ggtitle(expression(paste("Rate square ", sigma^2,sep=""))) + theme(plot.title = element_text(hjust = 0.5))+scale_fill_discrete(labels=c('clade','root'))+ scale_x_discrete(breaks=c("0","1"),labels=c("High", "Low"))+scale_fill_manual(legend_title,values=c("orange","blue"))+ coord_cartesian(ylim = ylim1)+scale_y_continuous(trans=scales::pseudo_log_trans(base = 10),breaks=c(1, 10, 100, 1000, 10000, 100000), labels = expression(1, 10, 10^2, 10^3, 10^4, 10^5),expand = c(0, 0))
## Scale for 'fill' is already present. Adding another scale for 'fill', which
## will replace the existing scale.
grid.arrange(p1.intercept,p1.slope,p1.sigma2,ncol=3)
## Warning in self$trans$transform(x): NaNs produced
## Warning: Transformation introduced infinite values in continuous y-axis
## Warning: Removed 2742 rows containing non-finite values (stat_boxplot).
## Warning in self$trans$transform(x): NaNs produced
## Warning: Transformation introduced infinite values in continuous y-axis
## Warning: Removed 8192 rows containing non-finite values (stat_boxplot).
####################################################
##aic so no difference between polytype and level###
####################################################
ggviolin(full.sim.output,y="aic",x="colbil",fill="coltype") +xlab("polytomy levels") +ggtitle("AIC across polytomy types and levels") + theme(plot.title = element_text(hjust = 0.5))+scale_fill_discrete(labels=c('clade','root'))+ scale_x_discrete(breaks=c("0","1"),labels=c("High", "Low"))+scale_fill_manual(legend_title,values=c("orange","red"))+ scale_y_continuous(trans='log10',labels=function(x) format(x, big.mark = ",", scientific = FALSE))
## Scale for 'fill' is already present. Adding another scale for 'fill', which
## will replace the existing scale.
###############################
#### This one for rsq##########
###############################
ggboxplot(full.sim.output,y="r.squared",x="colbil",fill="coltype") +xlab("polytomy levels")+ ylab(expression(r^2)) +ggtitle(expression(paste(r^2," across polytomy types and level ",sep=""))) + theme(plot.title = element_text(hjust = 0.5))+scale_fill_discrete(labels=c('clade','root'))+ scale_x_discrete(breaks=c("0","1"),labels=c("High", "Low"))+scale_fill_manual(legend_title,values=c("orange","red"))
## Scale for 'fill' is already present. Adding another scale for 'fill', which
## will replace the existing scale.
###############################
#### This one for mse##########
###############################
ggboxplot(full.sim.output,y="mse",x="colbil",fill="coltype") +xlab("polytomy levels")+ ylab("MSE") +ggtitle(expression(paste("MSE"," across polytomy types and level ",sep=""))) + theme(plot.title = element_text(hjust = 0.5))+scale_fill_discrete(labels=c('clade','root'))+ scale_x_discrete(breaks=c("0","1"),labels=c("High", "Low"))+scale_fill_manual(legend_title,values=c("orange","red"))+ylim(1,3)
## Scale for 'fill' is already present. Adding another scale for 'fill', which
## will replace the existing scale.
## Warning: Removed 76077 rows containing non-finite values (stat_boxplot).
ggviolin(full.sim.output,y="mse",x="colbil",fill="coltype", draw_quantiles = 0.5) +xlab("polytomy levels")+ ylab("MSE") +ggtitle(expression(paste("MSE"," across polytomy types and level ",sep=""))) + theme(plot.title = element_text(hjust = 0.5))+scale_fill_discrete(labels=c('clade','root'))+ scale_x_discrete(breaks=c("0","1"),labels=c("High", "Low"))+scale_fill_manual(legend_title,values=c("orange","red"))+ylim(1,3)
## Scale for 'fill' is already present. Adding another scale for 'fill', which
## will replace the existing scale.
## Warning: Removed 76077 rows containing non-finite values (stat_ydensity).
## Warning: Removed 136 rows containing missing values (geom_violin).
######################
#########tree########
#####################
full.sim.output$tree<-factor(full.sim.output$tree,levels=c("birthdeath","pbtree","BalancedTree", "rcoal"))
##aic
# this is due to the size incease with the likeluhood hence increase the likelihood
ggboxplot(full.sim.output,y="aic",x="coltype",fill="tree")+xlab("Tree") + ylab("AIC")+ggtitle(expression(paste("AIC"," across trees and polytomy types ",sep=""))) + theme(plot.title = element_text(hjust = 0.5))+scale_fill_discrete(labels=c('birth death','pure birth','balanced','coal'))+xlab("polytomy types")+ scale_y_continuous(trans='log10',labels=function(x) format(x, big.mark = ",", scientific = FALSE))
##rsq
# this is due to the size incease with the likeluhood hence increase the likelihood
ggboxplot(full.sim.output,y="r.squared",x="coltype",fill="tree")+xlab("Tree") + ylab(expression(r^2))+ggtitle(expression(paste(r^2," across trees and polytomy types ",sep=""))) + theme(plot.title = element_text(hjust = 0.5))+scale_fill_discrete(labels=c('birth death','pure birth','balanced','coal'))+xlab("polytomy types")
################################################
######## mse This one for poly level############
################################################
# compute lower and upper whiskers
ylim1 = boxplot.stats(full.sim.output$mse)$stats[c(1, 5)]
ggboxplot(full.sim.output,y="mse",x="colbil",fill="tree")+scale_y_continuous(trans='log10') + coord_cartesian(ylim = ylim1*1.10)+xlab("polytomy levels") +ggtitle("MSE across trees and polytomy levels") + theme(plot.title = element_text(hjust = 0.5))+scale_fill_discrete(labels=c('birth death','pure birth','balanced','coal'))+ scale_x_discrete(breaks=c("0","1"),labels=c("High", "Low"))
###############################
#### This one taking model ####
###############################
#############################
###############Intercept#####
#############################
ylim1 = boxplot.stats(full.sim.output$intercept)$stats[c(1, 5)]
full.sim.output$model<-factor(full.sim.output$model, levels=c("delta","EB","BM","kappa","OUrandomRoot","OUfixedRoot","lambda"))
p1.intercept.colbil<-ggboxplot(full.sim.output,y="intercept",x="colbil",fill="model")+scale_y_continuous(trans='log10')+ scale_x_discrete(breaks=c("0","1"),labels=c("High", "Low")) +xlab("polytomy levels")+ylab("values") +ggtitle("Intercept a") + theme(plot.title = element_text(hjust = 0.5))+ coord_cartesian(ylim = ylim1)
p1.intercept.colbil
## Warning in self$trans$transform(x): NaNs produced
## Warning: Transformation introduced infinite values in continuous y-axis
## Warning: Removed 2742 rows containing non-finite values (stat_boxplot).
p1.intercept.coltype<-ggboxplot(full.sim.output,y="intercept",x="coltype",fill="model")+scale_y_continuous(trans='log10') +xlab("polytomy types")+ylab("values") +ggtitle("Intercept a") + theme(plot.title = element_text(hjust = 0.5))+ coord_cartesian(ylim = ylim1)
p1.intercept.coltype
## Warning in self$trans$transform(x): NaNs produced
## Warning: Transformation introduced infinite values in continuous y-axis
## Warning: Removed 2742 rows containing non-finite values (stat_boxplot).
grid.arrange(p1.intercept.coltype,p1.intercept.colbil,ncol=2)
## Warning in self$trans$transform(x): NaNs produced
## Warning: Transformation introduced infinite values in continuous y-axis
## Warning: Removed 2742 rows containing non-finite values (stat_boxplot).
## Warning in self$trans$transform(x): NaNs produced
## Warning: Transformation introduced infinite values in continuous y-axis
## Warning: Removed 2742 rows containing non-finite values (stat_boxplot).
#############################
###############slope#########
#############################
ylim1 = boxplot.stats(full.sim.output$slope)$stats[c(1, 5)]
full.sim.output$model<-factor(full.sim.output$model, levels=c("delta","EB","BM","kappa","OUrandomRoot","OUfixedRoot","lambda"))
p1.slope.colbil<-ggboxplot(full.sim.output,y="slope",x="colbil",fill="model")+scale_y_continuous(trans='log10')+ scale_x_discrete(breaks=c("0","1"),labels=c("High", "Low")) +xlab("polytomy levels")+ylab("values") +ggtitle("Slope b") + theme(plot.title = element_text(hjust = 0.5))+ coord_cartesian(ylim = ylim1)
p1.slope.colbil
## Warning in self$trans$transform(x): NaNs produced
## Warning: Transformation introduced infinite values in continuous y-axis
## Warning: Removed 8192 rows containing non-finite values (stat_boxplot).
p1.slope.coltype<-ggboxplot(full.sim.output,y="slope",x="coltype",fill="model")+scale_y_continuous(trans='log10') +xlab("polytomy types")+ylab("values") +ggtitle("Slope b") + theme(plot.title = element_text(hjust = 0.5))+ coord_cartesian(ylim = ylim1)
p1.slope.coltype
## Warning in self$trans$transform(x): NaNs produced
## Warning: Transformation introduced infinite values in continuous y-axis
## Warning: Removed 8192 rows containing non-finite values (stat_boxplot).
grid.arrange(p1.slope.coltype,p1.slope.colbil,ncol=2)
## Warning in self$trans$transform(x): NaNs produced
## Warning: Transformation introduced infinite values in continuous y-axis
## Warning: Removed 8192 rows containing non-finite values (stat_boxplot).
## Warning in self$trans$transform(x): NaNs produced
## Warning: Transformation introduced infinite values in continuous y-axis
## Warning: Removed 8192 rows containing non-finite values (stat_boxplot).
#############################
###############sigma2#########
#############################
ylim1 = boxplot.stats(full.sim.output$sigma2)$stats[c(1, 5)]
full.sim.output$model<-factor(full.sim.output$model, levels=c("OUrandomRoot","OUfixedRoot", "EB","BM","delta","lambda","kappa"))
p1.sigma2.colbil<-ggboxplot(full.sim.output,y="sigma2",x="colbil",fill="model")+scale_y_continuous(trans='log10')+ scale_x_discrete(breaks=c("0","1"),labels=c("High", "Low")) +xlab("polytomy levels")+ylab("values") +ggtitle(expression(paste("rate ", sigma,sep=""))) + theme(plot.title = element_text(hjust = 0.5))+ coord_cartesian(ylim = ylim1)
p1.sigma2.colbil
p1.sigma2.coltype<-ggboxplot(full.sim.output,y="sigma2",x="coltype",fill="model")+scale_y_continuous(trans='log10') +xlab("polytomy types")+ylab("values") +ggtitle(expression(paste("rate ", sigma,sep=""))) + theme(plot.title = element_text(hjust = 0.5))+ coord_cartesian(ylim = ylim1)
p1.sigma2.coltype
grid.arrange(p1.sigma2.coltype,p1.sigma2.colbil,ncol=2)
#############################
###############aic#########
#############################
ylim1 = boxplot.stats(full.sim.output$aic)$stats[c(1, 5)]
full.sim.output$model<-factor(full.sim.output$model, levels=c("delta","EB","BM","kappa","OUrandomRoot","OUfixedRoot","lambda"))
p1.aic.colbil<-ggboxplot(full.sim.output,y="aic",x="colbil",fill="model")+scale_y_continuous(trans='log10')+ scale_x_discrete(breaks=c("0","1"),labels=c("High", "Low")) +xlab("polytomy levels")+ylab("values") +ggtitle("aic") + theme(plot.title = element_text(hjust = 0.5))+ coord_cartesian(ylim = ylim1)
p1.aic.colbil
p1.aic.coltype<-ggboxplot(full.sim.output,y="aic",x="coltype",fill="model")+scale_y_continuous(trans='log10') +xlab("polytomy types")+ylab("values") +ggtitle("aic") + theme(plot.title = element_text(hjust = 0.5))+ coord_cartesian(ylim = ylim1)
p1.aic.coltype
grid.arrange(p1.aic.coltype,p1.aic.colbil,ncol=2)
#############################
###############r.squared#########
#############################
full.sim.output$model<-factor(full.sim.output$model, levels=c("delta","EB","BM","kappa","OUrandomRoot","OUfixedRoot","lambda"))
p1.r.squared.colbil<-ggboxplot(full.sim.output,y="r.squared",x="colbil",fill="model")+ scale_x_discrete(breaks=c("0","1"),labels=c("High", "Low")) +xlab("polytomy levels")+ylab("values") +ggtitle(expression(r^2)) + theme(plot.title = element_text(hjust = 0.5))+ylim(0,1)
p1.r.squared.colbil
p1.r.squared.coltype<-ggboxplot(full.sim.output,y="r.squared",x="coltype",fill="model") +xlab("polytomy types")+ylab("values") +ggtitle(expression(r^2)) + theme(plot.title = element_text(hjust = 0.5))+ylim(0,1)
p1.r.squared.coltype
grid.arrange(p1.r.squared.coltype,p1.r.squared.colbil,ncol=2)
#############################
###############mse#########
#############################
ylim1 = boxplot.stats(full.sim.output$mse)$stats[c(1, 5)]
full.sim.output$model<-factor(full.sim.output$model, levels=c("delta","EB","BM","kappa","OUrandomRoot","OUfixedRoot","lambda"))
p1.mse.colbil<-ggboxplot(full.sim.output,y="mse",x="colbil",fill="model")+scale_y_continuous(trans='log10')+ scale_x_discrete(breaks=c("0","1"),labels=c("High", "Low")) +xlab("polytomy levels")+ylab("values") +ggtitle("mse") + theme(plot.title = element_text(hjust = 0.5))+ ylim(2,3)
## Scale for 'y' is already present. Adding another scale for 'y', which will
## replace the existing scale.
p1.mse.colbil
## Warning: Removed 136214 rows containing non-finite values (stat_boxplot).
p1.mse.coltype<-ggboxplot(full.sim.output,y="mse",x="coltype",fill="model")+scale_y_continuous(trans='log10') +xlab("polytomy types")+ylab("values") +ggtitle("mse") + theme(plot.title = element_text(hjust = 0.5))+ ylim(2,3)
## Scale for 'y' is already present. Adding another scale for 'y', which will
## replace the existing scale.
p1.mse.coltype
## Warning: Removed 136214 rows containing non-finite values (stat_boxplot).
grid.arrange(p1.mse.coltype,p1.mse.colbil,ncol=2)
## Warning: Removed 136214 rows containing non-finite values (stat_boxplot).
## Warning: Removed 136214 rows containing non-finite values (stat_boxplot).
