Msprime: Test growth rates in multiple merger coalescents.

[jeromekelleher]

1059 added support for growth rates in multiple merger coalescents, but it's untested. We should either add some tests to check it works or insist that growth_rate is always 0.

cc @JereKoskela, @eldonb, @TPPSellinger

See also #985

[TPPSellinger]

Sorry for the late answer. I could build simple tests to check for population size variation in multiple merger. Did @JereKoskela or @eldonb already start working on this issue ?

PS: I might not be very efficient in the next few weeks as I have to submit my thesis by end of summer.

[jeromekelleher]

I don't think anyone is working on this @TPPSellinger, we'd be very grateful for any help.

[JereKoskela]

I haven't made any progress on implementation thus far, but did realise that we know the expected merger times for a sample of size 2, which gives us something against which to compare. My intention was to implement that at some point this week.

How about I implement that in the next few days, and you can take over with further tests if they're needed once your thesis is handed in @TPPSellinger?

[TPPSellinger]

Okay! Let me know when you do the PR, I am curious on how to test this. When you mean sample 2, you mean 2 diploid or 2 haploid individual ?

[JereKoskela]

2 haploids. When there are just 2 lineages and no recombination, then every coalescent is Kingman, so we can get the exact expected coalescence time as a function of demographic parameters.

[JereKoskela]

1082 established that the growth rate being simulated matches the expected merger time. Do we need anything else? I guess changes in growth rate are the only other related piece of functionality, but those should work identically to the Hudson case, in that the same code is used to implement the change in both multiple merger and Hudson models. So I think we should be safe.

[TPPSellinger]

Should we put a note/warning for when beta tends to 2 ? Computation time could be very high for beta value higher than 1.99.

[jeromekelleher]

I think we can close the issue @JereKoskela - please reopen if you disagree.

[JereKoskela]

Should we put a note/warning for when beta tends to 2 ? Computation time could be very high for beta value higher than 1.99.

Is this something you've noticed in practice @TPPSellinger? If so, could you open a separate issue with a minimal example? The way the rejection sampling is implemented means that performance shouldn't deteriorate as the Beta-coalescent parameter tends to 2, and in fact it can even speed up.

[TPPSellinger]

Computation time was high in the previous implementation when beta was tending to 2. I didn't check with the new implementation but I will check it tomorrow and let you know if anything strange happens.

[TPPSellinger]

I did a quick check. msprime is very fast even for beta value tending to 2. Only one question to @JereKoskela. I have more SNPs for smaller beta value. I am afraid for some beta value very very close to 2, there might be no observed SNPs. Should we suggest users to use Kingman for beta value higher than 1.9999 (Although the beta coalescent does not tend to Kingman when beta tends to 2) ?

[JereKoskela]

I did a quick check. msprime is very fast even for beta value tending to 2. Only one question to @JereKoskela. I have more SNPs for smaller beta value. I am afraid for some beta value very very close to 2, there might be no observed SNPs. Should we suggest users to use Kingman for beta value higher than 1.9999 (Although the beta coalescent does not tend to Kingman when beta tends to 2) ?

That's a consequence of the way the Schweinsberg model scales again. As β -> 2, the timescale becomes proportional to N, as in the Kingman coalescent, but the constant of proportionality decays to zero because B(2-β, β) in the denominator diverges to infinity (for the form of the constant, see the docs here: https://msprime.readthedocs.io/en/latest/api.html#multiple-merger-coalescents). That means that you need a very large population size to counter the effect of the constant. We should probably add a warning to that effect in the docs, thanks for pointing it out @TPPSellinger.

[jeromekelleher]

Should we open an issue to track this @JereKoskela?

[JereKoskela]

Should we open an issue to track this @JereKoskela?

I don't think that's necessary; it's a minor tweak to the docs. I realise discussing it in a closed, unrelated issue hasn't been ideal, but I'll open a PR to fix it this afternoon and that should settle the matter.

[jeromekelleher]

Great, thanks @JereKoskela. Opening issues is a good way to make sure we don't forget about things, but if you're following it up straight away then there's no need.