Why Simulate Power?
For a simple comparison, a formula can sometimes approximate power from the donor count, effect size, and variance. In a multi-condition single-cell study, however, the planned comparison may depend on repeated measurements within donors, an interaction contrast, batch adjustment, or donor-specific responses. Simulation lets the investigator reproduce these design features and estimate how often the planned analysis would detect a specified effect.
The unit of biological replication is usually the donor, not the cell. A large number of cells from a few donors can improve the precision of each donor’s measurement, but it does not create the same information as recruiting additional independent donors. A simulation should therefore represent donor-level variation and the planned way of summarizing or modeling cell-level measurements. The examples below use a continuous donor-level expression summary, such as a suitably transformed pseudobulk value for a gene in a prespecified cell population. A count-based analysis requires a simulation and fitted test that also reflect its count distribution and mean–variance relationship.
Specify the Design and the Test Before Simulating
A power calculation is conditional on its assumptions. It does not discover an effect size or a variance structure from the design alone. Before generating data, state the biological unit, the conditions measured on each donor, the batches and their allocation, the effect of interest, the variability and correlation assumptions, and the exact test that will be applied.
A useful general model for a donor-level outcome is
Here, \(\mathbf{y}\) contains the simulated outcomes, \(X\) represents the fixed design terms (such as condition and batch), and \(\boldsymbol{\beta}\) contains their effects. \(Z\mathbf{u}\) represents random donor effects, and \(\boldsymbol{\varepsilon}\) represents residual variation. The random effects and residuals are assumed independent in this model. The matrices \(G\) and \(R\) specify the relevant variances and correlations. For example, a donor random intercept induces correlation among conditions measured on the same donor; a random slope allows donors to differ in their response to a condition.
The target is commonly a prespecified contrast \(c^\mathsf{T}\boldsymbol{\beta}\), such as a difference between two conditions or a factorial interaction. The null hypothesis is \(H_0:c^\mathsf{T}\boldsymbol{\beta}=0\). Each simulated data set must be analyzed using the model and test intended for the real study. If that test is a two-sided Wald test, for instance, reject when the simulated test statistic exceeds the appropriate critical value for the chosen \(\alpha\) and degrees of freedom. Do not simulate under one model and calculate power using a materially different analysis.
The model assumes that the specified design is estimable: the contrast must not be aliased with batch or another design term. It also assumes that the chosen distributions and variance components are reasonable approximations for the planned outcome. A simulation cannot repair complete condition–batch confounding, incorrect biological replication, or an unsuitable analysis model.
A Simulation Workflow
Make one record for each donor-condition measurement, with the actual condition, donor, and batch assignments. Include all repeated measurements specified by the design.
Specify the contrast effect and donor-level and residual variation. Use a covariance structure consistent with the repeated-measurement plan.
Draw random donor effects and residuals under the stated model, then add fixed condition and batch effects. Repeat to create one hypothetical study data set.
Fit the model, calculate the prespecified contrast and test statistic, and record whether the test rejects \(H_0\) at \(\alpha\).
For \(B\) simulated data sets, estimate power by the proportion of rejections. Report \(B\), the assumptions, and Monte Carlo uncertainty.
If \(I_b=1\) when simulation \(b\) rejects and \(I_b=0\) otherwise, then \(\widehat{\text{Power}}=\sum_{b=1}^{B} I_b/B\). Because this is a sample proportion from the simulation, its approximate Monte Carlo standard error is \(\sqrt{\widehat{\text{Power}}(1-\widehat{\text{Power}})/B}\). This quantifies simulation noise, not uncertainty about whether the assumed effect or variance is correct.
Worked Example: A Paired Condition Contrast
Suppose 16 hypothetical donors each contribute measurements under a reference condition and an active condition. The planned analysis is a two-sided paired \(t\)-test at \(\alpha=0.05\) on donor-level differences. Assume the marginal standard deviation under each condition is \(0.8\) log-expression units and the within-donor correlation is \(0.65\). The prespecified mean difference is \(0.5\) log-expression units.
For two equally variable measurements with correlation \(\rho\), the variance of their difference is \(2\sigma^2(1-\rho)\). Thus the assumed standard deviation of donor differences is
The standard error under these assumptions is \(0.6693/\sqrt{16}=0.1673\) log-expression units. The noncentrality parameter for the planned effect is \(0.5/0.1673=2.989\). A two-sided paired \(t\)-test has \(df=16-1=15\) and critical value about \(2.131\). These calculations provide a check on the expected power scale; the simulation still estimates power by generating paired outcomes, calculating each simulated difference, and applying that same \(t\)-test.
In a hypothetical run of \(B=10{,}000\) simulations, suppose 7,900 tests reject. The estimated power is \(7{,}900/10{,}000=0.790\). Its Monte Carlo standard error is \(\sqrt{0.790(0.210)/10{,}000}=0.00407\), or about 0.41 percentage points. This estimate describes power under the specified effect, correlation, and variance assumptions; it is not a guarantee of detection in an actual study.
Correlated Responses and Interaction Contrasts
The covariance structure matters because repeated measurements on one donor are not independent. A paired comparison can remove a donor-wide random intercept: if \(y_{ij}=\mu_j+u_i+\varepsilon_{ij}\), then \(y_{i1}-y_{i0}=(\mu_1-\mu_0)+(\varepsilon_{i1}-\varepsilon_{i0})\). The shared \(u_i\) cancels, but donor-specific differences in treatment response do not cancel if the model includes random slopes.
For a factorial design, the target may be an interaction rather than a condition’s average effect. Simulating the interaction requires generating all relevant measurements for each donor and applying the interaction contrast to each simulated data set. Simulating only a single group-level difference can give the wrong variance and therefore the wrong power.
Worked Example: Power for a Factorial Interaction
In a hypothetical \(2\times2\) design, 20 donors each contribute measurements in all four combinations of factors \(A\) and \(B\). Let the interaction contrast be \(y_{11}-y_{10}-y_{01}+y_{00}\), with a prespecified mean of \(0.6\) log-expression units. Suppose the model includes a donor-specific interaction deviation with standard deviation \(0.4\) and independent cell-condition residuals with standard deviation \(0.5\). Batch is recorded as a fixed effect, and each donor has the same batch for all four measurements.
The donor-specific interaction deviation contributes variance \(0.4^2=0.16\) to the contrast. Four independent residual terms, each with coefficient \(+1\) or \(-1\), contribute \(4(0.5^2)=1.00\). The variance of a donor’s interaction contrast is therefore \(1.16\) square log-expression units, and its standard error across 20 independent donors is
The noncentrality parameter is \(0.6/0.2408=2.491\). A paired-contrast \(t\)-test has \(df=19\); its two-sided 0.05 critical value is about \(2.093\). A simulation should generate all four outcomes per donor with the stated random interaction and residual terms, include the batch term, fit the prespecified model, and test the interaction. If a hypothetical 10,000-replicate run yields 6,400 rejections, estimated power is \(0.640\). The result is conditional on the assumed donor-specific interaction variability; omitting that component would understate the contrast variance and overstate power.
Batch Structure, Diagnostics, and Alternatives
Batch belongs in the simulation in the same way it belongs in the planned analysis. If condition and batch are crossed so that the contrast is estimable, simulate the actual allocation and include batch in the analysis. When every condition is measured for each donor and a donor stays in one batch, a batch fixed effect common to that donor’s measurements cancels in a within-donor contrast. In a design where conditions are distributed across batches, batch may instead affect the contrast precision directly.
Worked Example: Batch Allocation and Monte Carlo Precision
Consider 12 hypothetical donors assigned four each to three batches. Each donor is measured in both conditions, so the condition contrast is paired within donor. Assume the standard deviation of donor differences is \(0.6\) log-expression units and the target difference is \(0.5\). Batch has a fixed effect on the outcome, but because each donor contributes both conditions in the same batch, that fixed effect cancels in each paired difference.
The contrast standard error is \(0.6/\sqrt{12}=0.1732\) log-expression units, and the noncentrality parameter is \(0.5/0.1732=2.887\). The paired test uses \(df=11\), with a two-sided 0.05 critical value of about \(2.201\). A simulation that follows the allocation and includes batch should therefore produce power consistent with a paired test at this scale. If a hypothetical run of 10,000 replicates gives 7,200 rejections, estimated power is \(0.720\), with Monte Carlo standard error \(\sqrt{0.720(0.280)/10{,}000}=0.00449\), or 0.449 percentage points.
Now suppose instead that all reference samples are processed in one batch and all active samples in another, with no overlap. The condition effect and batch effect cannot be separated from these data. A simulation cannot make that contrast identifiable: the design matrix is rank deficient, so the analysis cannot estimate the condition effect independently of batch. The appropriate conclusion is that the design must be changed, not that a numerical power estimate resolves the confounding.
Simulation is useful when the planned design or test makes a closed-form power formula awkward. A formula based on a simplified paired \(t\)-test can be a quick benchmark, as in the examples, but may not capture random slopes, several conditions, or batch allocation. Simulation can represent those features, provided the data-generating model is credible and the same analysis is applied to every replicate. Resampling pilot data may preserve some observed complexity, but it cannot reliably represent effects or variance components that the pilot does not estimate well.
Common diagnostics include checking that simulated means, variances, correlations, and batch patterns match their specified values; confirming that the fitted model converges in nearly all replicates; and verifying that the contrast and rejection rule are exactly those planned. Under the null, the procedure should reject at approximately the intended \(\alpha\), subject to Monte Carlo error. If it rejects far too often, the test calibration, model assumptions, or implementation may be wrong. Also investigate whether the simulation accidentally treats cells as independent donors, whether all intended measurements were generated, and whether failures to fit are silently discarded. The estimated power applies to the specified model and analysis; it does not automatically account for uncertain assumptions, missing donors, or unmodeled sources of variation.
Check Your Understanding
Answer each question, then compare with the solutions below.
- Why does adding cells from the same donor not generally provide the same power gain as adding another independent donor?
- For a paired measurement model with a shared donor random intercept, what happens to that intercept in the within-donor difference?
- A simulation runs 5,000 replicates and rejects in 3,750. Calculate estimated power and its approximate Monte Carlo standard error.
- What should be concluded if condition is perfectly confounded with batch?
- Why must an interaction-power simulation generate all four measurements for each donor in a \(2\times2\) design?
Solutions. (1) Cells from one donor are clustered and share donor-level biology; they do not constitute independent biological replicates. (2) The shared intercept cancels in the difference, although residual variation and any donor-specific treatment response remain. (3) Estimated power is \(3{,}750/5{,}000=0.750\); the Monte Carlo standard error is \(\sqrt{0.750(0.250)/5{,}000}=0.00612\), or about 0.612 percentage points. (4) The condition effect is not separately estimable from batch, so redesign is needed before a meaningful power calculation for that effect. (5) The interaction is a contrast across all four condition combinations, and its variance depends on the joint measurements and their covariance.