Simulating complete phylogenies to test dating with constraints
Here we replicate simulation setting 1 from Zhang et al., “Total-Evidence Dating under the Fossilized Birth–Death Process”, Sys. Biol. 2016 (https://www.ncbi.nlm.nih.gov/pubmed/26493827).
setwd("~/data/datingWithConstraintsSimulations/")
#install.packages("TreeSim")
library(TreeSim)
#install.packages("BioGeoBEARS")
library(BioGeoBEARS)Now we perform a single simulation as in Zhang et al., 2016.
n<-100
lambda <- 0.3
mu <- 0.2
frac <-0.5
numbsim<-2
# Each extant species is included in final tree with probability frac
# (the tree has n species AFTER sampling):
tre<-sim.bd.taxa(n,numbsim,lambda,mu,frac,complete=TRUE,stochsampling=TRUE)
is.phylo(tre[[1]])[1] TRUEPlotting the complete phylogeny
plot(tre[[1]], show.tip.label=FALSE)
Plotting only extant species
extant<-drop.extinct(tre[[1]])
plot(extant, show.tip.label=FALSE)
Plotting only extinct species
extantLeaves <-setdiff( tre[[1]]$tip.label, is.extinct(tre[[1]]))
extinct<-drop.tip(tre[[1]], extantLeaves)
plot(extinct, show.tip.label=FALSE)
Getting fossil node ages
treProperties <- prt(tre[[1]])prt(t): root=602
[1] "prt(t): root=602\n"
node ord_ndname node_lvl node.type parent_br edge.length ancestor
1 1 1 27 tip 223 1.519101878 1113
2 2 2 31 tip 212 0.542249524 1162
3 3 3 28 tip 1098 0.357140978 1173
4 4 4 24 tip 827 1.732232014 1104
5 5 5 28 tip 191 2.281456752 1083
6 6 6 29 tip 590 0.709050061 1153
7 7 7 29 tip 913 0.490273573 1166
8 8 8 26 tip 660 1.008502658 1138
9 9 9 17 tip 45 6.272976210 931
10 10 10 26 tip 343 4.401445525 992
11 11 11 30 tip 987 3.380649923 1044
12 12 12 30 tip 725 0.457085748 1167
13 13 13 30 tip 592 0.452297561 1168
14 14 14 26 tip 1147 1.228578318 1128
15 15 15 20 tip 625 3.705995101 1025
16 16 16 33 tip 249 0.171280051 1183
17 17 17 28 tip 397 2.509901474 1070
18 18 18 26 tip 265 0.256861309 1177
19 19 19 24 tip 73 0.719083794 1152
20 20 20 27 tip 142 0.337844053 1174
21 21 21 29 tip 186 0.136633734 1186
22 22 22 23 tip 270 2.458336080 1073
23 23 23 26 tip 107 1.426805239 1117
24 24 24 24 tip 672 2.747510146 1063
25 25 25 27 tip 425 6.312253395 930
26 26 26 23 tip 823 1.623037257 1106
27 27 27 27 tip 659 0.082344117 1190
28 28 28 33 tip 407 1.279358835 1126
29 29 29 29 tip 423 2.006648051 1091
30 30 30 26 tip 600 0.449414185 1169
31 31 31 28 tip 183 1.564783243 1109
32 32 32 29 tip 187 0.136633734 1186
33 33 33 30 tip 1155 0.206264061 1179
34 34 34 25 tip 1134 4.431675729 991
35 35 35 26 tip 264 0.256861309 1177
36 36 36 25 tip 93 0.706286107 1154
37 37 37 19 tip 623 5.755383454 950
38 38 38 24 tip 147 3.789349096 1019
39 39 39 26 tip 77 0.050498393 1197
40 40 40 32 tip 419 0.045964439 1198
41 41 41 22 tip 606 10.894585478 804
42 42 42 25 tip 1041 1.184193615 1131
43 43 43 29 tip 259 2.317365447 1079
44 44 44 29 tip 566 0.123908185 1188
45 45 45 26 tip 138 1.047609991 1136
46 46 46 25 tip 943 3.009348011 1052
47 47 47 23 tip 271 2.458336080 1073
48 48 48 24 tip 110 2.201684093 1084
49 49 49 27 tip 222 1.519101878 1113
50 50 50 31 tip 442 0.638112218 1158
51 51 51 28 tip 154 0.501592400 1164
52 52 52 22 tip 96 3.814334533 1016
53 53 53 26 tip 328 1.451855759 1115
54 54 54 26 tip 279 1.765324282 1101
55 55 55 26 tip 106 1.426805239 1117
56 56 56 29 tip 710 0.067380003 1194
57 57 57 29 tip 953 0.145743433 1184
58 58 58 33 tip 414 1.533590009 1111
59 59 59 24 tip 281 4.012598972 1006
60 60 60 25 tip 1009 0.190978522 1180
61 61 61 28 tip 657 0.080618202 1191
62 62 62 23 tip 824 1.623037257 1106
63 63 63 24 tip 661 1.597009099 1108
64 64 64 26 tip 599 0.449414185 1169
65 65 65 28 tip 340 0.641848598 1157
66 66 66 25 tip 92 0.706286107 1154
67 67 67 26 tip 143 1.103440321 1134
68 68 68 23 tip 66 1.221890289 1129
69 69 69 26 tip 76 0.050498393 1197
70 70 70 32 tip 420 0.045964439 1198
71 71 71 27 tip 1105 1.421228760 1118
72 72 72 33 tip 543 0.022570324 1200
73 73 73 28 tip 155 0.501592400 1164
74 74 74 22 tip 456 1.851211062 1098
75 75 75 31 tip 421 1.876475893 1097
76 76 76 31 tip 974 3.297862412 1046
77 77 77 25 tip 597 3.530366770 1035
78 78 78 32 tip 250 0.692375059 1155
79 79 79 22 tip 455 1.851211062 1098
80 80 80 29 tip 915 0.790534784 1142
81 81 81 24 tip 828 1.732232014 1104
82 82 82 30 tip 403 3.932324866 1013
83 83 83 28 tip 379 2.085511189 1090
daughter_nds node_ht time_bp fossils label
1 45.795655 0.0000000 FALSE t316
2 45.795655 0.0000000 FALSE t471
3 45.795655 0.0000000 FALSE t294
4 45.795655 0.0000000 FALSE t329
5 45.795655 0.0000000 FALSE t126
6 45.795655 0.0000000 FALSE t571
7 45.795655 0.0000000 FALSE t204
8 45.795655 0.0000000 FALSE t237
9 45.795655 0.0000000 FALSE t418
10 45.795655 0.0000000 FALSE t125
11 45.795655 0.0000000 FALSE t12
12 45.795655 0.0000000 FALSE t266
13 45.795655 0.0000000 FALSE t389
14 45.795655 0.0000000 FALSE t185
15 45.795655 0.0000000 FALSE t224
16 45.795655 0.0000000 FALSE t402
17 45.795655 0.0000000 FALSE t353
18 45.795655 0.0000000 FALSE t325
19 45.795655 0.0000000 FALSE t450
20 45.795655 0.0000000 FALSE t470
21 45.795655 0.0000000 FALSE t235
22 45.795655 0.0000000 FALSE t249
23 45.795655 0.0000000 FALSE t392
24 45.795655 0.0000000 FALSE t349
25 45.795655 0.0000000 FALSE t458
26 45.795655 0.0000000 FALSE t371
27 45.795655 0.0000000 FALSE t22
28 45.795655 0.0000000 FALSE t152
29 45.795655 0.0000000 FALSE t420
30 45.795655 0.0000000 FALSE t307
31 45.795655 0.0000000 FALSE t137
32 45.795655 0.0000000 FALSE t186
33 45.795655 0.0000000 FALSE t108
34 45.795655 0.0000000 FALSE t344
35 45.795655 0.0000000 FALSE t410
36 45.795655 0.0000000 FALSE t306
37 45.795655 0.0000000 FALSE t8
38 45.795655 0.0000000 FALSE t457
39 45.795655 0.0000000 FALSE t112
40 45.795655 0.0000000 FALSE t136
41 45.795655 0.0000000 FALSE t360
42 45.795655 0.0000000 FALSE t254
43 45.795655 0.0000000 FALSE t231
44 45.795655 0.0000000 FALSE t305
45 45.795655 0.0000000 FALSE t408
46 45.795655 0.0000000 FALSE t566
47 45.795655 0.0000000 FALSE t478
48 45.795655 0.0000000 FALSE t497
49 45.795655 0.0000000 FALSE t587
50 45.795655 0.0000000 FALSE t144
51 45.795655 0.0000000 FALSE t354
52 45.795655 0.0000000 FALSE t81
53 45.795655 0.0000000 FALSE t273
54 45.795655 0.0000000 FALSE t568
55 45.795655 0.0000000 FALSE t430
56 45.795655 0.0000000 FALSE t393
57 45.795655 0.0000000 FALSE t59
58 45.795655 0.0000000 FALSE t381
59 45.795655 0.0000000 FALSE t422
60 45.795655 0.0000000 FALSE t403
61 45.795655 0.0000000 FALSE t352
62 45.795655 0.0000000 FALSE t193
63 45.795655 0.0000000 FALSE t584
64 45.795655 0.0000000 FALSE t74
65 45.795655 0.0000000 FALSE t278
66 45.795655 0.0000000 FALSE t374
67 45.795655 0.0000000 FALSE t593
68 45.795655 0.0000000 FALSE t34
69 45.795655 0.0000000 FALSE t42
70 45.795655 0.0000000 FALSE t85
71 45.795655 0.0000000 FALSE t88
72 45.795655 0.0000000 FALSE t330
73 45.795655 0.0000000 FALSE t251
74 45.795655 0.0000000 FALSE t397
75 45.795655 0.0000000 FALSE t54
76 45.795655 0.0000000 FALSE t445
77 45.795655 0.0000000 FALSE t347
78 45.795655 0.0000000 FALSE t172
79 45.795655 0.0000000 FALSE t301
80 45.795655 0.0000000 FALSE t508
81 45.795655 0.0000000 FALSE t303
82 45.795655 0.0000000 FALSE t105
83 45.795655 0.0000000 FALSE t548
[ reached getOption("max.print") -- omitted 1118 rows ]fossilAges <- treProperties[which(treProperties$node.type=="tip" & treProperties$time_bp >0),]
#print(fossilAges)
summary(fossilAges$time_bp) Min. 1st Qu. Median Mean 3rd Qu. Max.
0.01606 3.04353 7.02925 10.60962 14.90311 42.29132 summary(fossilAges$node_ht) Min. 1st Qu. Median Mean 3rd Qu. Max.
3.504 30.893 38.766 35.186 42.752 45.780 Getting extant node ages
extantAges <- treProperties[which(treProperties$node.type=="tip" & treProperties$time_bp ==0),]
#print(extantAges)
summary(fossilAges$time_bp) Min. 1st Qu. Median Mean 3rd Qu. Max.
0.01606 3.04353 7.02925 10.60962 14.90311 42.29132 summary(extantAges$node_ht) Min. 1st Qu. Median Mean 3rd Qu. Max.
45.8 45.8 45.8 45.8 45.8 45.8 LS0tCnRpdGxlOiAiU2ltdWxhdGluZyBjb21wbGV0ZSBwaHlsb2dlbmllcyB0byB0ZXN0IGRhdGluZyB3aXRoIGNvbnN0cmFpbnRzIgpvdXRwdXQ6IGh0bWxfbm90ZWJvb2sKLS0tCgpIZXJlIHdlIHJlcGxpY2F0ZSBzaW11bGF0aW9uIHNldHRpbmcgMSBmcm9tIFpoYW5nIGV0IGFsLiwgIlRvdGFsLUV2aWRlbmNlIERhdGluZyB1bmRlciB0aGUgRm9zc2lsaXplZCBCaXJ0aOKAk0RlYXRoIFByb2Nlc3MiLCBTeXMuIEJpb2wuIDIwMTYgKGh0dHBzOi8vd3d3Lm5jYmkubmxtLm5paC5nb3YvcHVibWVkLzI2NDkzODI3KS4KCmBgYHtyfQpzZXR3ZCgifi9kYXRhL2RhdGluZ1dpdGhDb25zdHJhaW50c1NpbXVsYXRpb25zLyIpCiNpbnN0YWxsLnBhY2thZ2VzKCJUcmVlU2ltIikKbGlicmFyeShUcmVlU2ltKQojaW5zdGFsbC5wYWNrYWdlcygiQmlvR2VvQkVBUlMiKQpsaWJyYXJ5KEJpb0dlb0JFQVJTKQoKYGBgCgoKYGBge3IsIGluY2x1ZGU9RkFMU0V9Cm48LTEwCmxhbWJkYSA8LSAyLjAKbXUgPC0gMC41CmZyYWMgPC0wLjYKbnVtYnNpbTwtMgoKIwojIFNpbXVsYXRpbmcgbnVtYnNpbSB0cmVlcyB3aXRoIG4gc3BlY2llcyB1bmRlciBhIGJpcnRoLWRlYXRoIHByb2Nlc3Mgd2l0aAojIHNwZWNpYXRpb24gcmF0ZSBsYW1iZGEgYW5kIGV4dGluY3Rpb24gcmF0ZSBtdToKc2ltLmJkLnRheGEobixudW1ic2ltLGxhbWJkYSxtdSkKIyBFYWNoIGV4dGFudCBzcGVjaWVzIGlzIGluY2x1ZGVkIGluIGZpbmFsIHRyZWUgd2l0aCBwcm9iYWJpbGl0eSBmcmFjCiMgKHRoZSB0cmVlIGhhcyBuIHNwZWNpZXMgQUZURVIgc2FtcGxpbmcpOgpzaW0uYmQudGF4YShuLG51bWJzaW0sbGFtYmRhLG11LGZyYWMsY29tcGxldGU9RkFMU0Usc3RvY2hzYW1wbGluZz1UUlVFKQojIEEgZnJhY3Rpb24gZnJhYyBvZiB0aGUgZXh0YW50IHNwZWNpZXMgaXMgaW5jbHVkZWQgaW50byB0aGUgZmluYWwgdHJlZQojICh0aGUgdHJlZSBoYXMgbiBzcGVjaWVzIEFGVEVSIHNhbXBsaW5nKToKc2ltLmJkLnRheGEobixudW1ic2ltLGxhbWJkYSxtdSxmcmFjLGNvbXBsZXRlPUZBTFNFLHN0b2Noc2FtcGxpbmc9RkFMU0UpCmBgYAoKIyBOb3cgd2UgcGVyZm9ybSBhIHNpbmdsZSBzaW11bGF0aW9uIGFzIGluIFpoYW5nIGV0IGFsLiwgMjAxNi4KCgpgYGB7cn0KCm48LTEwMApsYW1iZGEgPC0gMC4zCm11IDwtIDAuMgpmcmFjIDwtMC41Cm51bWJzaW08LTIKCiMgRWFjaCBleHRhbnQgc3BlY2llcyBpcyBpbmNsdWRlZCBpbiBmaW5hbCB0cmVlIHdpdGggcHJvYmFiaWxpdHkgZnJhYwojICh0aGUgdHJlZSBoYXMgbiBzcGVjaWVzIEFGVEVSIHNhbXBsaW5nKToKdHJlPC1zaW0uYmQudGF4YShuLG51bWJzaW0sbGFtYmRhLG11LGZyYWMsY29tcGxldGU9VFJVRSxzdG9jaHNhbXBsaW5nPVRSVUUpCmlzLnBoeWxvKHRyZVtbMV1dKQoKYGBgCiMjIFBsb3R0aW5nIHRoZSBjb21wbGV0ZSBwaHlsb2dlbnkKCmBgYHtyfQpwbG90KHRyZVtbMV1dLCBzaG93LnRpcC5sYWJlbD1GQUxTRSkKYGBgCgojIyBQbG90dGluZyBvbmx5IGV4dGFudCBzcGVjaWVzCgpgYGB7cn0KZXh0YW50PC1kcm9wLmV4dGluY3QodHJlW1sxXV0pCnBsb3QoZXh0YW50LCBzaG93LnRpcC5sYWJlbD1GQUxTRSkKYGBgCgojIyBQbG90dGluZyBvbmx5IGV4dGluY3Qgc3BlY2llcwoKCmBgYHtyfQpleHRhbnRMZWF2ZXMgPC1zZXRkaWZmKCB0cmVbWzFdXSR0aXAubGFiZWwsIGlzLmV4dGluY3QodHJlW1sxXV0pKQpleHRpbmN0PC1kcm9wLnRpcCh0cmVbWzFdXSwgZXh0YW50TGVhdmVzKQpwbG90KGV4dGluY3QsIHNob3cudGlwLmxhYmVsPUZBTFNFKQoKYGBgCgojIyBHZXR0aW5nIGZvc3NpbCBub2RlIGFnZXMKYGBge3J9CnRyZVByb3BlcnRpZXMgPC0gcHJ0KHRyZVtbMV1dKQpmb3NzaWxBZ2VzIDwtIHRyZVByb3BlcnRpZXNbd2hpY2godHJlUHJvcGVydGllcyRub2RlLnR5cGU9PSJ0aXAiICYgdHJlUHJvcGVydGllcyR0aW1lX2JwID4wKSxdCiNwcmludChmb3NzaWxBZ2VzKQpzdW1tYXJ5KGZvc3NpbEFnZXMkdGltZV9icCkKc3VtbWFyeShmb3NzaWxBZ2VzJG5vZGVfaHQpCgpgYGAKCgojIyBHZXR0aW5nIGV4dGFudCBub2RlIGFnZXMKYGBge3J9CmV4dGFudEFnZXMgPC0gdHJlUHJvcGVydGllc1t3aGljaCh0cmVQcm9wZXJ0aWVzJG5vZGUudHlwZT09InRpcCIgJiB0cmVQcm9wZXJ0aWVzJHRpbWVfYnAgPT0wKSxdCiNwcmludChleHRhbnRBZ2VzKQpzdW1tYXJ5KGZvc3NpbEFnZXMkdGltZV9icCkKc3VtbWFyeShleHRhbnRBZ2VzJG5vZGVfaHQpCgoKYGBgCg==