Two-stage analysis strategy for identifying the IgM quantitative trait locus

Association studies currently offer an exciting approach to mapping complex quantitative trait loci (QTLs). Wallace et al. [1] recently recommended a generalized Hotelling's T² test for QTL linkage disequilibrium (LD) mapping, which is uniformly the best test for normally distributed data. However, if the assumption of a normal distribution is not met, T² may lose considerable power. When the trait distribution is unclear, some nonparametric tests may be preferred because they are only slightly less powerful than T² when the trait is normally distributed, but much more powerful than T² in some cases of non-normality. In general, it is unclear what test is the best when the trait distribution is unknown. Some investigators only report the most significant result from several statistics, but the type I error rate cannot be properly controlled when this is done. It is also not wise in this situation to use an approach such as the Bonfferoni method to control the type I error rate because the various tests are usually highly correlated and therefore the result will be overly conservative.

To achieve satisfactory power in most cases, without compromising the overall significance level, we consider adopting a two-stage adaptive analysis strategy: several statistics are compared on a portion of the samples in the first stage and the statistic that is found to be most powerful is then used for the remaining samples. Previously, two-stage strategies have been adopted in genetic association studies to reduce the cost of genotyping [2, 3] or the penalty due to multiple testing when modeling gene × gene interactions [4]. Here, we apply this strategy for a different purpose: to select a powerful test for the data at hand and hence obtain good power overall. We evaluate this procedure of adaptively selecting the optimal test by mapping the IgM QTL with the simulated data of Problem 3 in Genetic Analysis Workshop 15 (GAW15).

The procedure examines the power of various statistics using a portion of the data in an exploratory first stage and then applies this most powerful test to the rest of the data in the second stage. The statistics from the two stages are combined to make full use of the information. This approach of combining the results of the two stages is equivalent to a more general method of combining p-values. For the procedure of combining these p-values to be valid, however, we need to specify before the analysis which statistic will be used to obtain the p-value (p₁) from the exploratory stage in the combination. The p-value from the second stage (p₂) is calculated based on the most powerful statistic found at the first stage. Under the null hypothesis, each p value is, at least asymptotically, distributed uniformly on U(0, 1). The final decision then depends on a combining function f(p₁, p₂). The most common such function may be Fisher's combination test [5], which is defined by

f(p₁, p₂) = -2log(p₁p₂),

where under the null hypothesis Fisher's statistic will be distributed as a χ² with 4 degrees of freedom. Another example is the weighted inverse normal method,

f(p₁, p₂) = 1 - Φ[w₁Φ^-1(1 - p₁) + w₂Φ^-1(1 - p₂)],

where Φ is the cumulative distribution function of a standard normal distribution, 0 <w _i < 1 and w₁² + w₂² = 1. This statistic will be distributed as a standard normal distribution.

To obtain p₂, we have to estimate the power of the various statistics at the exploratory first stage. Traditional power calculation methods require the trait distribution to be known, which is not the case here. A bootstrap method of using the data from the exploratory stage can be adopted to approximate the power [6, 7]. The bootstrap and permutation are two often used nonparametric procedures. It is often desired to obtain "exact" p-values by employing a permutation procedure to generate the null distribution of the statistic that is used for a test. Here, on the other hand, we want to estimate the power of a statistic, and for this we need the distribution of the statistics under the alternative hypothesis; a permutation procedure cannot be directly applied for this purpose. Let the trait values of individuals with genotype g be denoted x _g , where g = 0, 1, 2 for an additive SNP marker and g = 0, 1 for a recessive/dominant marker. For this example, we assume a dominant model for the rarer allele. We denote the sample size for each genotype n _g . We assume the distribution of trait values for different genotypes have similar shape, but the locations of the distributions are shifted by d _g . The hypothesis to detect association between a marker and the trait is then defined as H₀:d = 0. The power function of the statistic T for d = δ at the significance level α is then given by P(T; δ, α) The method of Collings and Hamilton [6] to approximate P(T; δ, α) by a nonparametric bootstrap procedure is as follows:

1. For each genotype group g, a random sample of $n={\displaystyle {\sum }_g{n}_g}$, n₂ <n₀, trait values is drawn with replacement. The sampled trait values are denoted X _g ^b = (x_1g ^b ,...,x _ng ^b ). A simulation sample of trait values, Y _g ^b , is then obtained by adding X _g ^b to (0, δ), where 0 is a row vector of n₀ elements each of which is 0 and δis row vector of n₁ + n₂ elements, each of which is δ. The corresponding genotype groups are set to be G ^b = (0, 1).

2. Different statistics are calculated on the simulated sample values Y _g ^b and G ^b , and the corresponding p-values (p _g ^b ) are recorded.

3. Steps 1 and 2 are repeated B times. The estimated power function of ${\widehat{P}}_g(T;\delta ,\alpha )$ is given by $\frac{{\displaystyle \sum I_{\{p_g^b<\alpha \}}}}{B}$.

4. Finally, we estimate the power of the different statistics using the weighted average estimates of the different genotype groups, given by $\frac{{\displaystyle {\sum }_g(n_g{\widehat{P}}_g)}}{n}$.

Two-stage designs have been applied to large-scale genetic association studies to substantially reduce genotyping cost while maintaining power. In addition to the knowledge of which markers are promising, we can obtain information about the distribution of the phenotype based on the data from the exploratory stage. This knowledge is useful for the choice of a statistic to use at the second stage and can therefore lead to a considerable gain in power. In our analysis, we evaluated this idea by considering just two statistics. Hotelling's T² has been proved to be a powerful statistic, even with sample selection. However, the advantage of T² depends on the trait distribution. On the other hand, although a nonparametric statistic is not the most powerful one when normality of the trait holds, it usually works well. So it is reasonable to consider combining the p-value of a nonparametric statistic from the exploratory stage with the p-value of the most powerful statistic for the second stage.

The idea of a two-stage analysis can be further generalized in genetic association studies. Because LD patterns vary greatly, it is often unclear whether a single-marker analysis or a multiple-marker analysis or a haplotype-based analysis is most powerful for a specific data set. Further work on developing a data-driven adaptive procedure to choose the type of analysis to perform on the second stage data would be potentially useful.

The adaptive two-stage procedure can lead to considerable gain in power by guiding the choice of a test based on the knowledge learned from an exploratory stage. At the same time, the type I error rate can be well controlled.