Calibrating QTL Tests With Parametric Bootstrap

Why Calibrate a QTL Test by Simulation?

A family QTL test asks whether a genetic locus explains trait covariance among relatives beyond the covariance allowed by a no-QTL model. A common test statistic is a likelihood ratio. Its usual large-sample reference distribution can be unreliable when the QTL variance is constrained to be nonnegative, when nuisance parameters are estimated, or when the pedigree and sample are not large enough for the asymptotic approximation to be accurate. The parametric bootstrap offers a model-based alternative: simulate data under the fitted no-QTL model, repeat the test on each simulated data set, and compare the observed statistic with the resulting empirical null distribution.

This calibration is conditional on the assumed model and the fitted null parameters. It does not make the test assumption-free, and it does not establish that a QTL is causal. It answers a narrower question: if the fitted no-QTL model generated family phenotypes with the observed family and marker structure, how often would the specified testing procedure produce a likelihood ratio at least as large as the one observed?

Test Statistic and Bootstrap Target

Definition: Let $\ell_0$ be the maximized log likelihood under the no-QTL model, and let $\ell_1$ be the maximized log likelihood under the alternative that includes a nonnegative QTL variance component. The likelihood-ratio statistic is $$ T_{\mathrm{obs}}=2(\ell_1-\ell_0). $$ Both models must use the same data and the same specification for their shared parameters. The alternative must include the null as a special case, so the maximized alternative likelihood cannot be smaller than the maximized null likelihood (apart from numerical optimization error).

For example, a null covariance model for a family phenotype vector $Y$ might be $V_0$, containing the non-QTL covariance components. At a specified locus, the alternative adds a nonnegative QTL component:

V_1(Z)=V_0+\tau^2 C(Z), \qquad \tau^2\geq 0.

Here $Z$ denotes the latent identity-by-descent (IBD) configuration, and $C(Z)$ is the QTL covariance contribution associated with that configuration. The mean specification and the other covariance terms are shared by the two models. If marker data $M$ leave IBD uncertain, the alternative likelihood integrates over configurations using their probabilities given the marker data:

L_1(Y\mid M)=\sum_z P(Z=z\mid M)\, \phi\!\left(Y;\mu,V_0+\tau^2C(z)\right).

The bootstrap should reproduce the likelihood actually used for the observed test. Thus, in a conditional bootstrap, keep each family’s pedigree structure and observed marker information fixed, including its IBD probabilities, and use the same integrated likelihood when refitting the alternative. Under the null, $\tau^2=0$, so the phenotype-generating covariance is $V_0$; the alternative is still refitted with $\tau^2$ constrained to be nonnegative.

Parametric Bootstrap Procedure

Fit the null model.
Estimate the parameters of the no-QTL model from the observed family phenotypes. These fitted values, together with the specified pedigree and marker structure, define the bootstrap-generating model.

Generate a replicate under the null.
For every family, simulate a phenotype vector from the fitted no-QTL distribution. In a conditional analysis, retain the observed family structures and marker information. Do not generate phenotypes using a positive QTL variance component.

Refit both models.
For the simulated data, maximize the null likelihood and the alternative likelihood anew, including re-estimation of shared nuisance parameters. Apply the same constraints, IBD integration, and optimization settings as for the observed data.

Save the replicate statistic.
Calculate $T_b=2(\ell_{1,b}-\ell_{0,b})$ for replicate $b$. Repeat the simulation and refitting $B$ times to obtain the empirical null distribution.

Compare with the observed statistic.
Count the replicates with $T_b\geq T_{\mathrm{obs}}$ and use the empirical tail proportion, with a finite-simulation correction, as the bootstrap p-value.

Refitting matters. The likelihood-ratio statistic reflects not only random phenotypes but also the ability of both models to estimate their parameters from a finite sample. Holding the observed parameter estimates fixed in every replicate would not reproduce that part of the test procedure. The statistic also depends on constraints: because $\tau^2$ cannot be negative, the alternative fit can land at $\widehat{\tau}^2=0$, producing a statistic of zero.

Worked Example: Calculating an Empirical Bootstrap P-Value

Suppose a hypothetical family study gives $\ell_0=-42.60$ and $\ell_1=-39.85$. The alternative has the same mean and non-QTL covariance specification as the null, and its QTL variance is constrained to be nonnegative. Calculate the observed statistic.

Observed statistic. $$ T_{\mathrm{obs}}=2\bigl[-39.85-(-42.60)\bigr] =2(2.75)=5.50. $$

A parametric bootstrap produces $B=999$ replicate statistics, of which $b=23$ are at least 5.50. The common finite-simulation estimate is

$$ \widehat{p}_{\mathrm{boot}}=\frac{b+1}{B+1} =\frac{23+1}{999+1}=0.0240. $$

The plus-one correction avoids reporting a simulated p-value of zero when no replicate exceeds the observed statistic. In this hypothetical analysis, the estimated tail probability is 0.0240: under the fitted null model and the specified testing procedure, about 2.4% of replicate statistics are at least as large as 5.50. This is not the probability that the null hypothesis is true.

The approximate Monte Carlo standard error of this estimate is $$ \sqrt{\frac{0.0240(1-0.0240)}{1000}}=0.00484. $$ Thus, the reported tail estimate has noticeable simulation uncertainty. More replicates would be useful if a more precise estimate were needed.

What the Simulation Preserves

The bootstrap-generating model must correspond to the analysis being calibrated. In a conditional analysis, families are not resampled as if they were independent unrelated individuals. Instead, each observed family’s phenotype vector is generated using its fitted null covariance, which retains the dependence implied by that family’s pedigree. The observed marker information and its IBD probabilities are also retained when the target is a test conditional on those data.

This distinction is especially important with uncertain IBD. Under the null, the phenotype simulation uses the fitted no-QTL covariance; the marker-derived IBD probabilities do not create a QTL covariance in the simulated phenotypes. They remain relevant when each replicate’s alternative likelihood is fitted, because that fit must integrate over IBD uncertainty in the same way as the observed analysis.

A different target may call for a different bootstrap. If the intended inference averages over new marker data or new pedigree configurations, those features must also be generated from an explicit model, and the analysis must recompute the corresponding IBD probabilities. That is not the same as a conditional bootstrap. The target, what is held fixed, and what is simulated should be stated before results are interpreted.

Worked Example: Simulating a Two-Relative Family Under the Null

Consider a hypothetical pair of relatives with a fitted null mean of 10 trait units and null covariance $$ V_0= \begin{pmatrix} 6&2\\ 2&6 \end{pmatrix}. $$ The diagonal variances are 6 and the covariance is 2. One way to generate a vector with this covariance is to write each phenotype as the mean plus a shared component, a relative-specific component, and residual variation: $$ Y_i=10+U+G_i+E_i, $$ where $U,G_1,G_2,E_1,E_2$ are independent, $U\sim N(0,2)$, $G_i\sim N(0,2)$, and $E_i\sim N(0,2)$. Each phenotype has variance $2+2+2=6$, and the two phenotypes share only $U$, so their covariance is 2.

For one illustrative draw, take $U=1$, $G_1=-0.5$, $G_2=0.5$, $E_1=0.25$, and $E_2=-0.25$. Then $$ Y_1=10+1-0.5+0.25=10.75,\qquad Y_2=10+1+0.5-0.25=11.25. $$ This is one possible null replicate, not a claim that these particular values are typical. In a bootstrap, the component values are drawn repeatedly from their fitted distributions. For every generated family data set, the null and QTL models are both refitted; the observed family’s marker information is retained for the alternative IBD-integrated likelihood.

Why the Empirical Reference Distribution Helps

When a variance component is tested at the boundary $\tau^2=0$, the usual regularity conditions for a standard chi-squared likelihood-ratio reference distribution do not all hold. In some simple settings a boundary mixture distribution is available, but its adequacy depends on the model and the testing setup. Small or irregular family samples, nuisance estimation, uncertain IBD, and a search across loci can further affect calibration. A parametric bootstrap can reproduce these features by generating data under the fitted null and running the same constrained fitting procedure on each replicate.

The justification for the empirical tail estimate is direct. Let $I_b$ equal 1 when bootstrap statistic $T_b$ is at least $T_{\mathrm{obs}}$, and 0 otherwise. Conditional on the fitted generating model, these indicators are Bernoulli variables with an exceedance probability $q$. The average of $B$ such indicators estimates $q$, and its variance is $q(1-q)/B$. The plus-one version is a finite-simulation convention that gives a nonzero estimate even when no simulated replicate exceeds the observed statistic.

\widehat{q}=\frac{1}{B}\sum_{b=1}^{B}I_b, \qquad \operatorname{SE}(\widehat{q})\approx \sqrt{\frac{\widehat{q}(1-\widehat{q})}{B}}.

This calculation describes Monte Carlo error only. It does not account for uncertainty in the fitted null parameters or for misspecification of the family covariance model. The bootstrap distribution is therefore an approximation under a fitted model, not an exact finite-sample distribution for every possible null parameter value.

Calibrating a Scan Rather Than One Locus

If the observed analysis tests one prespecified locus, each bootstrap replicate should reproduce that single-locus test. If the observed analysis searches across positions and reports the largest likelihood ratio, the replicate statistic must also be the largest across the same search. Otherwise, the bootstrap omits the extra opportunity for a large statistic created by searching.

Worked Example: Reproducing a Locus Scan

Suppose a hypothetical analysis evaluates 50 positions and reports the maximum likelihood-ratio statistic, $T_{\max,\mathrm{obs}}=7.20$. In each of 999 bootstrap replicates, the analyst simulates phenotypes under the fitted no-QTL model, refits the models at all 50 positions, and saves that replicate’s maximum. If 37 replicate maxima are at least 7.20, the scan-calibrated estimate is $$ \widehat{p}_{\mathrm{scan}}=\frac{37+1}{999+1}=0.0380. $$

This estimates the tail probability for obtaining a maximum of at least 7.20 under the fitted null and the specified 50-position search. It does not estimate the tail probability for a single prespecified position. The two procedures answer different questions, so their empirical distributions need not give the same p-value.

Common Pitfalls and Practical Checks

Simulating from the alternative. Generating phenotypes with a positive fitted QTL variance does not produce a null reference distribution. Generate from the fitted no-QTL model, with the QTL variance set to zero.
Failing to refit. Reusing the original fitted parameters in simulated data does not reproduce the likelihood-ratio procedure. Re-estimate both models for every replicate.
Changing the analysis between observed and simulated data. Use the same constraints, likelihood integration over uncertain IBD, locus grid, and optimization rules. If the observed statistic comes from a scan, save the scan maximum in every replicate.
Calling the bootstrap p-value an exact probability. It is an empirical estimate under a fitted parametric null. It has Monte Carlo error and can be inaccurate if the null model is misspecified.
Ignoring optimization failures. Check convergence and boundary estimates in both observed and replicate fits. A failed fit must not silently be treated as an ordinary statistic; investigate it and apply a consistent, defensible fitting rule.

A useful report identifies the null model and its fitted parameters, what was held fixed and what was simulated, how uncertain IBD was handled, the number of replicates, the exceedance count, and the empirical p-value. These details make clear which null distribution was approximated and whether the calibration matches the actual QTL testing procedure.

Key takeaway: A parametric bootstrap calibrates a QTL likelihood-ratio test by simulating phenotypes from the fitted no-QTL model and refitting the null and alternative models in every replicate. Preserve the analysis structure, including IBD integration and any locus search, and interpret the resulting tail estimate as conditional on the fitted model.

Check Your Understanding

Answer these questions using the definitions and procedures in this tutorial.

Why must a bootstrap replicate refit both the null and alternative models rather than reuse the estimates from the observed data?
In a conditional bootstrap with uncertain IBD, which information is retained, and how does IBD uncertainty enter the alternative fit?
For $B=499$ replicates, 14 of which have statistics at least as large as the observed statistic, calculate the plus-one empirical p-value.
Why should each bootstrap replicate save the maximum statistic when the observed analysis reports the maximum over a locus scan?
What does a parametric bootstrap p-value fail to account for when the fitted null model is misspecified?