Why Linkage Disequilibrium Matters for Multiple Testing
A genome-wide association study (GWAS) may test millions of SNPs for association with the same phenotype. Nearby SNPs are often in linkage disequilibrium (LD): their alleles, and therefore their genotype dosages, are correlated in the study population. When tests use the same people and phenotype, that genotype correlation can make their test statistics—and their p-values—dependent.
Dependence matters because the familiar Benjamini–Hochberg (BH) false discovery rate (FDR) guarantee has assumptions. BH controls FDR under independence of all p-values, and it also controls FDR under the positive regression dependence on a subset (PRDS) condition described in the earlier material on PRDS. LD is a biological and statistical source of dependence; it is not by itself a proof that the p-values satisfy PRDS. In particular, the conclusion depends on the test construction, the null hypotheses, and whether tests are one-sided or two-sided.
From Correlated Genotypes to Correlated Test Statistics
Consider a simple score-test setting. For participant \(i\), let \(Y_i\) be a quantitative phenotype and let \(G_{ij}\) be the dosage for SNP \(j\), usually coded as 0, 1, or 2 copies of a specified allele. Suppose the phenotype has been adjusted for any covariates, and write the residual under the null as \(\varepsilon_i\). For clarity, first consider the global null of no effect at either SNP and the model \(Y_i=\varepsilon_i\), where the residuals are independent, have mean 0, and have variance \(\sigma^2\). Assume the two genotype columns are fixed for this calculation. Define \(x_j\) as the centered genotype column for SNP \(j\), and use the standardized score statistic \(Z_j=x_j^\mathsf{T}\varepsilon/(\sigma\lVert x_j\rVert)\), where \(\varepsilon\) is the vector of residuals.
The zero means and common variance imply that the scores have mean 0 and variance 1. Their covariance is the correlation between the two centered genotype columns. These facts do not require Gaussian residuals. If the residuals are Gaussian, the score pair is jointly normal; otherwise, joint normality is not guaranteed, though a normal approximation may be reasonable under suitable large-sample conditions. The distinction is important: knowing a covariance does not establish a joint distribution or a PRDS property.
The second equality follows because independent residuals with variance \(\sigma^2\) have covariance matrix \(\sigma^2 I\), so \(\operatorname{Cov}(x_1^\mathsf{T}\varepsilon,x_2^\mathsf{T}\varepsilon) =\sigma^2 x_1^\mathsf{T}x_2\). The same calculation gives \(\operatorname{Var}(Z_j)=1\). This is a useful justification for why genotype LD can lead to test-statistic dependence under a shared null. It is not a universal model for all GWAS analyses: covariate adjustment, estimated variance, population structure, relatedness, and different samples can change the dependence structure.
Worked Example: Calculating Dosage Correlation in a Hypothetical Sample
Suppose eight hypothetical participants have the following paired binary genotype values at two SNPs. The values 0 and 1 stand in for two possible allele states, to keep the calculation small.
| Genotype pair | Count |
|---|---|
| (0, 0) | 3 |
| (0, 1) | 1 |
| (1, 0) | 1 |
| (1, 1) | 3 |
Each genotype column has four 1s, so each sample mean is \(4/8=0.5\). The sum of the products of paired genotype values is 3. Thus the centered cross-product is \(3-(4)(4)/8=1\). The centered sum of squares for either column is \(4-4^2/8=2\). Therefore the sample correlation is
In this hypothetical sample, the two genotype columns have positive correlation 0.5. Under the simple shared-null score model above, the score statistics also have covariance and correlation 0.5. The calculation demonstrates a possible route from LD to statistic dependence; it does not establish that the p-values from an arbitrary GWAS pipeline are PRDS.
What Positive Dependence Can—and Cannot—Guarantee
For BH, the relevant positive-dependence condition is PRDS on the set of true null hypotheses. Informally, as a true-null p-value becomes larger, the conditional probability of an increasing event in the full p-value vector does not decrease. An increasing event is one that remains true if any p-value is made larger. This condition formalizes a particular pattern of positive dependence; it is stronger and more specific than saying that two genotypes or two statistics have positive correlation.
When all p-values are valid and independent, BH controls FDR at no more than \((m_0/m)\alpha\), where \(m\) is the total number of tests, \(m_0\) is the number of true nulls, and \(\alpha\) is the target level. Under valid p-values and PRDS on the true nulls, BH controls FDR at no more than \(\alpha\). These are guarantees about the expected false discovery proportion, not statements that a particular rejection set has exactly that proportion of false discoveries.
In some Gaussian testing settings, nonnegative correlations among one-sided test statistics provide conditions under which the resulting p-values have the needed positive dependence. But this is a setting-specific result, not a rule that can be applied from an LD plot alone. Two-sided p-values discard the sign of the statistic, and their dependence must be assessed for the actual testing construction. Mixed null and nonnull effects, covariate adjustment, sample overlap, and population structure can also affect the joint behavior. Consequently, positive LD between SNPs is evidence that dependence is plausible, not sufficient evidence for a formal BH guarantee.
Worked Example: Interpreting a Score Correlation Under the Model
Suppose two centered genotype columns in a hypothetical study have squared lengths \(\lVert x_1\rVert^2=20\) and \(\lVert x_2\rVert^2=18\), and cross-product \(x_1^\mathsf{T}x_2=12\). Assume the shared-null model with independent residuals of mean 0 and variance \(\sigma^2\). The score correlation is
The scores have unit variance, so their covariance is also approximately 0.6325. If, in addition, the residual vector is Gaussian, the score pair is jointly normal with this correlation. Without that additional Gaussian assumption, the covariance result still holds but joint normality does not follow from the stated assumptions. In neither case does this calculation alone verify PRDS for a set of p-values; it only characterizes the pairwise score dependence in this specified model.
BH Decisions When Nearby SNPs Are Correlated
A BH calculation is performed on the p-values in the usual way, whether or not the tests are independent. Its rejections are determined by the ordered p-values and the BH thresholds. Dependence enters the question of what error-rate guarantee is justified, not the arithmetic of the rejection rule.
Worked Example: BH Rejections Are Not Independent Confirmations
Consider five hypothetical GWAS p-values at target level \(\alpha=0.05\), ordered as \(0.004,\ 0.012,\ 0.018,\ 0.20,\ 0.70\). The BH threshold at rank \(i\) is \(i\alpha/m\), with \(m=5\). The thresholds are \(0.01,\ 0.02,\ 0.03,\ 0.04,\ 0.05\), respectively.
At rank 1, \(0.004\leq0.01\); at rank 2, \(0.012\leq0.02\); and at rank 3, \(0.018\leq0.03\). At rank 4, \(0.20>0.04\), and at rank 5, \(0.70>0.05\). The largest passing rank is 3, so the BH step-up rule rejects the first three hypotheses.
Suppose two of these p-values correspond to SNPs in LD. The algorithm still rejects all three, but those correlated rejections should not be described as three independent confirmations of an association. They may reflect overlapping information about a shared genetic signal. Whether the nominal FDR guarantee applies depends on the joint dependence conditions, not simply on the fact that the BH arithmetic was carried out correctly.
Practical Interpretation, Alternatives, and Failure Modes
LD is useful for understanding why a discovery list may contain clusters of nearby SNPs. A cluster can arise because several correlated markers track one underlying signal, even when only one marker is causal—or when none is causal and the pattern is due to confounding or chance. Association at correlated markers does not by itself identify a causal variant. That question requires additional assumptions and analyses beyond FDR adjustment.
When the dependence structure is not known to meet independence or PRDS conditions, the Benjamini–Yekutieli procedure gives an FDR guarantee under arbitrary dependence, provided the individual p-values are valid under their nulls. It uses a more stringent threshold, multiplying the BH denominator by the harmonic sum \(H_m=\sum_{i=1}^{m}1/i\). This conservatism can be substantial when \(m\) is large. Choosing between BH and BY is therefore a choice about assumptions and power, not a claim that one procedure removes LD from the data.
A frequent mistake is to argue, “The SNPs have positive LD, so BH is guaranteed to control FDR.” The correction is that positive genotype correlation does not automatically imply PRDS for the p-values in the analysis. The test type, joint distribution, and null configuration matter. A second mistake is to count nearby BH discoveries as independent replications. Their shared genotype information makes that interpretation unjustified; independent evidence requires an appropriate independent data source or design.
Other warning signs include treating a pairwise correlation estimate as a complete description of dependence across millions of tests, ignoring overlap in participants, or assuming that valid marginal p-values are necessarily independent. LD estimates are also population- and sample-dependent. A dependence guarantee should be justified for the actual analysis or replaced by a procedure whose guarantee covers the dependence allowed. Diagnostics can describe observed patterns, but they do not by themselves prove PRDS.
Check Your Understanding
Answer each question, then compare your reasoning with the solutions.
- Under the shared-null score model, two centered genotype columns have correlation \(-0.4\). What is the correlation between their score statistics?
- Independent residuals have mean 0 and common variance. Does this alone imply that two score statistics are jointly normal? What additional assumption would imply joint normality in the model described here?
- A researcher says that positive LD between two SNPs proves PRDS for their two-sided GWAS p-values. What is the flaw in this argument?
- BH rejects four nearby SNP hypotheses in LD. Can those rejections be called four independent confirmations? Explain.
Solutions. (1) The score correlation is the genotype-column correlation, \(-0.4\), under the stated shared-null model. (2) No. The assumptions give the covariance structure but do not guarantee a normal distribution; Gaussian residuals imply joint normality here. (3) Genotype correlation alone does not establish PRDS for p-values. The result depends on the test construction and joint distribution, and two-sided tests require particular care. (4) No. The SNPs may carry overlapping information because of LD, so the rejections are not automatically independent confirmations.