[Help wanted] Reversible jump between species partitions in multispecies coalescent

[thauffe]

Hi,

I am trying to use reversible jump MCMC to switch between species partitions in an multispecies
MultispeciesCoalescent.zip
coalescent model.
The idea is similar to the program bpp and the STACEY module of BEAST: We have a species tree in which the gene trees are embedded and jump among species trees where nodes are collapsed, and thus the number of species varies. This helps with species delimitation where collapsed species are probably no good species.

Unfortunately, I cannot get any of the two reversible jump approaches from the script Model averaging of substitution models tutorial working.

The attached zip file contains scripts for three multispecies coalescent analyses.

1. mcmc_msc.Rev is recreating the *BEAST example to check if I use the multispecies coalescent (more or less) correct. No reversible jump here.
2. mcmc_msc_rj_v1.Rev tries to sample a model indicator variable (line 83). This is then used to select one of two taxon objects (with 7 or 8 species) and to generate the species tree (line 110). The problem is that the model indicator is somehow not a node in the model and I need to add to the model (line 198). If I include it in the model, it will be sampled but this does not affect the species tree. This always contains the same set of species (checking the posterior distribution of output/trees/msc_rj_v1.trees in e.g. FigTree).
3. mcmc_msc_rj_v1.Rev uses the two taxon maps to create two species trees and then tries to jump between the trees (line 118). Here the problem is that always the same species tree is sampled.

Any idea on how to implement reversible jump moves in the context of multispecies coalescent will be highly appreciated!

Thanks,
Torsten

[davidcerny]

This might not be easy to solve, unfortunately...

I'm not quite sure why your v1 script doesn't work as intended, but the fact that model_indicator isn't properly connected to the model graph unless you explicitly include it in the model() function might point to some underlying problems. It might be worth opening an issue for this, especially if it can be reproduced using the simpler example from the tutorial. I remember running into the same problem before, too.

Your v2 script probably doesn't work because of this part:

#moves.append( mvGibbsMixtureAllocation(species_tree, weight = 5) )
moves.append( mvTreeScale(species_tree, weight = 5) )

Here, it really is the mvGibbsMixtureAllocation that you want, but I understand why you commented it out – it doesn't accept objects of type TimeTree. Extending it so that it does might not be trivial. It's certainly not as easy as adding addType( new Move_GibbsMixtureAllocation<Tree>( ) ); to RbRegister_Move – I tried that, and got a segfault...

[bredelings]

I'll try and look into this on Thursday.

[thauffe]

Thanks so much for your feedback and for checking this!

[thauffe]

Sorry to bother, @bredelings but did you had time to look into this issue?

Thanks a lot!

[thauffe]

I made a bit of progress in my aim to jump among species partitions after coming across issue #533.

I create a species tree per species partition (lines 127-140 in the attached script) and then select one of them by an indicator variable with a user defined function.

function fnSelectTreeByIndicator(tree_list, indicator) {
    return tree_list[indicator]
}

This runs and jumps among three species partitions. Awesome!

However, this only works when I initialize the model_indicator with the 3rd species partition model_indicator.setValue(3) (line 102). There is a core dump when starting with partition 1 or 2 (although they are sampled when we start with partition 3).
It is odd thought that when I limit the potential partitions to 1 and 2 (i.e. n_species_partitions <- 2 in line 18), I can initialize with partition 2. Moreover, for these two partitions an likelihood is calculated for the first MCMC generation and the core dump happens just after Iter 0 (just comment out mymcmc.burnin).

I also encounter a core dump with my 4th species partition no matter with which partition I initialize. You can replicate this by changing in line 18 from 3 to 4 n_species_partitions <- 4

Any idea what is causing the core dumps or how I can find out?

As always,
thanks!

MultispeciesCoalescentRJ.zip

[thauffe]

Using the debug mode rb --setOption logMCMC=3 scripts/mcmc_msc_rj.Rev it turned out that some out-of-bounds indexing of the species partition caused a core dump.
This has been fixed thanks to Benjamin. Discussion PR# 709

However, I still get a core dump due to a segmentation fault. Unfortunately, there is no specific error message printed on screen.

The data and script are attached.
Interestingly, the core dump does not happen when changing line 14 toinitial_partition <- 3. However, all other combinations of n_species_partitions (in the range [1, 5]) and initial_partition (range [1, n_species_partitions]) do not work.

The script is a modification of the one from the previous post and allows easier selection of the species partitions by changing the order and adding/removing numbers in order_species_partitions.

As always, any input is highly appreciated!

MultispeciesCoalescentRJ.zip

[bredelings]

I'll take a look at this tomorrow.

[bredelings]

I made the coalescent code a bit more robust, and now it returns this error:

generation = 0    proposal = 14/421    RandomGeometricWalk(model_indicator)
    BEFORE:   model_indicator:  -1.0986122887
    BEFORE:   gene_tree[1]:  134.91608597
    BEFORE:   gene_tree[2]:  136.4970039
    BEFORE:   gene_tree[3]:  141.26630882

   Error:	No tree node for species 'Thomomys_umbrinus'
   Error:	Problem processing line 272 in file "scripts/mcmc_msc_rj.Rev"

[thauffe]

Thanks a lot for checking this!

Really odd that the species has no node. You mean node in the RevBayes DAG, not a node in the phylogeny, right?

Because the uploaded script sets as initial species phylogeny one that contains Thomomys umbrinus.
When I set the initial_partition to 3, all three species partitions are sampled (see tracer screenshot with the frequency of sampled partitions). The species trees of partition 2 and 3 do not contain Thomomys umbrinus because it's individuals are collapsed into another nominal species.

[thauffe]

@bredelings suggested in #721 to check my coalescent script for bugs:

Also, #608 is not fixed yet. Fixing x[i] is one more step towards the fix -- possibly there is a bug in the #608 script, and possibly there are more revbayes bugs.

I created a very simple script where we jump among three phylogenies with different sets of taxa. There is no coalescence process involved because we have one DNA alignment for each of these species partitions.

The script works when I initialize with the first set of species and the model indicator jumps from 1->2->3.

Iter        |      Posterior   |     Likelihood   |          Prior   |   model_indi..   |    elapsed   |        ETA   |
-----------------------------------------------------------------------------------------------------------------------
0           |       -1523.43   |       -1590.01   |        66.5771   |              1   |   00:00:00   |   --:--:--   |
50          |       -993.606   |       -1036.78   |         43.178   |              2   |   00:00:00   |   --:--:--   |
100         |       -992.943   |       -1038.84   |        45.8934   |              2   |   00:00:01   |   00:00:04   |
150         |       -990.818   |       -1037.79   |         46.975   |              2   |   00:00:01   |   00:00:02   |
200         |       -984.703   |       -1037.08   |        52.3741   |              3   |   00:00:01   |   00:00:01   |
250         |       -993.525   |       -1039.94   |        46.4111   |              3   |   00:00:01   |   00:00:01   |
300         |        -989.61   |       -1038.21   |        48.5997   |              3   |   00:00:02   |   00:00:01   |
350         |       -999.019   |       -1042.82   |        43.7992   |              3   |   00:00:02   |   00:00:00   |
400         |       -992.232   |       -1040.53   |         48.299   |              3   |   00:00:02   |   00:00:00   |
450         |       -999.065   |       -1037.67   |        38.6062   |              3   |   00:00:02   |   00:00:00   |
500         |        -998.38   |       -1037.09   |        38.7141   |              3   |   00:00:02   |   00:00:00   |

However, when I initialize with initial_partition <- 2 (or 3 in line 9), I still encounter a segmentation fault.
rb --setOption debugMCMC=1 --setOption logMCMC=3 scripts/mcmc_rj.Rev

generation = 0    proposal = 6/248    RandomIntegerWalk(model_indicator)
    BEFORE:   model_indicator:  -1.0986122887
    BEFORE:   phySeq:  -1560.1392457

Segmentation fault (core dumped)

Interestingly, when running under prior with mymodel.ignoreAllData(), it always works irrespective of the initial partition.

I used the refactor-phyloctmc-caching branch #729, which should include all previous fixes for reversible jump.

Thanks for any input!
RJ_DNA.zip

[bredelings]

I took at look at this using the development branch, which has merged all the proposed fixes so far.

A few things:

• you no longer need the fnSelect* functions -- you can just write array[index] directly.
• dnPhyloCTMC is currently assuming that the number of tree nodes does not change. The SEGFAULT is coming from the fact that the number of nodes is increasing to more than the amount of memory initially allocated.
• clamping deterministic nodes doesn't currently work.

Clamping deterministic nodes seems especially problematic:

• In general, it doesn't make sense for the data to differ depending on the model. The whole idea in model comparison is that you are comparing equivalent explanations for the SAME data.
• This script is looking at sequence evolution on a species tree. Why wouldn't you look at sequence evolution on a gene tree, and then ask which species tree best fits the gene tree?

Is there a reason why you have coded the script in this way?

[thauffe]

To explore further why we can jump among trees when ignoring the DNA sequence data, I experimented with a single continuous trait instead of the DNA alignment.

When initializing with initial_partition <- 1, I receive an error

generation = 0    proposal = 10/227    NodeTimeSlideUniform(psi[2])
    BEFORE:   psi[2]:  22.099200935
    BEFORE:   BM:  10.173187741

   Error:	Executing NodeTimeSlideUniform(psi[2]): PDFs not up-to-date before proposal!
       BM: 10.173187741 != 13.792849549    diff = -3.6196618074
   
   Error:	Problem processing line 171 in file "scripts/mcmc_rj.Rev"

Using 2 or 3 as initial_partition results again in a segmentation fault.
RJ_Morpho.zip