Calibrating the SKAT-O Parameter Search

Why the SKAT-O Search Needs Calibration

SKAT-O combines burden-like and variance-component evidence through a correlation parameter, often denoted by \(\rho\). A small \(\rho\) emphasizes the SKAT component; a value near one emphasizes the burden component. In practice, the analysis may evaluate several \(\rho\) values and report the most significant result. That search can make the reported p-value smaller than a p-value from any one prespecified test, even when the null hypothesis is true.

The relevant null distribution is therefore not just the distribution of one SKAT-O statistic. It is the distribution of the smallest marginal p-value across the searched values. The tests are generally dependent because they are constructed from the same genotype data and score vector. Correct calibration must account for that dependence rather than treating the marginal p-values as independent by default.

Definition: For a prespecified grid \(\rho_1,\ldots,\rho_m\), let \(p_{\rho_j}\) be the marginal null p-value for the SKAT-O statistic at \(\rho_j\). The search statistic is \(T_{\min}=\min_{1\le j\le m}p_{\rho_j}\). Its search-adjusted p-value is \(P_0(T_{\min}\le T_{\min,\mathrm{obs}})\), where \(P_0\) denotes probability under the null model. It measures how often a null analysis would produce a minimum marginal p-value at least as small as the observed minimum.

Null Model and Dependence Across the Grid

The calibration begins with the null score distribution used for SKAT. Suppose the covariate-adjusted score vector for the variants is \(U\), and under the null it is approximated by a mean-zero multivariate normal vector with covariance matrix \(V\): \(U\sim N(0,V)\). This is the usual score-test null approximation, with the nuisance parameters and covariate adjustment handled as in the fitted SKAT model. For each fixed \(\rho\), write the quadratic statistic as \(Q_\rho=U^\mathsf{T}A_\rho U\), where \(A_\rho\) is the symmetric matrix defining that SKAT-O statistic, including the variant weights.

Under this model, the marginal null distribution of \(Q_\rho\) is a weighted sum of independent one-degree-of-freedom chi-square variables. Its weights are the eigenvalues of \(V^{1/2}A_\rho V^{1/2}\), so \(Q_\rho\) has the distribution \(\sum_k\lambda_{\rho k}\chi^2_{1,k}\). The marginal p-value \(p_\rho\) is the upper-tail probability from that distribution, evaluated at the observed \(Q_\rho\). This is the mixture-chi-square calibration used for SKAT statistics.

Although each statistic has its own marginal distribution, all \(Q_{\rho_j}\) use the same \(U\). In particular, under the stated normal model,

\operatorname{Cov}_0(Q_{\rho_j},Q_{\rho_k}) =2\operatorname{tr}(A_{\rho_j}V A_{\rho_k}V).

The shared score vector induces dependence even when the matrices differ. Consequently, the probability of the minimum cannot generally be obtained by multiplying marginal probabilities. A null calibration must use the joint behavior of the statistics, or a valid bound that does not require knowing that joint behavior exactly.

Assumptions for the calibration: The grid of \(\rho\) values and the statistic at each value are specified before examining the association results; the null score distribution and covariance \(V\) adequately represent the data under the null; and each marginal p-value is computed from the corresponding null distribution. Simulation from \(U\sim N(0,V)\) gives a model-based calibration under these assumptions. If the null approximation or covariance model is inadequate, the resulting calibration may also be inadequate.

Deriving the Search-Adjusted Probability

For a threshold \(t\), the event that the minimum p-value is no greater than \(t\) occurs if at least one grid-specific p-value is no greater than \(t\). Equivalently, its complement is the event that every p-value exceeds \(t\). Thus the null cumulative distribution of \(T_{\min}\) is

P_0(T_{\min}\le t) =1-P_0(p_{\rho_1}>t,\ldots,p_{\rho_m}>t).

This identity is the basic calibration result. It does not assume independence. The joint probability on the right must be calculated from the actual dependence structure or estimated by a null simulation that preserves it.

For comparison, if all \(m\) marginal p-values were independent and exactly uniform under the null, then \(P_0(p_{\rho_j}>t)=1-t\) for each \(j\), and independence would give

P_0(p_{\rho_1}>t,\ldots,p_{\rho_m}>t)=(1-t)^m, \qquad P_0(T_{\min}\le t)=1-(1-t)^m.

SKAT-O p-values across \(\rho\) are not generally independent, so this formula is a reference calculation, not a default calibration. Positive dependence often makes the independence calculation too large, but dependence can be more complicated; one should not assume a direction without checking the joint null distribution.

Worked Example: The Search Event and the Independence Calculation

Suppose a hypothetical analysis searches five fixed parameter values and obtains \(T_{\min,\mathrm{obs}}=0.012\). First consider the simplifying assumption that the five marginal p-values are independent and uniform under the null. The probability that all five exceed 0.012 is

P_0(p_{\rho_1}>0.012,\ldots,p_{\rho_5}>0.012) =(1-0.012)^5=0.988^5\approx 0.9414.

Therefore, the probability that at least one is at most 0.012, and hence the independence-based search-adjusted p-value, is

P_0(T_{\min}\le 0.012) =1-0.9414=0.0586.

The unadjusted minimum is 0.012, while the independence-based adjusted value is about 0.0586. This calculation is correct only under the stated independence and uniformity assumptions. Because SKAT-O statistics across \(\rho\) share the same score vector, those assumptions usually do not describe the actual search.

Joint Null Simulation

A direct way to account for dependence is to generate null score vectors with the fitted null covariance structure, calculate every statistic on the grid for each vector, and record the minimum marginal p-value. For a positive semidefinite covariance matrix \(V\), one can write \(U=V^{1/2}g\), where the entries of \(g\) are independent standard normal variables. Each simulated vector then has covariance \(V\), and the same vector is used across all \(\rho\) values within a replicate. That shared vector is essential: simulating each grid point independently would destroy the dependence being calibrated.

For each replicate, compute \(Q_{\rho_j}=U^\mathsf{T}A_{\rho_j}U\) for every grid value and convert each statistic to its marginal p-value using that value’s null mixture distribution. Record the smallest of those p-values. Comparing the observed minimum with this simulated null distribution estimates the search-adjusted p-value.

Fix the search.
Specify the \(\rho\) grid and the statistic to be evaluated at each grid point. The search represented in the calibration must match the search used for the reported result.

Generate a null score vector.
For replicate \(b\), draw \(g^{(b)}\sim N(0,I)\) and set \(U^{(b)}=V^{1/2}g^{(b)}\), or use an equivalent null-generation method that preserves the fitted covariance structure.

Evaluate the whole grid together.
Use the same \(U^{(b)}\) to compute all \(Q_{\rho_j}^{(b)}\), convert them to their marginal p-values, and retain \(T_{\min}^{(b)}=\min_j p_{\rho_j}^{(b)}\).

Compare with the observed minimum.
Estimate the adjusted p-value by the fraction of null replicates whose minimum is no greater than \(T_{\min,\mathrm{obs}}\), with a finite-simulation correction when appropriate.

Worked Example: Calibrating by a Simulated Joint Null

Suppose a hypothetical null simulation uses \(B=100{,}000\) score vectors, each generated from the same fitted null covariance matrix and evaluated across the entire prespecified SKAT-O grid. The observed minimum marginal p-value is \(0.0040\), and 317 simulated minima are less than or equal to \(0.0040\). Using the plus-one Monte Carlo estimate gives

\widehat{p}_{\mathrm{adj}} =\frac{1+\#\{b:T_{\min}^{(b)}\le 0.0040\}}{B+1} =\frac{1+317}{100{,}000+1} =\frac{318}{100{,}001} \approx 0.003180.

Under the fitted null model, the calibrated probability of obtaining a grid minimum at least as small as the observed 0.0040 is estimated to be about 0.00318. This is a search-adjusted p-value, not a probability that the null hypothesis is true. The plus-one estimate avoids reporting a p-value of zero when no simulated replicate crosses the threshold.

The estimate also has Monte Carlo error. Treating the exceedance proportion as a binomial estimate, its approximate simulation standard error is

\sqrt{\frac{\widehat{p}_{\mathrm{adj}}(1-\widehat{p}_{\mathrm{adj}})}{B+1}} \approx \sqrt{\frac{0.003180(1-0.003180)}{100{,}001}} \approx 0.000178.

This is uncertainty from the finite number of null replicates, not uncertainty in the underlying association effect. More replicates are needed when a very small adjusted p-value must be estimated precisely.

Alternative Bounds and What They Do Not Use

When joint-null simulation or an analytic joint calculation is unavailable, the Bonferroni inequality gives a conservative bound. If each of the \(m\) marginal p-values is valid under the null, then for any \(t\),

P_0(T_{\min}\le t) =P_0\left(\bigcup_{j=1}^m\{p_{\rho_j}\le t\}\right) \le \sum_{j=1}^mP_0(p_{\rho_j}\le t) \le mt.

Consequently, \(\min(1,mT_{\min,\mathrm{obs}})\) is a valid, possibly conservative, search-adjusted p-value for the finite grid. This bound does not require independence or a model for the dependence, but it may be less precise than calibration using the correct joint null distribution. The effective cost of the search can be smaller than \(m\) when the tests are strongly dependent, but it should be estimated from the joint null rather than guessed.

Worked Example: A Dependence-Free Bonferroni Bound

For a hypothetical search of 12 fixed \(\rho\) values, suppose the smallest marginal p-value is \(0.0065\). The Bonferroni-adjusted value is

\min(1,12\times 0.0065)=0.0780.

Thus, without making a dependence assumption, the Bonferroni bound is 0.0780. It is a valid but potentially conservative adjustment if all 12 marginal p-values are valid. It is not an estimate of the exact joint-null probability unless the bound happens to be tight.

Interpretation and Common Pitfalls

The calibrated p-value refers to the entire parameter search that was included in the null calculation. If the analysis also selects among uncalibrated choices—for example, additional grids, variant sets, or testing procedures—the stated p-value does not automatically account for those choices. The calibration must reproduce the selection process relevant to the result being reported. Likewise, if a grid is changed after examining results, the original calibration no longer necessarily matches the search that produced the reported minimum.

A common mistake is to say that the chance that all \(m\) p-values exceed \(t\) is \(1-(1-t)^m\), then use that as the probability of obtaining at least one small p-value. Under independence and uniform marginal p-values, \((1-t)^m\) is the probability of the intersection that all p-values exceed \(t\). The probability of at least one p-value at or below \(t\) is its complement, \(1-(1-t)^m\). For SKAT-O, even this complement applies only under independence; joint simulation or another valid dependence-aware method is needed for the actual calibration.

Key takeaway: A minimum over SKAT-O correlation values is a selected statistic. Calibrate it by evaluating the same full search under the joint null distribution, preserving dependence across grid points. Report the calibrated probability for the observed minimum, and state whether it came from joint-null calculation, simulation, or a conservative bound.

Check Your Understanding

Answer these questions using the definitions and calibration logic above.

For a grid of \(m\) values, express \(P_0(T_{\min}\le t)\) using the probability that all marginal p-values exceed \(t\).
Under independence and uniform marginal p-values, what is the probability that all \(m\) p-values exceed \(t\), and what is the probability that at least one is at most \(t\)?
Why must a null simulation use the same score vector across all \(\rho\) values within a replicate?
A null simulation with \(B=20{,}000\) replicates has 39 minima no greater than the observed minimum. Calculate the plus-one Monte Carlo estimate.
What search choices are not accounted for if the calibration includes only a fixed grid of \(\rho\) values?