# Different effect sizes of mutations in different demes¶

New in version 0.7.0.

This section describes how to allow mutations to have different effect sizes in different demes, building on the material from Soft selection with discrete demes.

To allow for correlated effect sizes across demes, we used `fwdpy11.mvDES`

to model multivariate
distributions of effect sizes. This class inherits from `fwdpy11.Sregion`

(see Defining distributions of mutation effect sizes).
In general, there is no standard *general* method for drawing random deviates from arbitrary multivariate
distributions. The approach that we take is to use a multivariate Gaussian distribution as the underlying kernel.

First, we will describe the cases where the multivariate Gaussian kernel leads to output effect size distributions
that have straightforward interpretations. Then, we will move onto the more general case allowing you to construct
multivariate distributions from current `fwdpy11.Sregion`

types.

## The multivariate Gaussian¶

We may model Gaussian effect sizes using the existing `fwdpy11.MultivariateGaussianEffects`

in conjunction with `fwdpy11.mvDES`

. Using `fwdpy11.MultivariateGaussianEffects`

on its
own is used to model pleiotropic effects on a trait. Here, we are re-using this type to model correlated
effect sizes across demes.

At this time, it is probably best to look at an example. The following code models Gaussian stabilizing selection on a quantitative trait. The effects sizes within each deme are themselves given by Gaussian distributions and there is no correlation in the effect size in the two demes.

```
In [1]: pdict = {
...: "nregions": [],
...: "recregions": [],
...: "sregions": [
...: fwdpy11.mvDES(
...: fwdpy11.MultivariateGaussianEffects(0, 1, 1, np.identity(2)), np.zeros(2)
...: )
...: ],
...: "rates": (0, 1e-3, None),
...: "demography": fwdpy11.DiscreteDemography(
...: migmatrix=np.array([0.9, 0.1, 0.1, 0.9]).reshape((2, 2))
...: ),
...: "simlen": 100,
...: "gvalue": fwdpy11.Additive(
...: ndemes=2, scaling=2, gvalue_to_fitness=fwdpy11.GSS(optimum=0.0, VS=10.0)
...: ),
...: }
...:
```

Most of the above is standard. Let’s dissect the new bits:

An instance of

`fwdpy11.mvDES`

is our only region with selected mutations.This instance holds an instance of

`fwdpy11.MultivariateGaussianEffects`

that puts mutations on the interval \([0, 1)\) with weight 1 and an identity matrix specifies the correlation in effect sizes between demes 0 and 1. The identity matrix has the value zero for all off-diagonal elements, meaning no covariance in effect sizes across demes.The final constructor argument specifies the mean of each marginal Gaussian distribution. The means are both zero.

Our genetic value type accepts an ndemes parameter, telling it that it has to look for deme-specific effect sizes. This value must be set to the maximum number of demes that will exist during a simulation.

Let’s evolve the model now:

```
In [2]: params = fwdpy11.ModelParams(**pdict)
In [3]: pop = fwdpy11.DiploidPopulation([100, 100], 1.0)
In [4]: rng = fwdpy11.GSLrng(1010)
In [5]: fwdpy11.evolvets(rng, pop, params, 10)
```

Let’s extract the effect sizes from each deme:

```
In [6]: for i in pop.tables.mutations:
...: print(pop.mutations[i.key].esizes)
...:
[ 0.09742853 -0.79068629]
[-0.46314495 -0.31686816]
```

Let’s look at another example where effect sizes covary negatively across demes and raise the mutation rate a bit:

```
In [7]: vcv = np.array([1.0, -0.99, -0.99, 1.0]).reshape((2, 2))
In [8]: pdict["sregions"] = [
...: fwdpy11.mvDES(fwdpy11.MultivariateGaussianEffects(0, 1, 1, vcv), np.zeros(2))
...: ]
...:
In [9]: pdict["rates"] = (0, 5e-3, None)
In [10]: params = fwdpy11.ModelParams(**pdict)
In [11]: pop = fwdpy11.DiploidPopulation([100, 100], 1.0)
In [12]: fwdpy11.evolvets(rng, pop, params, 10)
In [13]: for i in pop.tables.mutations:
....: print(pop.mutations[i.key].esizes)
....:
[-0.65688902 0.67458945]
[ 1.44236014 -1.19184102]
[ 0.53515552 -0.49690106]
[0.11239562 0.03076654]
[ 0.0470034 -0.14227415]
[-0.44707712 0.39814511]
[-0.64860038 0.53073426]
[-0.12680578 -0.00521323]
[-1.702493 1.78230999]
[-0.12798166 -0.14205828]
[ 0.82686083 -0.68915806]
```

Now we see that the effect sizes often differ in sign between the two demes.

## The multivariate lognormal¶

If \(X\) is a multivariate Gaussian distribution, \(N(\mathbf{\mu}, \mathbf{\sum})\), where \(\mathbf{\mu}\) is a vector of mean values and \(\mathbf{\sum}\) is the covariance matrix, then \(Y = e^X\) is a multivariate lognormal random variable with mean \(E[Y]_i = e^{\mu_i + \frac{1}{2}\sum_{ii}}\) and covariance matrix \(Var[Y]_{i,j} = e^{\mu_i + \mu_j + \frac{1}{2}(\sum_{ii} + \sum_{jj})}(e^{\sum_{ij}}-1)\).

To specify a multivariate lognormal distribution of effect sizes, we use
the static class method `fwdpy11.LogNormalS.mv()`

. The following code
constructs a distribution of effect sizes such that -2Ns (where N is the
size of a single deme) is a multivariate lognormal with means zero and an
identity matrix as a covariance matrix used to specify the multivate
Gaussian kernel.

```
In [14]: mvdes = fwdpy11.mvDES(
....: fwdpy11.LogNormalS.mv(0, 1, 1, scaling=-200), np.zeros(2), np.identity(2)
....: )
....:
```

Note

The lognormal distribution returns deviates \(> 0\). To model deleterious mutations/effect sizes < 0, use the scaling parameter with a negative value like we just did!

Let’s put it in a simulation and run it:

```
In [15]: pdict = {
....: "nregions": [],
....: "recregions": [],
....: "sregions": [mvdes],
....: "rates": (0, 1e-3, None),
....: "demography": fwdpy11.DiscreteDemography(
....: migmatrix=np.array([0.9, 0.1, 0.1, 0.9]).reshape((2, 2))
....: ),
....: "simlen": 100,
....: "gvalue": fwdpy11.Multiplicative(ndemes=2, scaling=2),
....: }
....:
In [16]: params = fwdpy11.ModelParams(**pdict)
In [17]: pop = fwdpy11.DiploidPopulation([100, 100], 1.0)
In [18]: fwdpy11.evolvets(rng, pop, params, 10)
In [19]: for i in pop.tables.mutations:
....: print(pop.mutations[i.key].esizes)
....:
[-0.00183755 -0.0009758 ]
[-0.00110216 -0.00049425]
```

## “Custom” multivariate distributions¶

The previous two sections cover cases where the methods for generating deviates from a multivariate distribution are straightforward and agreed upon.

In order to simulate multivariate distributions of effect sizes based on
`fwdpy11.Sregion`

types, we follow a fairly intuitive approach
described in [Song2000]. Briefly, the multivariate Gaussian kernel is
used to produce deviates. Then, the quantiles from the cummulative distribution
of each marginal Gaussian are used to generate a deviate from the desired output distribution of interest.

For a simulation with n populations we need:

A

`list`

of n`fwdpy11.Sregion`

objectsAn array of n means for the multivariate Gaussian

An n-by-n covariance matrix for the multivariate Gaussian

The following generates exponentially distributed effect sizes in each deme with a high correlation across demes:

```
In [20]: mvdes = fwdpy11.mvDES(
....: [fwdpy11.ExpS(0, 1, 1, -0.5)] * 2,
....: np.zeros(2),
....: np.matrix([1, 0.9, 0.9, 1]).reshape((2, 2)),
....: )
....:
In [21]: pdict["sregions"] = [mvdes]
In [22]: params = fwdpy11.ModelParams(**pdict)
In [23]: pop = fwdpy11.DiploidPopulation([100, 100], 1.0)
In [24]: fwdpy11.evolvets(rng, pop, params, 10)
In [25]: for i in pop.tables.mutations:
....: print(pop.mutations[i.key].esizes)
....:
[-0.05297616 -0.05994487]
[-0.00600096 -0.00063 ]
```

We can mix and match our distributions. Here, the distribution of effect sizes in deme 0 is exponential and the distribution in deme 1 is gamma. The two distributions have means with opposite signs and the magnitudes of the marginal deviates negatively covary:

```
In [26]: mvdes = fwdpy11.mvDES(
....: [fwdpy11.ExpS(0, 1, 1, -0.5), fwdpy11.GammaS(0, 1, 1, mean=0.1, shape_parameter=1)],
....: np.zeros(2),
....: np.matrix([1, -0.9, -0.9, 1]).reshape((2, 2)),
....: )
....:
In [27]: pdict["sregions"] = [mvdes]
In [28]: params = fwdpy11.ModelParams(**pdict)
In [29]: pop = fwdpy11.DiploidPopulation([100, 100], 1.0)
In [30]: fwdpy11.evolvets(rng, pop, params, 10)
In [31]: for i in pop.tables.mutations:
....: print(pop.mutations[i.key].esizes)
....:
[-0.17037555 0.17324704]
[-0.22122016 0.22419496]
[-1.78012379 0.01713378]
[-0.36227095 0.10331541]
[-0.24379116 0.12145104]
```

The type `fwdpy11.ConstantS`

has intuitive behavior:

```
In [32]: mvdes = fwdpy11.mvDES(
....: [fwdpy11.ExpS(0, 1, 1, -0.5), fwdpy11.ConstantS(0, 1, 1, -0.1)],
....: np.zeros(2),
....: np.matrix([1, -0.9, -0.9, 1]).reshape((2, 2)),
....: )
....:
In [33]: pdict["rates"] = (0, 5e-3, None)
In [34]: pdict["sregions"] = [mvdes]
In [35]: params = fwdpy11.ModelParams(**pdict)
In [36]: pop = fwdpy11.DiploidPopulation([100, 100], 1.0)
In [37]: rng = fwdpy11.GSLrng(1010)
In [38]: fwdpy11.evolvets(rng, pop, params, 10)
In [39]: for i in pop.tables.mutations:
....: print(pop.mutations[i.key].esizes)
....:
[-0.99150137 -0.1 ]
[-0.22828726 -0.1 ]
[-0.10570464 -0.1 ]
[-0.06271967 -0.1 ]
```

## Recipes¶

### Different signs in different demes¶

Consider two demes. You want any beneficial mutation in one deme to be deleterious in the other and vice-versa.

For the multivariate Gaussian, use the covariance matrix as done above. Note
that this approach only generates a *tendency* to different signs in different demes.

With the multivariate lognormal, the best we can do is to use negative correlations such that deleterious mutations in deme 0 are less deleterious in deme 1, etc.:

```
In [40]: sregions = [
....: fwdpy11.mvDES(
....: fwdpy11.LogNormalS.mv(0, 1, 1, scaling=-200),
....: np.zeros(2),
....: np.matrix([1, -0.99, -0.99, 1]).reshape((2, 2)),
....: )
....: ]
....:
In [41]: sregions.append(
....: fwdpy11.mvDES(
....: fwdpy11.LogNormalS.mv(0, 1, 1, scaling=200),
....: np.zeros(2),
....: np.matrix([1, -0.99, -0.99, 1]).reshape((2, 2)),
....: )
....: )
....:
In [42]: pdict["sregions"] = sregions
In [43]: params = fwdpy11.ModelParams(**pdict)
In [44]: pop = fwdpy11.DiploidPopulation([100, 100], 1.0)
In [45]: rng = fwdpy11.GSLrng(1010)
In [46]: fwdpy11.evolvets(rng, pop, params, 10)
In [47]: for i in pop.tables.mutations:
....: print(pop.mutations[i.key].esizes)
....:
[0.00039308 0.0693669 ]
[-0.01488498 -0.00173671]
[-0.00514551 -0.00519719]
[-0.00377736 -0.00592563]
[0.00069116 0.0453049 ]
[0.00355528 0.00643812]
[-0.00208255 -0.01409883]
[-0.01590584 -0.00175813]
[-0.00045263 -0.04933448]
[0.00227164 0.01008629]
[0.01128436 0.00198865]
[0.01109052 0.00222656]
[-0.00256051 -0.01002711]
[-0.00152782 -0.01807099]
```

In the output, we see that an effect size in deme i has a corresponding effect size in deme j that is a about an order of magnitude smaller in absolute value.

For the general approach, simply create a `list`

of objects with the desired mean (or constant) effect sizes. For example:

```
In [48]: sregions = [
....: fwdpy11.mvDES(
....: [fwdpy11.ExpS(0, 1, 1, -0.5), fwdpy11.ExpS(0, 1, 1, 0.1)],
....: np.zeros(2),
....: np.identity(2),
....: )
....: ]
....:
In [49]: sregions.append(
....: fwdpy11.mvDES(
....: [fwdpy11.ExpS(0, 1, 1, 0.1), fwdpy11.ExpS(0, 1, 1, -0.5)],
....: np.zeros(2),
....: np.identity(2),
....: )
....: )
....:
In [50]: pdict["sregions"] = sregions
In [51]: params = fwdpy11.ModelParams(**pdict)
In [52]: pop = fwdpy11.DiploidPopulation([100, 100], 1.0)
In [53]: rng = fwdpy11.GSLrng(1010)
In [54]: fwdpy11.evolvets(rng, pop, params, 10)
In [55]: for i in pop.tables.mutations:
....: print(pop.mutations[i.key].esizes)
....:
[-0.99150137 0.08290456]
[ 0.04565745 -0.16002424]
[-0.10570464 0.21656547]
[ 0.41101789 -0.43686025]
[-0.88406416 0.16572598]
[ 0.25571088 -0.08508201]
[ 0.15472792 -0.292425 ]
[-0.14495801 0.08997944]
[-0.06271967 0.15351897]
```