Calibrating the SKAT-O Minimum P-Value

Why the Smallest Component P-Value Needs Calibration

SKAT-O combines burden-like and variance-component evidence by evaluating tests over a set of values of \(\rho\). Each value gives a component statistic and a corresponding p-value. Choosing the smallest of these p-values is a useful way to adapt to different patterns of variant effects, but the selected minimum is not automatically a valid omnibus p-value. The same data, and therefore the same score information, contribute to every component test. Their p-values are generally dependent.

The distinction is important: a component p-value measures evidence for one specified test, while an omnibus p-value measures evidence after considering the entire set of component tests. Calibration must account for the chance that at least one component appears extreme under the null.

The grid is part of the testing procedure. It must be fixed in advance or otherwise included in the calibration. Changing the grid changes the number and dependence structure of the component tests, and can change the distribution of \(T\). The observed minimum alone does not contain enough information to determine its null probability.

How the Joint Null Distribution Arises

Under the null model, SKAT-O is based on a score vector \(U\) for the variants in the set. In the usual large-sample score approximation, \(U\) is modeled as multivariate normal with mean zero and covariance matrix \(V\). The component statistic at \(\rho\) is a quadratic form in this same score vector, which can be written generally as \(Q_\rho=U^\mathsf{T}A_\rho U\), where \(A_\rho\) is determined by the test definition, variant weights, and \(\rho\).

For each fixed \(\rho\), the null distribution of \(Q_\rho\) is a weighted sum of independent chi-square random variables, with weights determined by the eigenvalues of the corresponding covariance-adjusted quadratic form. This yields the marginal component p-value. But the vector \((Q_{\rho_1},\ldots,Q_{\rho_m})\) is formed from one shared \(U\), not from independently generated score vectors. Its coordinates, and therefore the component p-values, are correlated.

Joint calibration can be performed by drawing null score vectors from the fitted null distribution, calculating all component statistics from each same draw, and converting each statistic to its marginal p-value. The simulated minimum then approximates the null distribution of \(T\). In implementations using a more detailed null-generation method, the essential requirement is the same: each replicate must preserve the joint dependence across the entire rho grid.

Why Marginal P-Values Do Not Settle the Problem

For any threshold \(t\), the event that the minimum is at most \(t\) is the union of the component events:

Thus its probability depends on the intersections among these events. If the component p-values are exactly uniform under the null, then \(P_0(p_j\leq t)=t\) for each \(j\), but this marginal fact does not determine the probability of their union. If p-values are only valid, meaning \(P_0(p_j\leq t)\leq t\), each marginal probability is at most \(t\), not necessarily equal to \(t\). Continuity and validity alone do not imply exact uniformity.

A useful bound follows from the union bound:

when each component p-value is valid. This bound motivates a Bonferroni adjustment, but it can be conservative when component tests are strongly dependent. Exact joint calibration uses the dependence structure rather than discarding it.

Worked Calibration for a SKAT-O Grid

The following hypothetical analysis uses a fixed grid of five rho values. Assume the null model is correctly specified, the score vector has a suitable multivariate-normal null approximation, and the fitted covariance matrix captures its dependence. Genotypes and covariates are held fixed while null score vectors are generated. For each replicate, all five component tests use the same generated score vector.

Justification of the Joint Simulation

The calibration procedure follows directly from the definition of the target probability. Under the null, let \(U^{(b)}\) be a simulated score vector for replicate \(b\). For each grid point, calculate \(Q_{\rho_j}^{(b)}=(U^{(b)})^\mathsf{T}A_{\rho_j}U^{(b)}\), then obtain its marginal upper-tail p-value \(p_j^{(b)}\) using that component's null distribution. Define \(T^{(b)}=\min_j p_j^{(b)}\). The estimated fraction of replicates with \(T^{(b)}\leq T_{\mathrm{obs}}\) estimates the probability in the definition.

The dependence-preserving step is essential. If a separate independent score vector were generated for each rho value, the simulated minima would reflect an artificial independence assumption and could have the wrong null distribution.

Alternatives, Diagnostics, and Limitations

Bonferroni calibration controls the family-wise probability of at least one false positive by comparing the minimum to \(\alpha/m\), where \(m\) is the number of component tests. It is simple and valid under arbitrary dependence when the component p-values are valid, but it does not exploit their dependence and can be conservative. Joint calibration targets the minimum's actual null distribution under the specified model. Reporting the unadjusted minimum is not an alternative calibration: it ignores the selection step.

Check that the rho grid used in null replicates exactly matches the observed grid and that the score covariance and marginal component distributions are computed consistently. A useful simulation diagnostic is to compare empirical null rejection rates of the calibrated omnibus p-values at several thresholds with those thresholds, using enough replicates to distinguish discrepancies from Monte Carlo noise. If the null generator omits dependence, the model is misspecified, or component tail probabilities are inaccurate, the resulting omnibus p-value may also be miscalibrated.

For binary phenotypes with sparse data, a normal score approximation may be inadequate; the calibration must use a null procedure appropriate to that setting. Relatedness, population structure, covariate adjustment, and estimated nuisance parameters also affect the score covariance or null generation. Joint calibration does not repair an incorrect phenotype model, poorly estimated covariance, or a variant set chosen after inspecting association results unless that selection is included in the testing procedure.