Differences Between WF and Hudson's Model Outputs

[Rhinezly]

I have encountered some unexpected results when comparing temporal statistical outputs between different simulation models. Specifically, using a constant population size of 30,000, the Standard Coalescent (Hudson’s model) produces nearly constant statistics over time, while the Discrete Time Wright-Fisher (DTWF) model diverges from Hudson's after a couple of generations. I then tested other population sizes as well and observed similar patterns, though it seems that as population size increases, the DTWF model results become more similar to those from Hudson's. The plots below are from 10 replicates with 95% confidence intervals.
statistics over time.pdf

Here's a sample script using the DTWF model. The only change made to run the Hudson model was replacing msprime.DiscreteTimeWrightFisher() with msprime.StandardCoalescent(). I’m unsure how to interpret the differences and would greatly appreciate any insight or guidance.

import msprime

Ne = 30000
seqlength = 1e6
theta = 0.0045
mu = theta / (4 * Ne)
recombination_rate = 10 * mu
timepoints = list(range(0, 301, 3))  # Time points from 0 to 300 in steps of 3
samples = [msprime.SampleSet(50, time=t) for t in timepoints]

ts = msprime.sim_ancestry(
    samples=samples,
    population_size=Ne,
    sequence_length=seqlength,
    recombination_rate=recombination_rate,
    model=[msprime.DiscreteTimeWrightFisher()]
)

mts = msprime.sim_mutations(
    ts,
    rate=mu
)

[jeromekelleher]

Thanks for this @Rhinezly, it's not obvious to me why taking samples at different time points should matter here. It think it's sufficient to look at the number of mutations plot. Would you mind sharing the code used to compute that, and placing the plots inline here as images, so it's easier for people to view them?

[Rhinezly]

Hi @jeromekelleher, I was going to explore how population decline will affect the statistics through time, thus the temporal samples. But the results seemed a bit weird so I tried to find out the reason and figured it should be caused by the models. The constant scenario is just a clearer way to show how the two models differ.

Here's the code for the statistics. I've output the tree files from repeated simulations before this.

import tskit
import pandas as pd
import glob

# Get the prefix of the tree file
def get_tree_prefixes():
    prefixes = []
    for f in glob.glob("*.trees"):
        prefixes.append(f.replace(".trees", ""))
    return prefixes

prefix = get_tree_prefixes()

# Calculate summary statistics
def oneway_stats(ts):
    stats = {
        "Number of trees": ts.num_trees,
        "Number of mutations": ts.num_mutations,
        "Density of segregating sites": ts.segregating_sites(),
        "Nucleotide diversity (pi)": ts.diversity(),
        "Tajima's D": ts.Tajimas_D()
    }
    return stats

def oneway_output(ts, prefix):
    results = {}
    for t in timepoints:
        # Extract samples for the current time point
        sample_ids = [n.id for n in ts.nodes() if n.time == t]
        sub_ts = ts.simplify(samples=sample_ids)

        # Calculate statistics
        stats = oneway_stats(sub_ts)
        results[t] = stats

    # Convert results to a DataFrame
    oneway_df = pd.DataFrame(results).T
    oneway_df.index.name = "t"
    oneway_df.to_csv(f"{prefix}.csv")

def process_trees_file():
    for prefix in get_tree_prefixes():
        ts = tskit.load(f"{prefix}.trees")
        oneway_output(ts, prefix)

if __name__ == "__main__":
    process_trees_file()

[petrelharp]

Perhaps the problem is here:

         sample_ids = [n.id for n in ts.nodes() if n.time == t]

In 'hudson', the only nodes whose times are exactly equal to the sampling times will indeed be the sample nodes, but in dtwf (since it's in discrete time), there will be other nodes whose time is exactly equal to those times. Could you swap out this line for

        sample_ids = ts.samples(time=t)

and see if that changes things?

[Rhinezly]

The results from two models are quite similar now, thank you so much!

[jeromekelleher]

Great spot, thanks @petrelharp!