-
Notifications
You must be signed in to change notification settings - Fork 0
/
GelmanRubinTest.R
103 lines (77 loc) · 3.2 KB
/
GelmanRubinTest.R
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
# First part of the code gets Gelman Rubin diagnostic for Mean behavior variables. Second part for the 144
# covariance terms of the Phylo matrix corresponding to the 12 most common behaviors. Third part
# the Gelman Rubin test for the covariance terms of the individual matrix corresponding to the same 12 most
#common behaviors.
library(coda,lib.loc="Rpackages")
ListPsrf=c()
for(k in 1:134) {
FileName=paste("~/ownCloud/flytree/data/Anc134BhModel0PrTrueAllMinusYak00",1,sep="")
load(FileName)
Temp=model0$Sol
ChainMatrix=Temp[,k]
for(j in 2:20) {
FileName=paste("~/ownCloud/flytree/data/Anc134BhModel0PrTrueAllMinusYak00",j,sep="")
load(FileName)
Temp=model0$Sol
ChainT=Temp[,k]
ChainMatrix=rbind(ChainT,ChainMatrix)
rm(Temp)
}
ChainsMCMC=as.mcmc.list(lapply(as.data.frame(ChainMatrix), mcmc))
PSRFTemp=gelman.diag(ChainsMCMC, confidence = 0.95, transform=FALSE, autoburnin=FALSE,
multivariate=TRUE)
ListPsrf=cbind(c(PSRFTemp$psrf[[1]], PSRFTemp$psrf[[2]]),ListPsrf)
print(ListPsrf)
rm(ChainT,ChainMatrix)
}
writeMat("~/ownCloud/flytree/data/MeanVariablesGelmanRubinTest.mat", ListPsrf=ListPsrf)
# Gelman Rubin for the 10% most common behaviors for Phylo
FileName="~/ownCloud/flytree/data/IndexPhyloHighBhForR.mat"
Index=readMat(FileName)
ListPsrfVCVPhylo=c()
for(k in 1:length(Index$InR)) {
Filename=paste("~/ownCloud/flytree/data/Anc134BhModel0PrTrueAllMinusYak00",1,sep="")
load(Filename)
Temp=model0$VCV
ChainMatrix=Temp[,Index$InR[[k]]]
for(j in 2:20) {
FileName=paste("~/ownCloud/flytree/data/Anc134BhModel0PrTrueAllMinusYak00",j,sep="")
load(FileName)
Temp=model0$VCV
ChainT=Temp[,Index$InR[[k]]]
ChainMatrix=rbind(ChainT,ChainMatrix)
rm(Temp)
}
ChainsMCMC=as.mcmc.list(lapply(as.data.frame(ChainMatrix), mcmc))
PSRFTempVCV=gelman.diag(ChainsMCMC, confidence = 0.95, transform=TRUE, autoburnin=FALSE,
multivariate=TRUE)
ListPsrfVCVPhylo=cbind(c(PSRFTempVCV$psrf[[1]], PSRFTempVCV$psrf[[2]]),ListPsrfVCVPhylo)
print(ListPsrfVCVPhylo)
rm(ChainT,ChainMatrix)
}
writeMat("~/ownCloud/flytree/data/VCVPhyloVariablesGelmanRubinTest.mat", ListPsrfVCVPhylo=ListPsrfVCVPhylo)
# Gelman Rubin for the 10% most common behaviors for Individual matrix
FileName="~/ownCloud/flytree/data/IndexIndiHighBhForR.mat"
Index=readMat(FileName)
ListPsrfVCVIndi=c()
for(k in 1:length(Index$InR)) {
Filename=paste("~/ownCloud/flytree/data/Anc134BhModel0PrTrueAllMinusYak00",1,sep="")
load(Filename)
Temp=model0$VCV
ChainMatrix=Temp[,Index$InR[[k]]]
for(j in 2:20) {
FileName=paste("~/ownCloud/flytree/data/Anc134BhModel0PrTrueAllMinusYak00",j,sep="")
load(FileName)
Temp=model0$VCV
ChainT=Temp[,Index$InR[[k]]]
ChainMatrix=rbind(ChainT,ChainMatrix)
rm(Temp)
}
ChainsMCMC=as.mcmc.list(lapply(as.data.frame(ChainMatrix), mcmc))
PSRFTempVCV=gelman.diag(ChainsMCMC, confidence = 0.95, transform=TRUE, autoburnin=FALSE,
multivariate=TRUE)
ListPsrfVCVIndi=cbind(c(PSRFTempVCV$psrf[[1]], PSRFTempVCV$psrf[[2]]),ListPsrfVCVIndi)
print(ListPsrfVCVIndi)
rm(ChainT,ChainMatrix)
}
writeMat("~/ownCloud/flytree/data/VCVIndiVariablesGelmanRubinTest.mat", ListPsrfVCVIndi=ListPsrfVCVIndi)