Msprime: Something odd with migration to pop 0 in DTWF

[nspope]

Hi, thank you for developing this wonderful tool.

When simulating from the DTWF model with island-with-migration population structure (where migration rate varies; e.g. with distance) and very low effective population sizes per deme, I notice something like an artefact in simulated tree topology with regards to tip order.

For example, (A) I generate a random migration matrix by simulating points in the unit square and mapping distance to migration rate with an exponential function; (B) permute rows and columns of the migration matrix; (C) set deme (diploid) Ne to 1; (D) simulate independent genealogies of a single haploid per deme, with deme Ne set to 1; (E) calculate the average pairwise TMRCA of a given tip, averaged over all other tips and loci; (F) repeat B through E ... then this is what the distribution of the average pairwise TMRCA looks like, across permutations:

I am using msprime 0.7.4 and can replicate with the latest version of msprime/master. This does not seem like intended behavior, because the first tip generally has a greater TMRCA (regardless of the spatial location of the first tip, which is really the only thing in the model that makes tips non-exchangeable)? Here is code to reproduce:

import msprime
import numpy
import scipy.spatial.distance as dist
import pandas as pd
import matplotlib.pyplot as plt

# expectation of TMRCA under DTWF island model
def island_model_TMRCA (migr_mat, reps = 100, seed = 1):
    migr_mat   = numpy.array(migr_mat, dtype=float)
    num_popul  = migr_mat.shape[0]
    pop_confg  = []
    for i in range(num_popul):
      pop_confg += [msprime.PopulationConfiguration(
        sample_size=1, initial_size=1.)]
    sim = msprime.simulate(
            population_configurations = pop_confg,
            migration_matrix = migr_mat,
            mutation_rate = 0.,
            Ne = 1.,
            num_replicates = reps,
            random_seed = seed,
            model = "dtwf")
    tmrca = numpy.zeros(migr_mat.shape)
    for locus in sim:
       for deme1 in range(num_popul):
         for deme2 in range(num_popul):
           tmrca[deme1,deme2] += locus.at(0.).tmrca(deme1, deme2) 
    return tmrca / reps

# random "spatial" migration matrix
def simulate_migration_matrix(seed, dim):
  numpy.random.seed(seed)
  coords = numpy.random.uniform(0, 1, (dim, 2))
  distance = dist.squareform(dist.pdist(coords))
  return numpy.exp(-distance**2) - numpy.eye(dim)

# get average of TMRCA involving a given tip, across permutations of spatial locations
def average_TMRCA_across_permutations (seed, migr_mat, num_perm, num_loci):
  indices = range(migr_mat.shape[0])
  avg_tmrca = numpy.zeros([num_perm, migr_mat.shape[0]])
  for perm in range(num_perm):
    numpy.random.seed(seed)
    permuted_indices = numpy.random.permutation(indices)
    permuted_migr_mat = migr_mat[:,permuted_indices]
    permuted_migr_mat = permuted_migr_mat[permuted_indices,:]
    avg_tmrca[perm,:] = numpy.mean(island_model_TMRCA(permuted_migr_mat, num_loci, seed), 1)
    seed += 10000
  return avg_tmrca

# density plots of average TMRCA with other tips, across permutations
migr_mat = simulate_migration_matrix(seed = 1, dim = 10)
df = pd.DataFrame(average_TMRCA_across_permutations(seed = 1, migr_mat = migr_mat, num_perm = 50, num_loci = 50))
fig, ax = plt.subplots(1,1)
for deme in df.columns:
  if deme == 0: col = 'red'
  else: col = 'black'
  df[deme].plot(kind='density', color=col, label="tip"+str(deme))

ax.set_xlabel("Pairwise TMRCA averaged over other tips")
ax.set_ylabel("Density across permutations of spatial location")
fig.legend(loc="best")

fig.savefig("tmrca_density.png", dpi=199)

I get similar behavior if I allow islands to merge at some point in the past, or sample multiple haploids per deme (in which case, all samples from the first deme have "elevated" TMRCA with other tips.)

[petrelharp]

Thanks for the report, @nspope - I can confirm that this happens only with dtwf, not with the standard coalescent, which sure seems like a dtwf-specific bug. (The problem shows up with many fewer reps; I've edited your report so the code runs faster.)

[petrelharp]

Here's a more minimal version of your (very nice) example, that also shows something else very strange:

import msprime
import numpy

def island_model_TMRCA (rate, reps, seed, num_popul):
    migr_mat = rate * (numpy.ones((num_popul, num_popul)) - numpy.eye(num_popul))
    pop_confg  = []
    for i in range(num_popul):
      pop_confg += [msprime.PopulationConfiguration(
        sample_size=1, initial_size=1.)]
    sim = msprime.simulate(
            population_configurations = pop_confg,
            migration_matrix = migr_mat,
            num_replicates = reps,
            random_seed = seed,
            model = "dtwf")
    tmrca = numpy.zeros(migr_mat.shape)
    for locus in sim:
       for deme1 in range(num_popul):
         for deme2 in range(num_popul):
           t = locus.first()
           tmrca[deme1,deme2] += t.tmrca(deme1, deme2) 
    return tmrca / reps

for r in [0.9, 0.99, 0.999]:
    avg_tmrcas = island_model_TMRCA(rate=r, reps=100, seed=2, num_popul=3)
    print(f"rate = {r}")
    print(avg_tmrcas)

This produces:

rate = 0.9
[[ 0.   14.86 13.96]
 [14.86  0.    5.55]
 [13.96  5.55  0.  ]]
rate = 0.99
[[  0.   100.01  95.64]
 [100.01   0.     7.24]
 [ 95.64   7.24   0.  ]]
rate = 0.999
[[   0.   1012.29 1009.97]
 [1012.29    0.      4.87]
 [1009.97    4.87    0.  ]]

Two things are odd: (a) the migration matrix is symmetric, so the numbers should all be roughly the same, and (b) they should not blow up as the migration rate approaches 1.0. (and increasing the rate to 0.9999 makes this take quite a long time).

I think this is definately a bug in the migration code for dtwf. @DomNelson, do you feel like having a look at this? Or know who might?

[petrelharp]

Looking at the code, I don't see anything obviously wrong, but changing the loop here to go in the opposite direction produces a segfault, which seems suspicious.

[jeromekelleher]

@ivan-krukov, have you any idea what might be happening here?

[nspope]

I have an idea of what is going on here ... It is foolery on my part (of course) but also highlights the need for an internal check.

The migration matrix should be right stochastic. In my example and the modification by @petrelharp, this is not going to be the case if the average migration rates per row are high. When making my example, I was thinking about the migration rate matrix as the infintesimal generator of a continuous time Markov process, but of course that is not actually how it is used.

What happens internally is that the empty diagonal is "filled in" with one minus the sum of offdiagonal entries per row:
https://github.com/tskit-dev/msprime/blob/389f4862bcb4fc08089e5f8a620363844255be10/lib/msprime.c#L4386-L4395

So we end up with negative numbers on the diagonal, of much larger magnitude than those on the offdiagonal. Whatever gsl_ran_multinomial does internally when an element of the probability vector is negative (I haven't bothered to check), it definitely does not generate variates from the intended multinomial: regardless of which index of the probability vector is negative, it seems the first element of the output vector always has a higher empirical weight ... (perhaps it is simulating using series of conditional binomials; and the first in the series gets messed up when there's a negative?)

[nspope]

... I was a bit off but it is clear that the problem is with gsl_ran_multinomial and negative "probabilities", so the conclusion still stands. Here is a minimal example:

// gsl_ran_multinom_issue.c
#include <stdio.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>

int main (void)
{
  const gsl_rng_type *rng_type;
  gsl_rng *rng;

  gsl_rng_env_setup();

  rng_type = gsl_rng_default;
  rng = gsl_rng_alloc (rng_type);

  double prob = 0.9;
  unsigned int draws = 1000;
  unsigned int dim = 3;
  unsigned int n[dim];
  double p[dim];

  for(unsigned i=0;i<dim;++i)
  {
    printf("pop %u,\tprobs:", i);
    for(unsigned j=0;j<dim;++j) 
    {
      p[j] = i==j ? 1.-prob*(dim-1) : prob;
      printf("\t% 5.3f",p[j]);
    }
    gsl_ran_multinomial(rng, dim, draws, p, n);
    printf(",\tempirical prop:");
    for(unsigned j=0;j<dim;++j) 
      printf("\t% 5.3f",float(n[j])/float(draws));
    printf("\n");
  }

  gsl_rng_free(rng);
}

and ...

g++ -Wall gsl_ran_multinom_issue.c -o gsl_ran_multinom_issue.out -lgsl -lblas && \
    GSL_RNG_SEED=1 ./gsl_ran_multinom_issue.out
# pop 0,  probs:  -0.800   0.900   0.900, empirical prop:  0.000   0.487   0.513
# pop 1,  probs:   0.900  -0.800   0.900, empirical prop:  0.901   0.000   0.099
# pop 2,  probs:   0.900   0.900  -0.800, empirical prop:  0.899   0.101   0.000

... and I still suspect the reason that the first index is consistently different, is that the GSL's implementation of the multinomial RNG uses a series of conditional binomials.

[petrelharp]

Whoo! Yes! Good call, @nspope - I inserted

            assert(mig_tmp[j] >= 0.0);

right after this line and ran my script above, and - indeed - the assert was triggered. So I guess we just need to catch this problem earlier.

What to do here is a little unclear, though: as you point out, the standard coalescent is a continuous-time Markov process, parameterized by rates, while the DTWF is discrete-time. This contradicts what [the docs}(https://msprime.readthedocs.io/en/latest/api.html#discrete-time-wright-fisher) say: "with all parameters denoting the same quantities in both models".

One way to deal with this would be to take the rate matrix G and use it to compute the probability matrix expm(G) (matrix exponential), and use those for the probabilities. I don't think anyone would really appreciate that, though, since they'd actually like to specify the usual WF migration parameter directly, and don't want migrants moving directly from i->j when G[i,j] = 0.

So, maybe we just check if the probabilies for any population sum to more than 1? Summing to anything more than 1/2 would be really wierd, but I don't see a reason to disallow it. What do you think, @nspope?

And, we should update the DTWF docs to reflect the difference.

[nspope]

@petrelharp Yeah, I definitely don't think it's intuitive to always interpret the migration rate matrix as an infinitesimal generator and convert it to transition probabilities.

I agree: just check that the input matrix is right stochastic, and update the docs accordingly. I don't think the docs are too far off. I mean, biologically the elements of the migration matrix are interpreted the same regardless of whether it's continuous time or discrete time ... That is, "the [j, k]th element giving the fraction of population j that consists of migrants from population k in each generation" (from the current API docs) applies just fine to discrete time. It's just that a limit is being taken in the continuous time case.

I actually don't think it's unreasonable to want to simulate from the DTWF with high overall migration ... for example, exploring behavior that is close to panmixia or trying to approximate continuous spatial structure on a lattice.

[jeromekelleher]

Thanks so much for following up on this @nspope, this is hugely helpful! A pull request implementing a check would be yet more helpful again :wink:

[jeromekelleher]

I should say also @nspope - this would be a significant bug fix, so would certainly qualify for authorship on the forthcoming 1.0 paper: https://github.com/tskit-dev/msprime-1.0-paper/issues/2

[nspope]

Hi @jeromekelleher, sure, I'd be happy to put together a PR implementing checks and updating documentation!

[nspope]

resolved with #1130