Likelihood surfaces

Note

Advanced. Most analyses never need these tasks — they exist to inspect the likelihood directly rather than to sample from the posterior. If you want an answer, use infer.

Four tasks evaluate the log-likelihood at fixed parameter values and write it out, instead of running an MCMC:

Task	Varies
calculateLL	nothing — one value per locus at the current parameters
calculateLLForS	the selection coefficient, over its whole grid
calculateLLForSH	selection and dominance, over the full 2-D grid
calculateLLForN	the population size, over a grid you specify

The commonest use by far is the last one: seeing where your data put N before choosing a prior for it.

What they share

All four take the same data and model arguments as infer — --counts, --meta, --mu_a_A, the grids, the priors, --linkage, and so on. Three things behave differently, and all three matter:

Important

They evaluate at each parameter’s initial value. No MCMC runs, so nothing is estimated. Whatever --N, --m, --s, --h, --betaMean … are set to is what the likelihood is evaluated at. Setting a starting value is normally a minor detail; here it is the whole input.

Each task writes its own file, named after it, so the four can share an --out without overwriting one another:

Task	Output
calculateLL	<out>_LL.txt
calculateLLForS	<out>_LLForS.txt
calculateLLForSH	<out>_LLForSH.txt
calculateLLForN	<out>_LLForN.txt

Note

The values are shifted. Emissions are binomial densities taken relative to their maximum, i.e. log10 P(c|n,x) − log10 P(c|n,c/n). Every reported log-likelihood therefore differs from the true one by a constant that depends on the locus. Differences within a locus — across s, h or N — are meaningful, which is all these tasks are for. Absolute values are not, and are not comparable to other software.

calculateLL

One log-likelihood per locus at the current parameter values, plus their sum in the log.

$ ./sweeplink calculateLL --counts counts.txt --meta meta.txt --mu_a_A 1e-8 --out ll

Output ll_LL.txt:

chr	pos	LL

--writeBeta

Also dump the forward/backward diffusion vectors at every intermediate time step to <out>_betas.txt. Debugging output — large, and written per locus.

calculateLLForS

For each locus, the log-likelihood at every value on the grid of selection coefficients, holding everything else fixed.

$ ./sweeplink calculateLLForS --counts counts.txt --meta meta.txt \
      --mu_a_A 1e-8 --numStatesAbsS 8 --out lls

One column per grid state, named with the value it stands for:

chr	pos	LL0_s_-0.0741	LL1_s_-0.0714	...	LL17_s_0.08

Remember the grid is symmetric and contains two zeros — see how it works.

--popID <n>

Which population’s s to vary. Default 0.

calculateLLForSH

The full two-dimensional surface: every combination of dominance and selection, per locus. Columns run with h outer and s inner, giving numStatesH × 2 × (numStatesAbsS + 1) of them.

chr	pos	LL0_s_-0.0741_h_0	LL1_s_-0.0714_h_0	...

--popID <n>

Which population’s s to vary. Default 0.

calculateLLForN

The log-likelihood over a grid of population sizes — the one worth reaching for routinely, because it shows where the data alone put N, independently of the prior.

$ ./sweeplink calculateLLForN --counts counts.txt --meta meta.txt --mu_a_A 1e-8 \
      --minN 100 --maxN 1e7 --numN 100 --logN --sumLL --out probeN

--minN <value> (default 1), --maxN <value> (default 1000000)

Ends of the grid.

--numN <n> (default 100)

Number of grid points. Each one rebuilds the diffusion for that population, so runtime scales with this.

--logN

Space the grid evenly on the log10 scale rather than linearly. Almost always what you want, since plausible values of N span orders of magnitude.

--sumLL

Write one row summed over all loci instead of one row per locus. For choosing N this is the row you want.

--popID <n>

Which population’s N to vary. Default 0.

With --sumLL the output is a single row whose first two fields are sumLL and -:

chr	pos	LL0_N_1.000000	LL1_N_1.149757	...
sumLL	-	-4954.27	-4953.52	...

The largest value marks the best-supported population size. Use it to sanity-check the prior you set with --log10N.alpha / --log10N.lambda — see the FAQ — not to replace it.

Note

--logN (the flag above) and log10N (the model parameter) are different things with unhelpfully similar names. The flag only controls the spacing of this task’s grid.