26 Statistics and Experimental Design: Concepts
26.1 The problem
Every part of this book before this one ends in a table: a variant call set, a count matrix, a peak list, an abundance table. Each has tens of thousands of rows and, attached to each row, a number that is supposed to say whether that row means anything.
This part is about that number. It is the shortest stretch of a pipeline and the one where the most claims are lost, because two things are true at once that feel contradictory. Almost none of the interesting decisions here are about which test to run — the tests are standard and the packages are correct. And almost all of the ways an analysis goes wrong at this stage were fixed long before the test was run: when samples were assigned to processing days, when a sample size was chosen by what the budget allowed, when a threshold was picked after looking at what it produced.
Two failures look identical in a results table and have nothing in common. In the first, the effect is real and you cannot see it: too few samples, too much variance, a threshold made severe by the number of tests. In the second, the effect is not real and you can see it: enough tests, or enough flexibility in how the analysis was run, that something was bound to come out. The first is a power problem and is fixed before the experiment exists. The second is a selection problem and is fixed by accounting for the selection. Neither remedy touches the other, and most published advice is about the second while most damage is done by the first.
26.2 What a p-value is a statement about
A p-value is the probability of observing data at least as extreme as yours if the null hypothesis were true. It is a statement about data under a hypothesis, not about a hypothesis given data, and the two are not interchangeable: the fraction of significant results that are wrong depends on how many of the hypotheses tested were null to begin with, which the p-value cannot know (Ioannidis 2005).
The property that makes p-values useful at genome scale is narrower and more practical. If the null is true and the test is correctly specified, p-values are uniformly distributed on the interval from zero to one. Every value is as likely as every other. That is not a technicality to be recited; it is the single most useful diagnostic in this chapter, because it makes an entire analysis checkable in one picture.
Read the shape before reading the numbers. A flat background with a spike near zero is the healthy case, and the height of the flat part is an estimate of how many of your tests were null — a quantity that turns out to be worth money, as the next section explains. A flat histogram with no spike means there is nothing to find, or nothing you were powered to find. Anything else is a finding about your test rather than about your biology: a slope rising toward one usually means a conservative or discrete test, a hump in the middle means the null is mis-specified, and a spike at one means something is badly wrong, often a one-sided test pointed the wrong way.
Every multiple testing procedure in the next section assumes this histogram is flat under the null. Applying one to a distribution that is visibly not flat produces adjusted values with no guarantee attached, and nothing downstream will complain. The plot costs one line and is the only check that sees the whole analysis at once.
26.3 Two guarantees, and they are not the same
Test 20,000 genes at a threshold of 0.05 with nothing going on, and about 1,000 of them come out significant. That much is arithmetic. The interesting question is what to do about it, and there are two answers that sound similar and protect completely different things.
Family-wise error rate control bounds the probability of making any false discovery. Bonferroni’s threshold — the target rate divided by the number of tests — achieves it, and Holm’s step-down procedure achieves it while rejecting at least as much, so there is no reason to prefer Bonferroni except that it is easier to explain.
False discovery rate control bounds the expected proportion of false discoveries among the things you called significant (Benjamini and Hochberg 1995). The Benjamini–Hochberg procedure sorts the p-values and compares each to a threshold that rises with its rank, which is why it can accept a p-value tens of times larger than Bonferroni would.
The difference is what you are protecting. Family-wise control protects a claim: if you report that something here is real, you want to be confident that at least one thing is. False discovery rate control protects a list: you accept that some entries are wrong in exchange for the list being long enough to be useful, and you fix what fraction you will tolerate.
Genomics almost always produces a list, which is why false discovery rate control is the default here and family-wise control is the default in a clinical trial with one endpoint. The choice follows from what happens next to the output, not from the size of the experiment.
Two refinements are worth knowing because they change the answer rather than the presentation. The Benjamini–Hochberg threshold divides by the total number of tests, which implicitly treats all of them as null and is conservative when many are not; estimating the true null proportion from the flat part of the histogram and dividing by that instead recovers the difference, which is what a q-value does (Storey 2002; Storey and Tibshirani 2003). And the quantity most people want when they look at a single row is not the false discovery rate of the whole list but the probability that this gene is null given its statistic — the local false discovery rate, which needs the shape of the distribution rather than just its tail (Efron 2004).
Dependence is the assumption people worry about most and it is usually the benign direction. Benjamini–Hochberg is valid under independence and under a form of positive dependence that co-expressed gene modules satisfy; the alternative that holds under arbitrary dependence costs a factor of the harmonic number and is rarely necessary (Benjamini and Yekutieli 2001).
A list of 400 genes at a false discovery rate of 5% is a list with about twenty wrong entries in it, and nothing tells you which twenty. That is a perfectly good basis for choosing what to follow up and a poor basis for a claim about any individual gene on it. If one row is going to become an experiment, the question you want answered is the local one.
26.4 The null you assumed is not always the null you have
Everything above rests on the theoretical null being the right null. In association studies it frequently is not: residual population structure, unmodelled cell composition, or an unrecorded batch shifts the entire distribution of test statistics, and the histogram in Figure 26.1 tilts. Correcting for multiple testing then controls an error rate against a reference that is wrong, which controls nothing.
The response is to estimate the null from the bulk of the statistics themselves, on the reasoning that most tests are null, so whatever the middle of the distribution looks like is the null (Efron 2004). It is a real remedy with a real cost, and the cost is not fixable: a genuinely widespread biological signal and a confounder both inflate the statistics, and no amount of cleverness applied to the statistics alone can tell them apart. Which is the argument for the next section, where the distinction is decided.
26.5 Batch effects and confounding
A batch effect is a systematic difference between groups of samples caused by how they were processed rather than by what they are — day, reagent lot, flow cell, operator, position on a plate. Confounding is what happens when that grouping lines up with the variable you care about.
The first is a nuisance and is estimable, which makes it a statistics problem. The second is not a statistics problem at all: it is a design that cannot answer the question, and no method recovers it, because the two effects are the same vector in the data (Leek et al. 2010).
Batch effects are not a minor correction. In a typical high-throughput experiment the processing batch is among the largest sources of variation in the matrix — often larger than the biological effect under study.
The picture is the same whether the design is sound or ruined, which is the point worth internalizing. In the simulation above, cases and controls are balanced across all three batches, so the batch term can be estimated and removed and the condition effect survives. Had every case been processed in batch one and every control in batch three, this figure would look approximately identical — and the analysis would be dead, with nothing in the data to say so.
Three responses to a batch effect, in the order to prefer them:
Put it in the model. If you recorded the batch, it is a covariate, and a linear model with a batch term is doing the correction as part of the inference rather than before it. This is the correct default, it costs one column in the design matrix, and it is what the differential expression tools of Section 17.3 are built to do.
Estimate it when you did not record it. Some unwanted variation is real and unlabelled — degradation, an undocumented process change. Surrogate variable analysis and related factor methods estimate such components from the data and hand them back as covariates, which then go into the same design matrix.
Produce a corrected matrix, but only for the things that need one: visualization, clustering, a machine learning input. Testing on a corrected matrix understates uncertainty, because the test is never told that the correction consumed degrees of freedom, and the resulting confidence intervals are too narrow and the p-values too small (Nygaard et al. 2016).
Randomize what you can, block what you cannot, and record everything you did not control. Processing samples in an order that mirrors the biology — cases first, then controls; one group per flow cell — is the single most expensive mistake available at this stage, and it is made before any data exist (Auer and Doerge 2010). The corresponding single-cell case is Section 18.7, where integration methods remove whatever variation the labels you supply are associated with, disease effect included.
26.6 Power is a property of the design
Power is the probability of detecting an effect that is really there. It depends on the effect size, the variance, the sample size and the significance threshold — and at genome scale the threshold is not 0.05, because the multiple testing burden of the previous section enters here as a direct cost in samples.
An underpowered experiment does not simply fail to find things. It distorts what it does find, in a direction that is easy to miss: when power is low, only the largest apparent effects clear the threshold, so every effect you report is overestimated, and a larger fraction of your significant results are false (Button et al. 2013; Ioannidis 2005). A field that publishes small studies therefore accumulates an effect size literature that is systematically wrong, and each new study is designed against those inflated numbers.
Post-hoc power — computing power from the effect you observed, after the fact — is not evidence about anything. It is a monotone function of the p-value you already have, so a non-significant result always yields low computed power and the exercise cannot distinguish a real absence from a missed effect. The useful question after a null result is not “was I powered?” but “what effect sizes can I now rule out?”, which is a confidence interval.
Doing it properly needs three numbers you do not have yet: the smallest effect worth detecting, the variance, and the threshold your multiple testing burden implies. The first is a scientific judgment and nobody else can supply it. The second comes from a pilot or a public dataset in the same tissue. The third follows from the first section of this chapter. A power calculation is therefore an argument with three stated assumptions, and its value is mostly in having had to state them.
Cells are not samples. Ten thousand cells from three donors carry roughly three donors’ worth of information about a population, not ten thousand, and a test that treats each cell as independent will report extraordinary p-values for nothing. The same applies to technical replicates, to multiple sections from one tissue block, and to repeated measures on one subject. The design variable in single-cell work is the number of donors; cells per donor and reads per cell are a budget split under it.
26.7 Inference after selection
The last failure mode is the one that survives every remedy above, because it does not look like a statistical decision at all.
A p-value is only interpretable if the hypothesis it tests was fixed before the data were seen. Use the data to choose the hypothesis and the guarantee is gone — not weakened, gone — and the resulting number can be arbitrarily small with nothing behind it. This has a name in every field that has been caught by it; in neuroimaging it is double dipping, and the canonical form is selecting voxels because they respond, then testing whether they respond (Kriegeskorte et al. 2009).
Genomics does this constantly, and usually without noticing:
- Cluster cells, then test for genes that differ between clusters. The clusters were built from those genes, so the test is asking whether a split made on the basis of a difference exhibits that difference (Zhang et al. 2019).
- Raise the clustering resolution until a population of interest appears, then report that it appears — the specific case Section 18.6 warns about.
- Filter to the features that look most different, then test them.
- Try several normalizations, batch corrections or differential expression tools, and report the one whose results look best.
The last is the widest and least visible: every analysis choice made after seeing results is a comparison you have made and not counted.
Not all data-dependent filtering is illegitimate, and the distinction is sharp. Filtering on a statistic that is independent of the p-value under the null — mean expression, total count — removes tests that had no power without touching the null distribution of what remains, and it increases power. It is the same principle that lets a covariate be used to weight hypotheses. Filtering on the test statistic itself, or on anything correlated with it under the null, is selection, and it invalidates the tests that survive.
The remedies all buy back independence by paying data for it. Splitting samples into a discovery half and a test half is the oldest and it is expensive at genomics sample sizes. Splitting the counts instead is the newer answer: a count can be thinned into two independent counts under the same model, so every sample contributes to both halves, and the selection can be done on one and the test on the other (Neufeld et al. 2024). Where the selection is a clustering, the alternative is to test the split itself rather than the genes downstream of it (Grabski et al. 2023). And where the selection is a regression model, the inference can sometimes be corrected exactly, conditional on the selection that was made (Tibshirani et al. 2016).
26.8 What to take forward
- A p-value histogram is the one diagnostic that sees the whole analysis; a distribution that is not flat under the null invalidates everything applied to it afterwards.
- Family-wise error control protects a claim, false discovery rate control protects a list. Genomics produces lists, which is why it uses the second.
- A false discovery rate is a statement about a list, not about any row in it. If one row is going to become an experiment, the local quantity is the one you want.
- A batch effect is estimable and a confounded design is not. The remedy for the first is a term in the model; there is no remedy for the second.
- Correcting a matrix and then testing on it produces intervals that are too narrow, so correct within the model wherever a model is going to be fitted.
- The multiple testing threshold is a cost in sample size, and it belongs in the power calculation you do before the experiment rather than the one you do afterwards, which is not evidence.
- Any hypothesis chosen using the data needs its selection accounted for, and the practical routes are splitting the data, splitting the counts, or testing the selection itself.
Section 27.1 names what to run.