---
title: "Comparison of pairwise models under the LRT test"
output: html_document
---

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


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


BM11.1$Lobs
BM11.1$Lconstrained

OU11.1$Lobs.ou
OU11.1$Lconstrained.ou






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

BM11.1$Prob

# Step1: screen reasonable AICs  (exclude -2e+10)
# Record those reasonable AICs () within each models

# output 1: check consistency between AIC(comprare two AICs) and LRT(p-value) 

# output 2: report AIC across all dataset in the aspect of model adequacy.  (BM support obs, does PMM also support obs)





# 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")

A=array(NA,c(2,3))
A[2,1]<-1
A[2,2]<-0
A

p.value.all<-pout<-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)
        pout[mIndex,dIndex,tIndex]<-ifelse(getobjtoload$Prob<0.001,1,0)
        p.value.all[mIndex,dIndex,tIndex]<-getobjtoload$Prob
        
        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)
        }
      }
    }
  }
}

table(pout[1,,])
table(aicout[1,,])


table(pout[2,,])
table(aicout[2,,])

summary(aic.value.all)





#modelset<-c("BM","PMM","OU","EB","ID")
BMpout<-pout[1,,]# 
PMMpout<-pout[2,,]
OUpout<-pout[3,,]
EBpout<-pout[4,,]
IDpout<-pout[5,,]


BMp.value.all<-p.value.all[1,,]# 
PMMp.value.all<-p.value.all[2,,]
OUp.value.all<-p.value.all[3,,]
EBp.value.all<-p.value.all[4,,]
IDp.value.all<-p.value.all[5,,]


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


# Assessing the performance of models in the aspect of the violation to the null hypothesis.

BMpout[1:5,1:10]
PMMpout[1:5,1:10]
OUpout[1:5,1:10]
EBpout[1:5,1:10]
IDpout[1:5,1:10]


#BMPMM
agree=0
notagree=0

sensbmpmm<-0
senspmmbm<-0

aicobsbmpmm<-0
aicobspmmbm<-0

aicconbmpmm<-0
aicconpmmbm<-0

for (i in 1:14){
  for(j in 1:100){
    if(!is.na(BMpout[i,j]) & !is.na(PMMpout[i,j]) ){
      if(BMpout[i,j]==PMMpout[i,j]){agree=agree+1}
      if(BMpout[i,j]!=PMMpout[i,j]){notagree=notagree+1}
      
      if(BMp.value.all[i,j]<PMMp.value.all[i,j]){sensbmpmm=sensbmpmm+1}else{senspmmbm=senspmmbm+1}
    }
    
    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]){
        aicobsbmpmm=aicobsbmpmm+1}else{aicobspmmbm=aicobspmmbm+1}
      
      if(BMaic.value.all.con[i,j] < PMMaic.value.all.con[i,j]){
        aicconbmpmm=aicconbmpmm+1}else{aicconpmmbm=aicconpmmbm+1}
    }
  }
}
print(agree)
print(notagree)

agreebmpmm<-agree
notagreebmpmm<-notagree
sensbmpmm
senspmmbm

sensbmpmm/(sensbmpmm+senspmmbm)
#  0.7236842,  bm is more violate to the null assumption


aicobsbmpmm

aicobspmmbm

aicobsbmpmm/(aicobsbmpmm+aicobspmmbm)

aicconbmpmm/(aicconbmpmm+aicconpmmbm)






#BMEB
agree=0
notagree=0

sensbmeb<-0
sensebbm<-0

for (i in 1:14){
  for(j in 1:100){
    if(!is.na(BMpout[i,j]) & !is.na(EBpout[i,j]) ){
      if(BMpout[i,j]==EBpout[i,j]){agree=agree+1}
      if(BMpout[i,j]!=EBpout[i,j]){notagree=notagree+1}
      
      if(BMp.value.all[i,j]<EBp.value.all[i,j]){sensbmeb=sensbmeb+1}else{sensebbm=sensebbm+1}
      
    }
  }
}
print(agree)
print(notagree)

agreebmeb<-agree
notagreebmeb<-notagree
sensbmeb
sensebbm

sensbmeb/(sensbmeb+sensebbm)
#   0.6020408,  bm is more violate to the null assumption


#BMOU
agree=0
notagree=0

sensbmou<-0
sensoubm<-0

for (i in 1:14){
  for(j in 1:100){
    if(!is.na(BMpout[i,j]) & !is.na(OUpout[i,j]) ){
      if(BMpout[i,j]==OUpout[i,j]){agree=agree+1}
      if(BMpout[i,j]!=OUpout[i,j]){notagree=notagree+1}
      
      if(BMp.value.all[i,j]<OUp.value.all[i,j]){sensbmou=sensbmou+1}else{sensoubm=sensoubm+1}
      
    }
  }
}
print(agree)
print(notagree)

agreebmou<-agree
notagreebmou<-notagree
sensbmou
sensoubm

sensbmou/(sensoubm+sensbmou)

#   0.6119403,  bm is more violate to the null assumption


#BMID
agree=0
notagree=0

sensbmid<-0
sensidbm<-0

for (i in 1:14){
  for(j in 1:100){
    if(!is.na(BMpout[i,j]) & !is.na(IDpout[i,j]) ){
      if(BMpout[i,j]==IDpout[i,j]){agree=agree+1}
      if(BMpout[i,j]!=IDpout[i,j]){notagree=notagree+1}
      
      if(BMp.value.all[i,j]<IDp.value.all[i,j]){sensbmid=sensbmid+1}else{sensidbm=sensidbm+1}
      
    }
  }
}
print(agree)
print(notagree)

agreebmid<-agree
notagreebmid<-notagree
sensbmid
sensidbm

sensbmid/(sensbmid+sensidbm)

#   0.2318841,  bm is more violate to the null assumption


#PMMEB
agree=0
notagree=0

senspmmeb<-0
sensebpmm<-0

for (i in 1:14){
  for(j in 1:100){
    if(!is.na(PMMpout[i,j]) & !is.na(EBpout[i,j]) ){
      if(PMMpout[i,j]==EBpout[i,j]){agree=agree+1}
      if(PMMpout[i,j]!=EBpout[i,j]){notagree=notagree+1}
      
      if(PMMp.value.all[i,j]<EBp.value.all[i,j]){senspmmeb=senspmmeb+1}else{sensebpmm=sensebpmm+1}
      
    }
  }
}
print(agree)
print(notagree)

agreepmmeb<-agree
notagreepmmeb<-notagree
senspmmeb
sensebpmm

senspmmeb/(senspmmeb+sensebpmm)

#PMMOU
agree=0
notagree=0

senspmmou<-0
sensoupmm<-0

for (i in 1:14){
  for(j in 1:100){
    if(!is.na(PMMpout[i,j]) & !is.na(OUpout[i,j]) ){
      if(PMMpout[i,j]==OUpout[i,j]){agree=agree+1}
      if(PMMpout[i,j]!=OUpout[i,j]){notagree=notagree+1}
      
      if(PMMp.value.all[i,j]<OUp.value.all[i,j]){senspmmou=senspmmou+1}else{sensoupmm=sensoupmm+1}
      
    }
  }
}
print(agree)
print(notagree)

agreepmmou<-agree
notagreepmmou<-notagree
senspmmou
sensoupmm

senspmmou/(sensoupmm+senspmmou)

#PMMID
agree=0
notagree=0

senspmmid<-0
sensidpmm<-0

for (i in 1:14){
  for(j in 1:100){
    if(!is.na(PMMpout[i,j]) & !is.na(IDpout[i,j]) ){
      if(PMMpout[i,j]==IDpout[i,j]){agree=agree+1}
      if(PMMpout[i,j]!=IDpout[i,j]){notagree=notagree+1}
      
      if(PMMp.value.all[i,j]<IDp.value.all[i,j]){senspmmid=senspmmid+1}else{sensidpmm=sensidpmm+1}
      
    }
  }
}
print(agree)
print(notagree)

agreepmmid<-agree
notagreepmmid<-notagree
senspmmid
sensidpmm

senspmmid/(senspmmid+sensidpmm)

#EBOU
agree=0
notagree=0

sensebou<-0
sensoueb<-0

for (i in 1:14){
  for(j in 1:100){
    if(!is.na(EBpout[i,j]) & !is.na(OUpout[i,j]) ){
      if(EBpout[i,j]==OUpout[i,j]){agree=agree+1}
      if(EBpout[i,j]!=OUpout[i,j]){notagree=notagree+1}
      
      if(EBp.value.all[i,j]<OUp.value.all[i,j]){sensebou=sensebou+1}else{sensoueb=sensoueb+1}
      
    }
  }
}
print(agree)
print(notagree)

agreeebou<-agree
notagreeebou<-notagree
sensebou
sensoueb


sensebou/(sensoueb+sensebou)

#EBID
agree=0
notagree=0

sensebid<-0
sensideb<-0

for (i in 1:14){
  for(j in 1:100){
    if(!is.na(EBpout[i,j]) & !is.na(IDpout[i,j]) ){
      if(EBpout[i,j]==IDpout[i,j]){agree=agree+1}
      if(EBpout[i,j]!=IDpout[i,j]){notagree=notagree+1}
      
      if(EBp.value.all[i,j]<IDp.value.all[i,j]){sensebid=sensebid+1}else{sensideb=sensideb+1}
      
    }
  }
}
print(agree)
print(notagree)

agreeebid<-agree
notagreeebid<-notagree
sensebid
sensideb

sensebid/(sensebid+sensideb)


#OUID
agree=0
notagree=0

sensouid<-0
sensidou<-0

for (i in 1:14){
  for(j in 1:100){
    if(!is.na(OUpout[i,j]) & !is.na(IDpout[i,j]) ){
      if(OUpout[i,j]==IDpout[i,j]){agree=agree+1}
      if(OUpout[i,j]!=IDpout[i,j]){notagree=notagree+1}
      
      if(OUp.value.all[i,j]<IDp.value.all[i,j]){sensouid=sensouid+1}else{sensidou=sensidou+1}
      
    }
  }
}
print(agree)
print(notagree)

agreeouid<-agree
notagreeouid<-notagree
sensouid
sensidou

sensouid/(sensouid+sensidou)




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

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("sens",models[i],models[j],sep="")),",",
                           get(paste("sens",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("sens",models[i],models[j],sep=""))/(get(paste("sens",models[i],models[j],sep=""))+get(paste("sens",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 LRT")+
  geom_text(aes(label=round(ratiotoPlot,2)), position=position_dodge(width=0.9), vjust=-0.25)


# Ensure dplyr is loaded
library(dplyr)

# Ensure dplyr and ggplot2 are loaded
library(dplyr)
library(ggplot2)

# Arrange the data.frame by 'ratiotoPlot' in descending order
df <- df %>% arrange(desc(ratiotoPlot))

# Reorder namesVar according to the order of ratiotoPlot
df$namesVar <- factor(df$namesVar, levels = df$namesVar)

# 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 LRT") +
  theme(plot.title = element_text(hjust = 0.5), # this centers the title
        legend.position = "none", # this removes the legend
        panel.background = element_rect(fill = "white")) + # this makes the background white
  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)) %>%
  arrange(desc(ratiotoPlot))

# Reorder namesVar according to the order of ratiotoPlot
df$namesVar <- factor(df$namesVar, levels = df$namesVar)

# 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 LRT") +
  theme(plot.title = element_text(hjust = 0.5), # this centers the title
        legend.position = "none", # this removes the legend
        panel.background = element_rect(fill = "white")) + # this makes the background white
  geom_text(aes(label=round(ratiotoPlot,2)), position=position_dodge(width=0.9), vjust=-0.25)


# A smaller p-value indicate stronger violation to the null hypothesis for homogeneous rate of evolution across traits 
# here we compare two models on testing the homogeneity of rates via their sensitivity to the violation to the null. 
# Since all models use the same datasets, 
# for each dataset, we say the model A has stronger evidence to violate the null hypothesis than the model B if the p-value of the likelihood ratio test under model A is smaller than the the p-value of the likelihood ratio test under model B.

```
