---
title: "comparison of pairwise models in violation of the null hypothesis using AIC"
output: html_document
---


```{r}
#############################################
# 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)

BM1.7$AICc.obs
BM1.7$AICc.constrained

BM11.1$AICc.obs
BM11.1$AICc.constrained



# 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,,]))
}

sumtable<-t(sumtable)
colnames(sumtable)<-modelset
rownames(sumtable)<-c("obs","con")
sumtable

library(xtable)
xtable(sumtable)


table(aicout[2,,])
summary(aic.value.all)


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
aic.obs.pmmbm

aic.obs.bmpmm/(aic.obs.bmpmm+aic.obs.pmmbm)

aic.con.bmpmm
aic.con.pmmbm

aic.con.bmpmm/(aic.con.bmpmm+aic.con.pmmbm)


#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
aic.obs.ebbm

aic.obs.bmeb/(aic.obs.bmeb+aic.obs.ebbm)

aic.con.bmeb
aic.con.ebbm

aic.con.bmeb/(aic.con.bmeb+aic.con.ebbm)


#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
aic.obs.oubm

aic.obs.bmou/(aic.obs.bmou+aic.obs.oubm)

aic.con.bmou
aic.con.oubm

aic.con.bmou/(aic.con.bmou+aic.con.oubm)


#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
aic.obs.idbm

aic.obs.bmid/(aic.obs.bmid+aic.obs.idbm)

aic.con.bmid
aic.con.idbm

aic.con.bmid/(aic.con.bmid+aic.con.idbm)


#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
aic.obs.ebpmm

aic.obs.pmmeb/(aic.obs.pmmeb+aic.obs.ebpmm)

aic.con.pmmeb
aic.con.ebpmm

aic.con.pmmeb/(aic.con.pmmeb+aic.con.ebpmm)


#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
aic.obs.oupmm

aic.obs.pmmou/(aic.obs.pmmou+aic.obs.oupmm)

aic.con.pmmou
aic.con.oupmm

aic.con.pmmou/(aic.con.pmmou+aic.con.oupmm)


#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
aic.obs.idpmm

aic.obs.pmmid/(aic.obs.pmmid+aic.obs.idpmm)

aic.con.pmmid
aic.con.idpmm

aic.con.pmmid/(aic.con.pmmid+aic.con.idpmm)


#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
aic.obs.oueb

aic.obs.ebou/(aic.obs.ebou+aic.obs.oueb)

aic.con.ebou
aic.con.oueb

aic.con.ebou/(aic.con.ebou+aic.con.oueb)


#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
aic.obs.ideb

aic.obs.ebid/(aic.obs.ebid+aic.obs.ideb)

aic.con.ebid
aic.con.ideb

aic.con.ebid/(aic.con.ebid+aic.con.ideb)


#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
aic.obs.idou

aic.obs.ouid/(aic.obs.ouid+aic.obs.idou)

aic.con.ouid
aic.con.idou

aic.con.ouid/(aic.con.ouid+aic.con.idou)



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

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


library(xtable)
xtable(outcount,digits=3)

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

```{r}
library(dplyr)
# 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

```{r}
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


library(xtable)
xtable(outcount,digits=3)

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))
```
