05 Model selection: marginal likelihoods and heating

Status: draft (2026-11-12), sections 1-9. Run with migrate-n 7.0.110 and 7.0.111 (7.0.111 fixes a crash in the warning list that hung some runs of the wrong models; the numbers are the same).

Which migration model fits the data better: one with migration in both directions, one with equal rates, or one with migration in one direction only? Migrate answers this with the marginal likelihood of each model and Bayes factors. This recipe shows where the marginal likelihoods are in the outfile, checks the four estimators against an exact value, and then asks how often they pick the true model in simulated data, and how many heated chains that needs. Read 00-basics and 01-migration first.

Contents

1. Marginal likelihoods and Bayes factors
2. Marginal likelihoods in the parmfile and the outfile
3. An exact test
4. Larger data: how many chains?
5. Four migration models, 5 loci
6. Four migration models, 20 loci
7. More chains, a narrower prior
8. Model probabilities
9. What next

1. Marginal likelihoods and Bayes factors

The marginal likelihood of a model is the probability of the data under the model, averaged over the prior of its parameters (and over all genealogies):

P(D | model) = integral of P(D | parameters) p(parameters) d parameters

Migrate reports its logarithm, log P(D | model), summed over loci. The model with the higher value fits better. The difference of two log marginal likelihoods is the log Bayes factor:

log BF = log P(D | model A) - log P(D | model B)

A model with more parameters can fit the data better, but it pays for each extra parameter through the prior: the average over a wide prior includes many values that fit badly. The marginal likelihood therefore depends on the prior, not only on the data (section 7).

A common reading of Bayes factors (Kass and Raftery 1995) uses 2 log BF: 0-2 not worth more than a mention, 2-6 positive, 6-10 strong, above 10 very strong evidence. In log units: a difference below 1 is a tie, 1-3 positive, 3-5 strong, above 5 very strong evidence.

2. Marginal likelihoods in the parmfile and the outfile

Migrate computes the marginal likelihood from the heated chains of a Bayesian run. Every run with heating=YES reports it; nothing else has to be switched on. The heated chains sample from P(D | parameters)^beta p(parameters) with beta = 1/temperature, from beta = 1 (the cold chain) to almost 0 (temperature 1,000,000). The more chains, the better the estimate. Two kinds of ladders:

POWER:N sets beta = x^3.33 for x evenly spaced between 1 and 0, which puts more chains at small beta, where the mean log likelihood changes fastest; LINEAR:N spaces beta itself evenly. The recipes use POWER.

Each chain costs about as much as the cold chain, so 16 chains take about four times as long as 4 (sections 4 and 7 give the times). Under MPI the chains of one locus run on one worker.

The table “Log-Probability of the data given the model” in the outfile has one line per locus and two lines for the sum (here for a 20-locus run of section 7):

Locus          TI(1a)       BTI(1b)         SS(2)         HS(3)
---------------------------------------------------------------
      1      -2283.57      -2283.56      -2283.10      -2248.55
      2      -2145.86      -2145.86      -2145.31      -2102.66
...
     20      -1973.26      -1973.26      -1972.89      -1928.46
---------------------------------------------------------------
  All       -42882.66     -42882.64     -42872.67     -41996.88   (independent parameters)
  All       -42879.11     -42879.08     -42869.11     -41993.32   (joint)

Use the “All … (joint)” line for model comparison: it adds the per-locus values with the joint parameters of the multi-locus combination (01-migration section 7), so that all loci share one set of parameters, as in the model. The “(independent parameters)” line is the older combination, which treats the parameters of the loci as independent; it can be biased (it gave a spurious preference for an extra parameter in our tests), use it only to compare with older runs.

3. An exact test

How good are the four estimators? For one population, one locus and few sequences the marginal likelihood can be computed directly: draw many coalescent genealogies, compute the data likelihood for every Theta on a fine grid, average, and integrate over the prior. scripts/exact_ml.py does this for data/small.txt (6 sequences, 80 sites, 7 segregating; simulated with Theta 0.02) under the Jukes-Cantor model and the default prior Theta ~ U(0, 0.1):

python3 scripts/exact_ml.py data/small.txt 200000
sequences 6, patterns 11, sites 80, segregating 7
log P(D) = -159.5613   (4 batches: -159.5288 -159.6101 -159.5195 -159.5900, sd 0.0448)
posterior Theta quantiles 2.5/25/50/75/97.5%: 0.0170 0.0360 0.0518 0.0708 0.0963

The posterior agrees with migrate without heating (parmfiles/parmfile.exact, 2.5/50/97.5%: 0.0169, 0.0510, 0.0962), so both compute the same thing. scripts/run_exact.sh runs migrate with five ladders and five random seeds (about 6 s per run with 4 chains, 66 s with 32). Mean over the seeds (the seeds differ by less than 0.05, HS by up to 1):

SS is right with every ladder. TI and BTI approach the exact value as the number of chains grows; with 4 chains TI is 4 units too low and BTI 4 too high. HS is 3 units too high whatever the ladder.

4. Larger data: how many chains?

With more data the log likelihood changes by much more between the hot and the cold chain, and the estimators need more chains. The dna-2pop data of 00-basics (two populations, 5 loci, 10 sequences each), full model, run with 6 MPI processes; there is no exact value, so the reference is the mean of TI, BTI and SS with 32 chains (-4624.5, two seeds within 0.5):

(Differences from the reference; negative = too low.) With the default four chains every estimator is off by hundreds of log units. From 8 chains on, SS is within about one unit; TI and BTI need 16 or more. The error does not cancel automatically between models: models with different parameters have different heated chains. Whether a short ladder still picks the right model is the question of the next sections.

5. Four migration models, 5 loci

study/make_study.py simulates 20 data sets with the dna-2pop truth (Theta_1 = 0.01, Theta_2 = 0.005, M_2->1 = 100, M_1->2 = 400; 5 loci of 500 bp; 10 sequences from North and South) and writes parmfiles for four models (rows are the receiving populations, see 01-migration) and three ladders:

python3 study/make_study.py study5        # from cookbooks/05-model-selection
# run every study5/<ladder>/<model>/repNN/parmfile, then
python3 scripts/compare.py study5

How often each model won (out of 20 data sets):

The true model rarely wins. With POWER:8 the three estimators agree, and all models are within a few log units: the full model loses to the best other model by 1.6 (median; range -4 to +0.5). Five loci with Theta 0.01 and 0.005 have few differences per locus; they do not tell the four models apart, and the model with fewer parameters, south-only, wins most often: the evidence for a migration rate of 100 into North is weak, and a simpler model that fits about as well is preferred. Section 7 shows that part of this preference comes from the ladder.

With 4 chains the estimators disagree: TI favours the restricted models by 26 units (median), SS scatters by up to 70 units in either direction. BTI with 4 chains is close to the POWER:8 values.

6. Four migration models, 20 loci

More informative data: the same migration rates, but Theta_1 = 0.02, Theta_2 = 0.01, 20 loci of 1000 bp:

python3 study/make_study.py study20 --loci 20 --sites 1000 --theta 0.02 0.01

The runs used 21 MPI processes (one per locus) and took 2-13 minutes each on a cluster.

North-only (no immigration into South, where the true rate is 400) now loses clearly, by about 9.5 units with POWER:8. But south-only still wins most data sets, although the full model estimates M_2->1 well: the joint posterior median lies between 39 and 372 (truth 100), and the lower bound of the 95% interval is above 0 in 17 of 20 data sets (but below 10 in six of them). With POWER:5, BTI prefers south-only in all 20 data sets by a margin of 11 units (median), which looks decisive.

7. More chains, a narrower prior

Two things can make the true model lose: the ladder (the estimators have not converged), and the prior (a wide prior on M_2->1 costs the full model more). Five of the data sets, full and south-only only, with 16 and 32 chains, and with the default prior on M (uniform 0-1000) or a narrower one (uniform 0-500):

python3 study/make_study.py check --loci 20 --sites 1000 --theta 0.02 0.01 \
    --reps 5 --ladders P8 P16 P32 --models full south-only
python3 study/make_study.py check-m500 --loci 20 --sites 1000 --theta 0.02 0.01 \
    --reps 5 --ladders P8 P16 P32 --models full south-only --m-prior 500
python3 scripts/compare.py check full south-only

Full minus south-only (BTI; positive = the true model wins):

Run time per data set and model with 21 MPI processes: 14 min with 16 chains, 29 min with 32; with the narrower prior 5 min (POWER:8), 11 min (POWER:16) and 12-22 min (POWER:32).

The ladder. The preference for south-only shrinks with every doubling of the chains, and with 32 chains the two models tie (mean -0.3). TI and BTI agree from POWER:8 on, and SS gives the same differences with 32 chains (within 0.4). The large margins of sections 5 and 6 were made by the short ladders. In these data the error of a short ladder favours the model with fewer parameters; that is a conservative error, but it is not a rule: in 03-growth (section 5) the short ladder favoured the growth model, the one with more parameters. The absolute values have still not converged with 32 chains (rep00: -42879 with 32 chains, -42914 with 16), only the differences between the models have.

The prior. With M uniform on 0-500 the full model wins all five data sets with 16 and with 32 chains (mean +3.8 and +4.5): more chains do not remove this preference. Halving the range of M_2->1 would add only log 2 = 0.7 to the full model. Most of the change comes from the other rate: without M_2->1, the south-only model moves M_1->2 up to make up for the missing migration (posterior medians 550-980 with the wide prior, truth 400). The narrower prior stops it at 500 (medians 490), and the south-only model fits worse. The full model’s M_1->2 (medians 320-490) is also pressed against the bound, but it needs less of it.

So with 20 loci and these migration rates, the data cannot tell whether M_2->1 is 100 or 0 once the estimators have converged; they can if the prior excludes the high M_1->2 values that the simpler model needs. Choose the prior from what is plausible for the organism, not to make a model win, and report it with the Bayes factors.

8. Model probabilities

With several models, the log marginal likelihoods give each model a probability (all models equally likely a priori):

P(model i | D) = exp(lml_i - lml_max) / sum over j of exp(lml_j - lml_max)

where lml_max is the largest value. For two models with a difference d = log P(D | A) - log P(D | B): P(A | D) = 1 / (1 + exp(-d)). For the five data sets of section 7 with 32 chains:

A difference of 3 log units gives the better model 95%, a difference of 5 more than 99%. Differences of 10 or more, as with the short ladders in sections 5 and 6, would look like certainty, so check that the ladder is long enough before reading them that way.

Practical rules from this recipe:

• Use BTI or SS from the “All … (joint)” line, never HS.
• Run the competing models with the same data, the same priors (as far as the models share parameters) and the same ladder.
• Repeat the comparison with twice the chains. If the order of the models stays the same and the margin is large, the answer holds; if the margin shrinks, keep doubling until it is stable. With 20 loci these data needed 16-32 chains; with 5 loci 8 were enough to show that the models tie.
• A margin of a few log units is not a decision; report the model probabilities.
• The prior is part of the model: a wider prior lowers the marginal likelihood of the model with the extra parameter, and a prior that truncates a parameter can penalize the model that needs it.

9. What next

• 06-assignment: assigning individuals to populations, which uses the same marginal likelihoods to compare assignments.
• 08-mittag-leffler: comparing the fractional coalescent with the Kingman coalescent.
• The constant-size versus growth comparison of 03-growth (section 5) and the number of skyline segments in 04-skyline (sections 5 and 7) are model choices of the same kind.