Run the program using the following strom.conf file:
datafile = rbcl10.nex
treefile = rbcl10.tre
ncateg = default:4
pinvar = default:0.1
seed = 13579
usedata = no
nchains = 4
burnin = 10000
niter = 100000
printfreq = 10000
samplefreq = 10
Note that no data will be used for this run, so the MCMC analysis will explore the prior. You’ll remember that we did this in the previous step, but now there is a convenient program option to allow easy switching between using data and ignoring data.
After the run has finished, load the resulting params.txt file into the program Tracer to summarize the results. You should find that the means and variances of all parameters are the same as in the previous step: adding multiple chains does not change the distribution being explored by the cold chain, it simply makes moving around that landscape more efficient.
Note that chain summaries are provided for all four chains.
Chain 0 (power 1.00000)
Updater Tuning Param. Accept % No. Updates
State Frequencies 1000.000 100.0 4203
Exchangeabilities 1000.000 100.0 4255
Gamma Rate Variance 10.568 29.6 4250
Proportion of Invariable Sites 1.000 100.0 4209
Tree Topol. and Edge Prop. 1000.000 100.0 78939
Tree Length 13.137 29.7 4144
Chain 1 (power 0.66667)
Updater Tuning Param. Accept % No. Updates
State Frequencies 1000.000 100.0 4173
Exchangeabilities 1000.000 100.0 4222
Gamma Rate Variance 13.741 30.3 4229
Proportion of Invariable Sites 1.000 100.0 4122
Tree Topol. and Edge Prop. 1000.000 100.0 79047
Tree Length 13.716 28.1 4207
Chain 2 (power 0.50000)
Updater Tuning Param. Accept % No. Updates
State Frequencies 1000.000 100.0 4165
Exchangeabilities 1000.000 100.0 4154
Gamma Rate Variance 10.863 30.6 4159
Proportion of Invariable Sites 1.000 100.0 4279
Tree Topol. and Edge Prop. 1000.000 100.0 79040
Tree Length 11.490 29.4 4203
Chain 3 (power 0.40000)
Updater Tuning Param. Accept % No. Updates
State Frequencies 1000.000 100.0 4204
Exchangeabilities 1000.000 100.0 4217
Gamma Rate Variance 13.100 30.6 4224
Proportion of Invariable Sites 1.000 100.0 4204
Tree Topol. and Edge Prop. 1000.000 100.0 79095
Tree Length 12.140 29.4 4056
The acceptance fractions are near 30% (which is what we specified as the target acceptance percentage in the Chain::createUpdaters function) for Gamma Rate Variance and Tree Length. For State Frequencies, Exchangeabilities, and Tree Topology and Edge Proportions, the tuning parameter has maxed out at a boldness of 1000 and the acceptance percentages are both 100%. For Proportion of Invariable Sites, the acceptance percentage is also 100% but the tuning parameter maxed out at 1.0 rather than 1000. Why is this? The parameters modified by these four updaters are the only parameters with truly flat priors, and, under a flat prior, all possible combinations have equal probability density. Because there are no downhill moves when updating state frequencies, exchangeabilities, tree topology, edge length proportions, or the proportion of invariable sites, all proposals are accepted. It is impossible to tune the samplers in this case because the proposal will be accepted regardless of how far away the proposed point is from the current point. You should not see this behavior when analyzing real data because sequence data contains information about these parameters that results in a posterior that is not perfectly flat. The fact that the Proportion of Invariable Sites tuning parameter maxes out at 1 rather than 1000 is due to the highlighted lines at the beginning of the PinvarUpdater::proposeNewState function:
inline void PinvarUpdater::proposeNewState() {
<span style="color:#0000ff"><strong>if (_lambda > 1.0)</strong></span>
<span style="color:#0000ff"><strong>_lambda = 1.0;</strong></span>
// Save copy of _curr_point in case revert is necessary.
_prev_point = getCurrentPoint();
// Propose new value using uniform window of width _lambda centered over _prev_point
double p = (_prev_point - _lambda/2.0) + _lambda*_lot->uniform();
if (p < 0.0)
p = -p;
else if (p > 1.0)
p = 1.0 - (p - 1.0);
_asrv->setPinvar(p);
_log_hastings_ratio = 0.0; //symmetric proposal
// This proposal invalidates all transition matrices and partials
_tree_manipulator->selectAllPartials();
_tree_manipulator->selectAllTMatrices();
}
Because the proposal method used for pinvar is a window of width _lambda, it is prudent to limit _lambda to the maximum possible window width, namely 1.0. Allowing _lambda to get larger than 1.0 risks having to reflect multiple times if the proposed new state lies outside the valid inteval 0.0 to 1.0.
The swap summary shows that more than 16000 swaps were attempted between each pair of chains, fewer swaps were successful between chains that differed substantially in their heating power. For example, 10012381/16377 = 76% of swaps were accepted between the two hottest chains whereas only 1004118/16876 = 24% were accepted between the cold chain and the hottest chain. This is explained by the fact that the hottest chain is often far away from a peak while the cold chain never strays far from a peak; the probability that the cold chain will accept a jump from the peak to the lowlands is small, so, despite the fact that the heated chain can accept the swap because it is uphill, fewer swaps succeed between the coldest and hottest chains than between chains visiting more similar landscapes. The reason that even 24% of such swaps succeed has a lot to do with the fact that the landscape is determined entirely by the prior in this case, so there is not a huge difference between the altitude of the peak compared to the lowlands.
Swap summary (upper triangle = no. attempted swaps; lower triangle = no. successful swaps):
0 1 2 3
----------------------------------------------------------------
0 --- 16520 16589 16824
1 11719 --- 16561 16908
2 8699 13130 --- 16598
3 6898 10922 13951 ---
----------------------------------------------------------------
Feel free to run the program again, this time setting usedata=yes (or leaving out the usedata program altogether, which will default to using the data).
You now have a program that can carry out a Bayesian phylogenetic analysis using either a GTR+I+G or a codon model, allowing data to be partitioned, and can use multiple chains during MCMC. With slight modification, your program can make use of the GPU version of BeagleLib when computing likelihoods, which greatly speeds up the codon model if GPUs are available.
In the next step in the tutorial, we will implement reversible-jump MCMC to make it possible to explore trees containing polytomies.