Tutorials › High-Dimensional Genomic Data Analysis › Calibration of Mixed-Model Association P-Values

Mixed-Model GWAS, Relatedness and Inflation Diagnostics · Tutorial 265 of 795

Calibration of Mixed-Model Association P-Values

Learn to distinguish genuine association signal from evidence that mixed-model p-values are miscalibrated, and to evaluate calibration using null simulations and summary diagnostics.

Advanced 12 min read

What You'll Learn

  • Define p-value calibration under a specified null hypothesis and test procedure.
  • Explain why valid continuous null p-values are uniformly distributed.
  • Estimate empirical type I error from null simulations and assess its uncertainty.
  • Compare observed p-value distributions with null expectations without mistaking real signals for miscalibration.
  • Interpret a genomic-control median summary and recognize what it cannot diagnose.
  • Identify model, simulation, and dependence failures that can distort calibration checks.

What Calibration Means

A mixed-model association scan produces a p-value for each variant, but a small p-value is informative only if the test is calibrated: when the null hypothesis is true, the procedure should reject at approximately the advertised rate. Calibration is therefore a property of the test procedure under a specified data-generating null, not simply a visual feature of a particular genome-wide result.

For a variant \(j\), a typical null model includes fixed covariates and a random effect accounting for relatedness, for example \(y=X\beta+u+\varepsilon\), where \(u\sim N(0,\sigma_g^2K)\) and \(\varepsilon\sim N(0,\sigma_e^2I)\). An association test adds a variant effect, such as \(g_j\alpha_j\), and tests \(H_0:\alpha_j=0\). Calibration checks concern the p-values from the chosen Wald, score, or likelihood-ratio procedure when that null is true. They do not establish that the model is correct, nor do they determine whether any particular small p-value represents a true association.

Definition: A p-value procedure is calibrated under a null model if, for every threshold \(a\) between 0 and 1, the probability of obtaining a p-value at most \(a\), when the null hypothesis is true and the stated procedure is applied, is \(a\) for a continuous exact test. In practice, asymptotic tests and estimated nuisance parameters may yield approximate rather than exact calibration.

A GWAS result usually contains both null and non-null variants. Under a valid test, the p-values for true null variants should follow the null reference distribution, but genuine associations can add an excess of small p-values. Consequently, comparing every observed p-value with a uniform distribution is not a clean test of calibration when the set of null variants is unknown. Null simulations, or a defensible set of known-null tests, are more direct.

Why Uniform P-Values Are the Reference

The uniform reference follows from the definition of a continuous p-value. Suppose a test statistic \(T\) has continuous null cumulative distribution function \(F_0\), and use the upper-tail p-value \(P=1-F_0(T)\). For any \(a\) between 0 and 1,

$$ P(P\leq a\mid H_0) =P\{T\geq F_0^{-1}(1-a)\mid H_0\} =a. $$

The equality uses the fact that \(T\) follows \(F_0\) under the null. Thus, if the test is exact and continuous, the null p-value has a uniform distribution on \([0,1]\). The corresponding statement holds for a two-sided p-value when it is correctly constructed from the null distribution. For discrete tests, p-values are often conservative: the probability of being at most \(a\) can be less than \(a\), because only a limited set of p-values is attainable.

In mixed-model GWAS, the reference distribution may be approximate. Wald, score, and likelihood-ratio statistics often rely on large-sample or model-based approximations. Variance components are estimated, and the validity of the calculation depends on the fitted covariance structure, covariate specification, genotype coding, and implementation of the test. The earlier material on fitting a null model and conditional association tests gives the model-based setting; calibration checks whether the resulting procedure behaves as intended.

Null Simulations and Empirical Error Rates

A direct calibration assessment generates data under a null model, applies the same association pipeline, and compares the observed rejection frequency with its nominal level. In a mixed model, simulations should reflect the design features that matter: the fixed covariates, the relationship matrix \(K\), the residual and random-effect variation, and the analysis procedure, including how nuisance parameters are estimated. A parametric simulation from a fitted null model is a useful diagnostic, but it is a plug-in assessment: if the fitted model is wrong, the simulation may reproduce the wrong assumptions rather than reveal them.

For \(R\) independent null test outcomes, let \(I_r(a)\) equal 1 if the \(r\)th p-value is at most \(a\), and 0 otherwise. The estimated rejection rate is \(\hat q(a)=\sum_{r=1}^R I_r(a)/R\). Under valid calibration, each indicator has probability \(a\), so the count has a binomial distribution with mean \(Ra\) and variance \(Ra(1-a)\), provided the outcomes are independent. Its estimated rate has standard error \(\sqrt{a(1-a)/R}\). In real scans, markers are correlated through linkage disequilibrium; treating every marker as an independent replicate can therefore understate uncertainty. Independent simulated phenotype replicates, or uncertainty calculations that account for dependence, are safer units for inference.

Worked Example: Comparing a Simulated Rejection Rate With Its Target

Suppose a hypothetical calibration study produces 10,000 independent null test outcomes from repeated simulations. At the nominal level \(a=0.05\), 540 p-values are at most 0.05. The simulation preserves the covariates and relatedness structure, and each replicate is analyzed with the intended mixed-model procedure.

The estimated rejection rate is \(540/10{,}000=0.0540\), or 5.40%. If the procedure is calibrated, the expected number of rejections is \(10{,}000(0.05)=500\). The null standard deviation of that count is

$$ \sqrt{Ra(1-a)} =\sqrt{10{,}000(0.05)(0.95)} =\sqrt{475} \approx 21.79. $$

The observed excess is \(540-500=40\), or \(40/21.79\approx1.84\) null standard deviations. This is not, by itself, strong evidence of a substantial calibration failure: ordinary simulation variability can produce differences of this size. The estimated rate is a diagnostic estimate, not proof that the true type I error is exactly 5.40%. If the outcomes were correlated markers from one scan rather than independent simulation units, this binomial calculation would not give a valid uncertainty assessment without accounting for that dependence.

A single threshold can miss important features. Checking several thresholds, such as 0.01, 0.05, and 0.10, reveals whether excess rejection is concentrated in the extreme tail or extends across the distribution. The thresholds should be chosen before examining results where possible; selecting only the threshold that looks most concerning exaggerates apparent evidence of failure. Very stringent thresholds also require many null simulations: with \(R=10{,}000\) and \(a=10^{-5}\), the expected number of rejections is only 0.1, so that simulation is too small to estimate the tail rate precisely.

Worked Example: Checking P-Value Bins in Null Simulations

In a second hypothetical assessment, 10,000 null p-values are divided into ten equal-width bins: \([0,0.1),[0.1,0.2),\ldots,[0.9,1]\). The counts are 970, 1,030, 980, 1,020, 1,005, 995, 1,010, 990, 1,020, and 980. Under a uniform null distribution, each bin has expected count \(10{,}000/10=1{,}000\).

For an approximate goodness-of-fit check, calculate the Pearson statistic by summing squared deviations divided by expected counts:

$$ X^2 =\frac{(-30)^2+(30)^2+(-20)^2+(20)^2+5^2+(-5)^2+10^2+(-10)^2+20^2+(-20)^2}{1000} =\frac{3650}{1000} =3.65. $$

With ten bins and no parameters estimated from these binned counts, the reference degrees of freedom are \(10-1=9\). A value of 3.65 is not unusually large for a chi-square distribution with 9 degrees of freedom, so these counts show no obvious broad departure from uniformity. This check is approximate and assumes independent p-values; the bin counts can be dependent if they come from correlated markers. It also does not rule out a localized tail problem that a threshold-specific check might detect.

Genomic-Control Summaries

A compact summary often reported alongside a p-value scan is based on the median of one-degree-of-freedom test statistics. If \(T_j\) are such statistics, the genomic-control ratio is the observed median divided by the theoretical median of a \(\chi^2_1\) null statistic, approximately 0.4549:

$$ \lambda_{\mathrm{GC}} =\frac{\operatorname{median}(T_j)}{0.4549}. $$

A value near 1 indicates that the median statistic is near its one-degree-of-freedom null reference. This is a summary, not a full calibration analysis: different departures in the tails can have similar medians, and a median near its reference does not guarantee correct p-values at stringent significance thresholds. The ratio is intended for comparable one-degree-of-freedom statistics; statistics with different degrees of freedom require a suitable null reference. Detailed interpretation of this summary is the subject of the following material on the genomic inflation factor.

Worked Example: Interpreting a Median-Based Summary

Suppose a hypothetical scan reports a median one-degree-of-freedom statistic of 0.500. Its genomic-control ratio is

$$ \lambda_{\mathrm{GC}} =\frac{0.500}{0.4549} \approx 1.099. $$

The median statistic is about 9.9% above the theoretical null median. This is a reason to inspect calibration further, not a diagnosis of its cause. The difference could reflect model misspecification, numerical or testing issues, or genuine non-null signal affecting the distribution. A median summary alone cannot distinguish these explanations and cannot show whether the extreme tail is calibrated. In particular, it is not a direct estimate of the false-positive rate at a chosen p-value threshold.

Diagnostics, Alternatives, and Failure Modes

A useful calibration assessment combines evidence rather than relying on a single number. Null simulations test the full implemented procedure under a chosen null model. Threshold-specific rejection rates focus on type I error at inferentially relevant cutoffs. P-value bins or an empirical cumulative distribution compare broader features with the uniform reference. A genomic-control median provides a compact central summary. These approaches answer related but non-identical questions, and agreement among them is more informative than any one result.

ApproachWhat it assessesImportant limitation
Null simulationBehavior of the full procedure under the simulated nullConclusions depend on the realism of the simulated model
Threshold-specific ratesEmpirical type I error at selected cutoffsRare thresholds require many independent null outcomes
Bins or empirical distributionBroad departures from the uniform null referenceDependence complicates formal goodness-of-fit calculations
Median-based genomic controlCentral shift in one-degree-of-freedom test statisticsDoes not locate tail errors or identify their cause

Simulation design is a frequent source of misleading reassurance. Simulating independent residuals when the analyzed cohort has related individuals does not test whether the relationship adjustment is working. Holding nuisance estimates fixed in every replicate may also fail to reflect uncertainty from refitting the null model. A stronger check reproduces the analysis pipeline, including variance-component estimation and test-statistic computation, and records numerical failures as well as p-values.

Permutation is an alternative way to construct a null reference, but unrestricted phenotype permutation is generally not justified when observations are related or have unequal structure: it can break the dependence represented by \(K\). A permutation scheme is valid only when its exchangeability assumptions match the design and null model. Residual or restricted permutation methods have their own conditions and should not be treated as automatic fixes. Likewise, agreement with a fitted-model simulation cannot validate assumptions that the simulation itself imposes.

The observed genome-wide distribution also has a key limitation as a calibration diagnostic: real associations can create an excess of small p-values even when null tests are valid. Conversely, a scan with few true associations may look close to uniform while a particular subset of tests is miscalibrated. Calibration is most directly assessed from known-null tests or carefully designed null simulations; the full scan supplies context, not a definitive reference.

Common Pitfall

Wrong reasoning: “There are more small p-values than expected under a uniform distribution, so the mixed model is under-correcting for relatedness.” This conclusion treats every tested variant as null and assumes that any excess is an error. Genuine associations, polygenic signal, and model or computational problems can all contribute to the observed pattern.

Correction: State which tests are known to be null or define the null data-generating simulation before interpreting an excess as miscalibration. Examine error rates across thresholds, the dependence structure, and the analysis pipeline. A diagnostic can indicate that results warrant investigation; it usually cannot identify the cause on its own.

Key takeaway: Calibration means that null p-values meet their advertised error rates under the model and procedure being evaluated. Use null simulations and complementary distribution summaries, account for dependence, and do not equate a surplus of small p-values in a real scan with proof of miscalibration.

Check Your Understanding

Answer each question, then compare with the solutions.

  1. Under an exact continuous null test, what is the probability that a p-value is at most 0.02?
  2. A null simulation has \(R=2{,}000\) independent test outcomes. At threshold 0.05, what is the expected rejection count and the standard deviation of that count?
  3. Why might 10,000 correlated markers not provide the same precision for an empirical rejection rate as 10,000 independent null outcomes?
  4. A real scan has many very small p-values. Does this alone demonstrate failed calibration? Explain.
  5. For one-degree-of-freedom statistics with observed median 0.400, calculate the genomic-control ratio using 0.4549 as the null median and give a cautious interpretation.

Solutions. (1) The probability is 0.02 by null calibration. (2) The expected count is \(Ra=2{,}000(0.05)=100\), and its standard deviation is \(\sqrt{Ra(1-a)}=\sqrt{2{,}000(0.05)(0.95)}=\sqrt{95}\approx9.747\). (3) Correlation reduces the effective amount of independent information, so an independence-based binomial standard error can be too small; use independent simulation replicates or account for dependence. (4) No. True associations can generate small p-values, so assess calibration using known-null tests or an appropriate null simulation. (5) The ratio is \(0.400/0.4549\approx0.8793\). The median is below its null reference, which merits further checking but does not establish tail calibration or identify a cause.

Clinical Biostats