Repository: tskit-dev/msprime (language: python)

Issue:
More details on defaults for TPM
We set p and m to some values estimated by Sainudiin et al, but it's not obvious to me where these values come from. Need more detail on this.
Documentation for EL2 contradicts signature
The documentation for EL2 talks about the ``p`` parameter, but this is not included in the signature and fixed at 0. The documentation is muddled also and needs another clarification pass.
Justify default values for lo and hi
Not obvious where these came from. In particular, why 2 and not 1 (or zero) for lo?

I've put a placeholder into the MicrosatMutationModel class.



Discussion / hints:
We could say something like:

Sainudiin et al (2004) estimated model parameters using 644 dinucleotide repeat loci identified as orthologous in human and chimpanzee genomes using maximum likelihood. TPM parameter estimates of `p` and `m` are taken from those MLEs.

is that what you are looking for @jeromekelleher?
Just a reference to which Table they come from would be enough. I was having trouble tracking the figures down quickly - we just want to make it easy for users to trace back to exactly what bit of the paper we're pulling these numbers from.

More fundamentally though - should we be putting in defaults for humans here? Do people use these models for nonhumans too? We don't put in defaults for particular organisms anywhere else in the API, so unless this model is *only* for humans, I wonder if we should just leave them unparametrised and instead give these values in some examples?
simple enough to add the table ref, but yes good questions. i would imagine that 

1) there aren't a lot of folks doing statistical estimation of these models any more (future project potential here for someone) 

2) that users when faced with a choice would use the defaults, which are specific and maybe not appropriate to all (many?) organisms / scenarios

so, sure, maybe we nuke the estimate default params and use these parameter values in an example, as you say?

these defaults weren't chosen for any good reason. `lo=2` makes sense as a default rather than 1 because these are repeats-- so repeat number > 1 by definition (usually)
Sound reasonable to me - we just need to document these choices. People will assume that the defaults are the Right Choice.

Fix patch:
diff --git a/docs/mutations.md b/docs/mutations.md
index bb5b4c685..892030790 100644
--- a/docs/mutations.md
+++ b/docs/mutations.md
@@ -93,9 +93,6 @@ for more information.
 {class}`.BinaryMutationModel` (Binary ancestral/derived)
 : Binary mutation model with two flip-flopping alleles: "0" and "1".
 
-{class}`.MatrixMutationModel` (General finite state model)
-: Superclass of mutation models with a finite set of states
-
 {class}`.SMM` (Microsatellites)
 : Stepwise mutation model for microsatellites
 
@@ -111,6 +108,9 @@ for more information.
 {class}`.SLiMMutationModel` (Integers)
 : An infinite-alleles model producing SLiM-style mutations
 
+{class}`.MatrixMutationModel` (General finite state model)
+: Superclass of mutation models with a finite set of states
+
 ---
 
 (sec_mutations_rate)=
@@ -501,7 +501,7 @@ in these models see {ref}`sec_mutations_matrix_mutation_models_details`.
 
 (sec_mutations_matrix_mutation_models_details)=
 
-### Mutation Matrix Models Details
+#### Mutation Matrix Models Details
 
 Mutation matrix models are specified by three things: an alphabet,
 a root distribution, and a transition matrix.
@@ -522,7 +522,7 @@ You can define your own, but you probably don't need to:
 there are several mutation matrix models already implemented in `msprime`,
 using binary (0/1), nucleotide, or amino acid alphabets:
 
-### Defining your own finite-sites model
+#### Defining your own finite-sites model
 
 If you want to define your own {class}`.MatrixMutationModel`, you have a good
 deal of freedom. For instance, here's a "decomposition/growth/disturbance"
@@ -560,7 +560,7 @@ for v in mts.variants():
 
 (sec_mutations_matrix_mutation_theory)=
 
-### Parameterisation of Matrix Mutation Models
+#### Parameterisation of Matrix Mutation Models
 
 Mutation matrix models are specified by three things: an alphabet,
 a root distribution, and a transition matrix.
@@ -668,7 +668,8 @@ print(f"Genetic diversity: {theta}.")
 
 That's pretty close to 0.001! The difference is within statistical error.
 
-(sec_mutations_mutation_infinite_alleles)=
+
+(sec_mutations_microsats)=
 
 ### Microsatellite Mutation Models
 
@@ -678,7 +679,8 @@ The basic idea here is that the number of copies of a given repeat is
 tracked, and its evolution over time subject to one of a number of potential
 biases.
 
-Consider the [Combined DNA Index System](https://en.wikipedia.org/wiki/Combined_DNA_Index_System) (CODIS), a US National database
+Consider the [Combined DNA Index System](https://en.wikipedia.org/wiki/Combined_DNA_Index_System) 
+(CODIS), a US National database
 of genotypes maintained by the justice system for use in forensic DNA analysis.
 The CODIS system relies on 20 core microsatellite loci, that are unlinked across
 the human genome. Let's simulate a sample of 5 individuals at 20 unlinked loci,
@@ -758,6 +760,52 @@ plt.ylabel('Variance in repeat number');
 So here we are seeing a clear relationship between the variance
 in repeat number and mutation rate, as expected. 🔥🔥🔥
 
+#### A more complicated microsatellite model
+
+In the 1990s a lot of work went in to describing patterns of
+mutation at microsatellite loci, and several models were
+put forward describing various biases in expansion vs contraction,
+mutistep mutations, and mutation rate biases. The 
+{class}`.MicrosatMutationModel` implements the general parameterization 
+of microsatellite mutation models developed in 
+[Sainudiin et al. (2004)](https://doi.org/10.1534/genetics.103.022665),
+allowing users fine grained control of microsatellite mutation.
+
+One concrete instance of this general model is the equal rate, linear biased, two-phase
+mutation model of [Garza et al. (1995)](https://doi.org/10.1093/oxfordjournals.molbev.a040239)
+that we have implemented in {class}`.EL2`. Let's simulate large
+samples under this model with different strengths of linear bias, 
+and compare the distribution of allele sizes we get back.
+
+```{code-cell} python
+from matplotlib import pyplot as plt
+from scipy import stats
+
+biases = np.logspace(-7, 1, num=9)
+ts = msprime.sim_ancestry(
+    1000, random_seed=2, sequence_length=1, population_size=100_000)
+for v in biases:
+    # these are values of m and u from Sainudiin et al. (2004)
+    model = msprime.EL2(m=0.01, u=0.68, v=v)
+    mts = msprime.sim_mutations(ts, rate=1e-3, model=model)
+    C = copy_number_matrix(mts)
+    kde = stats.gaussian_kde(C.flatten())
+    x = np.linspace(2, 51)
+    plt.plot(x, kde(x), label=f'v={v}')
+plt.xlabel('copy number')
+plt.ylabel('frequency')
+plt.legend()
+plt.show();
+```
+From these simulations we can see that the linear bias parameter `v`
+has a strong effect on the distribution of allele frequencies we
+might expect from our model. Care should be taken when choosing
+parameters for your own simulations. Potentially appropriate
+values for dinucleotide repeats in humans and chimp could be taken
+from Table 2 of [Sainudiin et al. (2004)](https://doi.org/10.1534/genetics.103.022665), 
+but we offer no guarantees and your mileage may vary. 
+
+(sec_mutations_mutation_infinite_alleles)=
 ### Infinite Alleles Mutation Models
 
 You can also use a model of *infinite alleles* mutation: where each new mutation produces a unique,
diff --git a/msprime/mutations.py b/msprime/mutations.py
index 98e3a7fa0..5cbb22c14 100644
--- a/msprime/mutations.py
+++ b/msprime/mutations.py
@@ -311,6 +311,9 @@ class MicrosatMutationModel(MatrixMutationModel):
     to the model of
     `Sainudiin et al. (2004) <https://doi.org/10.1534/genetics.103.022665>`_
 
+    .. seealso:: See the :ref:`sec_mutations_microsats` section for more details
+       and examples.
+
     Concretely, mutation rates under this model are computed as follows.
     Let ``qij`` be the transition rate from equation (3) in Sainudiin et al,
     and let ``M`` be the largest row sum of the matrix whose i,j-th entry
@@ -322,8 +325,12 @@ class MicrosatMutationModel(MatrixMutationModel):
     and are unrelated to repeat length
     (e.g., the repeat itself can be dinucleotide, tri-nucleotide, etcetera).
 
-    The default values for ``lo`` and ``hi`` are chosen based on
-    FIXME.
+    .. warning::
+        The default values for ``lo`` and ``hi`` are chosen arbitrarily.
+        The default of ``lo=2`` is chosed because repeat number need be greater
+        than for it to be considered a repeat. Please ensure
+        that you set ``lo`` and ``hi``
+        to the appropriate values for your organism / locus of interest.
 
     :param float s: strength of length dependence on mutation rate.
         Defaults to 0, i.e., no length dependence.
@@ -431,12 +438,10 @@ class TPM(MicrosatMutationModel):
     This is a :class:`.MicrosatMutationModel` with alleles ``[lo, .. , hi]``.
     For a precise definition of the mutation rates, see :class:`.MicrosatMutationModel`
     with ``s=0``, ``u=0.5``, and ``v=0``.
-    Default values for ``p`` and ``m``
-    are those estimated by
-    `Sainudiin et al. (2004) <https://doi.org/10.1534/genetics.103.022665>`_
+    See :class:`.MicrosatMutationModel` for further details of parameterization.
 
-    :param float p: Probability of a single step mutation. (Default=0.9)
-    :param float m: Success prob of truncated gamma distribution. (Default=0.93)
+    :param float p: Probability of a single step mutation.
+    :param float m: Success prob of truncated gamma distribution.
     :param int lo: Repeat number lower bound. See the documentation for
         :class:`.MicrosatMutationModel` for details.
     :param int hi: Repeat number upper bound. See the documentation for
@@ -445,9 +450,7 @@ class TPM(MicrosatMutationModel):
         :class:`.MicrosatMutationModel` for details.
     """
 
-    def __init__(self, *, p=None, m=None, lo=None, hi=None, root_distribution=None):
-        p = 0.9 if p is None else p
-        m = 0.93 if m is None else m
+    def __init__(self, *, p, m, lo=None, hi=None, root_distribution=None):
         if not (0 < p < 1.0):
             raise ValueError("p must be between 0 and 1")
         if not (0 < m < 1.0):
@@ -466,11 +469,17 @@ class EL2(MicrosatMutationModel):
 
     This models evolution of microsatellite repeat number in a manner
     that allows for multi-step mutations in copy number.
-    This is parameterized by `p` the probability of a single step mutation,
-    and `m` the success probability of the truncated gamma distribution
+    This is parameterized by `m` the success probability of
+    the truncated gamma distribution
     describing the distribution of longer steps, such that the
-    mean multi-step mutation is length `1/m`.
-    For this model both `p` and `m` need to be set to < 1.0.
+    mean multi-step mutation is length `1/m`. `u` the constant
+    bias parameter which determines bias (if any) in expansion or
+    contraction of the repeat number, and `v` the linear bias
+    parameter which determines mutation rate variation associated
+    with repeat number.
+    For this model `m` need to be set to < 1.0,
+    and `u` needs to be set to between 0 and 1.
+    and `v` can take any real, positive or negative value.
 
     See :class:`.MicrosatMutationModel` for details on the  constant
     bias parameter `u` and a linear bias parameter `v`.
@@ -478,12 +487,11 @@ class EL2(MicrosatMutationModel):
     This is a :class:`.MicrosatMutationModel` with alleles ``[lo, .. , hi]``
     For a precise definition of the mutation rates,
     see :class:`.MicrosatMutationModel`
-    with ``s=0``. Default values for ``m``, ``u`` and ``v``
-    are those estimated by Sainudiin et al. (2004).
+    with ``s=0``.
 
-    :param float m: Success prob of truncated gamma distribution (0.43).
-    :param float u: Constant bias parameter (Default=0.68).
-    :param float v: Linear bias parameter (Default=0.037).
+    :param float m: Success prob of truncated gamma distribution.
+    :param float u: Constant bias parameter.
+    :param float v: Linear bias parameter.
     :param int lo: Repeat number lower bound. See the documentation for
         :class:`.MicrosatMutationModel` for details.
     :param int hi: Repeat number upper bound. See the documentation for
@@ -492,12 +500,7 @@ class EL2(MicrosatMutationModel):
         :class:`.MicrosatMutationModel` for details.
     """
 
-    def __init__(
-        self, *, m=None, u=None, v=None, lo=None, hi=None, root_distribution=None
-    ):
-        m = 0.43 if m is None else m
-        u = 0.68 if u is None else u
-        v = 0.037 if v is None else v
+    def __init__(self, *, m, u, v, lo=None, hi=None, root_distribution=None):
         if not (0 < m < 1.0):
             raise ValueError("m must be between 0 and 1")
         if not (0 < u < 1.0):


Test patch:
diff --git a/tests/test_mutations.py b/tests/test_mutations.py
index 36326ccd8..bb64cc57e 100644
--- a/tests/test_mutations.py
+++ b/tests/test_mutations.py
@@ -381,7 +381,11 @@ def test_TPM(self):
                 assert np.isclose(exp, obs)
 
     def test_EL2(self):
-        model = msprime.EL2()
+        model = msprime.EL2(
+            m=0.43,
+            u=0.68,
+            v=0.037,
+        )
         self.validate_model(model)
         self.validate_stationary(model)
 
@@ -391,7 +395,7 @@ def test_EL2(self):
         u = 0.68
         v = 0.037
         p = 0
-        model = msprime.EL2(lo=lo, hi=hi)
+        model = msprime.EL2(lo=lo, hi=hi, m=m, u=u, v=v)
         for i in range(lo, hi + 1):
             if lo < i < hi:
                 ii = i - lo
@@ -469,24 +473,44 @@ def test_bad_params(self, m, u, v, p, s, lo, hi):
 
     def test_bad_TPM_p_value(self):
         with pytest.raises(ValueError, match="p must be"):
-            msprime.TPM(p=1.1)
+            msprime.TPM(m=0.9, p=1.1)
 
     def test_bad_TPM_m_value(self):
         with pytest.raises(ValueError, match="m must be"):
-            msprime.TPM(m=1.1)
+            msprime.TPM(m=1.1, p=0.9)
+
+    def test_bad_TPM_missing_m(self):
+        with pytest.raises(TypeError):
+            msprime.TPM(p=0.1)
+
+    def test_bad_TPM_missing_p(self):
+        with pytest.raises(TypeError):
+            msprime.TPM(m=0.1)
 
     def test_bad_EL2_m_value(self):
         with pytest.raises(ValueError, match="m must be"):
-            msprime.EL2(m=1.1)
+            msprime.EL2(m=1.1, u=0.6, v=0.3)
 
     def test_bad_EL2_u_value(self):
         with pytest.raises(ValueError, match="u must be"):
-            msprime.EL2(u=1.1)
+            msprime.EL2(m=0.9, u=1.1, v=0.6)
+
+    def test_bad_EL2_missing_m(self):
+        with pytest.raises(TypeError):
+            msprime.EL2(u=0.6, v=0.5)
+
+    def test_bad_EL2_missing_u(self):
+        with pytest.raises(TypeError):
+            msprime.EL2(m=0.6, v=0.5)
+
+    def test_bad_EL2_missing_v(self):
+        with pytest.raises(TypeError):
+            msprime.EL2(m=0.6, u=0.5)
 
     @pytest.mark.parametrize("v", [np.inf, np.nan])
     def test_bad_EL2_v_value(self, v):
         with pytest.raises(ValueError, match="finite"):
-            msprime.EL2(v=v)
+            msprime.EL2(v=v, m=0.9, u=0.6)
 
 
 class TestMutateRateMap: