Posterior Predictive Checks for Branching Models

Why Check a Branching Model Predictively?

A branching model can produce plausible lineage curves and still fail to reproduce important features of the observed cells. It might predict too few cells on one branch, imply a smooth transition where the data show a sharp expression difference, or miss the greater variability seen in one lineage. Posterior predictive checks ask whether data simulated from the fitted model resemble the data that motivated the fit, using features chosen to reveal these failures.

The check is only as meaningful as its replication scheme. A simulated data set must represent the same kind of observations as the real data: for example, the same number of cells, the same observed sampling scheme, and the same expression scale. If a model distinguishes biological branch probabilities from unequal cell sampling, the check must use its sampling model rather than treating the observed cell counts as direct fate frequencies.

Posterior Predictive Distributions

Let \(D\) denote the observed data, including cell-level expression measurements and any recorded trajectory information. Let \(\theta\) collect model parameters, such as branch probabilities, split time, trajectory coefficients, and residual variances. Let \(z_i\) be the branch label for cell \(i\), which may be latent. A fitted Bayesian model gives a posterior distribution \(p(\theta,z\mid D)\). A posterior predictive data set \(D^{\mathrm{rep}}\) is generated by first drawing parameters and latent quantities from their posterior distribution and then simulating new observations from the model.

For each draw \(m=1,\ldots,M\), one can draw \(\theta^{(m)}\) from \(p(\theta\mid D)\), generate \(D^{\mathrm{rep},(m)}\) under the observation model, and calculate both discrepancies using that draw. A posterior predictive tail proportion for an upper-tail check is

This proportion describes how often the chosen replicated discrepancy is at least as large as the observed discrepancy under the posterior predictive distribution. It is a diagnostic, not a probability that the model is true or false. A value near 0 or 1 can flag a mismatch in the direction measured by \(T\); a value near the middle does not establish that every aspect of the model is adequate. The choice of discrepancy matters.

What Must Be Replicated?

A typical branching expression model specifies a trajectory for each branch and a distribution around it. For a transformed expression measurement \(Y_{ig}\) for cell \(i\) and gene \(g\), one possible model is

Here \(t_i\) is the cell's trajectory position, \(z_i\) is its branch, \(f_{kg}\) is the branch-specific mean expression curve, and \(\sigma_{kg}^{2}\) describes residual variation. This example assumes that the stated conditional distribution is appropriate for the expression scale, and that the model includes the relevant dependence on trajectory position and branch. The normal model is not universally suitable for single-cell measurements; the check must use the likelihood and scale actually adopted for the data.

For branch occupancy, a simple model might generate \(N\) observed cells with branch probabilities \(\pi_1,\ldots,\pi_K\), so that the vector of branch counts is multinomial conditional on \(N\). This assumes the probabilities describe the branch labels in the sampled cells. If sampling differs by branch, the model must account for that difference. A check that ignores ascertainment or unequal sampling can flag the sampling design rather than the lineage model.

For soft lineage assignments, do not automatically treat posterior assignment weights as if they were observed integer counts. Replicated labels should be generated or integrated over according to the same latent-label model used in fitting. If branch labels are known from an external measurement, preserve that information in the replicated-data procedure. In either case, compare like with like.

Worked Example: Checking Branch Occupancy

Checking Expression Continuity Near a Split

A model that shares an ancestral trajectory and branches continuously at split time \(\tau\) implies that the branch mean curves meet there. For a gene \(g\), this constraint can be written \(f_{1g}(\tau;\theta)=f_{2g}(\tau;\theta)\). The equality concerns the modeled mean curves, not the individual expression values: cells near a split can differ because of residual variation, different positions relative to the split, or measurement noise.

A useful check therefore compares an observed-data discrepancy with its replicated distribution. For example, choose a narrow interval around \(\tau\), calculate the absolute difference between the branch-specific mean expression values in that interval, and calculate the same statistic for each replicated data set at the same trajectory positions. A large observed difference relative to replicated differences may indicate that the continuity constraint, the local curve shape, or the residual distribution is inadequate. The check does not by itself identify which part is responsible.

Checking Branch-Specific Variation

Mean curves alone do not determine whether a model reproduces the spread of expression around those curves. For each branch, residuals can be calculated as observed expression minus the corresponding fitted mean. A discrepancy such as the ratio of the two branch-specific residual variances can reveal whether a model with a shared variance is compatible with the observed variation. Replicates must use the same residual definition and variance calculation as the observed data.

Justification, Diagnostics, and Alternatives

The posterior predictive distribution follows directly from the model's two-stage generative structure. Given parameters \(\theta\), the model assigns a probability \(p(D^{\mathrm{rep}}\mid\theta)\) to replicated data. Because \(\theta\) is uncertain after observing \(D\), average this sampling distribution over its posterior:

In simulation, the integral is approximated by posterior draws and corresponding replicated data sets. For a discrepancy event \(A\), the proportion of draws for which \(D^{\mathrm{rep}}\in A\) estimates its posterior predictive probability. This explains why parameter uncertainty belongs in the check: simulating only at a single point estimate would generally produce a narrower predictive distribution than the fitted Bayesian model implies.

Useful diagnostics include checking whether the result changes substantially with a reasonable alternative discrepancy, inspecting replicated trajectories rather than only summary statistics, and increasing the number of simulations when the estimated tail proportion is noisy. Also verify that replicated data use the same cell count, trajectory positions, branch-label mechanism, and expression scale as the observed data. If replicated data systematically miss the observed pattern, investigate model structure and data processing rather than merely increasing simulation count.

Posterior predictive checks complement, rather than replace, other approaches. Residual plots can expose local patterns that a single discrepancy hides, but do not directly summarize how surprising those patterns are under the posterior predictive distribution. Predictive model comparisons can assess relative performance across candidate models, but a model can predict better than its competitors and still fail an important absolute-fit check. A check focused on one feature also cannot certify the rest of the model.

Posterior predictive tail proportions are not generally calibrated as uniform p-values under a true model, because the observed data contribute to the posterior used to generate the replicates. They may consequently be conservative. Do not apply a conventional significance threshold mechanically, and do not interpret a moderate tail proportion as confirmation. The check's role is to make discrepancies visible and guide model revision.

Common Pitfalls and Scope

• Using the wrong replication scheme. Simulating a fixed number of cells when the model predicts a random number, or ignoring unequal sampling, changes the question being checked. Match the fitted observation process.
• Confusing a continuous mean with identical cell values. Continuity constrains the branch mean curves at the split; it does not require noisy observations on the two branches to match exactly.
• Checking only the fitted mean. A model can match average expression while missing branch occupancy or residual spread. Choose discrepancies tied to distinct features of the model.
• Treating a tail proportion as a formal null-hypothesis p-value. It is a posterior predictive diagnostic, not the probability that the model is correct.
• Overlooking the scale and preprocessing. A normal residual model for transformed expression does not imply normality on the original count scale. Check predictions on the scale specified by the observation model.

These checks are limited by the model and the observed data. Sparse cells near the split can make continuity diagnostics weak; uncertain branch labels can blur occupancy and variance summaries; and a discrepancy selected after inspecting many patterns can exaggerate apparent evidence of misfit. State which feature each check targets, examine more than one scientifically motivated feature, and avoid claiming that a successful check validates every modeling assumption.