Re: The Use of Inferred Haplotypes in Downstream Analysis

To the Editor: In a letter published in the March 2007 issue of the Journal, Lin and Huang1 described some potential pitfalls in haplotype analysis when only unphased genotyped data are available. In particular, they relate the problems associated with picking the most probable haplotype from a distribution of potential haplotypes to use in risk analysis, replacing the unobserved haplotypes with their “best estimates.” This particular form of “single imputation” can indeed be dangerous. However, another form of single imputation exists that is far more reliable, as has been shown in a series of articles describing analytical results and a series of simulation experiments.2,3 This “correct” single-imputation method is the expectation-substitution method first described for haplotype analysis by Zaykin et al.,4 but which has analogues in many parts of biostatistical analysis dealing with measurement errors.5,6 Consider a model in which the number of copies carried of a certain haplotype, h, is related to the odds of disease in a log-linear fashion. The expectation-substitution method works as follows: since the number of copies of the risk haplotype carried, nh, is unknown for most individuals when only genotype data are available, one replaces this unknown quantity with its conditional expectation, E(nh|G), where G is the observed set of (unphased) markers, and the expectation is generally computed under the assumption that the haplotype frequencies are known and under the assumption of Hardy-Weinberg equilibrium. This expectation (which generally takes noninteger values 0–2 and thus is definitely not equivalent to using the most likely value of nh, which must take values 0, 1, or 2) is then used as if it were the true value of nh, with no other allowance for the uncertainty of the estimation of nh. In particular, we perform score tests and likelihood-ratio tests of the null hypothesis (that h is unrelated to risk) as if nh were known for all subjects. This method is not limited to log-linear penetrance models and can be extended to dominant, recessive, or codominant models by calculating the expected values of appropriate diplotype codings.2 It has been used in a large number of recent analyses of haplotype-specific risk in candidate-gene studies.7–10 The following issues arise in the application of this method: 1. Because haplotype frequencies must be estimated from a finite amount of genotype data, the calculation of the expectation is itself uncertain, and this uncertainty is being neglected. What effect does this have? 2. Because odds-ratio (OR) models are not strictly linear in nh, the expected OR as a function of nh is not equal to the OR function applied to E(nh). Is this an important problem? 3. Doesn’t case-control sampling distort haplotype frequencies? Should the expectation be applied separately to the cases and controls, or should the cases and controls be combined to compute the expectation? Regarding issue 1, it turns out11 that the score test arising from the expectation method is both valid and asymptotically fully efficient for testing the null hypothesis of no haplotype-specific risk associated with haplotype h against the alternative (that disease odds are log-linear in nh). Thus, for score tests under the null hypothesis, there is no need to account for the uncertainty of the estimates of the haplotype frequencies on which the calculation of E(nh) is based. Similarly, under the null hypothesis, case-control sampling does not distort the distribution of the risk haplotype, so the calculation of the expectation is not affected. This means that we do not need to account for case-control sampling either when testing the null hypothesis. Simulations under the alternative hypothesis2 tend to show that, under the alternative hypothesis (that disease is associated with haplotype count), the expectation-substitution method also gives quite reasonable estimates and confidence limits for the value of the risk parameter (log-OR per copy of h) in most practical settings. Here, we present the results of two simulation experiments. We considered estimating the association between five-SNP haplotypes and disease risk in two situations: low haplotype diversity (table 1) and high haplotype diversity (table 2). For the first situation, only seven haplotypes were present: 00000, 00001, 00010, 00011, 00100, 01000, and 10000 with frequencies 0.35, 0.15, 0.15, 0.05, 0.1, 0.1, and 0.1, respectively (R2h>0.7 for all haplotypes, where R2h is the squared correlation of expected versus true haplotype counts12). For the second situation, all 32 possible haplotypes were equally likely. This is something of a worst-case scenario for haplotype association analysis, both because of the high ambiguity in haplotypes given the genotypes (R2h≈0.43 for all haplotypes) and because all haplotypes are relatively uncommon (frequency <0.04). We compared the performance of logistic regression, using the true haplotypes, the most likely haplotype pair given the genotype data, the expectation-substitution method, and the prospective likelihood implemented in the haplo.glm function of the R package haplo.stats.13 This last approach is similar in spirit to that advocated by Lin and Huang,1 since it maximizes a likelihood that integrates over the missing haplotype-phase information. (We chose haplo.glm primarily for computational simplicity. This allowed us to perform simulations, analysis, archiving, and summary in one statistical package, R. Although HAPSTAT, the graphical user interface implementation of the maximum-likelihood methods developed by Lin et al.,14–16 is quite user-friendly and suitable for the analysis of single data sets, it does not appear to be easily automated for the analysis of hundreds of simulated data sets.) Table 1. Comparison of Haplotype-Association Methods under Low Haplotype Diversity[Note] Table 2. Comparison of Haplotype-Association Methods under High Haplotype Diversity[Note] As can be seen from table 1, under low haplotype diversity, the expectation-substitution and the maximum-likelihood method give indistinguishable results. In fact, the average euclidean distance between the vectors of haplotype log-ORs for the two methods is <0.007 for all the situations presented in table 1, and the average distance between the haplotype frequency vectors was <0.001. Neither method has noticeable bias, and both have appropriate coverage for ORs of 1.5 and 3.0. The most-likely-haplotype method is slightly biased toward the null. On the other hand, both the expectation-substitution and maximum-likelihood methods show evidence of bias under high haplotype diversity (table 2). The expectation-substitution method shows modest bias toward the null, whereas the maximum-likelihood method shows stronger bias away from the null. Moreover, the nominal 95% CIs for the expectation-substitution method have appropriate coverage, whereas those for the maximum-likelihood method are far too small, and the average estimated variance for the estimated log-OR is noticeably smaller than the observed variance in the maximum-likelihood estimates. We hypothesize that this is due to the relatively large number of highly colinear parameters (64) that must be jointly estimated from data from 600 cases and 600 controls. We emphasize that the high-diversity situation should not often arise in practice, because haplotype association analyses are generally restricted to regions of low haplotype diversity (appropriately so, in our opinion), but it is interesting that even in this situation the expectation-substitution method shows only modest bias and retains appropriate coverage, whereas the maximum-likelihood method performs poorly, perhaps because of numerical difficulties. Although we used the expectation-maximization algorithm to estimate haplotype frequencies and posterior haplotype probabilities conditional on individual genotypes, one could also use more-sophisticated algorithms, such as those implemented in PHASE.17 PHASE has been shown to provide more-accurate estimates of haplotype frequencies than does the expectation-maximization algorithm,18 in part because it models mutation and recombination processes, whereas the expectationmaximization algorithm only assumes Hardy-Weinberg equilibrium. The fact that the expectation-substitution approach allows the user to choose from a range of haplotype-frequency–estimation algorithms is a potential advantage over the maximum-likelihood approach, although, for most situations where haplotype association analysis is applied—small numbers of SNPs in high linkage disequilibrium over short distances—we anticipate that the difference between haplotype-frequency estimates from PHASE and from the expectation-maximization algorithm will be quite small. Our conclusions based on this type of simulation are that the expectation-substitution method provides very reliable inference (correct type I error rates under the null hypothesis), good power under alternatives, and little bias either in overall estimates or in confidence limits. It appears to be that only when the true ORs become extremely large do some problems occur with the method, and, frankly, from an epidemiological perspective, we should be so lucky as to have very many association studies with this problem!