Stdpopsim: Choose sane defaults for SLiM engine

[grahamgower]

The SLiM engine supports enabling/disabling recapitation. Time can be scaled, or not. At the moment, the defaults are to do recaptitation and scale time by a factor of 10. Are these defaults appropriate?

Also see #420, #444.

It would be great to get some wider feedback on this and #420 from @popsim-consortium/all

[petrelharp]

I strongly vote for a default scaling factor of 1.0. One of the biggest points of SLiM is doing selection, and rescaling does not give equivalent results for selection. We don't want to confuse people here.

I do think doing recapitation by default is a good idea!

[grahamgower]

Just so we're all informed, it takes 3 minutes to run a 250 kb OutOfAfrica_3G09 simulation with a scaling factor of 1.0, and under 2 seconds with a scaling factor of 10.0.

$ time stdpopsim -e slim --slim-scaling-factor 10 HomSap -o /dev/null -l 0.001 -q -d OutOfAfrica_3G09 50 50 50

real    0m1,835s
user    0m1,659s
sys     0m0,174s
$ time stdpopsim -e slim --slim-scaling-factor 1 HomSap -o /dev/null -l 0.001 -q -d OutOfAfrica_3G09 50 50 50

real    2m58,202s
user    2m56,219s
sys     0m1,821s

[grahamgower]

@petrelharp Can you elaborate on the differences with regard to selection? Is this all kinds of selection (positive, negative, single locus, QTL, etc)? What kind of non-equivalence? Do you have any references you can point me at?

[RyanGutenkunst]

I strongly vote that the default be 1.0, but that it be made very easily changeable and that extensively documented. The appropriate scaling is going to vary dramatically depending on organism and computational resources. For example, humans can probably be done scaling by 1.0, but nothing at scale could be done with Drosophila without major scaling.

In our experience, even if scaling isn't perfect, many many useful simulation results do scale appropriately, particularly for modest sample sizes. For example, I find it pretty embarrassing that we failed at doing neutral genome-scale simulations for Drosophila in the paper with msprime. In our experience, SLiM simulations at a scaling of 1000 would complete quite quickly, and at the least the SFS-based analyses would give identical results to simulations with less aggressive scaling. We find similar results for models with background selection, at least with the summaries we're looking at.

[petrelharp]

Can you elaborate on the differences with regard to selection? Is this all kinds of selection (positive, negative, single locus, QTL, etc)?

Yep - there are no population genetic equivalences under rescaling when it comes to any kind of selection, so far as I am aware. The SLiM manual does a pretty good job of mentioning the various issues. For instance, rescaling by Q=10 is like replacing ten loci with small effect with a single locus with ten times the effect; and clearly many small effect loci are different than one large effect locus.

In our experience, even if scaling isn't perfect, many many useful simulation results do scale appropriately,

That's good to know! I'm in favor of putting it in, but not by default and making clear that it's not equivalent - too many people are confused already on this point.

[RyanGutenkunst]

Can you elaborate on the differences with regard to selection? Is this all kinds of selection (positive, negative, single locus, QTL, etc)?

Yep - there are no population genetic equivalences under rescaling when it comes to any kind of selection, so far as I am aware. The SLiM manual does a pretty good job of mentioning the various issues. For instance, rescaling by Q=10 is like replacing ten loci with small effect with a single locus with ten times the effect; and clearly many small effect loci are different than one large effect locus.

My intuition is shaped by the diffusion approach, where for single loci and pairs of loci the rescaling with population size is mathematically exact. So in the regimes covered by those approaches, there isn't any difference with rescaling. The semi-empirical question is how much rescaling matters for your particular parameter regime and statistical approach.

In our experience, even if scaling isn't perfect, many many useful simulation results do scale appropriately,

That's good to know! I'm in favor of putting it in, but not by default and making clear that it's not equivalent - too many people are confused already on this point.

I agree, but I would add language emphasizing that it is a very useful, and often computationally essential, approach. Maybe we could even suggest a rule of thumb or two.

For example, in my group, if we're running with rescaling factor F=f for computational reasons, we'll compare results extensively with rescaling factor F=2f or F=10f to judge how sensitive our results are to the rescaling factor. (This is easy because those more aggressive scalings are even faster). This doesn't protect against effects that might occur between F=1 and F=f but that somehow saturate before F=2f, so it's something like an uncontrolled extrapolation -- definitely not ideal, but necessary in some cases.

[jeromekelleher]

Sounds to me like the only debate here is "how much documentation do we need to write" and the everyone is saying "lots"!

[petrelharp]

My intuition is shaped by the diffusion approach, where for single loci and pairs of loci the rescaling with population size is mathematically exact.

Oh, good point - so, for pairs of loci under with s of order 1/N that are not near 0 or 1, the approximation should be good. Even for these, the boundary behavior isn't guaranteed to be the same under rescaling, though.

[andrewkern]

For example, in my group, if we're running with rescaling factor F=f for computational reasons, we'll compare results extensively with rescaling factor F=2f or F=10f to judge how sensitive our results are to the rescaling factor. (This is easy because those more aggressive scalings are even faster). This doesn't protect against effects that might occur between F=1 and F=f but that somehow saturate before F=2f, so it's something like an uncontrolled extrapolation -- definitely not ideal, but necessary in some cases.

So i've been in contact with Graham and he is onboard for us to explore this as part of our revision. I'd like to see a few scaling factors for each model and a comparison of summary stats to msprime and runtime

[grahamgower]

Thanks for the feedback everyone! New defaults were merged with #450.