Msprime: let's talk about time scaling!

[andrewkern]

raising this issue because i'm trying to get the time scaling from the spiffy new sweep implementation to match the internal time scaling in msprime, and well, it's a bit tricky. I see that @fbaumdicker also ran in to this on #1106, so figured why not hash it out in a dedicated issue.

i guess the big question is this bit here:

https://github.com/tskit-dev/msprime/blob/dab983c4bebb188d10ce6ff2b2f07b0636d25fc0/lib/msprime.c#L3639-L3644

so we are scaling time to 4N?

[fbaumdicker]

Yes, I had to rescale the gene conversion rates to match the theoretical results. However, I did not notice that this is a temporary change! In that case, if this changes back, we will have to adapt the rates back again and I guess the same would be true for a bunch of other cases. It also means that the results of people using msprime could depend on the version used.

[jeromekelleher]

I had hoped to make the code a bit more understandable by taking a factor of two out of here, and making the coalescence rate n choose 2, rather than n(n - 1) as it is now.

@JereKoskela - you understand this much better than me. Is there a way we can make the rate calcs less confusing, or is what we have now about as good as it gets (and we just need to change the comments)?

[JereKoskela]

I had hoped to make the code a bit more understandable by taking a factor of two out of here, and making the coalescence rate n choose 2, rather than n(n - 1) as it is now.

@JereKoskela - you understand this much better than me. Is there a way we can make the rate calcs less confusing, or is what we have now about as good as it gets (and we just need to change the comments)?

I assume we still want the population size to be 4N haploids, right @jeromekelleher? If so, then the unscaled Wright-Fisher per-generation coalescence probability is 1/(4N) per pair, or (n choose 2)/(4N) in total when there are n lineages. Likewise, the recombination rate is just r per link per generation. Currently, the numbers used in the code are n*(n-1)/(4N) and r, and I'm not sure where one factor of 1/2 has gone in the former.

Is the intention to make pop->initial_size = 4N, and get rid of the factor of 4 that way? That would certainly make sense.

[jeromekelleher]

Thanks @JereKoskela - what I'm hoping to do is make it all a bit more obvious what's going on. Very much open to ideas on how we can do this!

[jeromekelleher]

But yes, 4N haploids is the basic assumption I think.

[andrewkern]

4N haploids? Don't we want either 2N haploids or N diploids?

[jeromekelleher]

:exploding_head: yes, surely it's 2N haploids.

[andrewkern]

btw both Hudson (1983) and Tajima (1983) use 2N haploids as the scaling so it would be fair to stick with that

[jeromekelleher]

Well, we have to stick with whatever it is that msprime does now, there's no way we can change the default behaviour. The question is documenting clearly what this is so that developers can understand it. I would also like to make the ploidy an option, like I said in #557. Sorting this mess out internally so we understand what's going on would be the ideal first step.

But, like I say, @JereKoskela seems to be the only one who actually understands what the code is doing! (For the record, my approach was "msprime does whatever ms does", which seemed reasonable at the time.)

[JereKoskela]

Ok, 2N haploids is the easiest picture in my opinion, since it matches nicely with the current behaviour.

I suggest we fix the convention

pop->initial_size = N,

use the common ancestor event rate

lambda = n*(n-1.0) / 2.0;

and then have msp_get_common_ancestor_waiting_time_from_rate() scale the merger time using population size 2 * pop->initial_size instead of the current 4 * pop->initial_size. That matches the current implementation, but moves a factor of 1/2 from the scaling to the rate, as @jeromekelleher described above. Then we can just replace the 2 in 2 * pop->initial_size by a variable called, say, ploidy, which can take value 1 or 2 to toggle between a haploid and diploid population.

If that sounds ok, I'd be happy to put it together an initial branch so we can see what the code actually looks like.

[andrewkern]

Ok, 2N haploids is the easiest picture in my opinion, since it matches nicely with the current behaviour.

I suggest we fix the convention

pop->initial_size = N,

use the common ancestor event rate

lambda = n*(n-1.0) / 2.0;

and then have msp_get_common_ancestor_waiting_time_from_rate() scale the merger time using population size 2 * pop->initial_size instead of the current 4 * pop->initial_size. That matches the current implementation, but moves a factor of 1/2 from the scaling to the rate, as @jeromekelleher described above.

I like this suggestion!

[jeromekelleher]

If that sounds ok, I'd be happy to put it together an initial branch so we can see what the code actually looks like.

Yes, this sounds great, that would be amazing @JereKoskela!

[jeromekelleher]

See #1052 for original issue flagging the problems.

[jeromekelleher]

I'm going to close this as resolved now that we have a ploidy variable in the C code and the calculations have been clarified. Please reopen if you think we need to discuss more though.

[andrewkern]

was code ever merged to reflect the changes that @JereKoskela suggested?

[JereKoskela]

was code ever merged to reflect the changes that @JereKoskela suggested?

Yes, it's all there now @andrewkern. I think #1110 is the right place to look.

[andrewkern]

awesome thanks.