Open BenjaminPeter opened 3 months ago
Hm, you're right. This was certainly not intentional. And, it's true for f4
as well:
import msprime as ms
ts = ms.sim_ancestry(sequence_length=1000, samples=4, ploidy=2, random_seed=1)
m = ms.sim_mutations(ts, rate=5e-3, random_seed=2)
a, b, c, d = (0, 1), (2, 3), (4, 5), (6, 7)
f12 = m.f4((a, b, c, d), span_normalise=False)
f21 = m.f4((c, d, a, b), span_normalise=False)
assert f12 == f21, f"f4: {f12} !={f21}"
Here's the issue: the definition we gave for f4 (and by extension, f3 and f2) is
The average density of sites at which a and c agree but differs from b and d, minus the average density of sites at which a and d agree but differs from b and c
This sounds symmetric under switching (a,b) with (c,d), and indeed it is for biallelic loci, since
a and c agree, but differ from b and d
is the same for biallelic loci as
a=c != b=d
and hence the same as
b and d agree, but differ from a and c
However, with more than two alleles, we can have b and d differing yet still neither agreeing with a and c.
So, I see two options:
The average density of sites at which a and c agree but differs from b and d, plus the average density of sites at which b and d agree but differs from a and c, minus the average density of sites at which a and d agree but differs from b and c and the average density of sites at which b and c agree but differs from a and d, all divided by two.
I'm in favor of (2) because: a. people can symmetrize it themselves if they want; b. the current version seems a bit more natural and perhaps informative than the symmetrized version c. the only downside of the current version I'm aware of is that it's surprising, so documenting should fix that.
However, if you know of a reason that the symmetrized version (or something else?) is a more natural estimator or something, then that's worth considering?
FYI, we've defined all the statistics in a way that makes sense with more-than-biallelic sites: the strategy taken is to treat each allele as a binary split between "the allele" and "all other alleles", and compute statistics in a way that sums a summary over these splits (see the paper for more discussion). I'm not aware of any other issues like this one that arises from multiallelic sites: for instance, ts.divergence( )
is still symmetric.
Hi Peter,
Thanks a lot for your thorough investigation and explanation of the issue.
Perhaps a lot of the difference in intuition (and why this behavior of tskit was surprising to me) is that I tend to think of $F_4$ as a sum of $F_2$ statistics, whereas your approach seems to start with $F_4$, and treats $F_2$ as a special case.
However, I think from the $F_2$ perspective, the symmetric version you describe as the more intuitive and practical extension of $f_2$ to multiallelic sites, here is my pitch: ;-)
In addition, most applications of $F$-statistics I have seen do not calculate all permutations of arguments. This is both true for studies that use $F$-statistics for individual tests, and tools that use matrices of $F$-statistics for more complex models, such as qpAdm.
In summary, I think the symmetric version is just the nicer, more convenient and, to me, more intuitive extension of $F$-stats to multiallelic sites. It retains most of the properties and interpretations from the biallelic case, whereas the current version does not.
Ben
Thanks very much - this is compelling (but I'm not yet sure). I think we considered the $\pi$-based definitions you present here, but used the definitions we ended up with because (a) they were equivalent for biallelic sites, and (b) they were easier to explain in words. But, do you have a good and concise way to explain $F2(a,b) = \pi{a,b} - \pi_a/2 - \pi_b/2$ in terms of what statistic it is of the samples?
However, I think that our current version, symmetrized (i.e., (ts.f4(a,b,c,d) + ts.f4(c,d,a,b))/2
) is not equal to the pairwise divergences version. So, I think what you've got is a third proposal for the definitions, in terms of divergence and diversity, and distinct from the "let's just symmetrize" proposal. Or, have I got this wrong?
Thanks for helping get this right.
p.s. I ran across the thread in which we worked out the multiallelic thing; it would have been nice to have your input at the time!
The implementation of the $f_2$-statistic for data with multiallelic sites is not symmetric (as would be expected from the definition of the statistic), I narrowed it down to the following minimal example:
Digging a bit deeper, I found the issue is very likely due to the presence of a multiallelic site; coming from
ms
I assumedmsprime
would generate infinite sites by default, when it does not. Usingfixes the problem
As F-statistics are not commonly run (or even defined for) multiallelic sites, this might not be a great practical problem, but I decided to post it anyways since it appears like the code for many statistics is shared.
Other info:
The data set generated is
Versions used: