Joint endpoint posteriors

Joint endpoint posterior probabilities

CSUBST 1.16.0 makes joint endpoint posterior probabilities (joint PP) the default for search/analyze, sites, inspect, benchmark and scan. This changes the estimator of substitution support, so counts, candidate rankings and downstream statistics can differ from earlier runs. It is a scientific-method change as well as an implementation change.

What joint PP means

For a branch from parent node p to child node c, joint PP is

J(a,b) = P(X_p=a, X_c=b | D, fitted model, rooted tree).

Here D contains the observations at all tips for the site. The probability describes the two endpoints together, retaining their dependence along the branch. CSUBST computes it by tree pruning and parent-to-child conditioning; it does not sample ancestral histories.

The previous marginal estimator multiplied separately inferred probabilities:

M(a,b) = P(X_p=a | D) × P(X_c=b | D).

That product loses the dependence between adjacent nodes. Consider a zero-length branch whose endpoints are each uncertain between A and B with probabilities 0.5 and 0.5:

Endpoint pair	Marginal product	Joint PP
A → A	0.25	0.50
A → B	0.25	0
B → A	0.25	0
B → B	0.25	0.50

The marginal product assigns total change support 0.5 despite the branch having zero length. Joint PP assigns zero change support. For nonzero branches, the difference depends on the observations, fitted model and branch length; joint PP does not imply a uniform reduction in every statistic.

How CSUBST uses it

Off-diagonal codon endpoint pairs are classified into synonymous (S) and nonsynonymous (N) changes. Amino-acid and recoded-state events sum the relevant codon-pair probabilities over both endpoint axes. Recoding does not require treating the reduced alphabet as a separate Markov process.

For discrete codon-rate categories, CSUBST mixes category-conditional results using P(category | D) for the complete site. It does not substitute a posterior-mean rate. event_pp means an endpoint-pair probability conditional on the data and fitted model. Event thresholds remove probability mass without renormalizing it.

With --expectation_method codon_model, the expected endpoint weight is

sum_r P(r | D) P(X_p=a | D,r) exp(Q × rate_r × fitted_length)[a,b].

This fitted conditional prediction retains the parent state in the transition calculation and uses the fitted model lengths. It is distinct from the observed joint PP. The urn expectation instead uses the new observed weights with its existing count-based null.

Running joint or marginal analyses

Joint is the default; making it explicit helps document a run:

csubst search --alignment_file alignment.fa --rooted_tree_file tree.nwk \
  --substitution_posterior joint --endpoint_block_size 64

To retain the previous endpoint estimator for a comparison:

csubst search --alignment_file alignment.fa --rooted_tree_file tree.nwk \
  --substitution_posterior marginal

For scan, --scan_observation overrides the common posterior option and also accepts the separate CTMC bridge method:

csubst scan --alignment_file alignment.fa --rooted_tree_file tree.nwk \
  --foreground foreground.txt --scan_observation joint

Joint/bridge scan resolves unspecified exposure and length options to endpoint and raw. Marginal scan resolves them to q_weighted and n_rescaled. Remove old explicit exposure/length settings when switching modes, or select a compatible pair; incompatible settings fail. See the scan CTMC guide for bridge support and calibration options.

Inputs and migration

• Joint inference recomputes ancestral uncertainty from tip observations and the fitted model. It requires a rooted tree and compatible IQ-TREE model files. --ml_anc yes is incompatible.
• Supported codon models include ECMK07, ECMrest, GY, MG and MGK with uniform rates or a supported discrete rate-category table. MG/MGK require the fitted F1X4/F3X4 nucleotide-frequency context. Mixtures of different Q matrices and ascertainment-corrected models are rejected, without silently switching to marginal estimation.
• Native 3Di joint inference uses its independent fitted uniform GTR context and requires direct ASR. Joint/bridge scan does not currently support native 3Di; select the supported marginal scan route for that use case.
• Missing observations contribute neutral likelihood during inference, but branches with an unobserved child subtree are excluded from reporting at that site. Tables distinguish unavailable observations (NA) from a real zero and report eligible observation counts.
• The inserted root is inferred from the model. Root-adjacent branches can become eligible even when no root row existed in the imported ASR file.
• *_endpoint_model.json records fitted-model provenance, generator hash, rates, lengths and storage choices. Retain it with the command, CSUBST version, inputs and result tables when comparing analyses.
• marginal restores the estimator, not every historical bug. In particular, MG/MGK frequency weighting is corrected in both modes. To reproduce a historical release exactly, retain that release and its input/model files.
• To evaluate exported true ancestral states from a simulation, use marginal mode. To evaluate joint inference on simulated tips, refit those tips with IQ-TREE; joint mode does not use supplied true internal states as truth.

Memory, runtime and higher arities

Inference works in site blocks controlled by --endpoint_block_size. Smaller blocks reduce the main inference workspace, but total RAM also includes state arrays, retained event information and downstream tables. Requested statistics and outputs determine how much information must remain available.

For supported arity-2 searches, CSUBST can accumulate branch-pair results directly and discard each block's detailed projections within a bounded workspace. That route is disabled when max_arity > 2. Higher-order searches still benefit from shared block processing and reuse of expected reducers, but retain the projections or full events needed for later arities. In particular, spe2spe requires full events outside the bounded pairwise route.

Joint scan can keep large event/category arrays in temporary site files; CSUBST_TMPDIR selects the temporary directory. Selected-branch sites analyses can retain detailed events only for requested branches while preserving global totals. These optimizations preserve probability mass and requested outputs; they do not approximate joint PP by dropping small probabilities. Disk-backed storage trades temporary disk space and I/O for lower RAM.

Measured speedups and RAM savings depend on the workload. Scan improvements must not be extrapolated to high-arity search. See the memory measurements and additional scan/sites measurements for before/after conditions, output-equivalence checks and limitations.

Interpretation limits

Joint PP is exact for a single edge under the supplied fitted model, subject to numerical precision. It is conditional on estimated parameters; it does not integrate uncertainty in the tree or model fit.

CB/CS/CBS still multiply per-branch scores. Those products are not a joint posterior across multiple branches. Endpoint changes also do not count all substitution histories: multiple hits and reversions can have unchanged endpoints. CTMC bridge counts and stochastic mapping answer different questions.

Joint PP alone does not establish calibrated omegaC P values or biological false-positive rates. Posterior thresholding and candidate discovery introduce selection; their calibration needs to be evaluated separately.

For detailed model contracts, equations, implementation and reproducible benchmarks, see the technical endpoint guide.