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