Gibbs Sampling for Joint Cluster Assignments

Why Gibbs Updates Must Respect the Shared Allocation

In a shared-allocation mixture, each sample has one cluster label used across all measured omics blocks. The cluster profiles can differ by modality, but the allocation is shared: a sample assigned to cluster \(k\) contributes its observed RNA, protein, methylation, or other measurements to the corresponding cluster-specific profiles. Gibbs sampling updates these unknown quantities by repeatedly drawing each one from its conditional distribution given the current values of the others.

The central calculation for a label update is a product of likelihood contributions across the sample’s observed modalities. This product combines evidence without treating a missing measurement as evidence for or against any cluster. The procedure here is uncollapsed: cluster profiles and prevalence parameters remain explicit in the chain, and labels are sampled conditional on their current values.

Definition: For \(n\) samples and \(M\) omics modalities, let \(z_i\in\{1,\ldots,K\}\) be sample \(i\)'s shared cluster label, and let \(\pi=(\pi_1,\ldots,\pi_K)\) be cluster prevalences, where \(\sum_{k=1}^K\pi_k=1\). Let \(x_{im}\) denote sample \(i\)'s data block for modality \(m\), and let \(\theta_{mk}\) be the profile parameters for modality \(m\) in cluster \(k\). Write \(O_{im}=1\) if that block is observed and \(O_{im}=0\) if it is missing. Conditional on \(z_i=k\) and \(\theta_{mk}\), an observed block has density or mass function \(f_m(x_{im}\mid\theta_{mk})\).

Assume samples are conditionally independent given their labels and model parameters. Within each sample, modalities are conditionally independent given its label and the corresponding modality-specific profiles. Let \(\pi\sim\operatorname{Dirichlet}(\alpha_1,\ldots,\alpha_K)\), with each \(\alpha_k>0\). Assume independent priors \(p(\theta_{mk})\) for the profiles, chosen to match their likelihoods when conjugate updates are desired. For the missingness treatment below, assume the observation pattern is ignorable for inference on the modeled data: missingness does not contribute an additional likelihood factor depending on \(z_i\), \(\pi\), or \(\theta\). These assumptions define the updates; a dependent missingness process would require its own model.

The Joint Distribution and Label Conditional

Let \(n_k=\sum_{i=1}^n I(z_i=k)\), where \(I(\cdot)\) is 1 when its condition holds and 0 otherwise. Ignoring constants that do not depend on the unknowns, the joint posterior density is proportional to the prior densities times the likelihood:

p(z,\pi,\theta\mid x_{\mathrm{obs}}) \propto \left\{\prod_{k=1}^K\pi_k^{\alpha_k-1}\right\} \left\{\prod_{i=1}^n \pi_{z_i} \prod_{m:O_{im}=1} f_m(x_{im}\mid\theta_{m,z_i})\right\} \left\{\prod_{m=1}^M\prod_{k=1}^K p(\theta_{mk})\right\}.

Here \(x_{\mathrm{obs}}\) denotes the observed blocks. A missing block is omitted from the product rather than assigned a value or a likelihood factor. To update \(z_i\), hold \(\pi\), all profiles, and the other labels \(z_{-i}\) fixed. Terms in the joint density that do not depend on the candidate label cancel, leaving:

P(z_i=k\mid z_{-i},\pi,\theta,x_{\mathrm{obs}}) = \frac{w_{ik}}{\sum_{\ell=1}^K w_{i\ell}}, \qquad w_{ik}=\pi_k\prod_{m:O_{im}=1} f_m(x_{im}\mid\theta_{mk}).

This is a categorical draw over the \(K\) candidate clusters. If every modality for sample \(i\) is missing, the empty product equals 1, so the update reduces to \(P(z_i=k\mid\cdots)=\pi_k\). The label is then informed by the current prevalence but not by a nonexistent data likelihood. For numerical stability, implementations commonly calculate \(\log w_{ik}=\log\pi_k+\sum_{m:O_{im}=1}\log f_m(x_{im}\mid\theta_{mk})\), subtract the largest log weight, exponentiate, and normalize.

Worked Example: Combining Two Observed Modalities

Suppose a hypothetical model has two clusters with current prevalences \(\pi_1=0.6\) and \(\pi_2=0.4\). For one sample, an observed Gaussian RNA block has likelihood 0.20 under cluster 1 and 0.50 under cluster 2. Its observed Bernoulli protein measurement has likelihood 0.80 under cluster 1 and 0.25 under cluster 2. Both blocks are observed, and the conditional-independence assumption applies. The unnormalized label weights are:

w_{i1}=0.6(0.20)(0.80)=0.096, \qquad w_{i2}=0.4(0.50)(0.25)=0.050.

The normalizing sum is \(0.096+0.050=0.146\), so \(P(z_i=1\mid\cdots)=0.096/0.146=0.6575\) and \(P(z_i=2\mid\cdots)=0.3425\). A Gibbs iteration samples from these probabilities; it does not automatically assign the sample to the larger-probability cluster. The calculation uses likelihood values conditional on the current profiles, not posterior probabilities from each modality treated separately.

The derivation follows directly from the joint density: when \(z_i=k\), its allocation prior contributes \(\pi_k\), and every observed block contributes its cluster-\(k\) likelihood. The remaining factors, including other samples’ terms and the profile priors, are unchanged across candidate values of \(z_i\), so they cancel in normalization. This cancellation explains why the update uses the current \(\pi\) rather than a prevalence estimate based on counts with sample \(i\) removed. Counts influence the label indirectly through the current draw of \(\pi\).

Updating Prevalence and Modality Profiles

With labels fixed, the prevalence update is conjugate. The likelihood contribution involving \(\pi\) is \(\prod_k\pi_k^{n_k}\); multiplying by its Dirichlet prior gives a Dirichlet conditional with each prior parameter increased by its cluster count:

\pi\mid z,\theta,x_{\mathrm{obs}} \sim \operatorname{Dirichlet}(\alpha_1+n_1,\ldots,\alpha_K+n_K).

For a profile parameter \(\theta_{mk}\), only samples assigned to \(k\) with an observed block in modality \(m\) contribute to its likelihood. Thus the general conditional is:

p(\theta_{mk}\mid z,\pi,\theta_{-(mk)},x_{\mathrm{obs}}) \propto p(\theta_{mk}) \prod_{i:z_i=k,\;O_{im}=1} f_m(x_{im}\mid\theta_{mk}).

The restriction \(O_{im}=1\) is essential: an assigned sample with a missing block in modality \(m\) contributes no likelihood factor for \(\theta_{mk}\). If no sample in cluster \(k\) has that modality observed, the product is empty and the conditional is just the prior \(p(\theta_{mk})\). For conjugate prior-likelihood pairs, this expression yields a standard-form distribution; otherwise it specifies a density from which an appropriate conditional sampler is needed.

Worked Example: Updating a Gaussian Cluster Mean

Consider a hypothetical continuous modality with known within-cluster variance \(\sigma^2=4\). For cluster \(k\), the mean has prior \(\mu_{mk}\sim N(0,9)\). Two samples assigned to this cluster have observed values 2 and 4 in this modality; a third assigned sample is missing this block. The observed-data likelihood uses only 2 and 4. For a normal likelihood with known variance and a normal prior, the conditional variance and mean are:

V_{\mathrm{post}}= \left(\frac{1}{9}+\frac{2}{4}\right)^{-1} =\left(\frac{11}{18}\right)^{-1} =\frac{18}{11}=1.636, \qquad \mu_{\mathrm{post}}= V_{\mathrm{post}}\left(\frac{0}{9}+\frac{2+4}{4}\right) =\frac{18}{11}(1.5)=\frac{27}{11}=2.455.

Therefore \(\mu_{mk}\mid\cdots\sim N(2.455,1.636)\), where the second argument is the variance. The missing sample adds neither a value nor an additional precision term; including it would be unjustified because its likelihood is undefined. The posterior mean is in the same measurement units as the observations, and the posterior variance is in squared units.

Worked Example: Updating Cluster Prevalences

Suppose a hypothetical three-cluster model uses \(\pi\sim\operatorname{Dirichlet}(1,1,1)\). At the current label configuration, the counts are \(n_1=3\), \(n_2=1\), and \(n_3=0\). The conditional distribution is:

\pi\mid z,\cdots\sim\operatorname{Dirichlet}(4,2,1).

The parameters sum to 7, so the conditional mean prevalence vector is \((4/7,2/7,1/7)=(0.5714,0.2857,0.1429)\). These means are summaries of the conditional distribution, not the update itself: the Gibbs step draws a random vector from \(\operatorname{Dirichlet}(4,2,1)\). The empty third cluster has no assigned observations, but its prevalence remains represented through its positive prior parameter, and its profile parameters are updated from their priors if it has no observed assigned blocks.

One Complete Uncollapsed Gibbs Sweep

A sweep updates all unknowns in a sequence of conditional draws. One valid ordering is to update prevalence, then profiles, then labels; each conditional must use the latest available values. The following outline makes the data contributions explicit.

Update prevalence.
Draw \(\pi\) from \(\operatorname{Dirichlet}(\alpha_1+n_1,\ldots,\alpha_K+n_K)\), using the current labels.

Update profiles.
For each modality and cluster, draw \(\theta_{mk}\) from its prior times the likelihood for samples assigned to \(k\) whose modality-\(m\) blocks are observed.

Update shared labels.
For each sample, calculate and normalize \(w_{ik}\) across clusters using the current prevalence and profiles, then draw one label. Use all and only that sample’s observed modalities.

This is an uncollapsed scheme because \(\pi\) and \(\theta\) are sampled rather than integrated out before drawing labels. A different ordering can also define a valid Gibbs sweep when each draw is from the full conditional for the variables being updated. What matters is that each conditional corresponds to the stated joint model and that the chain’s state is updated consistently.

Diagnostics, Alternatives, and Limitations

Basic checks can catch implementation errors before interpreting a chain. At every label update, all weights should be nonnegative and their normalized probabilities should sum to 1 within numerical tolerance. Verify that a sample missing modality \(m\) has no \(f_m\) factor in its label weights and that it is excluded from the profile update for \(\theta_{mk}\). Check that the prevalence parameters equal \(\alpha_k+n_k\), and that the cluster counts sum to \(n\). In log-scale calculations, handle impossible assignments with log weight \(-\infty\); if every candidate has zero likelihood, the model or computation has failed for that observation and should not be “fixed” by arbitrary normalization.

This method relies on its likelihood and prior assumptions. A single shared label may be inadequate when different modalities support genuinely different sample partitions; conditional independence may be inadequate when residual cross-modality dependence remains within clusters; and an ignorable missingness assumption fails if the probability of observing a block depends on unmodeled values or labels. Prior choices can also affect sparse and empty clusters. Gibbs sampling does not remove these modeling limitations, and a draw of a label is not a certainty about biological subtype membership.

The main alternative developed in collapsed sampling integrates conjugate profile parameters out of the label conditional, rather than conditioning on their current draws. This can change mixing behavior and the form of the assignment weights; it is not the update derived here. Uncollapsed updates have the practical advantage that profile parameters are directly sampled and readily available at each iteration, but strong dependence between labels and profiles can make movement through the posterior slow. Neither approach makes the conditional-independence or shared-allocation assumptions true by itself.

Common Pitfalls

Multiplying by a missing block’s likelihood. A missing measurement does not have a likelihood value to insert. Omit that factor from both the label weight and the profile conditional.
Using a single modality’s evidence as the joint assignment. The shared label conditional combines the prevalence and every observed modality likelihood, not a separate hard assignment from each block.
Choosing the most probable label instead of sampling. A Gibbs update draws from the normalized conditional probabilities. Replacing the draw by the mode changes the algorithm.
Confusing an empty cluster with an invalid state. If its prior supports positive prevalence and its profile prior is proper, a cluster may have zero current members; its profile conditional then equals its prior for modalities with no observed assigned data.

Key takeaway: In an uncollapsed shared-allocation Gibbs sampler, the label weight is prevalence times the product of likelihoods for that sample’s observed modalities. Profile updates likewise use only assigned samples with observed values in that modality; missing blocks contribute no likelihood factor.

Check Your Understanding

Answer each question, then compare with the solutions below.

A sample has an observed RNA block but a missing protein block. Which likelihood factors enter its label weight?
In a two-cluster model, \(\pi=(0.7,0.3)\). The only observed modality has likelihoods \(0.1\) and \(0.4\) under clusters 1 and 2. Calculate the normalized label probabilities.
For a modality-specific profile, what happens to its conditional posterior if no assigned sample has that modality observed?
With prior \(\operatorname{Dirichlet}(2,1)\) and current cluster counts \((2,4)\), give the prevalence conditional distribution.

Solutions. (1) Include the RNA likelihood and omit the protein factor; the weights are proportional to \(\pi_k f_{\mathrm{RNA}}(x_{i,\mathrm{RNA}}\mid\theta_{\mathrm{RNA},k})\). (2) The weights are \(0.7(0.1)=0.07\) and \(0.3(0.4)=0.12\), with sum \(0.19\). Thus the probabilities are \(0.07/0.19=0.3684\) and \(0.12/0.19=0.6316\). (3) Its conditional equals its prior because the likelihood product is empty. (4) The conditional is \(\operatorname{Dirichlet}(2+2,1+4)=\operatorname{Dirichlet}(4,5)\).