Stdpopsim: validation of selection simulations

[grahamgower]

The test suite doesn't (and can't) statistically validate simulations. When the SLiM engine was added to stdpopsim, we visually inspected QQ plots comparing a range of summary statistics to the msprime engine (see validation.py). Similarly, msprime has a large number of statistical verification tests, separate to unit testing, that compare simulations against analytical results and other simulation software.

Now that we're adding non-neutral simulations to the SLiM engine, we should add new tests in our validation.py script to assess the correctness of their implementation. Suggestions for specific tests that could be implemented would be most welcome.

[jeromekelleher]

I guess comparisons with another forward simulator like fwdpy11 on simple models would be a good baseline. This would probably shake out any simple off-by-constant errors and so on. I'd imagine it'd be quite a bit of work to actually do, though.

[grahamgower]

Hmm, yeah that sounds quite hard. It might be almost as much work to just implement a fwdpy11 engine!

[jeromekelleher]

Yeah, that's the problem. Comparisons with some analytical results then? I guess there must be a bunch for simple selective models.

[grahamgower]

How about we compare the observed SFS with values computed using dadi/momi/moments?

[apragsdale]

How about we compare the observed SFS with values computed using dadi/momi/moments?

That probably makes the most sense here, I think? Running a bunch of slim replicates with a given demography, mutation rate, selection coefficient would be straightforward to compare to dadi or moments (or analytic results under steady-state demography) - I'd be happy to help with those comparisons.

[mufernando]

cc @petrelharp

[andrewkern]

why not validate against BGS results? like table 1 from here: https://www.genetics.org/content/141/4/1605.long

[grahamgower]

Looks interesting @andrewkern! Do you have an idea of what exactly to simulate in order to best match the expected values they calculated? Their expected value of π seems to correspond with a neutral site located in the middle of the contig---would it be sufficient to calculate π from all variable sites in a simulated contig (inc. deleterious)? all neutral sites?

[andrewkern]

i'd follow Figure 1 here where we have an embedded neutral locus in the middle of m loci that have deleterious mutations. in the paper Hudson and Kaplan use 10^4 sites with the middlemost being neutral.

[mufernando]

I'm interested in matching those Hudson and Kaplan equations with simulations, so might have a go at implementing the tests soon. I would need help with the structure/design of the statistical test script, though.

Just waiting on my PR #644 to be merged.

[David-Peede]

Just talked to @mufernando and our plan is to first recreate Hudson's and Kaplan's results in a general SLiM script outside of stdpopsim and then to write a stdpopsim version for validation.py

[petrelharp]

I think that we don't need to verify that SLiM does selection properly; that's a job for SLiM's unit tests. And matching to analytic predictions for selection is a giant pain. Here's my proposal:

1. get SLiM to optionally print out to a log file periodically (even, every generation?):
2. population sizes
3. migration rates
4. mutation types and their distribution of selection coefficients

5. (other stuff?)

6. Code up some consistency checks by hand in a few simple-but-tricky situations.

Is there anything we can't check this way?

[grahamgower]

I'm not so worried that SLiM has bugs (as you say, that's a job for SLiM's unit tests). I'm more worried that stdpopsim has a bug because we're telling SLiM to do the wrong thing. If you think it's reasonable to test this by comparing log output, then that would be great!

[petrelharp]

It seems like asking SLiM to tell us what it's doing, and comparing that to an independent tabulation what we think SLiM ought to be doing would test that? Getting the "independent tabulation" would be tricky to do generically, but shouldn't be bad on some specific test cases.