Skip to content

Commit

Permalink
Add vignettes on varying effect sizes across demes (#950)
Browse files Browse the repository at this point in the history
Closes #753
Closes #754
  • Loading branch information
molpopgen committed Nov 14, 2022
1 parent c0065f1 commit a3bc180
Show file tree
Hide file tree
Showing 7 changed files with 601 additions and 0 deletions.
7 changes: 7 additions & 0 deletions doc/_toc.yml
Original file line number Diff line number Diff line change
Expand Up @@ -32,6 +32,13 @@ parts:
- file: short_vignettes/demography_vignettes_intro
- file: short_vignettes/demes_vignette
- file: short_vignettes/demography_debugger_vignette
- caption: Short vignettes - varying fitness conditions across demes
chapters:
- file: short_vignettes/varying_fitness_conditions_across_demes_intro
- file: short_vignettes/multivariate_gaussian_effect_sizes_across_demes
- file: short_vignettes/multivariate_lognormal_effect_sizes
- file: short_vignettes/custom_multivarate_distributions
- file: short_vignettes/special_considerations_demes_models
- caption: Short vignettes - conditional models
chapters:
- file: short_vignettes/trackmutationfates
Expand Down
131 changes: 131 additions & 0 deletions doc/short_vignettes/custom_multivarate_distributions.md
Original file line number Diff line number Diff line change
@@ -0,0 +1,131 @@
---
jupytext:
formats: md:myst
text_representation:
extension: .md
format_name: myst
kernelspec:
display_name: Python 3
language: python
name: python3
---

(custom_multivariate_distributions_vignette)=

# "Custom" multivariate distributions

```{code-cell} python
---
tags: ['hide-input']
---
import demes
import demesdraw
import fwdpy11
import numpy as np
```

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
{class}`fwdpy11.Sregion` types, we follow a fairly intuitive approach
described in {cite}`Xue-Kun_Song2000-qn`. 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 {class}`list` of `n` {class}`fwdpy11.Sregion` objects
* An array of `n` means for the multivariate Gaussian
* An `n-by-n` covariance matrix for the multivariate
Gaussian

We will use the same demographic model as in the previous vignettes.
The details can be viewed by expanding the "click to show" button below".

```{code-cell}python
---
tags: ['hide-input']
---
yaml = """
description: Island model forever
time_units: generations
demes:
- name: A
epochs:
- start_size: 100
- name: B
epochs:
- start_size: 100
migrations:
- demes: [A, B]
rate: 0.10
"""
g = demes.loads(yaml)
model = fwdpy11.discrete_demography.from_demes(g, burnin=1)
initial_sizes = [v for v in model.metadata["initial_sizes"].values()]
pdict = {
"nregions": [],
"recregions": [],
"sregions": None, # Will get filled in below
"rates": (0, 5e-3, None),
"demography": model,
"simlen": model.metadata["total_simulation_length"],
"gvalue": fwdpy11.Multiplicative(ndemes=2, scaling=2),
}
rng = fwdpy11.GSLrng(123512)
```

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

```{code-cell} python
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)),
)
pdict["sregions"] = [mvdes]
params = fwdpy11.ModelParams(**pdict)
pop = fwdpy11.DiploidPopulation(initial_sizes, 1.0)
fwdpy11.evolvets(rng, pop, params, 10)
for i in pop.tables.mutations:
print(pop.mutations[i.key].esizes)
```

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:

```{code-cell} python
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)),
)
pdict["sregions"] = [mvdes]
params = fwdpy11.ModelParams(**pdict)
pop = fwdpy11.DiploidPopulation(initial_sizes, 1.0)
fwdpy11.evolvets(rng, pop, params, 10)
for i in pop.tables.mutations:
print(pop.mutations[i.key].esizes)
```

The type {class}`fwdpy11.ConstantS` has intuitive behavior:

```{code-cell} python
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)),
)
pdict["sregions"] = [mvdes]
params = fwdpy11.ModelParams(**pdict)
pop = fwdpy11.DiploidPopulation(initial_sizes, 1.0)
rng = fwdpy11.GSLrng(1010)
fwdpy11.evolvets(rng, pop, params, 10)
for i in pop.tables.mutations:
print(pop.mutations[i.key].esizes)
```
23 changes: 23 additions & 0 deletions doc/short_vignettes/genetic_values_multi_deme_models.md
Original file line number Diff line number Diff line change
@@ -0,0 +1,23 @@
---
jupytext:
formats: md:myst
text_representation:
extension: .md
format_name: myst
kernelspec:
display_name: Python 3
language: python
name: python3
---

(genetic_values_multi_deme_models)=

# Setting up genetic value objects.

NOTE: this will get too long -- split!!

## Case 1: constant effect size, different fitness effects.

## Case 2: different effect sizes, same fitness effects.

## Case 3: different effect sizes, different fitness effects.
146 changes: 146 additions & 0 deletions doc/short_vignettes/multivariate_gaussian_effect_sizes_across_demes.md
Original file line number Diff line number Diff line change
@@ -0,0 +1,146 @@
---
jupytext:
formats: md:myst
text_representation:
extension: .md
format_name: myst
kernelspec:
display_name: Python 3
language: python
name: python3
---

(multivariate_gaussian_effect_sizes_across_demes_vignette)=

# The multivariate Gaussian distribution

```{code-cell}python
---
tags: ['hide-input']
---
import demes
import demesdraw
import fwdpy11
import numpy as np
```


We may model Gaussian effect sizes using the existing {class}`fwdpy11.MultivariateGaussianEffects`
in conjunction with {class}`fwdpy11.mvDES`.

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.

We will simulate a demographic model of migration happening into the infinite past of two equal-sized demes:

```{code-cell}python
---
tags: ['hide-input']
---
yaml = """
description: Island model forever
time_units: generations
demes:
- name: ancestor
epochs:
- start_size: 100
end_time: 100
- name: A
ancestors: [ancestor]
epochs:
- start_size: 100
- name: B
ancestors: [ancestor]
epochs:
- start_size: 100
migrations:
- demes: [A, B]
rate: 0.10
"""
g = demes.loads(yaml)
model = fwdpy11.discrete_demography.from_demes(g, burnin=1)
demesdraw.tubes(g);
```

Let's set up a dictionary to hold the parameters of our model:

```{code-cell} python
pdict = {
"nregions": [],
"recregions": [],
"sregions": [
fwdpy11.mvDES(
fwdpy11.MultivariateGaussianEffects(0, 1, 1, np.identity(3)), np.zeros(3)
)
],
"rates": (0, 2.5e-3, None),
"demography": model,
"simlen": model.metadata["total_simulation_length"],
"gvalue": fwdpy11.Additive(
ndemes=3, 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 {class}`fwdpy11.mvDES` is our only region with selected mutations.
* This instance holds an instance of {class}`fwdpy11.MultivariateGaussianEffects`
that puts mutations on the interval {math}`[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:

```{code-cell} python
params = fwdpy11.ModelParams(**pdict)
# TODO: update this once we have a function to pull the sizes
# automatically from demes-derived models:
initial_sizes = [v for v in model.metadata["initial_sizes"].values()]
pop = fwdpy11.DiploidPopulation(initial_sizes, 1.0)
rng = fwdpy11.GSLrng(1010)
fwdpy11.evolvets(rng, pop, params, 10)
```

Let's extract the effect sizes from each deme:

```{code-cell} python
assert len(pop.tables.mutations) > 0
for i in pop.tables.mutations:
print(pop.mutations[i.key].esizes)
```

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

```{code-cell} python
# Effect sizes across demes will
# have a correlation coefficient of r=1/2
cor_matrix = np.array([-0.5]*9).reshape(3,3)
np.fill_diagonal(cor_matrix, np.array([1.0]*3))
# Get our covariance matrix
sd = np.array([0.1]*3)
D = np.identity(3)
np.fill_diagonal(D, sd)
vcv_matrix = np.matmul(np.matmul(D, cor_matrix), D)
pdict["sregions"] = [
fwdpy11.mvDES(fwdpy11.MultivariateGaussianEffects(0, 1, 1, vcv_matrix), np.zeros(3))
]
params = fwdpy11.ModelParams(**pdict)
# TODO: update this once we have a function to pull the sizes
# automatically from demes-derived models:
initial_sizes = [v for v in model.metadata["initial_sizes"].values()]
pop = fwdpy11.DiploidPopulation(initial_sizes, 1.0)
fwdpy11.evolvets(rng, pop, params, 10)
for i in pop.tables.mutations:
print(pop.mutations[i.key].esizes)
```

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

0 comments on commit a3bc180

Please sign in to comment.