---
title: "Model Adequacy"
output: html_document
---

```{r,eval=FALSE}
rm(list=ls())
library(phangorn)
#install.packages("MBSP")
library(MBSP) 
#install.packages("matrixNormal")
library(mvtnorm)
library(matrixNormal)
library(parallel)
setwd("~/Dropbox/MingHanHsu/rcodeMHH")
source("ComparingBMv2.r")
source("ComparingIDv2.r")
source("ComparingPMMv2.r")
source("ComparingEBv2.r")
source("ComparingOUv2.r")

cmtxfile<-list.files("~/Dropbox/MingHanHsu/fulldataset/DataSets/Matrices0707/")
#traitfile <-list.files("~/Dropbox/FCU/Teaching/Mentoring/2020Fall/MingHanHsu/fulldataset/DataSets/NewMean")

traitfile <-list.files("~/Dropbox/MingHanHsu/fulldataset/DataSets/NewMean")

sims<-1000
for(cmtxIndex in 1:length(cmtxfile)){
#  cmtxIndex<-6
  setwd("~/Dropbox/MingHanHsu/fulldataset/DataSets/Matrices0707/")
  C<-read.table(cmtxfile[cmtxIndex])  
  C<-t(C)+C
  diag(C)<-1
  phy<-upgma(2*(max(C)-C))
#  plot(phy)
  setwd("~/Dropbox/MingHanHsu/fulldataset/DataSets/NewMean//")
#  for(traitIndex in 1:length(traitfile)){
#    traitIndex<-1
    #   traitset<-read.table("~/Dropbox/FCU/Teaching/Mentoring/2020Fall/MingHanHsu/fulldataset/DataSets/NewMean/EOV-59.txt")
    #  traitset<-read.table("~/Dropbox/MingHanHsu/fulldataset/DataSets/NewMean/EOV-59.txt")
    traitset<-read.table(traitfile[cmtxIndex])
    #traitfile
    possiblebarcomb <- combn(floor(ncol(traitset)/2),3)
    possiblesdcomb <- floor(ncol(traitset)/2) +  combn(floor(ncol(traitset)/2),3)
    
    for(dataIndex in 1:dim(possiblebarcomb)[2]){
#      dataIndex<-1
      traitbar1 <-traitset[,possiblebarcomb[1,dataIndex]]
      traitbar2 <-traitset[,possiblebarcomb[2,dataIndex]]
      traitbar3 <-traitset[,possiblebarcomb[3,dataIndex]]
      
      x<-cbind(traitbar1,traitbar2,traitbar3)
      rownames(x)<-phy$tip.label
      x<-log(x)
      ones<-matrix(rep(1,dim(x)[1]),ncol=1)
##########################################################      
      # BM
      result.bm<-NULL
      try(result.bm<-CompareRates.multTraitBM(phy=phy,x=x,TraitCov=T,ms.err=NULL,ms.cov=NULL))
      if(!is.null(result.bm)){
        C <- vcv(phy)
        R <- result.bm$Robs
        M <- c(t(ones)%*%solve(C)%*%x)/c(t(ones)%*%solve(C)%*%ones)      
        M <- ones%*%M
        sim.likearr.bm <- array(NA,c(sims))  
        result.bm$Lobs
        for(simIndex in 1:sims){ 
#          simIndex<-1
          sim.x<- matrix.normal(M=M, U=C, V=R)
#          matrixNormal::dmatnorm(sim.x, M, C, R)
          sim.likearr.bm[simIndex]<- matrixNormal::dmatnorm(sim.x, M, C, R)
        }
      }#is.null
      
##########################################################
      
      # PMM
      result.pmm<-NULL
      try(result.pmm<-CompareRates.multTraitPMM(phy=phy,x=x,TraitCov=T,ms.err=NULL,ms.cov=NULL))
      if(!is.null(result.pmm)){
        # result.pmm$Lobs.pmm
        # p<-3
        # for(Index in 1:p){
        #   h<-phylosig(tree=phy,x=x[,Index],method="lambda")$lambda
        #   h.array[Index]<-h
        #   Sh<-h^2*C + (1-h^2)*I
        # }        
        Robs.pmm <- result.pmm$Robs.pmm
        V<-result.pmm$RkronSh
        #V<-nearPD(V)
        #V<-as.matrix(V$mat)
        # need the vcv matrix for pmm
        spIndex<-splitIndices(dim(V)[1], 3)
        H1<-V[spIndex[[1]],spIndex[[1]]]
        solveH1 <- solve(H1)
        dim(solveH1)
        m1<-as.numeric(t(ones)%*%solveH1%*%x[,1])/as.numeric(t(ones)%*%solveH1%*%ones)
        m1<-m1*ones
        m1
        H2<-V[spIndex[[2]],spIndex[[2]]]
        solveH2 <- solve(H2)
        m2<-as.numeric(t(ones)%*%solveH2%*%x[,2])/as.numeric(t(ones)%*%solveH2%*%ones)
        m2<-m2*ones
        m2
        H3<-V[spIndex[[3]],spIndex[[3]]]
        solveH3 <- solve(H3)
        m3<-as.numeric(t(ones)%*%solveH3%*%x[,3])/as.numeric(t(ones)%*%solveH3%*%ones)
        m3<-m3*ones
        m3
        mvec<-rbind(m1,m2,m3)
        mvec
        V<-array(0,dim(V))
        V[spIndex[[1]],spIndex[[1]]]<-H1
        V[spIndex[[2]],spIndex[[2]]]<-H2
        V[spIndex[[3]],spIndex[[3]]]<-H3
        V
        sim.likearr.pmm <- array(NA,c(sims))  
        for(simIndex in 1:sims){
          # simIndex <-1
#          result.pmm$Lobs.pmm
          sim.x<-rmvnorm(n=1,mean=mvec,sigma=V)
#          dmvnorm(sim.x, mean=mvec,sigma=V,log = T)
          sim.likearr.pmm[simIndex]<- dmvnorm(sim.x, mean=mvec,sigma=V,log = T)
        }
      }#is.null
      
      
      
      # OU
      result.ou <- NULL
      try(result.ou<-CompareRates.multTraitOU(phy=phy,x=x,TraitCov=T,ms.err=NULL,ms.cov=NULL))
      
      if(!is.null(result.ou)){
        Robs.ou <- result.ou$Robs.ou
        V<-result.ou$RkronSa
        #V<-nearPD(V)
        #V<-as.matrix(V$mat)
        
        # need the vcv matrix for eb
        spIndex<-splitIndices(dim(V)[1], 3)
        O1<-V[spIndex[[1]],spIndex[[1]]]
        solveO1 <- solve(O1)
        dim(solveO1)
        m1<-as.numeric(t(ones)%*%solveO1%*%x[,1])/as.numeric(t(ones)%*%solveO1%*%ones)
        m1<-m1*ones
        m1
        O2<-V[spIndex[[2]],spIndex[[2]]]
        solveO2 <- solve(O2)
        m2<-as.numeric(t(ones)%*%solveO2%*%x[,2])/as.numeric(t(ones)%*%solveO2%*%ones)
        m2<-m2*ones
        m2
        O3<-V[spIndex[[3]],spIndex[[3]]]
        solveO3 <- solve(O3)
        m3<-as.numeric(t(ones)%*%solveO3%*%x[,3])/as.numeric(t(ones)%*%solveO3%*%ones)
        m3<-m3*ones
        m3
        mvec<-rbind(m1,m2,m3)
        mvec
        V <-array(0,dim(V))
        V[spIndex[[1]],spIndex[[1]]]<-O1
        V[spIndex[[2]],spIndex[[2]]]<-O2
        V[spIndex[[3]],spIndex[[3]]]<-O3
        V
        sim.likearr.ou <- array(NA,c(sims))  
        for(simIndex in 1:sims){
          # simIndex <-1
          result.ou$Lobs
          sim.x<-rmvnorm(n=1,mean=mvec,sigma=V)
          dmvnorm(sim.x, mean=mvec,sigma=V,log = T)
          sim.likearr.ou[simIndex]<- dmvnorm(sim.x, mean=mvec,sigma=V,log = T)
        }
      }#is.null
      
      # EB
      result.eb <- NULL
      try(result.eb<-CompareRates.multTraitEB(phy=phy,x=x,TraitCov=T,ms.err=NULL,ms.cov=NULL))
      
      if(!is.null(result.eb)){
        Robs.eb <- result.eb$Robs.eb
        V<-result.eb$RkronSr
        #V<-nearPD(V)
        #V<-as.matrix(V$mat)
        
        # need the vcv matrix for eb
        spIndex<-splitIndices(dim(V)[1], 3)
        E1<-V[spIndex[[1]],spIndex[[1]]]
        solveE1 <- solve(E1)
        dim(solveE1)
        m1<-as.numeric(t(ones)%*%solveE1%*%x[,1])/as.numeric(t(ones)%*%solveE1%*%ones)
        m1<-m1*ones
        m1
        E2<-V[spIndex[[2]],spIndex[[2]]]
        solveE2 <- solve(E2)
        m2<-as.numeric(t(ones)%*%solveE2%*%x[,2])/as.numeric(t(ones)%*%solveE2%*%ones)
        m2<-m2*ones
        m2
        E3<-V[spIndex[[3]],spIndex[[3]]]
        solveE3 <- solve(E3)
        m3<-as.numeric(t(ones)%*%solveE3%*%x[,3])/as.numeric(t(ones)%*%solveE3%*%ones)
        m3<-m3*ones
        m3
        mvec<-rbind(m1,m2,m3)
        mvec
        V <-array(0,dim(V))
        V[spIndex[[1]],spIndex[[1]]]<-E1
        V[spIndex[[2]],spIndex[[2]]]<-E2
        V[spIndex[[3]],spIndex[[3]]]<-E3
        V
        sim.likearr.eb <- array(NA,c(sims))  
        for(simIndex in 1:sims){
          # simIndex <-1
          result.eb$Lobs
          sim.x<-rmvnorm(n=1,mean=mvec,sigma=V)
#          dmvnorm(sim.x, mean=mvec,sigma=V,log = T)
          sim.likearr.eb[simIndex]<- dmvnorm(sim.x, mean=mvec,sigma=V,log = T)
        }
      }#is.null
      
      # ID
      result.id<-NULL
      try(result.id<-CompareRates.multTraitID(phy=phy,x=x,TraitCov=T,ms.err=NULL,ms.cov=NULL))
      
      if(!is.null(result.id)){
        I <- diag(1,dim(x)[1])
        R <- result.id$Robs
        M <- c(t(ones)%*%solve(I)%*%x)/c(t(ones)%*%solve(I)%*%ones)      
        M <- ones%*%M
        sim.likearr.id <- array(NA,c(sims))  
        for(simIndex in 1:sims){ 
#          result.id$Lobs
          sim.x<- matrix.normal(M=M, U=C, V=R)
#          matrixNormal::dmatnorm(sim.x, M, I, R)
          sim.likearr.id[simIndex]<- matrixNormal::dmatnorm(sim.x, M, I, R)
        }
      }#is.null
      
        save( result.bm,    
              result.pmm,    
              result.ou,
              result.eb,
              result.id,
              sim.likearr.bm,
              sim.likearr.pmm,   
              sim.likearr.ou,
              sim.likearr.eb,
              sim.likearr.id,
              file=paste("~/Dropbox/MingHanHsu/simulation/adequacy/tree",cmtxIndex,"dataset",dataIndex,"simIndex",simIndex,".rda",sep=""))
    }#dataIndex
}#matrixIndex



#setwd("~/Dropbox/FCU/Teaching/Mentoring/2020Fall/MingHanHsu/rcodeMHH/")
```
###############
### Summary ###
###############


```{r, warning=FALSE}

setwd("~/Dropbox/FCU/Teaching/Mentoring/2016Spring-2021Fall_AssociateProfessor/2020Fall/MingHanHsu/simulation/adequacy/")
system("ls")
filetoload<-list.files()

result.rawLobs.full<-NULL#array(NA,c(length(filetoload),5))
result.simLobs.full<-NULL
for(fileIndex in 1:length(filetoload)){
  #  fileIndex<-1
  load(filetoload[fileIndex])  
  rawLobs <- c(result.bm$Lobs,result.pmm$Lobs,result.ou$Lobs,result.eb$Lobs,result.id$Lobs)
  result.rawLobs.full<-rbind(result.rawLobs.full,rawLobs)
  simLobs <- cbind(sim.likearr.bm,sim.likearr.pmm,sim.likearr.ou,sim.likearr.eb,sim.likearr.id)
  result.simLobs.full[[fileIndex]]<-  simLobs
}
colnames(result.rawLobs.full)<-c("bm","pmm","ou","eb","id")
#result.rawLobs.full

quantile.full <- array(NA,dim(result.rawLobs.full))
colnames(quantile.full)<-c("bm","pmm","ou","eb","id")

for(fileIndex in 1:length(filetoload)){
  #  fileIndex<-1
  for(modelIndex in 1:5){
    #  modelIndex <-1
    dim(result.rawLobs.full)
    quantile.full[fileIndex,modelIndex]<- sum(result.rawLobs.full[fileIndex,modelIndex] < result.simLobs.full[[fileIndex]][,modelIndex])/1000
  }
}

#quantile.full



library(ggplot2)
library(gridExtra)

modelset <- c("BM", "PMM", "OU", "EB")

plot_list <- list()

for(modelIndex in 1:length(modelset)){
  df <- data.frame(Quantile = quantile.full[,modelIndex], Data_Index = 1:length(quantile.full[,modelIndex]))
  
  p <- ggplot(df, aes(x = Data_Index, y = Quantile)) + 
    geom_point(aes(color = ifelse(Quantile < 0.10 | Quantile > 0.90, "Outside", "Inside")), shape = 3) +
    geom_hline(yintercept = c(0.10, 0.90), color = "red",linetype=2) +
    labs(title = modelset[modelIndex], x = "Data Index", y = "Quantile") +
    scale_color_manual(values = c("Outside" = "gray", "Inside" = "black")) +
    theme_minimal() +
    theme(legend.position = "none")+
    theme(plot.title = element_text(hjust = 0.5))
  
  plot_list[[modelIndex]] <- p
}

do.call(grid.arrange, c(plot_list, ncol = 2))
```