Yet another additional challenge is the characterization of the nature of the relationship, which can be complicated and nonlinear. Recursive partitioning (RP) is a simple data mining tool that can, nevertheless, adequately handle these problems (Young et al., this issue). RP has been previously applied to genetic problems as a tool for dealing with the detection of interaction and identification of homogeneous subgroups (Rao 1998; Zhang et al., 2000; Czika et al., 2001; Province et al., 2001; Shannon et al., 2001; Costello et al., 2003). However, RP can be criticized for over-reliance on the presence of “marginal effects”.

When RP makes the first sweep through predictor variables (e.g., a set of genetic markers) it tries to identify those variables that result in a significant association with the trait. The data are then divided into subgroups corresponding to selected predictor categories, and the association test is repeated in the subgroups using the remaining predictors. At this stage the interaction can be detected. Suppose the first RP split divides the sample into two subgroups with individuals carrying the genotype AA vs. Aa + aa; the test among individuals in the AA group can detect interaction of AA with another marker. This interaction can only be detected if the first split is “significant” (i.e. there is substantial “marginal” effect at the marker A). One purpose of this article is to illustrate that it is reasonable to assume that most interactions will induce such marginal effects at one or more of the markers involved in the interacting set.

From the clinical perspective, a most useful characteristic of a predictive set of genetic polymorphisms is its ability to provide estimates of the probability of an event (condition) Y, such as an adverse event or disease. For a binary trait (event present/absent, denoted by Y, N), we also describe estimates of probabilities of developing the condition Y given that a random individual is classified into a particular node of the RP tree according to its multilocus genetic profile. Along with these estimates, interval estimates will also be provided. These are available by using a hold-out sample (i.e., part of the sample set aside for the validation purposes).

where F(p) is the p-value cumulative distribution function under the alternative hypothesis, which reflects the power. For example, with the χ² test, where Ψ_d,λ and Ψ_d,0 are the cumulative distribution functions of central and non-central χ², respectively, and d denotes the degrees of freedom.

Supposing that statistical tests and their p-values are used to rank the hypotheses and that there is a single “true” association, H_A, among L−1 null hypotheses, H₀, where L is the number of combinations, then the probability that a randomly picked hypothesis represents a false discovery is

If we let (FDR) = FDR/(1 − FDR), then this is an increasing function of FDR. Using Pr(H₀) = 1 − 1/L,

One interpretation of this formula is that the power gained by considering combinations must increase more sharply than the number of combinations, L, otherwise FDR is going to increase. If L is an exponential function of the number of markers, k, the condition that FDR should not increase is difficult to achieve. Practically, this means that when L is large, the chances that a claimed association turns out to be a false positive are high.

For example, assume the power to detect a single true combination at α = 0.05 is 80%. Our analysis depends on the definition of a “true combination” and for illustrative purposes we assume that there is only a single subset. One can use a different definition, whereby every subset that partially overlaps with the combination of interest does represent a true discovery, which will reduce L in half. If we start with k = 10 markers and examine all L = 2^k−1 = 1023 combinations, then among p-values ≤ 0.05, 98% are going to be false positives, with the assumption of a single true combination.

It is worthy to note that simple multiple testing correction alone does not easily get around this problem: suppose the testing is at the Bonferroni-corrected α = 0.05/1023 and that the proportion of false positives among p-values as small as 0.05/1023 ≈ 4.9 × 10⁻⁵ is examined. The proportion of false positives still remains high; about 32%. The situation worsens with a larger k. For example, when k = 15 the proportion of false positives among significant p-values after applying the Bonferroni correction (1.5 × 10⁻⁶ or smaller) is 69%.

This problem is very hard to circumvent, because L quickly increases with the number of markers, but the sample size does not increase correspondingly (see Witte et al., 2000, for requirements on that). On the other hand, RP has the implicit assumption that at least one of the individual markers in the “true” associated set, as well as the smaller subsets (combinations) of markers in this set, should show some degree of association. As a consequence, the number of combinations remains linear as a function of the number of markers, keeping FDR relatively low (e.g., a tree with three splits will make approximately 3 × k comparisons).

How reasonable is the assumption that at least one of the individual effects will be present? One might think it is easy to come up with examples of marker interactions that preclude “marginal” effects. Consider the simplest possible scenario: a binary case (e.g., the presence/absence of a disease) and the control phenotype (Y, N) and two diallelic loci (two binary predictors), x₁ = {0, 1} and x₂ = {0, 1}. If one considers a situation where allele ”1” has to be present at both markers in order for the disease risk to be high, then this is a clear case of interaction; however, the “case” category is likely to be enriched with the alleles of type ”1” at both markers compared to controls. A more subtle and seemingly “no-marginal effects” scheme is the following: ”1,1” or ”0,0” combinations are both associated with high risk, whereas ”0,1” and ”1,0” have low risk. Similar situations have been considered as justification for methods that go after “pure interactions”. Jannot et al. 2003 considered such interactions for the case of genotypic interactions where the penetrance array is 3 × 3 for two loci. They evaluated an approach wherein all 2^k−1 marker combinations were analyzed and p-values associated with the best combination were obtained by a permutation algorithm. Lin et al. 2004 took a similar approach using transmission disequilibrium tests. They state that because of the usage of the permutation approach to multiple tests “the gains from considering haplotypes in our exhaustive allelic method are not overshadowed by the penalty paid for doing far more tests”. These approaches assume that the accurate type-I error protection through permutations is balanced with the appropriate increase in power provided by considering the proper subset. Ritchie et al. 2001 proposed an algorithm that starts with a multi-dimensional table that allows for all possible interactions for a given subset of markers. The multi-dimensional table is represented by one of 2^k−1 combinations of k markers. The high-risk cells are then pooled into one group and the low-risk cells into another group, and a cross-validation is used to assess robustness of the reduced model. A combinatorial partitioning method (CPM) proposed by Nelson et al. 2001 evaluates all partitions of a cells (multilocus genotypes) into b groups. The number of partitions is given by:

which is substantially larger than the number of marker subsets. For example, with three diallelic markers there are 2³−1 = 7 marker subsets that could be tested separately. However, assuming we start with the binary classification at each locus and a = 2³ multilocus combinations, there are 127 ways to partition the genotypes into two groups, and 966 ways to partition them into three groups. Another combinatorial algorithm by Tahri-Daizadeh et al. 2003 considers increasingly complex models and evaluates proposed models on the basis of information criterion (Akaike, 1974). For quantitative traits, Culverhouse et al. 2004 proposed a modification of CPM that reduces this number of putative models by successfully merging cells with similar values of response. Such short-cut algorithms may consider a single model after the collapsing is finished. Significance levels are adjusted accordingly by a permutation algorithm. Nevertheless, the null distribution obtained via permutations essentially accounts for the number of looks through the data. This can be thought of as the adjustment by the ”effective” number of tests, L_e ≤ L.

While methods optimized for discovery of pure interactions are important, even the most contrived penetrance configurations allow for marginal effects – the feature exploited by the RP. The advantage of the RP method is in drastically reducing the number of potential models, thus lowering Pr(H₀) in (1) and, therefore, reducing the probability of a false discovery. Cordell and Clayton (2002) proposed a stepwise logistic regression procedure for evaluating the effects of the different polymorphisms within a gene. Curran (2003) compared such approach using the combination of forward and backward selection with the performance of RP. The research indicated a remarkable similarity between the two methods in terms of the predictive power and the effects and interactions entering the final model.

Returning to the “no-marginal effects” scheme (”1,1”/”0,0” vs. ”0,1”/”1,0”) it may seem that all alleles are going to be present at equal frequencies in both phenotype groups. However, this is generally not the case and it is possible to discover such haplotype-trait associations by looking at just a single marker. The single marker effect may be smaller than a multilocus effect; however the reduction in degrees of freedom or the number of tests may balance a reduction in the effect magnitude.

The four possible haplotypes formed by the two SNPs are X = {00, 01, 10, 11}, with frequencies P = {p_00, …, p₁₁}. Four haplotype penetrance values are:

With random pairing of haplotypes and multiplicativity, the population prevalence of Y is:

In terms of the two SNPs that constitute a haplotype, the values γ_ij define two-locus interactions with respect to the binary phenotype. For example, a high value of Pr(Y | X = 11) relative to γ indicates that both x₁ = 1 and x₂ = 1 are required for this probability to be above the population average. Suppose we have allele information at one of the markers, x₁. Then the association can be detected if the penetrance of one of the alleles, e.g. γ_1·. = Pr(Y | x₁ = 1) is different from γ. This marginal allele penetrance value is

using the allele probabilities among the case and the control groups:

In turn, these are calculated from the frequencies of haplotypes and the penetrance values as

Equation (4) can also be written as

When γ_1· ≠ γ, the marginal effect of a SNP can be detected. Thus, in general, this effect depends not only on the haplotype penetrances but on the population frequencies of the haplotypes as well. Looking at the “no marginal effect” difference,

it can be seen that for a plausible “pure interaction” penetrance configuration, , the additional requirement on the population frequencies of haplotypes is such that the two pairs of frequencies should match:

As an example let the haplotype penetrance array {γ₀₀, γ₀₁, γ₁₀, γ₁₁} be Γ = {0.9, 0.2, 0.2, 0.9}, which may seem to imply that an SNP effect is unlikely, because of the ”orthogonal” structure of Γ. As often happens in reality, one of the four haplotypes can be of relatively high population frequency. In such a case, the marginal (SNP) effects can be quite pronounced, as can be seen from two following examples.

Example 1. Consider a situation when one of the low penetrance haplotype has a high frequency, P = {0.05, 0.05, 0.85, 0.05}. In this case, γ = 0.27 and γ_1· = 0.55. The frequencies among the cases and the controls are Pr(x₁ = 1 | Y) = 0.204 and Pr(x₁ = 1 | N) = 0.062 corresponding to the log_e of relative risk of 1.191.

Example 2. Now suppose that one of the high penetrance haplotypes has a high frequency, P = {0.05, 0.05, 0.05, 0.85}. The values become γ = 0.83 and γ_1· = 0.55. The frequencies among the cases and the controls are Pr(x₁ = 1 | Y) = 0.066 and Pr(x₁ = 1 | N ) = 0.265 corresponding to the log relative risk of −1.39. Such value provides 99% power (at α = 0.05) with samples of 160 cases and 160 controls.

Allelic penetrance equation (5) can be easily generalized to multiple genetic loci. Let ”j” index all possible haplotypes that contain allele ”1” at the locus of interest. Then the marginal effect associated with allele ”1” is

where the denominator gives the frequency of ”1”. The condition of no marginal effect associated with allele ”1” is again γ_1· = γ.

Example 3. As a three-locus illustration using RP, consider three diallelic loci, x₁, x₂, x₃, with eight possible haplotype combinations, {000, 001, ..., 111}, with penetrances and population frequencies sampled from the Dirichlet(1/4, ..., 1/4) distribution, as given in Table 1.

In this example there are two high risk haplotypes (000 and 111) carrying different alleles. By looking at the penetrance values alone, a technique that relies on the marginal effects associated with one of the three SNPs may appear unlikely to be successful. For simplicity, we considered a haploid population and sampled 250 case and 250 control haplotypes. A sample RP tree using HelixTree^® is shown in Figure 1. HelixTree is a software system that uses recursive partitioning algorithms customized to the field of genetics (www.goldenhelix.com).

The first split is on x₃ with a p-value 4.6 ×10⁻⁴, signifying a substantial effect associated with this SNP. Allele ”0” at this SNP appears predictive of the cases, with the sample prevalence increasing from 0.5 in the root node to 0.86 at the first split. Following splits yield even smaller p-values. It appears that as the interaction is successfully uncovered, the p-values can still get progressively smaller even though the sample size is decreased. The right branch involving all three SNPs is significant with a p-value of 1.1 × 10⁻¹⁷, indicating that there is substantial interaction among SNPs, i.e. “haplotype effect”. Therefore, even in this odd situation, RP found one haplotype perfectly (x₁ = x₂ = x₃ = 1). It also found two of the three alleles defining the second high penetrance haplotype (x₁ = x₂ = x₃ = 0).

To examine the “no marginal effect” condition in more detail, the haplotype classes need to be reordered so that the first i = 1, ..., k of the total of m haplotypes contain the SNP of interest. For example, suppose we are looking at the effect of allele ”0” at the first SNP in Table 1. Then i = 1, ..., 3. The marginal penetrance of allele ”0” is

The marginal penetrance of allele ”1” (or all remaining alleles) is

and the population prevalence is

The condition of ”no marginal effect” can be expressed as

Consider some specific situations.

where and . The following results are based on 150,000 simulations.

Examples and simulation results presented in this section suggest that even in situations where the penetrance configuration is in favor of no marginal SNP effect, there needs to be a very specific set of population frequencies of haplotypes in order for the marginal effect to be absent. We considered highly interacting models so that the effect of a high-penetrance joint set of predictors (e.g., “haplotype”) is expected to be larger than that of a SNP. Nevertheless, one would need to balance the potential increase in the effect size with the number of looks through potential models. In this light, an analyst using the RP algorithm makes a very reasonable “bet” that marginal effects are going to be present among the markers that form the best predictive combination. The multiple testing implication of this sequence of order k searches is the reduction in the false discovery rate.

where p = Pr(AA | Y), and q = Pr(AA | N). An estimate is obtained by plugging in estimates of p, q, using sample frequencies of AA in cases and controls,. If the genotype AA is predictive of Y, then is the positive predictive value (PPV) of AA.

Quite usefully, the very first split of a RP tree allows a direct calculation of . Figure 6 shows a two-level tree where the first split is made based on the presence of a predictive marker M.

We can classify Figure 6 observations at the first split as follows:

(For multi-way RP splits, the counts in the second ”not-AA”() column are obtained by adding up the counts from all nodes except ”AA”). The frequency estimates in cases and controls are

Now we can estimate susceptibilities at the first two ”AA” and ”Aa + aa” nodes as

Note that (PPV at the right node) is one minus the quantity referred to as the “negative predictive value”, NPV, (Pepe, 2003) at the left (”AA”) node.

Susceptibility estimation can proceed to subsequent RP splits using the same approach. Consider , which is the probability that an individual with the genotype ”Aa or aa” and ”BB” will develop the condition Y,

The estimate is

where

Note that equation (16) is similar to (14), except that it is based on the estimates conditional on the previous (”AA” vs. ”Aa + aa”) split. For example, in (16) is the estimated prevalence among individuals carrying ”Aa + aa” genotypes, whereas in (14) is the unconditional population prevalence. Thus, node-specific genetic susceptibility can be calculated recursively by using counts and proportions of cases at the node that give estimates, as well as the susceptibility at the previous split,.

A more direct estimate of is obtained by recursively substituting the definition of into , etc. up to . After simplification, this leads to the following general formula:

with products giving joint multilocus (defined by the RP branch above the node) frequencies in cases and controls. This expression assumes that the population prevalence, ₀, is known or well estimated from the external data. When interval statements are made about , the uncertainty in ₀ may become an issue. For clarity, we will assume here that ₀ is known, however it is possible to incorporate uncertainty in ₀ into the interval estimate of susceptibility (Zaykin et al., 2004).

The right hand side of (17) simplifies on the logit scale:

where are estimated frequency of the multilocus genetic profile defining the ith node in cases and controls respectively. Referring to Figure 1, node ”N12”, containing 5 observations with the proportion of cases 0.4, the frequency of the profile (x₃ = 0, x₁ = 1) in the cases is and the frequency of that profile in the controls is .

Define . Assuming ₀ is constant, the variance of the logit is

using first-order Taylor series approximation. Assuming binomial sampling,

Then the estimated variance in terms of the node counts is

The binomial assumption in (20) implies random sampling from the case and the control population as well as that the tree is not built in “exploratory” mode. Typically, a tree is built using either the most significant predictor for a split or is “guided” by the scientist. Both situations can result in biased phenotypic proportions in the nodes. In this exploratory mode, the RP splits are usually chosen so that the association test statistic value is one of the largest among potential split variables (genetic markers). The distribution of n_iu_i will follow the distribution of one of the first order statistics from the binomial distribution. Thus, the variance calculated by (21) may be underestimated. If there is a hold-out sample, it can be used to validate a particular RP tree built in the exploratory mode and will provide unbiased estimates. Using such a hold-out sample, it is possible to construct confidence intervals for estimates. The asymptotic normal-based confidence interval for the logit is calculated as:

Although RP does not provide a perfect solution, its appealing feature is the remarkably reduced number of putative models. Given that the distribution of haplotype frequencies is often skewed and is unlikely to match the corresponding susceptibilities in the sense of satisfying the condition (12), we suggest that marginal effects at individual markers are likely to be expected. Moreover, the indices in the sums of that condition need to be rearranged for every different SNP that is considered, greatly increasing chances of marginal effects associated with at least one of the SNPs.