Repository: popsim-consortium/stdpopsim (language: python)

Issue:
log-normal DFE doesn't allow negative values
Over in #1292 I turned on the test that actually simualted from all of the defined DFEs. The [`DroMel.LognormalPlusPositive_R16`](https://github.com/popsim-consortium/stdpopsim/blob/main/stdpopsim/catalog/DroMel/dfes.py#L53) model uses a lognormal distribution for deleterious mutations. However, the simulation run errored because the selection coefficients that it produced were positive, and sometimes very very large. The problem is that we're trying to use something like `rlnorm(meanlog=-3, sdlog=2)` to get negative numbers. That's in fact how many distributions work in eidos: passing a negative mean for a positive distribution gets you numbers drawn from that distribution multipied by -1. However, the arguments to `rlnorm( )` are the mean and SD on a log scale, i.e., `mean(log( rlnorm(10000, -2, 4) ))` will be close to -2. (and so `mean(rlnorm(10000, -2, 4))` will be close to `exp(-2)` which is `>0`)

So, we need to fix up our underlying log-normal implementation somehow to allow this possibility.

I've also commented out the DFE so it doesn't error (in that PR). I think there are other issues with it (for instance, in the code it calls the meanlog and sdlog parameters the "shape" and "scale", which they are not).

As for the fix: recall that for most distributions we can use SLiM's built-in distribution types, e.g. `"e"` for exponential; however, it doesn't have a built-in log-normal, so we provide that one [in eidos code](https://github.com/popsim-consortium/stdpopsim/blob/main/stdpopsim/dfe.py#L111).  I propose *removing" the `"l"` distribution type (which is what we used for "log-normal", and replacing it with "lp" for "positive log-normal" and "ln" for "negative lognormal".


Discussion / hints:
> As for the fix: recall that for most distributions we can use SLiM's built-in distribution types, e.g. `"e"` for exponential; however, it doesn't have a built-in log-normal, so we provide that one [in eidos code](https://github.com/popsim-consortium/stdpopsim/blob/main/stdpopsim/dfe.py#L111). I propose *removing" the `"l"` distribution type (which is what we used for "log-normal", and replacing it with "lp" for "positive log-normal" and "ln" for "negative lognormal".

seems like a perfectly fine fix to me. in practice, are people fitting models of this sort?

Fix patch:
diff --git a/stdpopsim/catalog/DroMel/dfes.py b/stdpopsim/catalog/DroMel/dfes.py
index f68a4e344..c7008d6ec 100644
--- a/stdpopsim/catalog/DroMel/dfes.py
+++ b/stdpopsim/catalog/DroMel/dfes.py
@@ -89,6 +89,9 @@ def _RagsdaleDFE():
     # selection strengths
     S_positive = 39.9
     s_positive = S_positive / 2 / Ne
+    # TODO: these are *not* scale and shape parameters:
+    # the arguments to lognorm are meanlog and sdlog.
+    # Rename these after checking!
     scale_negative = 5.42 / 2 / Ne
     shape_negative = 3.36
     h = 0.5
@@ -96,8 +99,8 @@ def _RagsdaleDFE():
     synonymous = stdpopsim.MutationType()
     negative = stdpopsim.MutationType(
         dominance_coeff=h,
-        distribution_type="l",  # log-normal distribution
-        distribution_args=[-scale_negative, shape_negative],
+        distribution_type="ln",  # negative log-normal distribution
+        distribution_args=[scale_negative, shape_negative],
     )
     positive = stdpopsim.MutationType(
         dominance_coeff=h,
@@ -115,6 +118,4 @@ def _RagsdaleDFE():
     )
 
 
-# This DFE is NOT correct due to a bug/deficiency: the log-normal implementation
-# does not allow negative values.
-# _species.add_dfe(_RagsdaleDFE())
+_species.add_dfe(_RagsdaleDFE())
diff --git a/stdpopsim/dfe.py b/stdpopsim/dfe.py
index 4ab530a68..78ac8e7c1 100644
--- a/stdpopsim/dfe.py
+++ b/stdpopsim/dfe.py
@@ -23,7 +23,13 @@ class MutationType(object):
     - ``g``: gamma, two parameters (mean, shape)
     - ``n``: normal, two parameters (mean, SD)
     - ``w``: Weibull, two parameters (scale, shape)
-    - ``l``: logNormal, two parameters (mean and sd on log scale; see rlnorm)
+    - ``lp``: positive logNormal, two parameters (mean and sd on log scale; see rlnorm)
+    - ``ln``: negative logNormal, two parameters (mean and sd on log scale; see rlnorm)
+
+    Type "lp" is always positive, and type "ln" is always negative: both use
+    the same log-normal distribution, but "ln" is multiplied by -1.  For
+    exponential and gamma, a negative mean can be provided, obtaining always
+    negative values.
 
     :ivar distribution_type: A one-letter abbreviation for the distribution of
         fitness effects that each new mutation of this type draws from (see below).
@@ -117,22 +123,26 @@ def __attrs_post_init__(self):
                 raise ValueError("The scale parameter must be positive.")
             if self.distribution_args[1] <= 0:
                 raise ValueError("The shape parameter must be positive.")
-        elif self.distribution_type == "l":
-            # An lognormally-distributed fitness effect (meanlog, sdlog).
+        elif self.distribution_type in ("lp", "ln"):
+            # A lognormally-distributed fitness effect (meanlog, sdlog),
+            # either positive or negative.
             # See Eidos documentation for rlnorm().
             if len(self.distribution_args) != 2:
                 raise ValueError(
-                    "Lognormally-distributed sel. coefs. (distribution_type='l') "
+                    "Lognormally-distributed sel. coefs. (distribution_type='lp'/'ln') "
                     "use a (meanlog, sdlog) parameterisation, requiring sdlog > 0."
                 )
             if self.distribution_args[1] < 0:
                 raise ValueError("The sdlog parameter must be nonnegative.")
-            self.distribution_type = "s"
             # dealing with lognormal distribution
             # (adding instead of multiplying the mean):
             logmean = self.distribution_args[0]
             logsd = self.distribution_args[1]
-            self.distribution_args = [f"return rlnorm(1, {logmean} + log(Q), {logsd});"]
+            sign = "" if self.distribution_type == "lp" else "-1 *"
+            self.distribution_args = [
+                f"return {sign}rlnorm(1, {logmean} + log(Q), {logsd});"
+            ]
+            self.distribution_type = "s"
         else:
             raise ValueError(
                 f"{self.distribution_type} is not a supported distribution type."


Test patch:
diff --git a/tests/test_dfes.py b/tests/test_dfes.py
index f0b860bd2..21e98a76a 100644
--- a/tests/test_dfes.py
+++ b/tests/test_dfes.py
@@ -30,7 +30,8 @@ def test_Q_scaled_index(self):
             "e": ([1], [0]),
             "g": ([0.014, 0.19], [0]),
             "n": ([0.5, 1], [0, 1]),
-            "l": ([0.5, 1], []),
+            "ln": ([0.5, 1], []),
+            "lp": ([0.5, 1], []),
         }
         for t in mut_params:
             print(f"{t}\t{mut_params[t]}")
@@ -157,19 +158,22 @@ def test_create_bad_mutation_type_message(self):
             )
 
         # Lognormally-distributed selection coefficients
-        with pytest.raises(ValueError, match="use a .meanlog, sdlog. parameterisation"):
-            dfe.MutationType(
-                distribution_type="l",
-                distribution_args=[1, 2, 3, 4],
-            )
+        for dt in ["lp", "ln"]:
+            with pytest.raises(
+                ValueError, match="use a .meanlog, sdlog. parameterisation"
+            ):
+                dfe.MutationType(
+                    distribution_type=dt,
+                    distribution_args=[1, 2, 3, 4],
+                )
 
-        with pytest.raises(
-            ValueError, match="The sdlog parameter must be nonnegative."
-        ):
-            dfe.MutationType(
-                distribution_type="l",
-                distribution_args=[1, -2],
-            )
+            with pytest.raises(
+                ValueError, match="The sdlog parameter must be nonnegative."
+            ):
+                dfe.MutationType(
+                    distribution_type=dt,
+                    distribution_args=[1, -2],
+                )
 
     def test_mutation_type_is_neutral(self):
         mt = dfe.MutationType()
@@ -200,12 +204,13 @@ def test_mutation_types(self):
             "e": ([0.1], [10], [5000], [0]),
             "n": ([-0.1, 0.2], [0.1, 0.1], [50, 50]),
             "w": ([0.1, 0.2], [0.1, 0.1], [50, 50]),
-            "l": ([-0.1, 0.2], [0.1, 0.1], [50, 50]),
+            "lp": ([-0.1, 0.2], [0.1, 0.1], [50, 50]),
+            "ln": ([-0.1, 0.2], [0.1, 0.1], [50, 50]),
         }
         for t in mut_params:
             for p in mut_params[t]:
                 mt = dfe.MutationType(distribution_type=t, distribution_args=p)
-                if t == "l":
+                if t in ("lp", "ln"):
                     assert mt.distribution_type == "s"
                 else:
                     assert mt.distribution_type == t
@@ -220,7 +225,8 @@ def test_bad_mutation_types(self):
             "e": ([], [0, 1], [0.1, 0.4, 0.5], [np.inf]),
             "n": ([], [0.1, -1], [0.1, 0.4, 0.5], [0.1], [0.3, np.inf]),
             "w": ([], [-0.1, 1], [0.1, -1], [0.1, 0.4, 0.5], [0.1], [np.inf, 2.3]),
-            "l": ([], [0.1, -1], [0.1, 0.4, 0.5], [0.1], [0.1, np.inf]),
+            "lp": ([], [0.1, -1], [0.1, 0.4, 0.5], [0.1], [0.1, np.inf]),
+            "ln": ([], [0.1, -1], [0.1, 0.4, 0.5], [0.1], [0.1, np.inf]),
         }
         for t in bad_mut_params:
             for p in bad_mut_params[t]:
@@ -662,7 +668,7 @@ def test_simulation_runs(self):
         samples = model.get_samples(2)
         engine = stdpopsim.get_engine("slim")
         ts = engine.simulate(
-            model, contig, samples, slim_scaling_factor=10, slim_burn_in=10
+            model, contig, samples, slim_scaling_factor=10, slim_burn_in=10, seed=42
         )
 
         mut_info = {}
diff --git a/tests/test_slim_engine.py b/tests/test_slim_engine.py
index 4d4b5a9ec..f73ecbcee 100644
--- a/tests/test_slim_engine.py
+++ b/tests/test_slim_engine.py
@@ -891,14 +891,15 @@ class TestGenomicElementTypes(PiecewiseConstantSizeMixin):
         "e": ([0.1],),
         "n": ([-0.1, 0.2],),
         "w": ([0.1, 0.4],),
-        "l": ([-0.1, 0.2],),
+        "lp": ([-5, 0.2],),
+        "ln": ([-5, 0.2],),
     }
-    example_mut_types = [stdpopsim.MutationType()] + [
-        stdpopsim.MutationType(distribution_type=t, distribution_args=p)
+    example_mut_types = [("f", stdpopsim.MutationType())] + [
+        (t, stdpopsim.MutationType(distribution_type=t, distribution_args=p))
         for t, params in mut_params.items()
         for p in params
     ]
-    assert len(example_mut_types) == 9
+    assert len(example_mut_types) == 10
 
     def get_example_dfes(self):
         # this is in a function because scoping is weird
@@ -908,14 +909,14 @@ def get_example_dfes(self):
                 id="simple",
                 description="just one mut type",
                 long_description="🐺 🦊 🐕 🦝 🦮 🐩",
-                mutation_types=[self.example_mut_types[4]],
+                mutation_types=[self.example_mut_types[4][1]],
             ),
             stdpopsim.DFE(
                 id="less_simple",
                 description="two mutation types",
                 long_description="🦕 🦖 🐳 🐋 🐬 🦭🐠 🐡 🦈 🐙",
                 proportions=[0.5, 0.5],
-                mutation_types=self.example_mut_types[2:4],
+                mutation_types=[m for _, m in self.example_mut_types[2:4]],
             ),
             stdpopsim.DFE(
                 id="everything",
@@ -927,14 +928,16 @@ def get_example_dfes(self):
                     / (len(self.example_mut_types) * (len(self.example_mut_types) - 1))
                     for j in range(len(self.example_mut_types))
                 ],
-                mutation_types=self.example_mut_types,
+                mutation_types=[m for _, m in self.example_mut_types],
             ),
         ]
         return example_dfes
 
     def test_slim_simulates_dfes(self):
         contig = get_test_contig()
-        contig.mutation_rate = 10 * contig.mutation_rate
+        # make theta = 40
+        contig.mutation_rate = 10 / (1000 * contig.recombination_map.sequence_length)
+        print(contig.mutation_rate)
         dfes = [
             stdpopsim.DFE(
                 id=str(j),
@@ -942,7 +945,7 @@ def test_slim_simulates_dfes(self):
                 long_description=f"mut_{j}",
                 mutation_types=[mt],
             )
-            for j, mt in enumerate(self.example_mut_types)
+            for j, (_, mt) in enumerate(self.example_mut_types)
         ]
         n = len(dfes)
         breaks = np.linspace(0, contig.length, n + 1)
@@ -954,16 +957,17 @@ def test_slim_simulates_dfes(self):
             demographic_model=self.model,
             contig=contig,
             samples=self.samples,
-            seed=123,
+            seed=124,
             verbosity=3,
         )
-        for j, mt in enumerate(self.example_mut_types):
+        for j, (t, mt) in enumerate(self.example_mut_types):
             sites = np.where(
                 np.logical_and(
                     ts.tables.sites.position >= breaks[j],
                     ts.tables.sites.position < breaks[j + 1],
                 )
             )[0]
+            print(mt)
             assert len(sites) > 0
             for k in sites:
                 s = ts.site(k)
@@ -971,6 +975,10 @@ def test_slim_simulates_dfes(self):
                     for md in mut.metadata["mutation_list"]:
                         if not mt.is_neutral:
                             assert md["selection_coeff"] != 0
+                    if t == "lp":
+                        assert md["selection_coeff"] > 0
+                    elif t == "ln":
+                        assert md["selection_coeff"] < 0
 
     def slim_metadata_key0(self, metadata, key):
         # Everything in SLiM is a vector, so you can't just put a singleton
@@ -1172,14 +1180,14 @@ def test_same_dfes(self):
             description="the first one",
             long_description="hello world",
             proportions=[1.0],
-            mutation_types=[self.example_mut_types[5]],
+            mutation_types=[self.example_mut_types[5][1]],
         )
         dfe1 = stdpopsim.DFE(
             id="dfe",
             description="the second one",
             long_description="I'm different but have the same name! =( =(",
             proportions=[0.2, 0.8],
-            mutation_types=self.example_mut_types[7:9],
+            mutation_types=[m for _, m in self.example_mut_types[7:9]],
         )
         contig.add_dfe(np.array([[0, 0.5 * L]], dtype="int"), dfe0)
         contig.add_dfe(np.array([[0.2 * L, 0.4 * L]], dtype="int"), dfe0)
@@ -1207,7 +1215,7 @@ def test_slim_produces_mutations(self):
             mutation_types=[
                 stdpopsim.MutationType(),
                 stdpopsim.MutationType(
-                    distribution_type="l", distribution_args=[0.01, 0.2]
+                    distribution_type="lp", distribution_args=[0.01, 0.2]
                 ),
             ],
             proportions=[0.5, 0.5],