-
Notifications
You must be signed in to change notification settings - Fork 0
/
PosteriorPredictive_TreeSummary.Rev
93 lines (66 loc) · 3.88 KB
/
PosteriorPredictive_TreeSummary.Rev
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
################################################################################
#
# RevBayes Example: Bayesian model testing using posterior predictive simulation
# This file calculates the Tree Summary Statistics
#
# authors: Lyndon M. Coghill, Sebastian Hoehna and Jeremy M. Brown
#
################################################################################
outfileName = "results_" + model_name + "/simulated_inference_" + analysis_name + ".csv"
write(file=outfileName, "simID", "mean_rf", "quantile25", "quantile50", "quantile75", "quantile99", "quantile999", "mean_tl", "var_tl", "entropy", sep=",", append=FALSE)
write(file=outfileName, "\n", sep="\t", append=TRUE)
#num_post_sims = 10
################### calculate the pps stats here #########################
## Iterate through all of the posterior tree files from the simulation analyses
for ( i in 1:num_post_sims) {
inFileName = "output_" + model_name + "/posterior_predictive_sim_" + i + "/" + analysis_name + "_posterior.trees"
## read in the trace
sim_tree_trace = readTreeTrace(inFileName,treetype="non-clock")
## Calculate the pairwise RF distances between all trees in a single posterior
rf_dists <- sim_tree_trace.computePairwiseRFDistances(credibleTreeSetSize=1.0,verbose=FALSE)
## This collects the vector of tree lengths needed for mean tree length and tree length variance
tree_length <- sim_tree_trace.computeTreeLengths()
## This calculates the entropy statistic
entropy <- sim_tree_trace.computeEntropy(credibleTreeSetSize=1.0,verbose=FALSE)
## calculate the stuff we care about for a single pps posterior
mean_rf <- mean(rf_dists)
quantile25 = quantile(rf_dists, 0.25)
quantile50 = quantile(rf_dists, 0.50)
quantile75 = quantile(rf_dists, 0.75)
quantile99 = quantile(rf_dists, 0.99)
quantile999 = quantile(rf_dists, 0.999)
mean_tl = mean(tree_length)
var_tl = var(tree_length)
## write it to a file
write(file=outfileName, i, mean_rf, quantile25, quantile50, quantile75, quantile99, quantile999, mean_tl, var_tl, entropy, sep=",", append=TRUE)
write(file=outfileName, "\n", sep=",", append=TRUE)
clear(rf_dists)
}
################### end of pps calculations ####################################
################### calculate the empirical stats here ###################
## set the empirical output file
outfileName = "results_" + model_name + "/empirical_inference_" + analysis_name + ".csv"
write(file=outfileName, "mean_rf", "quantile25", "quantile50", "quantile75", "quantile99", "quantile999", "mean_tl", "var_tl", "entropy", sep=",", append=FALSE)
write(file=outfileName, "\n", sep=",", append=TRUE)
## find the empirical file
empFileName = "output_" + model_name + "/" + analysis_name + "_posterior.trees"
## read in the trace
emp_tree_trace = readTreeTrace(empFileName,treetype="non-clock")
## Calculate the pairwise RF distances between all trees in a single posterior
rf_dists <- emp_tree_trace.computePairwiseRFDistances(credibleTreeSetSize=1.0,verbose=FALSE)
## This collects the vector of tree lengths needed for mean tree length and tree length variance
tree_length <- emp_tree_trace.computeTreeLengths()
## This calculates the entropy statistic
entropy <- emp_tree_trace.computeEntropy(credibleTreeSetSize=1.0,verbose=FALSE)
## calculate the stuff we care about for a single posterior
mean_rf = mean(rf_dists)
quantile25 = quantile(rf_dists, 0.25)
quantile50 = quantile(rf_dists, 0.50)
quantile75 = quantile(rf_dists, 0.75)
quantile99 = quantile(rf_dists, 0.99)
quantile999 = quantile(rf_dists, 0.999)
mean_tl = mean(tree_length)
var_tl = var(tree_length)
write(file=outfileName, mean_rf, quantile25, quantile50, quantile75, quantile99, quantile999, mean_tl, var_tl, entropy, sep=",", append=TRUE)
write(file=outfileName, "\n", sep=",", append=TRUE)
################### end of empirical calculations ##############################