Worked GWAS Error-Control and Calibration Analysis

From Association Results to Defensible Discoveries

A genome-wide association study can produce thousands or millions of p-values. A small p-value is evidence against a variant-specific null hypothesis only to the extent that the test is calibrated; even well-calibrated p-values require an explicit approach to multiplicity. A QQ plot, a summary inflation statistic, and an error-control procedure answer related but different questions. Taken together, they help describe a results set, but they do not establish that every reported association is free of bias or will replicate.

This tutorial follows one hypothetical GWAS results set through those checks. The workflow is to inspect the overall p-value pattern, assess whether a chosen calibration adjustment is justified, apply a multiplicity procedure to the resulting p-values, and state what the discoveries do and do not show. The numerical example uses one-degree-of-freedom chi-square association tests so the calibration calculation is explicit.

Definition: A QQ plot compares ordered observed test statistics or p-values with the values expected under a specified null distribution. The genomic inflation factor, $\lambda$, is the ratio of the observed median association-test chi-square statistic to its null median. A calibration adjustment transforms test statistics or p-values to better match a reference null distribution; a multiplicity procedure then controls a stated error rate across the family of tests.

Inspecting the QQ Plot and Inflation Summary

For $m$ tests with ascending p-values $p_{(1)}\leq \cdots \leq p_{(m)}$, a common p-value QQ plot places expected values $-\log_{10}((i-0.5)/m)$ against observed values $-\log_{10}(p_{(i)})$, for ranks $i=1,\ldots,m$. The half-rank offset avoids plotting an expected p-value of exactly zero. Under a correctly specified null model with calibrated tests, the bulk of the points should roughly follow the null reference line, subject to sampling variation. Points rising above that line indicate smaller observed p-values than expected at those ranks.

The shape matters. A broad upward displacement through the bulk of the plot may indicate test-statistic inflation, but it can also reflect many genuine small effects. An upward bend confined to the most extreme ranks can be consistent with a limited number of strong associations. A QQ plot describes the pattern; it does not identify its cause. The genomic inflation factor gives a compact summary of the median chi-square statistics, but it also cannot by itself distinguish confounding from widespread polygenic signal.

For one-degree-of-freedom chi-square statistics, the null median is approximately $0.4549$. Thus, if the observed median is $0.4913$, the genomic inflation estimate is $0.4913/0.4549$. Interpretation depends on the test model, the sample and phenotype, and the broader QQ-plot pattern. It is not a universal pass/fail threshold.

Worked Example: Reading a Hypothetical QQ Plot

Suppose a hypothetical GWAS tests 10,000 variants using one-degree-of-freedom chi-square tests. Its QQ plot follows the reference line through most of the distribution, lies slightly above the line across the middle ranks, and rises more sharply for the most extreme ranks. The median observed chi-square statistic is $0.4913$.

The inflation estimate is

\lambda=\frac{\operatorname{median}(T_1,\ldots,T_{10000})}{0.4549} =\frac{0.4913}{0.4549}=1.080

Here $T_i$ denotes the association-test chi-square statistic for variant $i$. The estimate summarizes a median that is about 8.0% above the one-degree-of-freedom null median. In context, the broad but modest QQ-plot displacement and $\lambda=1.080$ warrant a calibration check; neither proves that the results are confounded. The sharper extreme tail is compatible with some strong associations, although the plot alone cannot establish that those associations are genuine.

The interpretation assumes the statistics are intended to follow a one-degree-of-freedom chi-square null distribution. If the tests have different degrees of freedom or use a different null reference, the corresponding null median must be used instead.

When and How to Apply Genomic Control

A common genomic-control adjustment divides each one-degree-of-freedom chi-square statistic by an estimated $\lambda$, then recomputes its upper-tail p-value under the one-degree-of-freedom chi-square distribution. The transformation reduces the test statistic when $\lambda>1$ and increases it when $\lambda<1$. It changes the test evidence, not the estimated genetic effect itself.

Calibration formula: For an observed one-degree-of-freedom chi-square statistic $T_i$ and an inflation estimate $\lambda>0$, the genomic-control adjusted p-value is $$ p_{i,\mathrm{GC}}=P\left(\chi^2_1\geq \frac{T_i}{\lambda}\right). $$ This formula presumes that scaling the test statistic by the estimated factor is an appropriate calibration model. It is not a general correction for every source of bias or model misspecification.

Using $\lambda$ estimated from the same association results can be problematic when many true effects contribute to the distribution: the estimate may partly reflect real polygenic signal, and correcting by it can make genuine evidence less significant. Conversely, a modest $\lambda$ does not rule out localized problems, such as a poorly calibrated subset of tests. The QQ plot, study design, ancestry and relatedness handling, and choice of test model all inform whether this simple adjustment is defensible. Empirical-null calibration is another possible approach, but it requires its own model and checks.

Worked Example: Adjusting One Association Statistic

Suppose one hypothetical variant has an observed one-degree-of-freedom chi-square statistic $T=25.00$, and the study uses the $\lambda=1.080$ estimate above. Assume for this calculation that genomic-control scaling is considered appropriate for these tests.

First scale the statistic:

T_{\mathrm{GC}}=\frac{25.00}{1.080}=23.15

Then calculate the upper-tail probability under a one-degree-of-freedom chi-square null:

p_{\mathrm{GC}}=P(\chi^2_1\geq 23.15)\approx 1.5\times 10^{-6}.

The unadjusted statistic $25.00$ corresponds to a p-value of approximately $5.7\times10^{-7}$; after scaling, the p-value is larger, about $1.5\times10^{-6}$. These values quantify evidence against the variant-specific null under the stated test model. Neither is the probability that the null hypothesis is true. Also, genomic control has not changed the variant’s estimated effect or supplied evidence that the association is causal.

Applying Error Control to a Calibrated Results Set

Once the p-values to be analyzed have been selected, the multiplicity procedure must match the intended error criterion. Family-wise error rate (FWER) is the probability of making at least one false rejection in the family. The false discovery rate (FDR) is the expected proportion of false rejections among all rejections, with that proportion defined as zero if there are no rejections. FWER control is more stringent in many settings; FDR control can yield more discoveries while allowing a controlled expected false-discovery proportion.

Holm’s step-down procedure controls FWER for valid marginal p-values under arbitrary dependence. Order the p-values and compare $p_{(i)}$ in turn with $\alpha/(m-i+1)$, beginning at $i=1$. Reject sequentially until the first comparison that fails; that hypothesis and all later ones are not rejected. The sequential stopping rule is essential: a later p-value cannot be rejected after an earlier step has failed.

The Benjamini–Hochberg (BH) procedure controls FDR under independence and under certain positive-dependence conditions. Order the p-values, find the largest rank $k$ satisfying $p_{(k)}\leq k\alpha/m$, and reject the hypotheses at ranks $1,\ldots,k$. In a GWAS, linkage disequilibrium means tests are dependent; BH is appropriate only when the dependence conditions for its use are justified. The specific conditions and alternatives are not inferred from the p-values alone.

Worked Example: Holm and BH on the Same Hypothetical GWAS

Consider a hypothetical set of $m=10{,}000$ variant tests with $\alpha=0.05$. Suppose the p-values have already been calibrated using a justified procedure, and the first five ordered values are $8.0\times10^{-8}$, $1.2\times10^{-6}$, $3.8\times10^{-6}$, $4.9\times10^{-6}$, and $6.0\times10^{-6}$. Assume the remaining ordered p-values are such that the BH inequality holds through rank 10 and fails at every rank above 10. Also assume the dependence conditions for BH are satisfied in this hypothetical analysis.

For Holm’s procedure, the first four comparisons are:

\begin{aligned} p_{(1)}&=8.0\times10^{-8}\leq \frac{0.05}{10000}=5.000\times10^{-6},\\ p_{(2)}&=1.2\times10^{-6}\leq \frac{0.05}{9999}\approx5.001\times10^{-6},\\ p_{(3)}&=3.8\times10^{-6}\leq \frac{0.05}{9998}\approx5.001\times10^{-6},\\ p_{(4)}&=4.9\times10^{-6}\leq \frac{0.05}{9997}\approx5.002\times10^{-6}. \end{aligned}

The fifth comparison fails:

6.0\times10^{-6}> \frac{0.05}{9996}\approx5.002\times10^{-6}.

Holm therefore rejects the first four null hypotheses and stops at rank five. The result is four FWER-controlled rejections at the 0.05 level for this family of 10,000 tests.

For BH, the rank-10 threshold is

\frac{10(0.05)}{10000}=5.0\times10^{-5}.

Given the stated hypothetical ordering, rank 10 satisfies its threshold and no higher rank does, so the largest qualifying rank is $k=10$. BH rejects the first 10 null hypotheses. This is an FDR-controlled set under the assumed dependence conditions, not a claim that exactly 5% of these 10 discoveries are false. The difference from Holm is not a contradiction: the procedures target different error rates.

Justification for the Holm Stopping Rule

Holm’s guarantee can be understood through the Bonferroni bound. Suppose there are $m_0$ true null hypotheses among the $m$ tests. If Holm makes any false rejection, consider the first true null reached in the ordered sequence. At that step, all earlier rejected hypotheses are false nulls, so the number of hypotheses remaining is at least $m_0$. The true-null p-value at this step must pass a threshold no larger than $\alpha/m_0$. For any one true null with a valid p-value, the probability of a p-value at most $\alpha/m_0$ is at most $\alpha/m_0$. Summing this bound over the $m_0$ true nulls gives an upper bound of $\alpha$ on the probability of any false rejection. This argument does not require independence between tests.

The reasoning depends on valid marginal p-values: under each true null, the chance of a p-value at most $u$ must be no greater than $u$. If a calibration failure makes p-values systematically too small, a multiplicity adjustment cannot by itself restore this validity. Calibration and error control are linked, but they are not interchangeable steps.

Reporting the Analysis and Its Limits

A useful report makes the analysis auditable: identify the test family and test statistic, describe what the QQ plot and $\lambda$ indicate, state whether calibration was applied and how, give the error-control method and $\alpha$, and report the number of discoveries under that method. If the result depends on an assumption about dependence or calibration, state that assumption rather than treating it as a property established by the output.

For the hypothetical results above, a suitably limited conclusion would be: “The one-degree-of-freedom test statistics showed modest median inflation ($\lambda=1.080$) and a QQ-plot displacement that was most pronounced in the extreme tail. After the stated calibration, Holm rejected four of 10,000 variant null hypotheses at family-wise level 0.05; BH rejected 10 under the assumed dependence conditions and an FDR target of 0.05. These are statistical association findings under the specified models, not evidence by themselves of causal effects or replication in another cohort.”

The common mistake is to see $\lambda>1$, divide every statistic by $\lambda$, and conclude that the adjusted results are automatically unbiased. That reasoning treats the inflation estimate as a diagnosis of cause and the transformation as a universal remedy. The correction is to inspect the QQ-plot shape and study context, justify the calibration model, and then apply an error-control procedure to valid p-values. Likewise, a small adjusted p-value is not a probability that a null is true, and FDR control does not label each individual discovery as true or false.

Key takeaway: Calibration assesses whether the test statistics follow a credible null reference; multiplicity procedures control a specified error rate across tests. Report both steps, their assumptions, and the limits of the resulting association claims.

Check Your Understanding

Use the definitions and procedures above to answer these questions.

In a one-degree-of-freedom GWAS analysis, the observed median chi-square statistic is $0.4549$. What is $\lambda$, and what does that value summarize?
A QQ plot follows the reference line through most ranks but rises above it at the extreme tail. Give one interpretation consistent with that pattern and explain why the plot does not prove that the strongest results are genuine.
Why does Holm stop rejecting after the first ordered p-value that fails its step-down comparison?
In the worked GWAS, why do Holm and BH produce different numbers of discoveries without one procedure necessarily being incorrectly applied?
Why does genomic-control adjustment not establish that an association is causal or that its p-value is correctly calibrated under every possible source of bias?