Designing Simulations For Subtype Reproducibility

Why simulate subtype reproducibility?

A stability score is useful only if its behavior is understood in settings where the underlying structure is known. In real data, an unstable partition might reflect weak subtype separation, too few members of a rare subtype, noisy features, or batch effects. A stable partition can also be wrong: an algorithm may repeatedly recover a batch grouping rather than the biological subtypes of interest. Simulation allows these conditions to be varied deliberately and the resulting diagnostics to be compared with known subtype labels.

The goal is not to create a synthetic data set that perfectly resembles every aspect of a study. It is to isolate important mechanisms, test whether a diagnostic responds as expected, and discover conditions in which it can mislead. The earlier material on consensus clustering, subgroup-level stability, and end-to-end stability provides the diagnostics that can be evaluated here; this tutorial focuses on designing the simulation that evaluates them.

A multi-omics data-generating model

Let \(n\) be the number of subjects, \(M\) the number of omics modalities, and \(p_m\) the number of features in modality \(m\). Each subject \(i\) has a known subtype label \(z_i\in\{1,\ldots,K\}\), drawn with subtype probabilities \(\pi_1,\ldots,\pi_K\), where \(\sum_{k=1}^K\pi_k=1\). Each subject also has a batch label \(b_i\). A useful starting model for feature \(j\) in modality \(m\) is

Here \(\mu_{mj}\) is a baseline, \(a_{mjk}\) is the subtype effect, \(c_{mj,b_i}\) is the batch effect, and \(\varepsilon_{imj}\) is residual noise. In practice, the notation \(a_{mjk}\big|_{k=z_i}\) means to use the effect for subject \(i\)'s subtype. A convenient constraint is \(\sum_k\pi_k a_{mjk}=0\), so the subtype effects are centered at the population-average effect. Residual features may be correlated: within modality \(m\), one can specify \(\boldsymbol{\varepsilon}_{im}\sim N(\mathbf{0},\Sigma_m)\). If modalities are intended to be conditionally independent given subtype and batch, generate their residuals independently; if cross-modal residual correlation is plausible, specify it rather than allowing it to arise accidentally.

This model separates the factors the simulation should vary. Sample size \(n\) controls how much information is available. The probabilities \(\pi_k\) control subtype prevalence. Subtype effects determine separation, while the scale and covariance of residual noise determine how obscured the effects are. Batch effects and the allocation of subjects to batches determine whether technical structure competes with, or is confounded with, subtype structure. The analysis pipeline must also be specified: preprocessing, feature selection, integration, clustering, and any resampling steps should be repeated as they would be in the intended analysis.

Build scenarios that isolate the factors

Begin with a reference scenario, then change one factor at a time to make its effect interpretable. A second, factorial set of scenarios can assess interactions—for example, whether increasing sample size helps less when a subtype is rare and feature noise is high. Choose factor levels that span scientifically plausible conditions and include deliberately difficult cases. State the levels before examining results; choosing them after seeing which conditions perform well can make the evaluation overly reassuring.

Separation can be specified through a standardized contrast. For two subtype means \(\boldsymbol{\mu}_1\) and \(\boldsymbol{\mu}_2\) with common within-subtype covariance \(\Sigma\), their Mahalanobis separation is

This measures mean difference relative to residual variability and feature correlation. Larger \(\Delta\) generally makes the subtypes easier to distinguish, though performance also depends on dimension, prevalence, and the method. For multiple modalities, specify which modalities carry signal: a subtype visible in one modality but not another presents a different integration problem from a subtype with coordinated effects across modalities.

For batch structure, vary both the size of batch effects and the relationship between batch and subtype. Balanced allocation makes each batch contain similar subtype proportions. Confounded allocation makes subtype composition differ across batches. These are not interchangeable: a large batch effect may be manageable when allocation is balanced, while even a moderate effect may mimic subtype structure under strong confounding.

Run the full analysis and define what counts as success

In every replicate, retain the known labels and batch labels, generate all modalities, and run the complete pipeline. If stability is estimated by resampling, repeat the same resampling procedure within each replicate. Compare partitions with a label-invariant measure such as the adjusted Rand index (ARI), which compares pairs of subjects while adjusting for chance agreement. Use the same declared measure consistently when comparing a fitted partition with truth and when comparing two fitted partitions.

Define success before computing results. For example, a replicate might count as successful recovery if the ARI between the fitted partition and the known labels is at least \(0.70\). A stability diagnostic might flag a data set as stable if its prespecified stability score is at least \(0.75\). These thresholds are illustrative, not universal. Use a separate set of simulation replicates, or a prespecified rationale, to avoid selecting thresholds just because they make the diagnostic look effective.

Across many generated data sets, compare the diagnostic flag with the defined recovery outcome. Sensitivity is the proportion of replicates meeting the recovery criterion that the diagnostic flags as stable; specificity is the proportion not meeting the recovery criterion that it does not flag as stable. Also report the underlying continuous stability and recovery summaries. A threshold can conceal useful distinctions, and neither sensitivity nor specificity describes performance outside the simulated scenarios.

Monte Carlo precision, diagnostics, and alternatives

A simulation estimate is itself noisy because it is based on a finite number of replicates. If \(R\) independent replicates are generated under one condition and \(r\) meet a binary criterion, the estimated proportion is \(\hat{q}=r/R\). Its approximate Monte Carlo standard error is \(\sqrt{\hat{q}(1-\hat{q})/R}\). Increasing \(R\) reduces this simulation uncertainty; it does not correct a misspecified data-generating model.

Useful diagnostics include checking that generated subtype counts match their intended distribution on average, verifying feature means and variances against the specified model, and confirming the batch-by-subtype table for each allocation condition. Inspect whether preprocessing leaks information across resamples, whether a feature-selection step uses truth labels, and whether failures concentrate in replicates with very small subtype counts. If methods are stochastic, record random seeds and separate variation due to their random initialization from variation due to generated data.

A one-factor-at-a-time design is easy to interpret but may miss interactions. A factorial design can estimate interactions but may require many more scenarios. An empirical simulation that resamples observed feature vectors can preserve complicated dependence, but it may not provide control over known subtype effects or represent conditions absent from the observed data. A parametric model offers direct control over prevalence, separation, noise, and batch effects, but can be unrealistic if its distributional assumptions are too simple. These approaches answer different questions; a useful study may combine a controlled parametric design with sensitivity analyses using alternative covariance or noise structures.

Common pitfalls and limits

• Equating stability with truth. Repeatedly obtaining the same partition does not establish recovery. Compare assignments with the known labels and check whether they align with batch.
• Changing several factors at once without a plan. A difference between scenarios cannot be attributed to prevalence if sample size, separation, and noise also changed. Use controlled contrasts or an explicit factorial design.
• Using only easy scenarios. A diagnostic that performs well under clear, balanced, abundant subtypes may fail for rare or noisy subtypes. Include difficult conditions and report them separately.
• Ignoring the analysis pipeline. Simulating only the final clustering step can overstate performance if real analyses also estimate preprocessing, select features, or integrate modalities. Repeat those steps inside each replicate where appropriate.
• Overgeneralizing results. Simulation conclusions apply to the specified data-generating mechanisms and pipeline. They do not guarantee reproducibility in a new population or prove biological validity.

A simulation is most informative when its assumptions are explicit and its failures are treated as results rather than discarded. Report the factors and levels, the number of replicates, the full pipeline, the success criteria, and the uncertainty of estimated performance. The central question is not whether a diagnostic ever returns a high score, but whether it distinguishes the conditions where recovery is reliable from those where it is not.