Posterior Predictive Checks for Compositional Models

Why Check Replicated Compositions?

A fitted Bayesian composition model can produce estimates for cell-type proportions while still failing to reproduce important features of the data. For example, it may predict average cell-type proportions reasonably well but generate too few donors with no cells of a rare type, or much less donor-to-donor variability than is observed. Posterior predictive checks make these failures visible by asking a direct question: if the fitted model generated new data resembling this study, would those data look like the observed compositions in the features we care about?

Let $y_{ik}$ be the observed number of cells of type $k$ in donor $i$, for donors $i=1,\ldots,m$ and cell types $k=1,\ldots,K$. The observed total for donor $i$ is $n_i=\sum_{k=1}^K y_{ik}$. Let $\theta$ denote the model parameters, including any donor-level or group-level parameters. The relevant posterior predictive distribution is obtained by drawing parameters from their posterior distribution and then generating replicated counts from the model at those parameter values.

Definition: The posterior predictive distribution of replicated data $y^{\mathrm{rep}}$, given observed data $y$, is $$ p(y^{\mathrm{rep}}\mid y)=\int p(y^{\mathrm{rep}}\mid\theta)\,p(\theta\mid y)\,d\theta. $$ A posterior predictive check compares a feature, or discrepancy, $T$ of the observed and replicated data. A common tail-area summary is $$ p_B=P\{T(y^{\mathrm{rep}},\theta)\geq T(y,\theta)\mid y\}. $$

The notation allows $T$ to depend on $\theta$, which is useful when the discrepancy is defined relative to model-predicted counts. To approximate $p_B$, draw $\theta^{(s)}$ from the posterior, draw $y^{\mathrm{rep},(s)}$ from $p(y^{\mathrm{rep}}\mid\theta^{(s)})$, and compare the two discrepancy values using the same $\theta^{(s)}$. The estimated tail area is the fraction of replicated values at least as large as their corresponding observed values. For a Monte Carlo estimate based on $S$ draws, a simple finite-simulation correction is $(r+1)/(S+1)$, where $r$ is the number of exceedances.

A posterior predictive tail area is not a classical p-value with a uniform null distribution, nor is it the probability that the model is true. The data have already informed the posterior used to generate replicates, so these checks can be conservative. Their primary value is diagnostic: a small or large tail area can flag a feature that the model tends to underproduce or overproduce. The discrepancy itself and the predictive distribution should also be inspected; a single number is not a verdict on the whole model.

Choosing Discrepancies for Cell-Type Compositions

A check is informative only if it targets a feature that could reveal a relevant failure. For single-cell compositions, useful discrepancies include the number of donors with zero observed cells of a cell type, the largest donor-specific proportion, the distribution of counts across donors, or a Pearson-type measure of deviation from model-predicted counts. A check of the overall cell-type totals may detect an incorrect mean composition, while donor-specific checks can reveal heterogeneity that pooled totals conceal.

The data-generating model determines what can be checked. In a multinomial model conditional on $n_i$, the donor total is fixed by construction. Replicated compositions can assess how those totals are divided among cell types, but cannot test whether the observed totals themselves are plausible. Checking totals requires a model that also generates totals. Similarly, a model that imposes one shared composition across donors may be challenged by donor heterogeneity, whereas a model with donor-specific variation should be checked for whether it captures the observed spread without being excessively variable.

For a discrepancy such as a Pearson statistic, the expected counts and parameter estimates must be defined clearly. If $\pi_{ik}(\theta)$ is the model-predicted proportion for donor $i$ and type $k$, then the expected count is $n_i\pi_{ik}(\theta)$. A typical discrepancy is

T(y,\theta)=\sum_{i=1}^m\sum_{k=1}^K \frac{\{y_{ik}-n_i\pi_{ik}(\theta)\}^2}{n_i\pi_{ik}(\theta)}.

This statistic summarizes deviations, but small expected counts can make its components unstable and its reference distribution need not follow a simple chi-squared approximation. In a posterior predictive check, the safer comparison is usually to calculate the same discrepancy for each replicated data set using its paired posterior draw, rather than relying on a textbook reference distribution.

Worked Example: Checking How Often a Rare Type Is Absent

In a hypothetical study, 20 donors each contribute a composition of fixed size. The observed number of donors with zero cells of a particular rare type is 11. Suppose a simplified fitted posterior predictive model assigns probability $q=9/29$ to a zero count for each donor, independently, at every posterior draw. This point-concentrated predictive setup makes the count of zero donors exactly binomial; it is a simplifying assumption for the calculation, not a general property of posterior predictive distributions.

Let $Z$ be the number of donors with zero rare-type cells in a replicated data set. Under these stated assumptions, $Z\sim\operatorname{Binomial}(20,9/29)$. Its expectation is

E(Z)=20\left(\frac{9}{29}\right)=6.207.

The observed count, 11, is therefore substantially above the model's expected count. The exact upper-tail probability is

P(Z\geq 11)=\sum_{z=11}^{20} \binom{20}{z}\left(\frac{9}{29}\right)^z \left(\frac{20}{29}\right)^{20-z} \approx 0.02221.

Thus, under this specified predictive model, about 2.221% of replicated studies would have at least 11 zero-count donors. This is evidence that the model tends to generate fewer such donors than observed, for this diagnostic feature. It does not establish why the discrepancy occurs, and the tail area is not a probability that the model is false.

If the probability of a zero count varied across posterior draws, or if donors shared uncertain parameters, $Z$ would generally not have this simple binomial distribution after integrating over the posterior. The appropriate calculation would then draw replicated data across posterior draws and compare the replicated zero count with the observed count.

Checking Donor Heterogeneity

Pooling counts across donors can hide an important kind of lack of fit. Two data sets can have the same total count for every cell type while one has similar compositions in all donors and the other has substantial between-donor variation. A posterior predictive check should therefore retain the donor structure whenever the fitted model makes donor-level predictions.

In a Dirichlet-multinomial model, for example, donor compositions vary around a population composition, with a concentration parameter controlling the amount of variation. A multinomial model with one shared composition has less flexibility: conditional on its probabilities, donors vary through multinomial sampling alone. These models may reproduce similar pooled totals but generate different distributions of donor-level discrepancies. The goal of the check is not to prefer complexity automatically, but to see whether the model's implied variation resembles the observed variation.

Worked Example: A Donor-Level Pearson Check

Consider four hypothetical donors, each with 20 cells classified into three types. The observed count rows are $(15,4,1)$, $(14,5,1)$, $(5,8,7)$, and $(14,5,1)$. Suppose a fitted shared-composition model has posterior predictive proportions $(0.60,0.275,0.125)$, treated as fixed here to keep the arithmetic transparent. The expected counts per donor are therefore $(12,5.5,2.5)$.

For donor 1, the Pearson contribution is

\frac{(15-12)^2}{12}+ \frac{(4-5.5)^2}{5.5}+ \frac{(1-2.5)^2}{2.5} =0.7500+0.4091+0.9000=2.059.

For donor 2, the contribution is

\frac{(14-12)^2}{12}+ \frac{(5-5.5)^2}{5.5}+ \frac{(1-2.5)^2}{2.5} =0.3333+0.04545+0.9000=1.279.

Donor 4 has the same counts and contribution as donor 2. Donor 3 contributes

\frac{(5-12)^2}{12}+ \frac{(8-5.5)^2}{5.5}+ \frac{(7-2.5)^2}{2.5} =4.083+1.136+8.100=13.32.

The total observed discrepancy is $T(y)=2.059+1.279+13.32+1.279=17.94$, to rounding. It is driven especially by donor 3, whose rare-type count is much higher than the shared-composition model expects.

Now suppose 1,000 hypothetical posterior predictive replicates, each with four rows of 20 counts generated under this fitted model, produce a Pearson discrepancy at least as large as 17.94 in 43 replicates. The finite-simulation estimate is

\widehat p_B=\frac{43+1}{1000+1}=0.04396.

This check flags the observed degree of donor-level deviation as relatively unusual under the specified model. It does not prove that donor 3 is erroneous, nor does it identify a unique alternative model. A model with donor-specific composition variation might reproduce such patterns more often, but that possibility should be assessed by fitting and checking that model rather than inferred from this tail area alone.

A Short Justification for Posterior Predictive Checking

The construction follows directly from uncertainty about the parameters. If $\theta$ were known, a replicate would be generated from $p(y^{\mathrm{rep}}\mid\theta)$. Since $\theta$ is not known, averaging those conditional distributions over plausible values gives the posterior predictive distribution. For any event $A$ about a replicated data set, the law of total probability gives

P(y^{\mathrm{rep}}\in A\mid y) =\int P(y^{\mathrm{rep}}\in A\mid\theta)\,p(\theta\mid y)\,d\theta.

In simulation, the integral is approximated by averaging over posterior draws. Specifically, if $I_s$ is 1 when $T(y^{\mathrm{rep},(s)},\theta^{(s)})\geq T(y,\theta^{(s)})$ and 0 otherwise, then

\widehat p_B=\frac{1}{S}\sum_{s=1}^{S}I_s.

This estimate is a Monte Carlo approximation to the chosen posterior predictive tail area. Increasing $S$ reduces simulation error, but does not correct a poorly chosen discrepancy or a flawed model. For example, a discrepancy based only on pooled totals may miss substantial donor-level lack of fit even if its estimated tail area is stable.

Worked Example: Checking an Unusually Large Rare-Type Count

Suppose one hypothetical donor has 20 cells, of which 6 are a rare type. A simplified fitted model predicts a rare-type probability of $0.10$, with that probability fixed for this calculation. Conditional on the donor total, the rare-type count $X$ has distribution $\operatorname{Binomial}(20,0.10)$. For a one-sided check of an unusually large count,

P(X\geq 6)=1-P(X\leq 5) =1-\sum_{x=0}^{5}\binom{20}{x}(0.10)^x(0.90)^{20-x} \approx 0.01125.

The observed count of 6 or more is rare under this particular predictive distribution: its upper-tail probability is about 1.125%. This calculation checks one donor-level feature. A full posterior predictive analysis would repeat the comparison across parameter draws, and would typically examine the distribution of such counts over all donors rather than treating one selected donor as the sole diagnostic.

Diagnostics, Alternatives, and Limitations

A useful check combines a numerical summary with a visual comparison. Plot observed statistics alongside their replicated distributions, such as zero-donor counts, donor-specific proportions, or the maximum Pearson contribution. Inspect whether the observed value lies in the center, in a tail, or outside the range of most replicates. Also examine simulation diagnostics: posterior draws should represent the fitted posterior adequately, replicated counts must use the correct donor totals and model structure, and estimates based on a small number of exceedances should be reported with their Monte Carlo uncertainty.

For a tail-area estimate based on $S$ independent replicate comparisons, its approximate Monte Carlo standard error is $\sqrt{\widehat p_B(1-\widehat p_B)/S}$. For the estimate $44/1001$ in the donor example, using $S=1000$ gives an approximate standard error of $\sqrt{0.04396(1-0.04396)/1000}\approx 0.00648$. This quantifies simulation noise, not uncertainty about whether the model is scientifically adequate.

Posterior predictive checks complement, rather than replace, other diagnostics. Residual and overdispersion diagnostics can summarize systematic deviations in fitted abundance models; posterior predictive checks instead compare data features with distributions generated from the full fitted Bayesian model. Posterior summaries describe uncertainty about parameters, but a narrow posterior for mean proportions does not guarantee that donor-level variation or zero frequencies are well reproduced. Predictive checks can also be compared across candidate models, provided the same scientifically relevant discrepancies are used; a model should not be selected solely because it fits one chosen statistic.

Several limitations matter in compositional settings. First, the data are constrained: counts sum to $n_i$, so increasing one proportion necessarily changes the others. Component-wise checks should be interpreted in that joint context. Second, because the posterior is conditioned on the observed data, posterior predictive tail areas are often concentrated away from extreme values and should not be read as calibrated hypothesis-test p-values. Third, many possible discrepancies can be examined. Searching many checks and reporting only the most alarming one exaggerates the apparent evidence of misfit; it is better to specify a small set of checks tied to the model's assumptions and scientific aims.

Finally, a check can only detect failures expressed in the summaries it examines. A model may reproduce the pooled composition and still miss rare zeros, donor clusters, or a tail of unusually large cell-type proportions. Conversely, a discrepancy may be unusual because of one influential donor or a small expected count. Follow an alarm by inspecting the relevant donor-level data and by checking whether the model assumptions and sampling design support the proposed replication scheme. A model conditioned on cell totals cannot diagnose total-count variation, and a model with no donor-level structure cannot explain patterns of persistent donor heterogeneity merely by fitting pooled proportions.

Common Mistakes

Treating a posterior predictive tail area as a classical p-value. It is a model-based diagnostic probability, not a uniformly calibrated test or the probability that the null or model is true.
Checking only pooled counts. Similar pooled totals can conceal substantial donor-to-donor variation. Include donor-level discrepancies when donor heterogeneity matters.
Testing a feature the model fixes by conditioning. A multinomial model conditional on donor totals cannot assess whether those totals are plausible.
Assuming all donor zero indicators are independent after integrating over parameters. Shared posterior uncertainty can induce dependence; use replicated data generated from the full posterior predictive procedure.
Calling one discrepancy a complete assessment of fit. A check diagnoses only the feature it measures, so use a limited, prespecified set of complementary summaries.

Key takeaway: A posterior predictive check compares observed compositional features with the same features in data generated from posterior draws and the fitted model. Choose checks that target rare-type zeros, donor heterogeneity, or other relevant features, and interpret each discrepancy as a focused model diagnostic—not as a definitive test of the entire model.

Check Your Understanding

Answer each question, then compare your reasoning with the solutions.

A multinomial composition model conditions on each donor's observed cell total. Can a posterior predictive check from this model establish whether the donor totals are plausible? Why?
For a specified replicated-data discrepancy, 28 of 799 posterior predictive draws are at least as extreme as the observed discrepancy. Calculate the finite-simulation tail-area estimate.
Why might pooled cell-type totals fail to reveal poor fit in donor-level compositions?
Explain why a posterior predictive tail area of 0.02 does not mean there is a 2% probability that the fitted model is true.
In a model with uncertain shared cell-type probabilities, why may the number of donors with zero cells of a type fail to follow a binomial distribution after posterior integration?

Solutions. (1) No. Replicates hold those totals fixed, so the check concerns allocation among cell types conditional on the totals, not the plausibility of the totals themselves. (2) The estimate is $(28+1)/(799+1)=29/800=0.03625$. (3) Pooling discards how counts are distributed across donors; markedly heterogeneous donor profiles can sum to the same pooled composition as relatively similar profiles. (4) The tail area is the posterior predictive probability of a discrepancy at least as extreme as observed under the specified procedure. It is not a posterior probability assigned to the model being true. (5) Integrating over shared uncertain probabilities can make donor outcomes dependent: donors with a shared probability draw are more likely to have zero counts together when that draw implies a low expected abundance. The resulting mixture need not be binomial.