-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathcorHMM_Diet.r
159 lines (129 loc) · 11.3 KB
/
corHMM_Diet.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
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
library("corHMM")
library("mclust")
library("stringr")
library("qpcR")
# Loading data
phy<-read.nexus("16FC_16C_374_sp.tree")
df<-read.table("table_data_diet.tsv", header = TRUE)
df$Species<-str_replace((df$Species), " ", "_")
states<-cbind(c(str_replace((df$Species), " ", "_"), setdiff(phy$tip.label, str_replace((df$Species), " ", "_"))), c(df$traits, rep("?", length(setdiff(phy$tip.label, str_replace((df$Species), " ", "_"))))))
states_traits<-as.data.frame(states[!states[,1] %in% setdiff(states[,1], phy$tip.label),])
colnames(states_traits)<-c("Species", "Diet")
states_traits<-states_traits[match(phy$tip.label,states_traits[,1]),]
states_traits[is.na(states_traits)]<-"?"
states_traits<-states_traits[,c(1,2)]
LegendAndRateMat <- getStateMat4Dat(states_traits)
RateMat <- LegendAndRateMat$rate.mat
RateMat_trans <- dropStateMatPars(RateMat, c(1,3))
pars2equal <- list(c(1,3), c(2,4))
RateMat_trans_sym <- equateStateMatPars(RateMat_trans, pars2equal)
pars2equal <- list(c(1:4))
RateMat_trans_eq <- equateStateMatPars(RateMat_trans, pars2equal)
diet_1_eq<-corHMM(phy, states_traits, rate.cat = 1, rate.mat=NULL, model = "ER", node.states = "marginal",
fixed.nodes=FALSE, p=NULL, root.p="yang", ip=NULL, nstarts=0, n.cores=5,
get.tip.states = FALSE, lewis.asc.bias = FALSE, collapse = TRUE, lower.bound = 1e-9,
upper.bound = 100, opts=NULL)
diet_1_sym<-corHMM(phy, states_traits, rate.cat = 1, rate.mat=NULL, model = "SYM", node.states = "marginal",
fixed.nodes=FALSE, p=NULL, root.p="yang", ip=NULL, nstarts=0, n.cores=5,
get.tip.states = FALSE, lewis.asc.bias = FALSE, collapse = TRUE, lower.bound = 1e-9,
upper.bound = 100, opts=NULL)
diet_1_ard<-corHMM(phy, states_traits, rate.cat = 1, rate.mat=NULL, model = "ARD", node.states = "marginal",
fixed.nodes=FALSE, p=NULL, root.p="yang", ip=NULL, nstarts=0, n.cores=5,
get.tip.states = FALSE, lewis.asc.bias = FALSE, collapse = TRUE, lower.bound = 1e-9,
upper.bound = 100, opts=NULL)
diet_1_tran_I_M<-corHMM(phy, states_traits, rate.cat = 1, rate.mat=RateMat_trans, node.states = "marginal",
fixed.nodes=FALSE, p=NULL, root.p="yang", ip=NULL, nstarts=0, n.cores=5,
get.tip.states = FALSE, lewis.asc.bias = FALSE, collapse = TRUE, lower.bound = 1e-9,
upper.bound = 100, opts=NULL)
diet_1_tran_I_M_sym<-corHMM(phy, states_traits, rate.cat = 1, rate.mat=RateMat_trans_sym, node.states = "marginal",
fixed.nodes=FALSE, p=NULL, root.p="yang", ip=NULL, nstarts=0, n.cores=5,
get.tip.states = FALSE, lewis.asc.bias = FALSE, collapse = TRUE, lower.bound = 1e-9,
upper.bound = 100, opts=NULL)
diet_1_tran_I_M_eq<-corHMM(phy, states_traits, rate.cat = 1, rate.mat=RateMat_trans_eq, node.states = "marginal",
fixed.nodes=FALSE, p=NULL, root.p="yang", ip=NULL, nstarts=0, n.cores=5,
get.tip.states = FALSE, lewis.asc.bias = FALSE, collapse = TRUE, lower.bound = 1e-9,
upper.bound = 100, opts=NULL)
RateMat_trans_2_rates<-list(RateMat_trans, RateMat_trans)
RateMat_trans_sym_2_rates<-list(RateMat_trans_sym, RateMat_trans_sym)
RateMat_trans_eq_2_rates<-list(RateMat_trans_eq, RateMat_trans_eq)
RateClassMat <- getRateCatMat(2)
RateMat_trans_2_rates <- getFullMat(RateMat_trans_2_rates, RateClassMat)
RateMat_trans_sym_2_rates <- getFullMat(RateMat_trans_sym_2_rates, RateClassMat)
RateMat_trans_eq_2_rates <- getFullMat(RateMat_trans_eq_2_rates, RateClassMat)
diet_2_eq<-corHMM(phy, states_traits, rate.cat = 2, rate.mat=NULL, model = "ER", node.states = "marginal",
fixed.nodes=FALSE, p=NULL, root.p="yang", ip=NULL, nstarts=0, n.cores=5,
get.tip.states = FALSE, lewis.asc.bias = FALSE, collapse = TRUE, lower.bound = 1e-9,
upper.bound = 100, opts=NULL)
diet_2_sym<-corHMM(phy, states_traits, rate.cat = 2, rate.mat=NULL, model = "SYM", node.states = "marginal",
fixed.nodes=FALSE, p=NULL, root.p="yang", ip=NULL, nstarts=0, n.cores=5,
get.tip.states = FALSE, lewis.asc.bias = FALSE, collapse = TRUE, lower.bound = 1e-9,
upper.bound = 100, opts=NULL)
diet_2_ard<-corHMM(phy, states_traits, rate.cat = 2, rate.mat=NULL, model = "ARD", node.states = "marginal",
fixed.nodes=FALSE, p=NULL, root.p="yang", ip=NULL, nstarts=0, n.cores=5,
get.tip.states = FALSE, lewis.asc.bias = FALSE, collapse = TRUE, lower.bound = 1e-9,
upper.bound = 100, opts=NULL)
diet_2_tran_I_M<-corHMM(phy, states_traits, rate.cat = 2, rate.mat= RateMat_trans_2_rates, node.states = "marginal",
fixed.nodes=FALSE, p=NULL, root.p="yang", ip=NULL, nstarts=0, n.cores=5,
get.tip.states = FALSE, lewis.asc.bias = FALSE, collapse = TRUE, lower.bound = 1e-9,
upper.bound = 100, opts=NULL)
diet_2_tran_I_M_sym<-corHMM(phy, states_traits, rate.cat = 2, rate.mat=RateMat_trans_sym_2_rates, node.states = "marginal",
fixed.nodes=FALSE, p=NULL, root.p="yang", ip=NULL, nstarts=0, n.cores=5,
get.tip.states = FALSE, lewis.asc.bias = FALSE, collapse = TRUE, lower.bound = 1e-9,
upper.bound = 100, opts=NULL)
diet_2_tran_I_M_eq<-corHMM(phy, states_traits, rate.cat = 2, rate.mat=RateMat_trans_eq_2_rates, node.states = "marginal",
fixed.nodes=FALSE, p=NULL, root.p="yang", ip=NULL, nstarts=0, n.cores=5,
get.tip.states = FALSE, lewis.asc.bias = FALSE, collapse = TRUE, lower.bound = 1e-9,
upper.bound = 100, opts=NULL)
RateMat_trans_3_rates<-list(RateMat_trans, RateMat_trans, RateMat_trans)
RateMat_trans_sym_3_rates<-list(RateMat_trans_sym, RateMat_trans_sym, RateMat_trans_sym)
RateMat_trans_eq_3_rates<-list(RateMat_trans_eq, RateMat_trans_eq, RateMat_trans_eq)
RateClassMat <- getRateCatMat(3)
RateMat_trans_3_rates <- getFullMat(RateMat_trans_3_rates, RateClassMat)
RateMat_trans_sym_3_rates <- getFullMat(RateMat_trans_sym_3_rates, RateClassMat)
RateMat_trans_eq_3_rates <- getFullMat(RateMat_trans_eq_3_rates, RateClassMat)
diet_3_eq<-corHMM(phy, states_traits, rate.cat = 3, rate.mat=NULL, model = "ER", node.states = "marginal",
fixed.nodes=FALSE, p=NULL, root.p="yang", ip=NULL, nstarts=0, n.cores=5,
get.tip.states = FALSE, lewis.asc.bias = FALSE, collapse = TRUE, lower.bound = 1e-9,
upper.bound = 100, opts=NULL)
diet_3_sym<-corHMM(phy, states_traits, rate.cat = 3, rate.mat=NULL, model = "SYM", node.states = "marginal",
fixed.nodes=FALSE, p=NULL, root.p="yang", ip=NULL, nstarts=0, n.cores=5,
get.tip.states = FALSE, lewis.asc.bias = FALSE, collapse = TRUE, lower.bound = 1e-9,
upper.bound = 100, opts=NULL)
diet_3_ard<-corHMM(phy, states_traits, rate.cat = 3, rate.mat=NULL, model = "ARD", node.states = "marginal",
fixed.nodes=FALSE, p=NULL, root.p="yang", ip=NULL, nstarts=0, n.cores=5,
get.tip.states = FALSE, lewis.asc.bias = FALSE, collapse = TRUE, lower.bound = 1e-9,
upper.bound = 100, opts=NULL)
diet_3_tran_I_M<-corHMM(phy, states_traits, rate.cat = 3, rate.mat=RateMat_trans_3_rates, node.states = "marginal",
fixed.nodes=FALSE, p=NULL, root.p="yang", ip=NULL, nstarts=0, n.cores=5,
get.tip.states = FALSE, lewis.asc.bias = FALSE, collapse = TRUE, lower.bound = 1e-9,
upper.bound = 100, opts=NULL)
diet_3_tran_I_M_sym<-corHMM(phy, states_traits, rate.cat = 3, rate.mat=RateMat_trans_sym_3_rates, node.states = "marginal",
fixed.nodes=FALSE, p=NULL, root.p="yang", ip=NULL, nstarts=0, n.cores=5,
get.tip.states = FALSE, lewis.asc.bias = FALSE, collapse = TRUE, lower.bound = 1e-9,
upper.bound = 100, opts=NULL)
diet_3_tran_I_M_eq<-corHMM(phy, states_traits, rate.cat = 3, rate.mat=RateMat_trans_eq_3_rates, node.states = "marginal",
fixed.nodes=FALSE, p=NULL, root.p="yang", ip=NULL, nstarts=0, n.cores=5,
get.tip.states = FALSE, lewis.asc.bias = FALSE, collapse = TRUE, lower.bound = 1e-9,
upper.bound = 100, opts=NULL)
df_diet<-data.frame(cbind(c(diet_1_eq$loglik, diet_1_sym$loglik, diet_1_ard$loglik, diet_1_tran_I_M$loglik, diet_1_tran_I_M_sym$loglik, diet_1_tran_I_M_eq$loglik,
diet_2_eq$loglik, diet_2_sym$loglik, diet_2_ard$loglik, diet_2_tran_I_M$loglik, diet_2_tran_I_M_sym$loglik, diet_2_tran_I_M_eq$loglik,
diet_3_eq$loglik, diet_3_sym$loglik, diet_3_ard$loglik, diet_3_tran_I_M$loglik, diet_3_tran_I_M_sym$loglik, diet_3_tran_I_M_eq$loglik),
c(diet_1_eq$AICc, diet_1_sym$AICc, diet_1_ard$AICc, diet_1_tran_I_M$AICc, diet_1_tran_I_M_sym$AICc, diet_1_tran_I_M_eq$AICc,
diet_2_eq$AICc, diet_2_sym$AICc, diet_2_ard$AICc, diet_2_tran_I_M$AICc, diet_2_tran_I_M_sym$AICc, diet_2_tran_I_M_eq$AICc,
diet_3_eq$AICc, diet_3_sym$AICc, diet_3_ard$AICc, diet_3_tran_I_M$AICc, diet_3_tran_I_M_sym$AICc, diet_3_tran_I_M_eq$AICc),
akaike.weights(c(diet_1_eq$AICc, diet_1_sym$AICc, diet_1_ard$AICc, diet_1_tran_I_M$AICc, diet_1_tran_I_M_sym$AICc, diet_1_tran_I_M_eq$AICc,
diet_2_eq$AICc, diet_2_sym$AICc, diet_2_ard$AICc, diet_2_tran_I_M$AICc, diet_2_tran_I_M_sym$AICc, diet_2_tran_I_M_eq$AICc,
diet_3_eq$AICc, diet_3_sym$AICc, diet_3_ard$AICc, diet_3_tran_I_M$AICc, diet_3_tran_I_M_sym$AICc, diet_3_tran_I_M_eq$AICc))$deltaAIC,
akaike.weights(c(diet_1_eq$AICc, diet_1_sym$AICc, diet_1_ard$AICc, diet_1_tran_I_M$AICc, diet_1_tran_I_M_sym$AICc, diet_1_tran_I_M_eq$AICc,
diet_2_eq$AICc, diet_2_sym$AICc, diet_2_ard$AICc, diet_2_tran_I_M$AICc, diet_2_tran_I_M_sym$AICc, diet_2_tran_I_M_eq$AICc,
diet_3_eq$AICc, diet_3_sym$AICc, diet_3_ard$AICc, diet_3_tran_I_M$AICc, diet_3_tran_I_M_sym$AICc, diet_3_tran_I_M_eq$AICc))$weights,
c((max(as.vector(diet_1_eq$index.mat)[!is.na(as.vector(diet_1_eq$index.mat))])), (max(as.vector(diet_1_sym$index.mat)[!is.na(as.vector(diet_1_sym$index.mat))])), (max(as.vector(diet_1_ard$index.mat)[!is.na(as.vector(diet_1_ard$index.mat))])), (max(as.vector(diet_1_tran_I_M$index.mat)[!is.na(as.vector(diet_1_tran_I_M$index.mat))])), (max(as.vector(diet_1_tran_I_M_sym$index.mat)[!is.na(as.vector(diet_1_tran_I_M_sym$index.mat))])), (max(as.vector(diet_1_tran_I_M_eq$index.mat)[!is.na(as.vector(diet_1_tran_I_M_eq$index.mat))])),
(max(as.vector(diet_2_eq$index.mat)[!is.na(as.vector(diet_2_eq$index.mat))])), (max(as.vector(diet_2_sym$index.mat)[!is.na(as.vector(diet_2_sym$index.mat))])), (max(as.vector(diet_2_ard$index.mat)[!is.na(as.vector(diet_2_ard$index.mat))])), (max(as.vector(diet_2_tran_I_M$index.mat)[!is.na(as.vector(diet_2_tran_I_M$index.mat))])), (max(as.vector(diet_2_tran_I_M_sym$index.mat)[!is.na(as.vector(diet_2_tran_I_M_sym$index.mat))])), (max(as.vector(diet_2_tran_I_M_eq$index.mat)[!is.na(as.vector(diet_2_tran_I_M_eq$index.mat))])),
(max(as.vector(diet_3_eq$index.mat)[!is.na(as.vector(diet_3_eq$index.mat))])), (max(as.vector(diet_3_sym$index.mat)[!is.na(as.vector(diet_3_sym$index.mat))])), (max(as.vector(diet_3_ard$index.mat)[!is.na(as.vector(diet_3_ard$index.mat))])), (max(as.vector(diet_3_tran_I_M$index.mat)[!is.na(as.vector(diet_3_tran_I_M$index.mat))])), (max(as.vector(diet_3_tran_I_M_sym$index.mat)[!is.na(as.vector(diet_3_tran_I_M_sym$index.mat))])), (max(as.vector(diet_3_tran_I_M_eq$index.mat)[!is.na(as.vector(diet_3_tran_I_M_eq$index.mat))])))
))
rownames(df_diet)<-c("diet_1_eq", "diet_1_sym", "diet_1_ard", "diet_1_tran_I_M", "diet_1_tran_I_M_sym", "diet_1_tran_I_M_eq",
"diet_2_eq", "diet_2_sym", "diet_2_ard", "diet_2_tran_I_M", "diet_2_tran_I_M_sym", "diet_2_tran_I_M_eq",
"diet_3_eq", "diet_3_sym", "diet_3_ard", "diet_3_tran_I_M", "diet_3_tran_I_M_sym", "diet_3_tran_I_M_eq")
colnames(df_diet)<-c("loglik", "AICc", "Delta_AICc", "AICcWt", "K_rates")
write.table(df_diet, "Results/df_diet.tsv", sep ="\t")
saveRDS(eval(parse(text = rownames(df_diet)[which.min(df_diet$AICc)])), "ASE/diet_corHMM.rds")