-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathcorHMM.r
445 lines (365 loc) · 31.7 KB
/
corHMM.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
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
library("corHMM")
library("mclust")
library("stringr")
library("qpcR")
args<-commandArgs(trailingOnly = TRUE)
# Loading data
df<-read.csv("Species_data.tsv", sep="\t") # omit sep ="\t" for .csv files
phy<-read.nexus("16FC_16C_374_sp.tree")
mb1 = Mclust(as.numeric(df$Body.size))
summary(mb, parameters = TRUE)
df$Body.size <- as.factor(mb1$classification)
df$Species<-str_replace(df$Species, " ", "_")
setdiff(phy$tip.label, df$Species)
states<-cbind(df$Species, df$Body.size, df$Reproduction, df$Habitat)
states_traits<-states[!states[,1] %in% setdiff(states[,1], phy$tip.label),]
states_traits[is.na(states_traits)]<-"?"
states_traits<-states_traits[match(phy$tip.label,states_traits[,1]),]
states_size<-states_traits[,c(1,2)]
LegendAndRateMat <- getStateMat4Dat(states_size)
RateMat <- LegendAndRateMat$rate.mat
RateMat_trans <- dropStateMatPars(RateMat, c(2,5))
pars2equal <- list(c(1,2), c(3,4))
RateMat_trans_sym <- equateStateMatPars(RateMat_trans, pars2equal)
pars2equal <- list(c(1:4))
RateMat_trans_eq <- equateStateMatPars(RateMat_trans, pars2equal)
bds_1_eq<-corHMM(phy, states_size, 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)
bds_1_sym<-corHMM(phy, states_size, 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)
bds_1_ard<-corHMM(phy, states_size, 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)
bds_1_tran_s_l<-corHMM(phy, states_size, 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)
bds_1_tran_s_l_sym<-corHMM(phy, states_size, 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)
bds_1_tran_s_l_eq<-corHMM(phy, states_size, 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)
bds_2_eq<-corHMM(phy, states_size, 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)
bds_2_sym<-corHMM(phy, states_size, 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)
bds_2_ard<-corHMM(phy, states_size, 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)
bds_2_tran_s_l<-corHMM(phy, states_size, 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)
bds_2_tran_s_l_sym<-corHMM(phy, states_size, 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)
bds_2_tran_s_l_eq<-corHMM(phy, states_size, 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)
bds_3_eq<-corHMM(phy, states_size, 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)
bds_3_sym<-corHMM(phy, states_size, 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)
bds_3_ard<-corHMM(phy, states_size, 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)
bds_3_tran_s_l<-corHMM(phy, states_size, 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)
bds_3_tran_s_l_sym<-corHMM(phy, states_size, 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)
bds_3_tran_s_l_eq<-corHMM(phy, states_size, 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_size<-data.frame(cbind(c(bds_1_eq$loglik, bds_1_sym$loglik, bds_1_ard$loglik, bds_1_tran_s_l$loglik, bds_1_tran_s_l_sym$loglik, bds_1_tran_s_l_eq$loglik,
bds_2_eq$loglik, bds_2_sym$loglik, bds_2_ard$loglik, bds_2_tran_s_l$loglik, bds_2_tran_s_l_sym$loglik, bds_2_tran_s_l_eq$loglik,
bds_3_eq$loglik, bds_3_sym$loglik, bds_3_ard$loglik, bds_3_tran_s_l$loglik, bds_3_tran_s_l_sym$loglik, bds_3_tran_s_l_eq$loglik),
c(bds_1_eq$AICc, bds_1_sym$AICc, bds_1_ard$AICc, bds_1_tran_s_l$AICc, bds_1_tran_s_l_sym$AICc, bds_1_tran_s_l_eq$AICc,
bds_2_eq$AICc, bds_2_sym$AICc, bds_2_ard$AICc, bds_2_tran_s_l$AICc, bds_2_tran_s_l_sym$AICc, bds_2_tran_s_l_eq$AICc,
bds_3_eq$AICc, bds_3_sym$AICc, bds_3_ard$AICc, bds_3_tran_s_l$AICc, bds_3_tran_s_l_sym$AICc, bds_3_tran_s_l_eq$AICc),
akaike.weights(c(bds_1_eq$AICc, bds_1_sym$AICc, bds_1_ard$AICc, bds_1_tran_s_l$AICc, bds_1_tran_s_l_sym$AICc, bds_1_tran_s_l_eq$AICc,
bds_2_eq$AICc, bds_2_sym$AICc, bds_2_ard$AICc, bds_2_tran_s_l$AICc, bds_2_tran_s_l_sym$AICc, bds_2_tran_s_l_eq$AICc,
bds_3_eq$AICc, bds_3_sym$AICc, bds_3_ard$AICc, bds_3_tran_s_l$AICc, bds_3_tran_s_l_sym$AICc, bds_3_tran_s_l_eq$AICc))$deltaAIC,
akaike.weights(c(bds_1_eq$AICc, bds_1_sym$AICc, bds_1_ard$AICc, bds_1_tran_s_l$AICc, bds_1_tran_s_l_sym$AICc, bds_1_tran_s_l_eq$AICc,
bds_2_eq$AICc, bds_2_sym$AICc, bds_2_ard$AICc, bds_2_tran_s_l$AICc, bds_2_tran_s_l_sym$AICc, bds_2_tran_s_l_eq$AICc,
bds_3_eq$AICc, bds_3_sym$AICc, bds_3_ard$AICc, bds_3_tran_s_l$AICc, bds_3_tran_s_l_sym$AICc, bds_3_tran_s_l_eq$AICc))$weights,
c((max(as.vector(bds_1_eq$index.mat)[!is.na(as.vector(bds_1_eq$index.mat))])), (max(as.vector(bds_1_sym$index.mat)[!is.na(as.vector(bds_1_sym$index.mat))])), (max(as.vector(bds_1_ard$index.mat)[!is.na(as.vector(bds_1_ard$index.mat))])), (max(as.vector(bds_1_tran_s_l$index.mat)[!is.na(as.vector(bds_1_tran_s_l$index.mat))])), (max(as.vector(bds_1_tran_s_l_sym$index.mat)[!is.na(as.vector(bds_1_tran_s_l_sym$index.mat))])), (max(as.vector(bds_1_tran_s_l_eq$index.mat)[!is.na(as.vector(bds_1_tran_s_l_eq$index.mat))])),
(max(as.vector(bds_2_eq$index.mat)[!is.na(as.vector(bds_2_eq$index.mat))])), (max(as.vector(bds_2_sym$index.mat)[!is.na(as.vector(bds_2_sym$index.mat))])), (max(as.vector(bds_2_ard$index.mat)[!is.na(as.vector(bds_2_ard$index.mat))])), (max(as.vector(bds_2_tran_s_l$index.mat)[!is.na(as.vector(bds_2_tran_s_l$index.mat))])), (max(as.vector(bds_2_tran_s_l_sym$index.mat)[!is.na(as.vector(bds_2_tran_s_l_sym$index.mat))])), (max(as.vector(bds_2_tran_s_l_eq$index.mat)[!is.na(as.vector(bds_2_tran_s_l_eq$index.mat))])),
(max(as.vector(bds_3_eq$index.mat)[!is.na(as.vector(bds_3_eq$index.mat))])), (max(as.vector(bds_3_sym$index.mat)[!is.na(as.vector(bds_3_sym$index.mat))])), (max(as.vector(bds_3_ard$index.mat)[!is.na(as.vector(bds_3_ard$index.mat))])), (max(as.vector(bds_3_tran_s_l$index.mat)[!is.na(as.vector(bds_3_tran_s_l$index.mat))])), (max(as.vector(bds_3_tran_s_l_sym$index.mat)[!is.na(as.vector(bds_3_tran_s_l_sym$index.mat))])), (max(as.vector(bds_3_tran_s_l_eq$index.mat)[!is.na(as.vector(bds_3_tran_s_l_eq$index.mat))])))
))
rownames(df_size)<-c("bds_1_eq", "bds_1_sym", "bds_1_ard", "bds_1_tran_s_l", "bds_1_tran_s_l_sym", "bds_1_tran_s_l_eq",
"bds_2_eq", "bds_2_sym", "bds_2_ard", "bds_2_tran_s_l", "bds_2_tran_s_l_sym", "bds_2_tran_s_l_eq",
"bds_3_eq", "bds_3_sym", "bds_3_ard", "bds_3_tran_s_l", "bds_3_tran_s_l_sym", "bds_3_tran_s_l_eq")
colnames(df_size)<-c("loglik", "AICc", "Delta_AICc", "AICcWt", "K_rates")
write.table(df_size, "ASE/df_size.tsv", sep ="\t")
saveRDS(eval(parse(text = rownames(df_size)[which.min(df_size$AICc)])), "ASE/bds_corHMM.rds")
states_reproduction<-states_traits[,c(1,3)]
LegendAndRateMat <- getStateMat4Dat(states_reproduction)
RateMat <- LegendAndRateMat$rate.mat
RateMat_trans_o <- dropStateMatPars(RateMat, c(1,2,4,7))
pars2equal <- list(c(1,6), c(2,4), c(3,7), c(5,8))
RateMat_trans_o_sym <- equateStateMatPars(RateMat_trans_o, pars2equal)
pars2equal <- list(c(1:8))
RateMat_trans_o_eq <- equateStateMatPars(RateMat_trans_o, pars2equal)
rp_1_eq<-corHMM(phy, states_reproduction, 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)
rp_1_sym<-corHMM(phy, states_reproduction, 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)
rp_1_ard<-corHMM(phy, states_reproduction, 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)
rp_1_trans_o<-corHMM(phy, states_reproduction, rate.cat = 1, rate.mat= RateMat_trans_o, 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)
rp_1_trans_o_sym<-corHMM(phy, states_reproduction, rate.cat = 1, rate.mat= RateMat_trans_o_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)
rp_1_trans_o_eq<-corHMM(phy, states_reproduction, rate.cat = 1, rate.mat= RateMat_trans_o_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_o_2_rates<-list(RateMat_trans_o, RateMat_trans_o)
RateMat_trans_o_sym_2_rates<-list(RateMat_trans_o_sym, RateMat_trans_o_sym)
RateMat_trans_o_eq_2_rates<-list(RateMat_trans_o_eq, RateMat_trans_o_eq)
RateClassMat <- getRateCatMat(2)
RateMat_trans_o_2_rates <- getFullMat(RateMat_trans_o_2_rates, RateClassMat)
RateMat_trans_o_sym_2_rates <- getFullMat(RateMat_trans_o_sym_2_rates, RateClassMat)
RateMat_trans_o_eq_2_rates <- getFullMat(RateMat_trans_o_eq_2_rates, RateClassMat)
rp_2_eq<-corHMM(phy, states_reproduction, 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)
rp_2_sym<-corHMM(phy, states_reproduction, 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)
rp_2_ard<-corHMM(phy, states_reproduction, 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)
rp_2_trans_o<-corHMM(phy, states_reproduction, rate.cat = 2, rate.mat= RateMat_trans_o_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)
rp_2_trans_o_sym<-corHMM(phy, states_reproduction, rate.cat = 2, rate.mat= RateMat_trans_o_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)
rp_2_trans_o_eq<-corHMM(phy, states_reproduction, rate.cat = 2, rate.mat= RateMat_trans_o_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_o_3_rates<-list(RateMat_trans_o, RateMat_trans_o, RateMat_trans_o)
RateMat_trans_o_sym_3_rates<-list(RateMat_trans_o_sym, RateMat_trans_o_sym, RateMat_trans_o_sym)
RateMat_trans_o_eq_3_rates<-list(RateMat_trans_o_eq, RateMat_trans_o_eq, RateMat_trans_o_eq)
RateClassMat <- getRateCatMat(3)
RateMat_trans_o_3_rates <- getFullMat(RateMat_trans_o_3_rates, RateClassMat)
RateMat_trans_o_sym_3_rates <- getFullMat(RateMat_trans_o_sym_3_rates, RateClassMat)
RateMat_trans_o_eq_3_rates <- getFullMat(RateMat_trans_o_eq_3_rates, RateClassMat)
rp_3_eq<-corHMM(phy, states_reproduction, 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)
rp_3_sym<-corHMM(phy, states_reproduction, 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)
rp_3_ard<-corHMM(phy, states_reproduction, 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)
rp_3_trans_o<-corHMM(phy, states_reproduction, rate.cat = 3, rate.mat= RateMat_trans_o_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)
rp_3_trans_o_sym<-corHMM(phy, states_reproduction, rate.cat = 3, rate.mat= RateMat_trans_o_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)
rp_3_trans_o_eq<-corHMM(phy, states_reproduction, rate.cat = 3, rate.mat= RateMat_trans_o_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_reproduction<-data.frame(cbind(c(rp_1_eq$loglik, rp_1_sym$loglik, rp_1_ard$loglik, rp_1_trans_o$loglik, rp_1_trans_o_sym$loglik, rp_1_trans_o_eq$loglik,
rp_2_eq$loglik, rp_2_sym$loglik, rp_2_ard$loglik, rp_2_trans_o$loglik, rp_2_trans_o_sym$loglik, rp_2_trans_o_eq$loglik,
rp_3_eq$loglik, rp_3_sym$loglik, rp_3_ard$loglik, rp_3_trans_o$loglik, rp_3_trans_o_sym$loglik, rp_3_trans_o_eq$loglik),
c(rp_1_eq$AICc, rp_1_sym$AICc, rp_1_ard$AICc, rp_1_trans_o$AICc, rp_1_trans_o_sym$AICc, rp_1_trans_o_eq$AICc,
rp_2_eq$AICc, rp_2_sym$AICc, rp_2_ard$AICc, rp_2_trans_o$AICc, rp_2_trans_o_sym$AICc, rp_2_trans_o_eq$AICc,
rp_3_eq$AICc, rp_3_sym$AICc, rp_3_ard$AICc, rp_3_trans_o$AICc, rp_3_trans_o_sym$AICc, rp_3_trans_o_eq$AICc),
akaike.weights(c(rp_1_eq$AICc, rp_1_sym$AICc, rp_1_ard$AICc, rp_1_trans_o$AICc, rp_1_trans_o_sym$AICc, rp_1_trans_o_eq$AICc,
rp_2_eq$AICc, rp_2_sym$AICc, rp_2_ard$AICc, rp_2_trans_o$AICc, rp_2_trans_o_sym$AICc, rp_2_trans_o_eq$AICc,
rp_3_eq$AICc, rp_3_sym$AICc, rp_3_ard$AICc, rp_3_trans_o$AICc, rp_3_trans_o_sym$AICc, rp_3_trans_o_eq$AICc))$deltaAIC,
akaike.weights(c(rp_1_eq$AICc, rp_1_sym$AICc, rp_1_ard$AICc, rp_1_trans_o$AICc, rp_1_trans_o_sym$AICc, rp_1_trans_o_eq$AICc,
rp_2_eq$AICc, rp_2_sym$AICc, rp_2_ard$AICc, rp_2_trans_o$AICc, rp_2_trans_o_sym$AICc, rp_2_trans_o_eq$AICc,
rp_3_eq$AICc, rp_3_sym$AICc, rp_3_ard$AICc, rp_3_trans_o$AICc, rp_3_trans_o_sym$AICc, rp_3_trans_o_eq$AICc))$weights,
c((max(as.vector(rp_1_eq$index.mat)[!is.na(as.vector(rp_1_eq$index.mat))])), (max(as.vector(rp_1_sym$index.mat)[!is.na(as.vector(rp_1_sym$index.mat))])), (max(as.vector(rp_1_ard$index.mat)[!is.na(as.vector(rp_1_ard$index.mat))])), (max(as.vector(rp_1_trans_o$index.mat)[!is.na(as.vector(rp_1_trans_o$index.mat))])), (max(as.vector(rp_1_trans_o_sym$index.mat)[!is.na(as.vector(rp_1_trans_o_sym$index.mat))])), (max(as.vector(rp_1_trans_o_eq$index.mat)[!is.na(as.vector(rp_1_trans_o_eq$index.mat))])),
(max(as.vector(rp_2_eq$index.mat)[!is.na(as.vector(rp_2_eq$index.mat))])), (max(as.vector(rp_2_sym$index.mat)[!is.na(as.vector(rp_2_sym$index.mat))])), (max(as.vector(rp_2_ard$index.mat)[!is.na(as.vector(rp_2_ard$index.mat))])), (max(as.vector(rp_2_trans_o$index.mat)[!is.na(as.vector(rp_2_trans_o$index.mat))])), (max(as.vector(rp_2_trans_o_sym$index.mat)[!is.na(as.vector(rp_2_trans_o_sym$index.mat))])), (max(as.vector(rp_2_trans_o_eq$index.mat)[!is.na(as.vector(rp_2_trans_o_eq$index.mat))])),
(max(as.vector(rp_3_eq$index.mat)[!is.na(as.vector(rp_3_eq$index.mat))])), (max(as.vector(rp_3_sym$index.mat)[!is.na(as.vector(rp_3_sym$index.mat))])), (max(as.vector(rp_3_ard$index.mat)[!is.na(as.vector(rp_3_ard$index.mat))])), (max(as.vector(rp_3_trans_o$index.mat)[!is.na(as.vector(rp_3_trans_o$index.mat))])), (max(as.vector(rp_3_trans_o_sym$index.mat)[!is.na(as.vector(rp_3_trans_o_sym$index.mat))])), (max(as.vector(rp_3_trans_o_eq$index.mat)[!is.na(as.vector(rp_3_trans_o_eq$index.mat))])))
))
rownames(df_reproduction)<-c("rp_1_eq", "rp_1_sym", "rp_1_ard", "rp_1_trans_o", "rp_1_trans_o_sym", "rp_1_trans_o_eq",
"rp_2_eq", "rp_2_sym", "rp_2_ard", "rp_2_trans_o", "rp_2_trans_o_sym", "rp_2_trans_o_eq",
"rp_3_eq", "rp_3_sym", "rp_3_ard", "rp_3_trans_o", "rp_3_trans_o_sym", "rp_3_trans_o_eq")
colnames(df_reproduction)<-c("logLik","AICc", "Delta_AICc", "AICcWt", "K_rates")
write.table(df_reproduction, "ASE/df_reproduction.tsv", sep ="\t")
saveRDS(eval(parse(text = rownames(df_reproduction)[which.min(df_reproduction$AICc)])), "ASE/rp_corHMM.rds")
states_habitat<-states_traits[,c(1,4)]
LegendAndRateMat <- getStateMat4Dat(states_habitat)
RateMat <- LegendAndRateMat$rate.mat
RateMat_trans_d <- dropStateMatPars(RateMat, c(5,9,17,1,2,4))
pars2equal <- list(c(1,8), c(2,5), c(3,9), c(4,12), c(6,10), c(13,7), c(11,14))
RateMat_trans_d_sym <- equateStateMatPars(RateMat_trans_d, pars2equal)
pars2equal <- list(c(1:14))
RateMat_trans_d_eq <- equateStateMatPars(RateMat_trans_d, pars2equal)
ht_1_eq<-corHMM(phy, states_habitat, 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)
ht_1_sym<-corHMM(phy, states_habitat, 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)
ht_1_ard<-corHMM(phy, states_habitat, 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)
ht_1_d<-corHMM(phy, states_habitat, rate.cat = 1, rate.mat=RateMat_trans_d, 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)
ht_1_d_sym<-corHMM(phy, states_habitat, rate.cat = 1, rate.mat=RateMat_trans_d_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)
ht_1_d_eq<-corHMM(phy, states_habitat, rate.cat = 1, rate.mat=RateMat_trans_d_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_d_2_rates<-list(RateMat_trans_d, RateMat_trans_d)
RateMat_trans_d_sym_2_rates<-list(RateMat_trans_d_sym, RateMat_trans_d_sym)
RateMat_trans_d_eq_2_rates<-list(RateMat_trans_d_eq, RateMat_trans_d_eq)
RateClassMat <- getRateCatMat(2)
RateMat_trans_d_2_rates <- getFullMat(RateMat_trans_d_2_rates, RateClassMat)
RateMat_trans_d_sym_2_rates <- getFullMat(RateMat_trans_d_sym_2_rates, RateClassMat)
RateMat_trans_d_eq_2_rates <- getFullMat(RateMat_trans_d_eq_2_rates, RateClassMat)
ht_2_eq<-corHMM(phy, states_habitat, 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)
ht_2_sym<-corHMM(phy, states_habitat, 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)
ht_2_ard<-corHMM(phy, states_habitat, 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)
ht_2_d<-corHMM(phy, states_habitat, rate.cat = 1, rate.mat=RateMat_trans_d_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)
ht_2_d_sym<-corHMM(phy, states_habitat, rate.cat = 1, rate.mat=RateMat_trans_d_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)
ht_2_d_eq<-corHMM(phy, states_habitat, rate.cat = 1, rate.mat=RateMat_trans_d_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_d_3_rates<-list(RateMat_trans_d, RateMat_trans_d, RateMat_trans_d)
RateMat_trans_d_sym_3_rates<-list(RateMat_trans_d_sym, RateMat_trans_d_sym, RateMat_trans_d_sym)
RateMat_trans_d_eq_3_rates<-list(RateMat_trans_d_eq, RateMat_trans_d_eq, RateMat_trans_d_eq)
RateClassMat <- getRateCatMat(3)
RateMat_trans_d_3_rates <- getFullMat(RateMat_trans_d_3_rates, RateClassMat)
RateMat_trans_d_sym_3_rates <- getFullMat(RateMat_trans_d_sym_3_rates, RateClassMat)
RateMat_trans_d_eq_3_rates <- getFullMat(RateMat_trans_d_eq_3_rates, RateClassMat)
ht_3_eq<-corHMM(phy, states_habitat, 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)
ht_3_sym<-corHMM(phy, states_habitat, 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)
ht_3_ard<-corHMM(phy, states_habitat, 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)
ht_3_d<-corHMM(phy, states_habitat, rate.cat = 1, rate.mat=RateMat_trans_d_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)
ht_3_d_sym<-corHMM(phy, states_habitat, rate.cat = 1, rate.mat=RateMat_trans_d_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)
ht_3_d_eq<-corHMM(phy, states_habitat, rate.cat = 1, rate.mat=RateMat_trans_d_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_habitat<-data.frame(cbind(c(ht_1_eq$loglik, ht_1_sym$loglik, ht_1_ard$loglik, ht_1_d$loglik, ht_1_d_sym$loglik, ht_1_d_eq$loglik,
ht_2_eq$loglik, ht_2_sym$loglik, ht_2_ard$loglik, ht_2_d$loglik, ht_2_d_sym$loglik, ht_2_d_eq$loglik,
ht_3_eq$loglik, ht_3_sym$loglik, ht_3_ard$loglik, ht_3_d$loglik, ht_3_d_sym$loglik, ht_3_d_eq$loglik),
c(ht_1_eq$AICc, ht_1_sym$AICc, ht_1_ard$AICc, ht_1_d$AICc, ht_1_d_sym$AICc, ht_1_d_eq$AICc,
ht_2_eq$AICc, ht_2_sym$AICc, ht_2_ard$AICc, ht_2_d$AICc, ht_2_d_sym$AICc, ht_2_d_eq$AICc,
ht_3_eq$AICc, ht_3_sym$AICc, ht_3_ard$AICc, ht_3_d$AICc, ht_3_d_sym$AICc, ht_3_d_eq$AICc),
akaike.weights(c(ht_1_eq$AICc, ht_1_sym$AICc, ht_1_ard$AICc, ht_1_d$AICc, ht_1_d_sym$AICc, ht_1_d_eq$AICc,
ht_2_eq$AICc, ht_2_sym$AICc, ht_2_ard$AICc, ht_2_d$AICc, ht_2_d_sym$AICc, ht_2_d_eq$AICc,
ht_3_eq$AICc, ht_3_sym$AICc, ht_3_ard$AICc, ht_3_d$AICc, ht_3_d_sym$AICc, ht_3_d_eq$AICc))$deltaAIC,
akaike.weights(c(ht_1_eq$AICc, ht_1_sym$AICc, ht_1_ard$AICc, ht_1_d$AICc, ht_1_d_sym$AICc, ht_1_d_eq$AICc,
ht_2_eq$AICc, ht_2_sym$AICc, ht_2_ard$AICc, ht_2_d$AICc, ht_2_d_sym$AICc, ht_2_d_eq$AICc,
ht_3_eq$AICc, ht_3_sym$AICc, ht_3_ard$AICc, ht_3_d$AICc, ht_3_d_sym$AICc, ht_3_d_eq$AICc))$weights,
c((max(as.vector(ht_1_eq$index.mat)[!is.na(as.vector(ht_1_eq$index.mat))])), (max(as.vector(ht_1_sym$index.mat)[!is.na(as.vector(ht_1_sym$index.mat))])), (max(as.vector(ht_1_ard$index.mat)[!is.na(as.vector(ht_1_ard$index.mat))])), (max(as.vector(ht_1_d$index.mat)[!is.na(as.vector(ht_1_d$index.mat))])), (max(as.vector(ht_1_d_sym$index.mat)[!is.na(as.vector(ht_1_d_sym$index.mat))])), (max(as.vector(ht_1_d_eq$index.mat)[!is.na(as.vector(ht_1_d_eq$index.mat))])),
(max(as.vector(ht_2_eq$index.mat)[!is.na(as.vector(ht_2_eq$index.mat))])), (max(as.vector(ht_2_sym$index.mat)[!is.na(as.vector(ht_2_sym$index.mat))])), (max(as.vector(ht_2_ard$index.mat)[!is.na(as.vector(ht_2_ard$index.mat))])), (max(as.vector(ht_2_d$index.mat)[!is.na(as.vector(ht_2_d$index.mat))])), (max(as.vector(ht_2_d_sym$index.mat)[!is.na(as.vector(ht_2_d_sym$index.mat))])), (max(as.vector(ht_2_d_eq$index.mat)[!is.na(as.vector(ht_2_d_eq$index.mat))])),
(max(as.vector(ht_3_eq$index.mat)[!is.na(as.vector(ht_3_eq$index.mat))])), (max(as.vector(ht_3_sym$index.mat)[!is.na(as.vector(ht_3_sym$index.mat))])), (max(as.vector(ht_3_ard$index.mat)[!is.na(as.vector(ht_3_ard$index.mat))])), (max(as.vector(ht_3_d$index.mat)[!is.na(as.vector(ht_3_d$index.mat))])), (max(as.vector(ht_3_d_sym$index.mat)[!is.na(as.vector(ht_3_d_sym$index.mat))])), (max(as.vector(ht_3_d_eq$index.mat)[!is.na(as.vector(ht_3_d_eq$index.mat))])))
))
rownames(df_habitat)<-c("ht_1_eq", "ht_1_sym", "ht_1_ard", "ht_1_d", "ht_1_d_sym", "ht_1_d_eq",
"ht_2_eq", "ht_2_sym", "ht_2_ard", "ht_2_d", "ht_2_d_sym", "ht_2_d_eq",
"ht_3_eq", "ht_3_sym", "ht_3_ard", "ht_3_d", "ht_3_d_sym", "ht_3_d_eq")
colnames(df_habitat)<-c("logLik","AICc", "Delta_AICc", "AICcWt", "K_rates")
write.table(df_habitat, "ASE/df_habitat.tsv", sep ="\t")
saveRDS(eval(parse(text = rownames(df_habitat)[which.min(df_habitat$AICc)])), "ASE/ht_corHMM.rds")