#############################################
# Step1: screen reasonable AICs (exclude -2e+10)
# Record those reasonable AICs () within each models
# output 1: report AIC across all dataset in the aspect of model adequacy. (BM support obs, does PMM also support obs)
# output 2: check consistency between AIC(comprare two AICs) and LRT(p-value)
#############################################
rm(list = ls())
#load("~/Dropbox/MingHanHsu/fulldataset/DataSets/results/empR1.rda")
load("~/Dropbox/FCU/Teaching/Mentoring/2016Spring-2021Fall_AssociateProfessor/2020Fall/MingHanHsu/fulldataset/DataSets/results/empR1.rda")
# 5 models
# two AICs in each model
# 14 trees
# few dataset within each tree
AICc<-function(aic=aic,n=n,k=k){
return( aic+2*k*(k+1)/(n-k-1))
}
AICc(aic=237.2,n=30,k=6)
## [1] 240.8522
BM1.7$AICc.obs
## [1] -168.4176
BM1.7$AICc.constrained
## [1] -2e+10
BM11.1$AICc.obs
## [1] 43.26753
BM11.1$AICc.constrained
## [1] 83.28728
# ifelse(BM2.1$Prob<0.05,1,0)
#
# log10(BM2.1$Prob)
# log10(EB2.1$Prob)
# log10(PMM2.1$Prob)
# log10(OU2.1$Prob)
# log10(ID2.1$Prob)
### Assessing the consistency on estimating the rates in multiple trait space.
### H_0: rate are the same $sigma_i$ are equals
### H_a: not all $sigma_i$ are equals
## by LRT
## across all datasets ---
modelset<-c("BM","PMM","OU","EB","ID")
aicout<-array(NA,c(5,14,100))
aic.value.all<-array(NA,c(2,5,14,100))
for(mIndex in 1:length(modelset)){
#mIndex<-1
for (dIndex in 1:length(matrixsetfile)){
#dIndex<-11
for (tIndex in 1:100){
# tIndex<-1
# print(paste(modelset[mIndex],dIndex,".",tIndex,sep=""))
objtoload<-paste(modelset[mIndex],dIndex,".",tIndex,sep="")
if(exists(objtoload)){
getobjtoload<-get(objtoload)
if( getobjtoload$AICc.obs>-2e+7 && getobjtoload$AICc.constrained>-2e+7){
aic.value.all[1,mIndex,dIndex,tIndex]<-getobjtoload$AICc.obs
aic.value.all[2,mIndex,dIndex,tIndex]<-getobjtoload$AICc.constrained
aicout[mIndex,dIndex,tIndex]<-ifelse(getobjtoload$AICc.constrained< getobjtoload$AICc.obs,1,0)
}
}
}
}
}
sumtable<-NULL
for(i in 1:5){
print(table(aicout[i,,]))
sumtable<-rbind(sumtable,table(aicout[i,,]))
}
##
## 0 1
## 153 3
##
## 0 1
## 12 114
##
## 0 1
## 72 130
##
## 0 1
## 77 63
##
## 0 1
## 104 4
sumtable<-t(sumtable)
colnames(sumtable)<-modelset
rownames(sumtable)<-c("obs","con")
sumtable
## BM PMM OU EB ID
## obs 153 12 72 77 104
## con 3 114 130 63 4
library(xtable)
xtable(sumtable)
## % latex table generated in R 4.2.3 by xtable 1.8-4 package
## % Sun Jun 25 09:02:41 2023
## \begin{table}[ht]
## \centering
## \begin{tabular}{rrrrrr}
## \hline
## & BM & PMM & OU & EB & ID \\
## \hline
## obs & 153 & 12 & 72 & 77 & 104 \\
## con & 3 & 114 & 130 & 63 & 4 \\
## \hline
## \end{tabular}
## \end{table}
table(aicout[2,,])
##
## 0 1
## 12 114
summary(aic.value.all)
## Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
## -285670.11 50.35 72.35 -1035.00 113.15 3575.20 12536
BMaic.value.all.obs<-aic.value.all[1,1,,]# observation
BMaic.value.all.con<-aic.value.all[2,1,,]# constraint
PMMaic.value.all.obs<-aic.value.all[1,2,,]# observation
PMMaic.value.all.con<-aic.value.all[2,2,,]# constraint
OUaic.value.all.obs<-aic.value.all[1,3,,]# observation
OUaic.value.all.con<-aic.value.all[2,3,,]# constraint
EBaic.value.all.obs<-aic.value.all[1,4,,]# observation
EBaic.value.all.con<-aic.value.all[2,4,,]# constraint
IDaic.value.all.obs<-aic.value.all[1,5,,]# observation
IDaic.value.all.con<-aic.value.all[2,5,,]# constraint
#BMPMM
aic.obs.bmpmm<-0
aic.obs.pmmbm<-0
aic.con.bmpmm<-0
aic.con.pmmbm<-0
for (i in 1:14){
for(j in 1:100){
if(!is.na(BMaic.value.all.obs[i,j]) & !is.na(PMMaic.value.all.obs[i,j]) ){
if(BMaic.value.all.obs[i,j] < PMMaic.value.all.obs[i,j]){
aic.obs.bmpmm=aic.obs.bmpmm+1}else{aic.obs.pmmbm=aic.obs.pmmbm+1}
if(BMaic.value.all.con[i,j] < PMMaic.value.all.con[i,j]){
aic.con.bmpmm=aic.con.bmpmm+1}else{aic.con.pmmbm=aic.con.pmmbm+1}
}
}
}
aic.obs.bmpmm
## [1] 41
aic.obs.pmmbm
## [1] 15
aic.obs.bmpmm/(aic.obs.bmpmm+aic.obs.pmmbm)
## [1] 0.7321429
aic.con.bmpmm
## [1] 3
aic.con.pmmbm
## [1] 53
aic.con.bmpmm/(aic.con.bmpmm+aic.con.pmmbm)
## [1] 0.05357143
#BMEB
aic.obs.bmeb<-0
aic.obs.ebbm<-0
aic.con.bmeb<-0
aic.con.ebbm<-0
for (i in 1:14){
for(j in 1:100){
if(!is.na(BMaic.value.all.obs[i,j]) & !is.na(EBaic.value.all.obs[i,j]) ){
if(BMaic.value.all.obs[i,j] < EBaic.value.all.obs[i,j]){
aic.obs.bmeb=aic.obs.bmeb+1}else{aic.obs.ebbm=aic.obs.ebbm+1}
if(BMaic.value.all.con[i,j] < EBaic.value.all.con[i,j]){
aic.con.bmeb=aic.con.bmeb+1}else{aic.con.ebbm=aic.con.ebbm+1}
}
}
}
aic.obs.bmeb
## [1] 68
aic.obs.ebbm
## [1] 16
aic.obs.bmeb/(aic.obs.bmeb+aic.obs.ebbm)
## [1] 0.8095238
aic.con.bmeb
## [1] 69
aic.con.ebbm
## [1] 15
aic.con.bmeb/(aic.con.bmeb+aic.con.ebbm)
## [1] 0.8214286
#BMOU
aic.obs.bmou<-0
aic.obs.oubm<-0
aic.con.bmou<-0
aic.con.oubm<-0
for (i in 1:14){
for(j in 1:100){
if(!is.na(BMaic.value.all.obs[i,j]) & !is.na(OUaic.value.all.obs[i,j]) ){
if(BMaic.value.all.obs[i,j] < OUaic.value.all.obs[i,j]){
aic.obs.bmou=aic.obs.bmou+1}else{aic.obs.oubm=aic.obs.oubm+1}
if(BMaic.value.all.con[i,j] < OUaic.value.all.con[i,j]){
aic.con.bmou=aic.con.bmou+1}else{aic.con.oubm=aic.con.oubm+1}
}
}
}
aic.obs.bmou
## [1] 64
aic.obs.oubm
## [1] 38
aic.obs.bmou/(aic.obs.bmou+aic.obs.oubm)
## [1] 0.627451
aic.con.bmou
## [1] 13
aic.con.oubm
## [1] 89
aic.con.bmou/(aic.con.bmou+aic.con.oubm)
## [1] 0.127451
#BMID
aic.obs.bmid<-0
aic.obs.idbm<-0
aic.con.bmid<-0
aic.con.idbm<-0
for (i in 1:14){
for(j in 1:100){
if(!is.na(BMaic.value.all.obs[i,j]) & !is.na(IDaic.value.all.obs[i,j]) ){
if(BMaic.value.all.obs[i,j] < IDaic.value.all.obs[i,j]){
aic.obs.bmid=aic.obs.bmid+1}else{aic.obs.idbm=aic.obs.idbm+1}
if(BMaic.value.all.con[i,j] < IDaic.value.all.con[i,j]){
aic.con.bmid=aic.con.bmid+1}else{aic.con.idbm=aic.con.idbm+1}
}
}
}
aic.obs.bmid
## [1] 34
aic.obs.idbm
## [1] 67
aic.obs.bmid/(aic.obs.bmid+aic.obs.idbm)
## [1] 0.3366337
aic.con.bmid
## [1] 79
aic.con.idbm
## [1] 22
aic.con.bmid/(aic.con.bmid+aic.con.idbm)
## [1] 0.7821782
#PMMEB
aic.obs.pmmeb<-0
aic.obs.ebpmm<-0
aic.con.pmmeb<-0
aic.con.ebpmm<-0
for (i in 1:14){
for(j in 1:100){
if(!is.na(PMMaic.value.all.obs[i,j]) & !is.na(EBaic.value.all.obs[i,j]) ){
if(PMMaic.value.all.obs[i,j] < EBaic.value.all.obs[i,j]){
aic.obs.pmmeb=aic.obs.pmmeb+1}else{aic.obs.ebpmm=aic.obs.ebpmm+1}
if(PMMaic.value.all.con[i,j] < EBaic.value.all.con[i,j]){
aic.con.pmmeb=aic.con.pmmeb+1}else{aic.con.ebpmm=aic.con.ebpmm+1}
}
}
}
aic.obs.pmmeb
## [1] 37
aic.obs.ebpmm
## [1] 20
aic.obs.pmmeb/(aic.obs.pmmeb+aic.obs.ebpmm)
## [1] 0.6491228
aic.con.pmmeb
## [1] 43
aic.con.ebpmm
## [1] 14
aic.con.pmmeb/(aic.con.pmmeb+aic.con.ebpmm)
## [1] 0.754386
#PMMOU
aic.obs.pmmou<-0
aic.obs.oupmm<-0
aic.con.pmmou<-0
aic.con.oupmm<-0
for (i in 1:14){
for(j in 1:100){
if(!is.na(PMMaic.value.all.obs[i,j]) & !is.na(OUaic.value.all.obs[i,j]) ){
if(PMMaic.value.all.obs[i,j] < OUaic.value.all.obs[i,j]){
aic.obs.pmmou=aic.obs.pmmou+1}else{aic.obs.oupmm=aic.obs.oupmm+1}
if(PMMaic.value.all.con[i,j] < OUaic.value.all.con[i,j]){
aic.con.pmmou=aic.con.pmmou+1}else{aic.con.oupmm=aic.con.oupmm+1}
}
}
}
aic.obs.pmmou
## [1] 61
aic.obs.oupmm
## [1] 43
aic.obs.pmmou/(aic.obs.pmmou+aic.obs.oupmm)
## [1] 0.5865385
aic.con.pmmou
## [1] 69
aic.con.oupmm
## [1] 35
aic.con.pmmou/(aic.con.pmmou+aic.con.oupmm)
## [1] 0.6634615
#PMMID
aic.obs.pmmid<-0
aic.obs.idpmm<-0
aic.con.pmmid<-0
aic.con.idpmm<-0
for (i in 1:14){
for(j in 1:100){
if(!is.na(PMMaic.value.all.obs[i,j]) & !is.na(IDaic.value.all.obs[i,j]) ){
if(PMMaic.value.all.obs[i,j] < IDaic.value.all.obs[i,j]){
aic.obs.pmmid=aic.obs.pmmid+1}else{aic.obs.idpmm=aic.obs.idpmm+1}
if(PMMaic.value.all.con[i,j] < IDaic.value.all.con[i,j]){
aic.con.pmmid=aic.con.pmmid+1}else{aic.con.idpmm=aic.con.idpmm+1}
}
}
}
aic.obs.pmmid
## [1] 8
aic.obs.idpmm
## [1] 34
aic.obs.pmmid/(aic.obs.pmmid+aic.obs.idpmm)
## [1] 0.1904762
aic.con.pmmid
## [1] 40
aic.con.idpmm
## [1] 2
aic.con.pmmid/(aic.con.pmmid+aic.con.idpmm)
## [1] 0.952381
#EBOU
aic.obs.ebou<-0
aic.obs.oueb<-0
aic.con.ebou<-0
aic.con.oueb<-0
for (i in 1:14){
for(j in 1:100){
if(!is.na(EBaic.value.all.obs[i,j]) & !is.na(OUaic.value.all.obs[i,j]) ){
if(EBaic.value.all.obs[i,j] < OUaic.value.all.obs[i,j]){
aic.obs.ebou=aic.obs.ebou+1}else{aic.obs.oueb=aic.obs.oueb+1}
if(EBaic.value.all.con[i,j] < OUaic.value.all.con[i,j]){
aic.con.ebou=aic.con.ebou+1}else{aic.con.oueb=aic.con.oueb+1}
}
}
}
aic.obs.ebou
## [1] 30
aic.obs.oueb
## [1] 61
aic.obs.ebou/(aic.obs.ebou+aic.obs.oueb)
## [1] 0.3296703
aic.con.ebou
## [1] 30
aic.con.oueb
## [1] 61
aic.con.ebou/(aic.con.ebou+aic.con.oueb)
## [1] 0.3296703
#EBID
aic.obs.ebid<-0
aic.obs.ideb<-0
aic.con.ebid<-0
aic.con.ideb<-0
for (i in 1:14){
for(j in 1:100){
if(!is.na(EBaic.value.all.obs[i,j]) & !is.na(IDaic.value.all.obs[i,j]) ){
if(EBaic.value.all.obs[i,j] < IDaic.value.all.obs[i,j]){
aic.obs.ebid=aic.obs.ebid+1}else{aic.obs.ideb=aic.obs.ideb+1}
if(EBaic.value.all.con[i,j] < IDaic.value.all.con[i,j]){
aic.con.ebid=aic.con.ebid+1}else{aic.con.ideb=aic.con.ideb+1}
}
}
}
aic.obs.ebid
## [1] 13
aic.obs.ideb
## [1] 50
aic.obs.ebid/(aic.obs.ebid+aic.obs.ideb)
## [1] 0.2063492
aic.con.ebid
## [1] 31
aic.con.ideb
## [1] 32
aic.con.ebid/(aic.con.ebid+aic.con.ideb)
## [1] 0.4920635
#OUID
aic.obs.ouid<-0
aic.obs.idou<-0
aic.con.ouid<-0
aic.con.idou<-0
for (i in 1:14){
for(j in 1:100){
if(!is.na(OUaic.value.all.obs[i,j]) & !is.na(IDaic.value.all.obs[i,j]) ){
if(OUaic.value.all.obs[i,j] < IDaic.value.all.obs[i,j]){
aic.obs.ouid=aic.obs.ouid+1}else{aic.obs.idou=aic.obs.idou+1}
if(OUaic.value.all.con[i,j] < IDaic.value.all.con[i,j]){
aic.con.ouid=aic.con.ouid+1}else{aic.con.idou=aic.con.idou+1}
}
}
}
aic.obs.ouid
## [1] 30
aic.obs.idou
## [1] 46
aic.obs.ouid/(aic.obs.ouid+aic.obs.idou)
## [1] 0.3947368
aic.con.ouid
## [1] 69
aic.con.idou
## [1] 7
aic.con.ouid/(aic.con.ouid+aic.con.idou)
## [1] 0.9078947
############
models<-c("bm", "pmm", "ou", "eb", "id" )
outcount<-array(NA,c(length(models),length(models)))
rownames(outcount)<-colnames(outcount)<-models
for(i in 1:length(models)){
for(j in 1:length(models)){
if(i!=j){
outcount[i,j]<-paste("(", get(paste("aic.con.",models[i],models[j],sep="")),",",
get(paste("aic.con.",models[j],models[i],sep="")),")",sep="")
}
}
}
outcount
## bm pmm ou eb id
## bm NA "(3,53)" "(13,89)" "(69,15)" "(79,22)"
## pmm "(53,3)" NA "(69,35)" "(43,14)" "(40,2)"
## ou "(89,13)" "(35,69)" NA "(61,30)" "(69,7)"
## eb "(15,69)" "(14,43)" "(30,61)" NA "(31,32)"
## id "(22,79)" "(2,40)" "(7,69)" "(32,31)" NA
library(xtable)
xtable(outcount,digits=3)
## % latex table generated in R 4.2.3 by xtable 1.8-4 package
## % Sun Jun 25 09:02:41 2023
## \begin{table}[ht]
## \centering
## \begin{tabular}{rlllll}
## \hline
## & bm & pmm & ou & eb & id \\
## \hline
## bm & & (3,53) & (13,89) & (69,15) & (79,22) \\
## pmm & (53,3) & & (69,35) & (43,14) & (40,2) \\
## ou & (89,13) & (35,69) & & (61,30) & (69,7) \\
## eb & (15,69) & (14,43) & (30,61) & & (31,32) \\
## id & (22,79) & (2,40) & (7,69) & (32,31) & \\
## \hline
## \end{tabular}
## \end{table}
namesVar<-array(NA,c(10))
ratiotoPlot<-array(NA,c(10))
k=0
for(i in 1:4){
for(j in (i+1):5){
k=k+1
ratiotoPlot[k]<-get(paste("aic.con.",models[i],models[j],sep=""))/(get(paste("aic.con.",models[i],models[j],sep=""))+get(paste("aic.con.",models[j],models[i],sep="")))
namesVar[k]<-paste(modelset[i],"-",modelset[j],sep="")
# print(outcount[i,j])
}
}
names(ratiotoPlot)<-namesVar
barplot(ratiotoPlot)
library(ggplot2)
df=data.frame(namesVar=namesVar,ratiotoPlot=ratiotoPlot)
ggplot(df,aes(x=namesVar,y=ratiotoPlot,fill=ratiotoPlot))+geom_bar(stat="identity")+ geom_hline(yintercept=0.5, linetype="dashed", color = "red") + xlab("Pair of Models")+ ylab("Ratio")+ggtitle("Compared by AICc: Homogenerous Rate")+
geom_text(aes(label=round(ratiotoPlot,2)), position=position_dodge(width=0.9), vjust=-0.25)
# Reorder factor levels by ratiotoPlot
df$namesVar = with(df, reorder(namesVar, ratiotoPlot))
# Plot
ggplot(df, aes(x = namesVar, y = ratiotoPlot, fill = ratiotoPlot)) +
geom_bar(stat = "identity") +
geom_hline(yintercept = 0.5, linetype = "dashed", color = "red") +
xlab("Pair of Models") +
ylab("Ratio") +
ggtitle("Compared by AICc: Homogeneous Rate") +
geom_text(aes(label = round(ratiotoPlot, 2)), position = position_dodge(width = 0.9), vjust = -0.25)
Homogeneous Rate
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
# Filter out rows with 'ID' in the namesVar
df <- df %>%
filter(!grepl('ID', namesVar))
# Reorder factor levels by ratiotoPlot
df$namesVar <- with(df, reorder(namesVar, ratiotoPlot))
# Plot
ggplot(df, aes(x = namesVar, y = ratiotoPlot, fill = ratiotoPlot)) +
geom_bar(stat = "identity") +
geom_hline(yintercept = 0.5, linetype = "dashed", color = "red") +
xlab("Pair of Models") +
ylab("Ratio") +
ggtitle("Compared by AICc: Homogeneous Rate") +
theme(plot.title = element_text(hjust = 0.5)) + # this centers the title
geom_text(aes(label = round(ratiotoPlot, 2)), position = position_dodge(width = 0.9), vjust = -0.25) + theme_classic()+
theme(legend.position="none") +
theme(plot.title = element_text(hjust = 0.5))
Heterogeneous Rate
models<-c("bm", "pmm", "ou", "eb", "id" )
outcount<-array(NA,c(length(models),length(models)))
rownames(outcount)<-colnames(outcount)<-models
for(i in 1:length(models)){
for(j in 1:length(models)){
if(i!=j){
outcount[i,j]<-paste("(", get(paste("aic.obs.",models[i],models[j],sep="")),",",
get(paste("aic.obs.",models[j],models[i],sep="")),")",sep="")
}
}
}
outcount
## bm pmm ou eb id
## bm NA "(41,15)" "(64,38)" "(68,16)" "(34,67)"
## pmm "(15,41)" NA "(61,43)" "(37,20)" "(8,34)"
## ou "(38,64)" "(43,61)" NA "(61,30)" "(30,46)"
## eb "(16,68)" "(20,37)" "(30,61)" NA "(13,50)"
## id "(67,34)" "(34,8)" "(46,30)" "(50,13)" NA
library(xtable)
xtable(outcount,digits=3)
## % latex table generated in R 4.2.3 by xtable 1.8-4 package
## % Sun Jun 25 09:02:43 2023
## \begin{table}[ht]
## \centering
## \begin{tabular}{rlllll}
## \hline
## & bm & pmm & ou & eb & id \\
## \hline
## bm & & (41,15) & (64,38) & (68,16) & (34,67) \\
## pmm & (15,41) & & (61,43) & (37,20) & (8,34) \\
## ou & (38,64) & (43,61) & & (61,30) & (30,46) \\
## eb & (16,68) & (20,37) & (30,61) & & (13,50) \\
## id & (67,34) & (34,8) & (46,30) & (50,13) & \\
## \hline
## \end{tabular}
## \end{table}
namesVar<-array(NA,c(10))
ratiotoPlot<-array(NA,c(10))
k=0
for(i in 1:4){
for(j in (i+1):5){
k=k+1
ratiotoPlot[k]<-get(paste("aic.obs.",models[i],models[j],sep=""))/(get(paste("aic.obs.",models[i],models[j],sep=""))+get(paste("aic.obs.",models[j],models[i],sep="")))
namesVar[k]<-paste(modelset[i],"-",modelset[j],sep="")
# print(outcount[i,j])
}
}
names(ratiotoPlot)<-namesVar
barplot(ratiotoPlot)
library(ggplot2)
df=data.frame(namesVar=namesVar,ratiotoPlot=ratiotoPlot)
ggplot(df,aes(x=namesVar,y=ratiotoPlot,fill=ratiotoPlot))+geom_bar(stat="identity")+ geom_hline(yintercept=0.5, linetype="dashed", color = "red") + xlab("Pair of Models")+ ylab("Ratio")+ggtitle("Compared by AICc: Heterogeneous Rate")+
geom_text(aes(label=round(ratiotoPlot,2)), position=position_dodge(width=0.9), vjust=-0.25)
# Reorder factor levels by ratiotoPlot
df$namesVar = with(df, reorder(namesVar, ratiotoPlot))
# Plot
ggplot(df, aes(x = namesVar, y = ratiotoPlot, fill = ratiotoPlot)) +
geom_bar(stat = "identity") +
geom_hline(yintercept = 0.5, linetype = "dashed", color = "red") +
xlab("Pair of Models") +
ylab("Ratio") +
ggtitle("Compared by AICc: Heterogeneous Rate") +
geom_text(aes(label = round(ratiotoPlot, 2)), position = position_dodge(width = 0.9), vjust = -0.25)
# Filter out rows with 'ID' in the namesVar
df <- df %>%
filter(!grepl('ID', namesVar))
# Reorder factor levels by ratiotoPlot
df$namesVar <- with(df, reorder(namesVar, ratiotoPlot))
# Plot
ggplot(df, aes(x = namesVar, y = ratiotoPlot, fill = ratiotoPlot)) +
geom_bar(stat = "identity") +
geom_hline(yintercept = 0.5, linetype = "dashed", color = "red") +
xlab("Pair of Models") +
ylab("Ratio") +
ggtitle("Compared by AICc: Heterogeneous Rate") +
theme(plot.title = element_text(hjust = 0.5)) + # this centers the title
geom_text(aes(label = round(ratiotoPlot, 2)), position = position_dodge(width = 0.9), vjust = -0.25) + theme_classic()+
theme(legend.position="none") +
theme(plot.title = element_text(hjust = 0.5))
