Why Simulate Pseudoreplication?
A statistical test can report strong evidence for a difference even when the biological groups have the same donor-population mean. One way this happens is pseudoreplication: cells from the same donor share biological features, but an analysis treats every cell as an independent replicate. Simulation makes this failure visible. We generate data under a null hypothesis that is true, preserve donor-level correlation, and repeatedly compare a cell-independent analysis with one that recognizes donors as the independent units.
The aim is not to simulate every feature of a single-cell experiment. It is to isolate one mechanism: shared donor variation makes cells from one donor more alike than cells from different donors. The introductory material on pseudoreplication describes how this dependence distorts uncertainty; here we make the mechanism operational and measure its consequences over repeated null data sets.
A Null Model With Donor-Level Correlation
Suppose there are two groups, indexed by \(g\in\{0,1\}\). Group \(g\) contains \(D\) independent donors, and donor \(d\) contributes \(m\) cells. Let \(Y_{gdi}\) be the expression value for cell \(i\) from donor \(d\) in group \(g\). A simple random-intercept model is
Assume donor effects are independent across donors and groups; residuals are independent across cells and independent of all donor effects. Here \(\mu_g\) is the group mean, \(b_{gd}\) is a shared donor deviation, and \(\varepsilon_{gdi}\) is cell-specific variation. Under the null hypothesis of no group difference, \(H_0:\mu_1-\mu_0=0\). The simulation must generate data under this null rather than merely resampling observed values whose group difference may not be zero.
The model assumes a common number of cells per donor, common variance components across groups, normally distributed donor effects and residuals, and independent donors. These are simplifying assumptions, not facts guaranteed by a data set. The normality assumption makes the donor means normally distributed and allows a standard donor-level \(t\)-test in the examples. The central issue—correlated cells—does not depend on normality, though the exact calibration of tests can.
Why the Two Analyses Behave Differently
For one donor, the mean of its \(m\) cells is \(\bar{Y}_{gd\cdot}\). Under the model, it equals \(\mu_g+b_{gd}\) plus the average residual. Its variance is \(\tau^2+\sigma^2/m\): collecting more cells reduces the residual contribution but does not eliminate the donor effect. If the target is a difference in group means across donors, a donor-level analysis uses the independent donor means as observations.
A cell-independent analysis instead counts all \(Dm\) cells in each group as independent. It estimates the cell-level marginal variance, \(\tau^2+\sigma^2\), but its usual standard error does not account for the shared donor effects. With equal group sizes, the actual variance of the difference in cell means is
The bracketed factor is the design effect. It is the ratio of the actual variance to the variance calculated as if the cells were independent. Thus the standard deviation of the null contrast is larger than the cell-independent standard error by a factor of \(\sqrt{1+(m-1)\rho}\). The donor-level variance formula gives the same result when the donors are equally weighted and each contributes \(m\) cells.
Under the normal model, independent donor means support a two-sample \(t\)-test at the donor level. If the group variances are equal, the pooled test uses \(2D-2\) degrees of freedom; if they are unequal, a Welch test can be used. In contrast, the usual cell-level test's independent-observation assumptions are false. Its nominal \(p\)-value is therefore not calibrated by the usual independent-cell \(t\) reference distribution. The design-effect calculation predicts the direction and approximate scale of the problem; simulation measures the behavior of the chosen tests under specified settings.
Worked Example: Quantifying Dependence Before Simulating
Suppose a hypothetical experiment has 8 donors, with 5 cells from each donor. Assume the null model has \(\tau^2=0.36\) and \(\sigma^2=0.64\), so a cell's marginal variance is 1 expression-unit squared. Calculate the within-donor correlation, the design effect, and the rough independent-cell equivalent of the 40 observations.
The within-donor correlation is $$ \rho=\frac{0.36}{0.36+0.64}=0.36. $$ With \(m=5\) cells per donor, the design effect is $$ 1+(m-1)\rho=1+4(0.36)=2.44. $$ Thus the cell-independent calculation understates the variance of a group mean by a factor of 2.44, or its standard error by \(\sqrt{2.44}\approx1.562\), relative to the model-based variance.
The familiar effective-size approximation is \(40/2.44\approx16.39\) independent-cell equivalents. This does not mean the experiment contains 16.39 independent donors; it summarizes the loss of precision in this particular equal-cluster setting. The actual independent biological units remain the 8 donors.
Designing a Null Simulation
A useful simulation states its data-generating assumptions, analysis rules, and performance measure in advance. For each repetition, draw independent donor effects, then draw cell residuals conditional on those effects. Generate both groups with the same \(\mu\), so every simulated data set satisfies the null. Apply both methods to the same generated data set; this paired comparison makes differences in decisions attributable to the analysis rather than to separate random data sets.
Choose donors per group \(D\), cells per donor \(m\), and variance components \(\tau^2,\sigma^2\). State that the group means are equal under the simulated null.
For each donor, draw one shared \(b_{gd}\); for each cell, draw its own \(\varepsilon_{gdi}\). Form \(Y_{gdi}=\mu+b_{gd}+\varepsilon_{gdi}\) in both groups.
Record whether the cell-independent test rejects at the prespecified \(\alpha\), and whether the donor-level test rejects at that same \(\alpha\).
Repeat \(R\) times. For each method, divide its number of rejections by \(R\) to estimate its null rejection probability.
For a test whose assumptions and calibration are correct, this null rejection probability is the Type I error rate and should be near \(\alpha\) over repeated data sets. A simulated rate will not equal \(\alpha\) exactly: it has Monte Carlo variation. If \(\hat{q}=K/R\) is the observed rejection proportion from \(K\) rejections, its approximate Monte Carlo standard error is \(\sqrt{\hat{q}(1-\hat{q})/R}\). This uncertainty describes simulation noise, not uncertainty about a biological effect.
Worked Example: A Donor-Level Test From Hypothetical Data
Consider hypothetical donor-level means in expression units. Group 0 has donor means \(0.8, 1.0, 1.2, 1.0\); group 1 has \(0.4, 0.6, 0.8, 0.6\). Assume independent donors and, for this illustration, normally distributed donor means with equal variances. The observed group means are \(\bar{x}_0=1.0\) and \(\bar{x}_1=0.6\), a difference of \(0.4\) expression units.
In each group, the deviations from its mean are \(-0.2,0,0.2,0\), whose squared deviations sum to \(0.08\). Therefore each sample variance is \(0.08/(4-1)=0.02667\) expression-units squared. The pooled variance is also \(s_p^2=0.02667\). The standard error for the difference is $$ SE=s_p\sqrt{\frac{1}{4}+\frac{1}{4}} =\sqrt{0.02667(0.5)} \approx0.1155\ \text{expression units}. $$
The test statistic is $$ t=\frac{0.6-1.0}{0.1155}\approx-3.464, \qquad df=4+4-2=6. $$ The two-sided \(p\)-value is approximately \(0.013\). Under the stated test assumptions, this is evidence against equal donor-population means at the 5% level. It is not evidence that the group difference is caused by a treatment; the example describes an association between group and donor-level expression.
Worked Example: Comparing Simulated False-Positive Rates
Now suppose a hypothetical simulation uses \(D=12\) donors per group and \(m=10\) cells per donor, with \(\tau^2=0.30\), \(\sigma^2=0.70\), and equal group means. Then the marginal cell variance is 1 and \(\rho=0.30\). The design effect is $$ 1+(10-1)(0.30)=3.70. $$ The cell-independent standard error, using the marginal variance and 120 cells per group, is \(\sqrt{2/120}\approx0.1291\) expression units. The actual null standard deviation of the difference in group means is $$ \sqrt{\frac{2(1)(3.70)}{120}}\approx0.2483\ \text{expression units}. $$
The ratio is \(0.2483/0.1291\approx1.923\). Thus a contrast that is roughly 1.923 times as variable as the cell-independent calculation assumes is tested against a reference distribution that treats its standard error as adequate. A normal approximation suggests a false-positive rate near $$ 2\Phi\left(-\frac{1.96}{\sqrt{3.70}}\right)\approx0.308, $$ rather than 0.05. This is an approximation, not an exact finite-sample result for the cell-level \(t\)-test.
For illustration, suppose one hypothetical run of \(R=10{,}000\) null data sets yielded 3,084 cell-independent rejections and 497 donor-level rejections at \(\alpha=0.05\). The estimated rates are \(3{,}084/10{,}000=0.3084\) and \(497/10{,}000=0.0497\). Their Monte Carlo standard errors are approximately $$ \sqrt{\frac{0.3084(1-0.3084)}{10{,}000}}\approx0.00462, \qquad \sqrt{\frac{0.0497(1-0.0497)}{10{,}000}}\approx0.00217. $$
The donor-level rate is compatible with a nominal 5% rate given simulation noise. The much higher cell-independent rate illustrates substantial null miscalibration under this design. The exact numbers from a simulation depend on its random draws; the model-based design effect explains why the cell-independent result is expected to be inflated.
Diagnostics, Alternatives, and Limits
Simulation is informative only if the generated data and the analysis match the question. Check that the simulated marginal variance and within-donor correlation are close to their targets, and confirm that each group truly has equal means under the null. Plot or tabulate rejection rates across several values of \(D\), \(m\), and \(\rho\), changing one feature at a time where possible. A large rejection rate that falls as \(\rho\) approaches zero is consistent with dependence-driven pseudoreplication. Also repeat the simulation with a larger \(R\) if Monte Carlo uncertainty is too large to distinguish plausible rates.
The principal comparison here is between treating cells as independent and analyzing independent donor summaries. A donor-level analysis makes the replication unit explicit, but summary-based analyses can discard information or require choices about how to summarize cells. Other cluster-aware models can represent donor dependence directly; their performance depends on their own assumptions and implementation. No method repairs confounding, biased donor selection, or an estimand that does not match the scientific question.
- Too few donors. Increasing cells per donor cannot substitute for independent donors when the target is a donor-population comparison. With few donors, donor-level estimates and tests can be unstable.
- Unequal cell counts. The simple design-effect formula assumes equal \(m\). Unequal cluster sizes can change the variance and the influence of donors, so simulate the actual sampling pattern if it is important.
- Misspecified correlation. A single common \(\rho\) may not describe all genes, cell types, or donors. Sensitivity analyses across plausible variance components are more informative than one chosen setting.
- Wrong interpretation. A simulation estimates the behavior of a specified procedure under a specified data-generating model. It does not estimate the probability that a biological null hypothesis is true.
Check Your Understanding
Answer each question, then compare with the solutions.
- If \(\tau^2=0.20\) and \(\sigma^2=0.80\), what is the within-donor correlation?
- With 6 cells per donor and \(\rho=0.25\), calculate the design effect.
- Why must both groups have the same mean in a simulation intended to estimate a null false-positive rate?
- A simulation has 1,000 null repetitions and 82 rejections. What is the estimated rejection rate?
- What is the main difference between the independent units in a cell-level test and in a donor-level test?
Solutions. (1) \(\rho=0.20/(0.20+0.80)=0.20\). (2) \(1+(6-1)(0.25)=2.25\). (3) Equal means ensure that the null hypothesis is true in each generated data set, so rejection frequency measures null behavior rather than power against a nonzero difference. (4) \(82/1{,}000=0.082\), or 8.2%. This exceeds 5%, though its Monte Carlo uncertainty should be considered before interpreting the discrepancy. (5) The cell-level test counts cells as independent observations; a donor-level test uses independent donors as observations, with cells contributing to each donor's summary.