PMCCPMCCPMCC

Search tips
Search criteria 

Advanced

 
Logo of nihpaAbout Author manuscriptsSubmit a manuscriptHHS Public Access; Author Manuscript; Accepted for publication in peer reviewed journal;
 
Biostatistics. Author manuscript; available in PMC 2009 May 17.
Published in final edited form as:
PMCID: PMC2683243
NIHMSID: NIHMS108120

Retrospective analysis of haplotype-based case–control studies under a flexible model for gene–environment association

Summary

Genetic epidemiologic studies often involve investigation of the association of a disease with a genomic region in terms of the underlying haplotypes, that is the combination of alleles at multiple loci along homologous chromosomes. In this article, we consider the problem of estimating haplotype–environment interactions from case–control studies when some of the environmental exposures themselves may be influenced by genetic susceptibility. We specify the distribution of the diplotypes (haplotype pair) given environmental exposures for the underlying population based on a novel semiparametric model that allows haplotypes to be potentially related with environmental exposures, while allowing the marginal distribution of the diplotypes to maintain certain population genetics constraints such as Hardy–Weinberg equilibrium. The marginal distribution of the environmental exposures is allowed to remain completely nonparametric. We develop a semiparametric estimating equation methodology and related asymptotic theory for estimation of the disease odds ratios associated with the haplotypes, environmental exposures, and their interactions, parameters that characterize haplotype–environment associations and the marginal haplotype frequencies. The problem of phase ambiguity of genotype data is handled using a suitable expectation–maximization algorithm. We study the finite-sample performance of the proposed methodology using simulated data. An application of the methodology is illustrated using a case–control study of colorectal adenoma, designed to investigate how the smoking-related risk of colorectal adenoma can be modified by “NAT2,” a smoking-metabolism gene that may potentially influence susceptibility to smoking itself.

Keywords: Case-control studies, EM algorithm, Gene-environment interactions, Haplotype, Semiparametric methods

1. Introduction

Genetic epidemiologic studies often involve investigation of the association between a disease and a candidate genomic region of biologic interest. Typically, in such studies, genotype information is obtained on multiple loci that are known to harbor genetic variations within the region of interest. An increasingly popular approach for analysis of such multilocus genetic data are haplotype-based regression methods, where the effect of a genomic region on disease risk is modeled through “haplotypes,” the combinations of alleles (gene variants) at multiple loci along individual homologous chromosomes. It is believed that association analysis based on haplotypes, which can efficiently capture inter-loci interactions as well as “indirect association” due to “linkage disequilibrium” of the haplotypes with unobserved causal variant(s), can be more powerful than more traditional locus-by-locus methods (Schaid, 2004).

A technical problem for haplotype-based regression analysis is that in traditional epidemiologic studies, the haplotype information for the study subjects is not directly observable. Instead, locus-specific genotype data are observed, which contain information on the pair of alleles a subject carries on his/her pair of homologous chromosomes at each of the individual loci but does not provide the “phase information,” that is which combinations of alleles appear across multiple loci along the individual chromosomes. In general, the genotype data of a subject will be phase ambiguous whenever the subject is heterozygous at 2 or more loci. Statistically, the lack of phase information can be viewed as a special missing data problem.

Recently, a variety of methods have been developed for haplotype-based analysis of case–control data using the logistic regression model (Zhao and others, 2003; Lake and others, 2003; Epstein and Satten, 2003; Satten and Epstein, 2004; Spinka and others, 2005; Lin and Zeng, 2006; Chatterjee and others, 2006). Two classes of methods, namely, “prospective” and “retrospective” have evolved. Prospective methods ignore the retrospective nature of the case–control design. In the classical setting, without any missing data, justification of prospective analysis of case–control data relies on the well-known result about the equivalence of prospective and retrospective likelihoods under a semiparametric model that allows the distribution of the underlying covariates to remain completely nonparametric (Andersen, 1970; Prentice and Pyke, 1979). Even with missing data, the equivalence of the prospective and retrospective likelihood may hold, provided the covariate distribution is allowed to remain unrestricted (Roeder and others, 1996). For haplotype-based genetic analysis, however, complete nonparametric treatment of the covariates, including haplotypes, may not be possible due to intrinsic identifiability issues for the phase-ambiguous genotype data (Epstein and Satten, 2003). Thus, in this setting, the proper retrospective analysis of case–control data requires special attention.

An attractive feature of the retrospective likelihood is that it can enhance efficiency of case–control analysis by directly incorporating certain type of covariate distributional constraints that are natural for genetic epidemiologic studies. The assumptions of Hardy–Weinberg equilibrium (HWE) and gene–environment independence are 2 prime examples of such constraints. The HWE model, which specifies simple relationships between “allele” and “genotype” frequencies at a given chromosomal locus or between haplotype and diplotype (pair of haplotypes on homologous chromosomes) frequencies across multiple loci, is a natural law for a random mating large stable population. Often, it is also natural to assume that a subject’s genetic susceptibility, a factor which is determined at birth, is independent of his/her subsequent environmental exposures. However, if these assumptions are violated in some situations, then retrospective methods can produce serious bias in odds ratio estimates (see, e.g. Satten and Epstein, 2004; Chatterjee and Carroll, 2005; Spinka and others, 2005). Thus, there is a need for alternative flexible models for specifying the joint distribution of genetic and environmental covariates that could be used to assess the sensitivity of the retrospective methods to underlying assumptions as well as to develop alternative robust methods.

Both Satten and Epstein (2004) and Lin and Zeng (2006) have described retrospective maximum likelihood analysis of case–control data under flexible population genetics models that can relax the HWE assumption. Moreover, Lin and Zeng considered a model that allows the conditional distribution of environmental exposure given unphased genotypes to remain completely nonparametric, but they assumed conditional independence between haplotypes and the environmental factors given the unphased genotypes. If, however, haplotypes are the underlying biologic units through which a mechanism of gene is determined, then it is more natural to allow for direct association between haplotypes and environmental exposures. Moreover, if such association could exist, then quantifying the association between haplotypes and certain type of environmental exposures, such as lifestyle and behaviorial factors, would be of scientific interest.

In this article, we propose methods for retrospective analysis of case–control data using a novel model for the gene–environment distribution that can account for direct association between haplotypes and environmental exposures. The model is developed in Section 2. We assume a standard logistic regression model to specify the disease risk conditional on diplotypes and environmental exposures. In addition, we assume a polytomous logistic regression model for specifying the population distribution of the diplotypes conditional on the environmental exposures, with the intercept parameters of the model specified in such a way that the “marginal” distribution of the diplotypes can follow certain population genetic constraints such as HWE. Moreover, by exploiting the equivalence of prospective and retrospective odds ratios under the polytomous regression model, we further incorporate certain constraints on the diplotype-exposure odds ratio parameters that could reflect specific “mode of effects” for the haplotypes. We allow the marginal distribution of the environmental exposure to remain completely nonparametric.

Under the proposed modeling framework, we then describe in Section 3 a “semiparametric” estimating equation method for inference about the finite-dimensional parameters of interest, namely the disease odds ratios, haplotype frequencies, and haplotype-exposure odds ratios. We develop a suitable expectation-maximization (EM) algorithm to account for the phase-ambiguity problem. We study asymptotic theory of the proposed estimator under the underlying semiparametric setting.

In Section 4, we assess the finite-sample performance of the proposed estimator based on case–control data that were simulated utilizing haplotype patterns and frequencies obtained from a real study. In Section 5, we apply the proposed methodology to a case–control study of colorectal adenoma to investigate whether certain haplotypes in the smoking metabolism gene, NAT2, could modify smoking-related risk of colorectal adenoma and whether the same haplotypes could influence an individual’s susceptibility to smoking as well. Section 6 contains concluding remarks. All technical details are in an appendix. A SAS macro is available from the Web site http://www.stat.sinica.edu.tw/yhchen/download.htm.

2. Notations and proposed model

For haplotype-based studies, the underlying genetic covariate for a subject is defined by “diplotypes,” that is, the 2 haplotypes the individual carries in his/her pair of homologous chromosomes, where each haplotype is the combination of alleles at the loci of interest along an individual chromosome. Following the notation developed in Spinka and others (2005), let the diplotype status for a subject be Hdi = (H1, H2), where H1 and H2 denote the constituent haplotypes. We assume that there are J possible haplotypes indexed by hj for j = 1, …, J. The diplotypes are then indexed by equation M1, j1 = 1, … j1, j2 = 1, …, j2. The diplotype data, however, is not directly observable. Instead, for each subject, the multilocus genotype data G is observed, which contains information on the pair of alleles the individual carries at each individual locus but does not provide the phase information, that is which combination of alleles appears along each of the individual chromosomes. Thus, the same genotype data G could be consistent with multiple diplotypes. We will denote An external file that holds a picture, illustration, etc.
Object name is nihms108120ig1.jpg (G) to be the set of all possible diplotypes that are consistent with the genotype data G.

Given the diplotype data Hdi and a set of environmental covariate X, we assume that the risk of the disease is given by the logistic regression model

equation M2
(2.1)

for some known function m(·, β1). Often one further imposes structural assumptions on the odds ratio parameters β1 by modeling the effect of the diplotypes through constituent haplotypes according to a “dominant,” “additive,” or “recessive” mode of effect (Wallenstein and others, 1998). For example, a logistic regression model which assumes an additive effect for each copy of a haplotype corresponds to

equation M3
(2.2)

where βX is the main effect of X, βhjk is the main effect of haplotype hjk, k = 1, 2, and βhjk X is the interaction effect of X with haplotype hjk, k = 1, 2. Such modeling may be necessary due to identifiability considerations (Epstein and Satten, 2003) and is desirable when the effects of the haplotypes themselves are of direct scientific interest.

Unlike Spinka and others (2005), who assumed independence of Hdi and X, we assume a general polytomous logistic regression for the conditional distribution of Hdi given X:

equation M4
(2.3)

where equation M5 is a chosen reference diplotype. Observe that model (2.3) allows association between Hdi and X through the regression parameters γ1j1j2. Let γ0 and γ1 denote the vectorized forms for the parameters γ0 j1j2 and γ1j1j2. Let qhap(hdi|x, γ0, γ1) denote pr(Hdi = hdi|X = x) as defined by model (2.3). We allow the marginal distribution of X, denoted by F(x), to remain completely unspecified. If Hdi were directly observable, then, in principle, no further assumptions are necessary, and one can estimate γ0 and γ1 together with the odds ratio parameters of the disease risk using the profile likelihood approach developed by Chatterjee and Carroll (2005). In the presence of phase ambiguity, however, the diplotypes being not directly observable, further constraints on the parameters γ0 and γ1 are needed for the purpose of identifiability. In the following, we show how certain natural genetic models can be used to impose these constraints.

Given that genetic susceptibility may influence environmental exposures and not vice versa, for causal interpretation of parameters it is more natural to consider a model for the environmental exposures given the diplotypes. However, the odds ratios associated with the distributions [X|H] and [H|X] being the same, the parameters in γ1 can be interpreted as measures of “diplotype effects” on the distribution of exposure. Thus, it is natural to specify the γ1 parameters according to certain mode of effects of the underlying haplotypes. For example, assuming an additive effect for the haplotypes, one can write γ1 j1j2 = γ1, j1 + γ1, j2, which allows the diplotype effects to be determined by a reduced set of “haplotype effect” parameters γ1, j; in this case, γ1 would denote the vectorized form for the parameters γ1, j. Similarly, other commonly used models, such as dominant or recessive models, could be used to impose natural constraints on the γ1 parameters in model (2.3). We also observe that the parametric model (2.3), combined with the non-parametric distribution F(x), imposes a semiparametric model on the distribution of [X|H] with a density

equation M6

This class of semiparametric models includes the parametric submodel where X|Hdi = hdi follows a multivariate normal distribution with mean μhdi and common variance–covariance matrix Σ. In this case, it is easy to see that equation M7, which is a measure of the shift in the mean of the distribution of X due to differences in the diplotypes.

The parameter γ0 in model (2.3) defines the population diplotype frequencies for a baseline value of the exposure X. It is common to use population genetics models, such as HWE, to specify a relationship between diplotype and haplotype frequencies. However, observe that if the diplotypes can influence certain environmental exposures, then the frequencies of the diplotypes within exposure categories may not follow the HWE constraints although the underlying population, as a whole, may be in HWE. Thus, the population-level marginal haplotype-pair distribution is assumed to follow HWE and is characterized by the parameters θ = (θ2, …, θJ) so that

equation M8
(2.4)

where h1 denotes the chosen reference haplotype and θ1 = 0. Let

equation M9

be the marginal frequency for the diplotype hdi. Recall that in the proposed model, γ0 is defined as an implicit function of γ1, θ, and F(x) through the relationship

equation M10
(2.5)

Note that F is left unspecified, and hence the model propoised is semiparametric.

3. Semiparametric estimating equation inference

3.1 Estimation with known haplotypes

In what follows, where there can be no confusion, we will write h for hdi.

Let H(x) = exp(x)/{1+ exp(x)} be the logistic distribution function. Write the risk model probability as

equation M11

Recall that qhap(h|X, γ0, γ1) = pr(Hdi = h|X; γ0, γ1) is the conditional model of Hdi given X that is specified as in (2.3).

To start with, consider the ideal case that the phase information is known so that Hdi is observed. Since F is treated nonparametrically, assume that F is discrete and has mass δk at xk, k = 1, …, K, where {x1, …, xK} are the distinct values of X that are observed in the case–control sample. Let ndkh be the number of subjects in the sample with (D = d, X = xk, Hdi = h). Ignoring the dependence of γ0 on F tentatively, the log-likelihood of the case–control data can then be written as

equation M12

Maximizing l with respect to δ for fixed values of ω = (β, γ0, γ1) then leads to

equation M13
(3.1)

and the profile log-likelihood

equation M14
(3.2)

where

equation M15

and

equation M16

with B = (β0, β1, κ)T and κ = β0 + log(n1/n0) − log{pr(D = 1)/pr(D = 0)}. The calculation is similar to that in Chatterjee and Carroll (2005).

As noted by Chatterjee and Carroll (2005), the parameter β0 is separable from κ and hence is theoretically identifiable. In practice, however, there is usually little information about β0 available in the observed data, and hence the information matrix is nearly singular. One way to bypass this problem is to use external information on the disease prevalence pr(D = 1), while another way is to use the rare-disease approximation when the disease is rare. The estimation method described below can be applied to both the 2 cases of pr(D = 1) being known and the rare-disease approximation being made, with suitable definitions on B and S(d, h, x, B, γ0, γ1). When pr(D = 1) is known, κ depends on β0 only, hence here we define B = (β0, β1)T. When the disease is rare so that

equation M17
(3.3)

we have

equation M18

Note that β0 does not appear in this expression, and hence we define B = (κ, β1)T.

Our goal is to estimate the parameters (B, θ, γ1) based on the profile log-likelihood (3.2), where γ0 is defined as an implicit function of (θ, γ1, F) through (2.5), and we write γ0 = An external file that holds a picture, illustration, etc.
Object name is nihms108120ig2.jpg (θ, γ1, F). Let Ω = (B, γ1, γ0), Ω* = {B, γ1, An external file that holds a picture, illustration, etc.
Object name is nihms108120ig2.jpg (θ, γ1, F)}, and Φ = (B, γ1, θ). Let LΩ(·) and LΦ(·) be, respectively, the derivatives of L(·) with respect to Ω and Φ, and An external file that holds a picture, illustration, etc.
Object name is nihms108120ig2.jpg θ and An external file that holds a picture, illustration, etc.
Object name is nihms108120ig2.jpg γ1 the derivatives of An external file that holds a picture, illustration, etc.
Object name is nihms108120ig2.jpg(·)with respect to and θ and γ1 We then have

equation M19

Explicit expressions for An external file that holds a picture, illustration, etc.
Object name is nihms108120ig2.jpg γ1 and An external file that holds a picture, illustration, etc.
Object name is nihms108120ig2.jpg θ; are given in Appendix C. Also, the information matrix is given by

equation M20

where x2110ΩΩ = E(−LΩΩ), with LΩΩ the second derivative of L with respect to Ω; note that the terms involving second derivatives of Ω* do not appear in the information matrix because E(LΩ) = 0, which is a direct consequence of the Lemmas A.1 and A.2 in Appendix A. We propose to obtain the estimate of Φ by solving the estimating equation

equation M21
(3.4)

where we have substituted an estimate F for F in An external file that holds a picture, illustration, etc.
Object name is nihms108120ig2.jpg(·); that is, for each fixed value of (θ, γ1), we solve γ0 from

equation M22
(3.5)

One convenient choice of F is the empirical estimate Femp, which is given by

equation M23

for the case where pr(D = 1) is known, where Femp,1(x) and Femp,0(x) are the empirical distributions of X in the case and control samples, and is given by Femp(x) = Femp,0(x) for the case where the rare-disease assumption can be made. An alternative choice of F (x) would be the profile likelihood estimate (3.1). Numerical calculations not given here show that the latter choice requires more computational efforts while yielding results very similar to those given by the empirical estimate Femp.

3.2 Estimation with ambiguous haplotype data

Now, we turn to the more practical case where the haplotype data cannot be observed directly and must be inferred from the unphased genotype data, that is, the haplotype information may be subject to ambiguity. In this case, we apply an EM-like algorithm to the “complete data” estimating equation (3.4). Let Gi denote the observed unphased genotype of subject i and An external file that holds a picture, illustration, etc.
Object name is nihms108120ig1.jpg (Gi) the set of diplotypes that are consistent with Gi. When only Gi instead of equation M24 is observed for each subject, we propose to obtain the estimate [Phi w/ hat] for Φ = (B, γ1, θ) as the solution of the weighted version of (3.4):

equation M25
(3.6)

where using the short-hand notation that [gamma with circumflex]0 = An external file that holds a picture, illustration, etc.
Object name is nihms108120ig2.jpg(θ, γ1, Femp), the weights are given by

equation M26
(3.7)

The limiting version of the weights is given as

equation M27
(3.8)

Solving the estimating equation (3.6) can be implemented simply by an EM-like algorithm as follows: starting with an initial value for Φ and hence an initial value for γ0, we

  1. calculate the weights {[omega with circumflex]i} from (3.7) and
  2. solve (3.6) to obtain an updated estimate of Φ using the weights {[omega with circumflex]i} given in (i); note that within this step we also need to solve (3.5) to obtain updated value of γ0.

The algorithm is iterated between the 2 steps until convergence. Note that the weights {[omega with circumflex]i} are only used in solving Φ from (3.6) and are not required in solving γ0 from (3.5).

3.3 Asymptotic theory

Make the following series of definitions. Expectations denoted as Ecc(·) are taken under the case–control sampling design, that is, for any random vector Y, equation M28, where μd = plim nd/n, d = 0, 1. Then, define

equation M29
(3.9)

Note that the second derivative of Ω* does not appear in I since E(LΦ) = E(LΦ) = 0, and the last identity in (3.9) is given by Lemma A.3 in Appendix A.

Define [p with hat]emp(Di) to be the mass of Femp(Xi), which is equal to equation M30 if πd = pr (D = d) is known and is equal to I (Di = 0)/n0 when the rare-disease approximation is used. Let qhap(X, γ1, γ0) = {qhap(h|X, γ1, γ0)} be the vector collection over h of qhap(h|X, γ1, γ0) for all diplotypes except the reference diplotype, and let qHWE(θ) be defined similarly. Define

equation M31

where x2110Ωγ0 is the obvious submatrix of x2110ΩΩ.

Theorem 3.1

Let

equation M32

Suppose that Ecc{E(·) ET (·)} exists and the matrix I is invertible. Then, n1/2([Phi w/ hat] − Φ) is asymptotically normal with mean zero and covariance matrix

equation M33
(3.10)

Remark 3.2

The asymptotic variance [Gamma with macron] can be readily estimated by replacing each component matrix with its empirical counterpart. Lemma A.3 gives useful expressions to facilitate this computation.

Remark 3.3

In our numerical experiments, the estimated covariance based on formula (3.10) is very close to that based on the “naive” covariance estimate obtained by naively treating the estimating equation (3.6) as a genuine score equation; namely, treating the Femp plugged into An external file that holds a picture, illustration, etc.
Object name is nihms108120ig2.jpg(·) as the true covariate distribution F. In this case, by applying Proposition 1(ii) in Chatterjee and Carroll (2005), the naive co-variance estimate can be obtained simply as the empirical counterpart of the matrix I−1I−1Ψx2110−1, where equation M34. Whether this naive estimate performs well in general is unknown, and we suggest using the estimate based on (3.10).

4. Simulations

4.1 Finite-sample performance under correct model

In this section, we study the finite-sample performance of the proposed estimator using simulated data generated under the proposed modeling framework. We simulated haplotypes following published data (Epstein and Satten, 2003) on haplotype patterns and frequencies for 5 single-nucleotide polymorphisms (SNPs) in a putative susceptibility gene for diabetes (see Table 1). The simulations involved a single environmental covariate X, assumed to follow a standard normal distribution in the population. Given X, the diplotypes (haplotype pair) for an individual were generated from a polytomous logistic regression of the form (2.3), where the diplotype-specific odds ratios were further specified according to an additive model of the form γ1j1j2= γ1, j1 + γ1, j2, where j1 and j2 denote the index for 14 haplotypes shown in Table 1. We assume γ1,4 = γ1,5 = −0.4 and γ1,12 = 0.4, and all the other γ1, j = 0. The parameters γ0j1j2 ’s in model (2.3) are then specified in such a way that the marginal diplotype distribution follow HWE with haplotype frequencies given in Table 1.

Table 1
Marginal haplotype frequencies and grouped haplotypes in the simulation. Here, j is index for the original haplotypes while j′ is index for the grouped haplotypes

For generating disease outcome, we chose the haplotype “01100” (j = 5) to be causal and used the logistic model

equation M35

where Z(5) denotes the number of the copies of the causal haplotype contained in Hdi. The true value of the parameter vector (β0, βH, βX, βH X) was set to (−3.0, 0.2, 0.1, 0.3). A case–control sample with 600 controls and 600 cases was then sampled. The results were based upon 1000 simulated data sets.

When analyzing the data, we only used the unphased genotype information. We did not assume the causal haplotype to be known. Thus, in both the disease-risk model (2.1) and the diplotype-frequency model (2.3), we choose the most common haplotype “10011” as a reference and estimated a separate regression parameter for each of the non-referent haplotypes. Since rare haplotypes may lead to unreliable estimates of the associated regression parameters, when estimating β and γ1, rare haplotypes with frequency <1% are grouped into the reference haplotype. The resulting 8-grouped haplotypes are labeled as h j, j′ = 2, …, 8; see Table 1 for details about how the haplotypes are grouped.

In each simulation, we obtain 2 sets of estimates from the proposed method, one using the rare-disease approximation (3.3) and the other using the known value of the population disease prevalence. Results shown in Table 2 show that both sets of estimates are essentially unbiased. Also, the standard error estimates are in close agreement with the true values, and the coverage probabilities are close to the nominal value (95%). As expected, the estimates for θ and γ1 are generally more efficient using external information on the disease prevalence than when using the rare-disease approximation, but no such efficiency gain is observed for the parameters β in the disease-risk model. Similar conclusion can be drawn from the simulations with a Bernoulli covariate (success probability = 0.5), showing the applicability of the proposed method to the categorical covariate. Detailed results for this latter set of simulations are included in the supplementary material available at Biostatistics online.

Table 2
Simulation results for the case with an additive genetic law. Here, mean is the mean over 1000 simulated data sets, equation M38 is the standard deviation of the estimates, SE is the mean of the estimated standard deviation of the parameter estimates, CP is the ...

4.2 Model robustness

Here, we consider a simulation study where we generate the data in such a way that the polytomous model for diplotype frequencies may not exactly hold. The main goal is to give an indication of the robustness of the estimate of the association parameters (β) from the proposed method when the model for [Hdi|X] is misspecified.

Following the argument of causality in Section 2, if [X|Hdi] follows a normal distribution with constant variance, then the polytomous model is exact. So, to show a modest violation of the polytomous model, for given diplotype we generate the data on X from

equation M36

where the diplotype data are again generated from the distribution in Table 1, ε is a t-distribution (d.f. = 3) truncated at ±5, λ4 = λ5 = −1.2, λ12 = 1.2, and all the other λj are 0. The disease status data are generated from the same logistic model as in the previous simulation. The simulated data on 600 cases and 600 controls are then analyzed with the proposed method, where the analysis models for the disease risk, [Hdi|X], and the marginal diplotype distribution are specified the same as those in the previous simulation. As a comparison, we also fit to the simulated data a model with the haplotype–environment (H-X) independence assumption, i.e. pr(Hdi|X) = pr(Hdi) = qHWE(Hdi), using the method proposed by Spinka and others (2005). The rare-disease approximation is made when applying both the 2 methods.

The results shown in Table 3 reveal that, for the estimation of the association parameters β, the proposed method may be quite robust to modest misspecification of the model for [Hdi|X]. On the other hand, the estimates from the H-X independence method does result in substantial bias, especially for parameters corresponding to haplotypes for which [X|Hdi] have nonzero mean; for example, the estimate for the interaction parameter between h5 and X is severely biased with the H-X independence method. The estimates for the marginal haplotype-frequency parameters θ seem to be robust to misspecification of [Hdi|X] for both the 2 methods.

Table 3
Simulation results for the case of a misspecified conditional diplotype distribution given covariates. Here, mean is the mean over 1000 simulated data sets, SE is the standard deviation of the estimates, equation M41 is the mean of the estimated standard deviation ...

5. Case–control study of colorectal adenoma study, NAT2 haplotype, and smoking

We illustrate the proposed modeling and estimating methodologies with an application to a case–control study of colorectal adenoma, a precursor of colorectal cancer. The study involved 628 prevalent advanced adenoma cases and 635 gender-matched controls, selected from the screening arm of the Prostate, Lung, Colorectal and Ovarian Cancer Screening Trial at the National Cancer Institute, USA (Gohagan and others, 2000; Moslehi and others, 2006). One of the main objectives of this study is to assess whether smoking-related risk of colorectal adenoma may be modified by certain haplotypes in NAT2, a gene known to be important in metabolism of smoking-related carcinogens. In addition, since NAT2 is involved in the smoking metabolism pathway, potentially it can influence an individual’s addiction to smoking. Thus, it was also of interest to identify potential haplotypes that could influence an individual’s susceptibility to smoking.

Genotype data were available on 6 SNPs. We initially applied the EM algorithm proposed by Li and others (2003) for haplotype-frequency estimation to derive 7 common haplotypes with estimated frequency greater than 0.5%, which are then included in our association analysis with the most frequent haplotype served as the reference haplotype. Subjects were categorized as “never,” “former,” or “current” smokers. We fit a logistic regression model (2.1) assuming an additive effect for each haplotype other than the reference one; see (2.2). The haplotype–environment interaction terms include only those for the haplotype “101010” with “Smk1” and “Smk2,” 2 dummy variables for former and never smokers, because they are the only promising interactions according to preliminary analysis. The disease-risk model was further adjusted for “age,” recorded in years, and “gender.” A polytomous logistic regression (2.3) is specified for the conditional distribution of diplotypes given the environmental covariates Smk1 and Smk2 with the marginal diplotype distribution being specified by the HWE constraints. The main parameters of interest include the disease-haplotype odds ratio parameters β1, the haplotype–environment odds ratio parameters γ1, and the marginal haplotype frequencies in the whole population. The marginal distribution for the environmental covariates is left unspecified. For estimation of regression parameters β and γ1, we grouped haplotypes with frequency less than 2% into the reference haplotype. The rare-disease approximation was made in deriving the estimating equation, and the EM algorithm proposed in Section 3 is utilized to accommodate the unphased genotype data.

Results from this application are displayed in Table 4. It is clear that current smokers can have significantly elevated risk for colorectal adenoma relative to nonsmokers, adjusting for gender and age. Relative to the reference haplotype “001100,” all the other haplotypes are associated with reduced risk for colorectal adenoma, but the statistical evidence is not significant. However, the significance of the interaction 101010 × Smk2 suggests that smoking-related risk of adenoma was much reduced for carriers of the haplotype 101010 than non-carriers. The finding is consistent with previous laboratory and epidemiologic studies that have identified the haplotype 101010, known as “NAT2*4,” as a rapid metabolizer for smoking-related carcinogens. The estimates for the parameter γ1 for the conditional diplotype distribution reveal that the susceptibility to smoking seems not to be influenced by any haplotypes we considered. Finally, the estimates for the marginal haplotype frequencies derived from the estimates of θ are quite close to those obtained by the EM algorithm of Li and others (2003) applied to the genotype data of the controls.

Table 4
Results of the colorectal adenoma study. Haplotypes 001010 and 001110 are grouped into the reference haplotype in the disease-risk and conditional diplotype distribution models

To check if the analysis is sensitive to model specification for the conditional distribution of diplotypes given the environmental covariates, we further fit the model (2.3) with various choices of the environmental covariates. The results (not shown) for the association parameters β and the marginal haplotype frequencies are fairly consistent across the analyses.

6. Concluding remarks

The model we have proposed for gene–environment association is suitable when the underlying haplotypes of a genomic region may causally influence the environmental exposure(s) under study. The model, however, requires special treatment for environmental factors, such as ethnicity or geographic region(s), which may be associated with the genomic region under study, not because of any causal relationship but merely due to population stratification. Suppose, in addition to the main environmental exposure X, there is a set of environmental factors S which could be used to divide the underlying population into K strata that are likely to be genetically heterogenous. In such a situation, a natural model for describing the association between diplotypes Hdi and environmental factors W = (X, S) is given by

equation M37
(6.1)

where the stratum-specific intercept parameters γ0j1j2 (S) should be specified in such a way that the diplotype frequencies, after marginalized over X, follow population genetics constraints, such as HWE, within each stratum defined by S. The disease-risk model could be also extended to include S as a risk factor. The proposed estimating equation methodology can be easily modified to estimate the gene–environment interaction and association parameters of interest under these extended models.

Acknowledgments

Chen’s research was supported by the National Science Council of the People’s Republic of China (NSC 95-2118-M-001-022-MY3). Chatterjee’s research was supported by the Intramural Research Program of the National Cancer Institute. Carroll’s research was supported by a grant from the National Cancer Institute (CA57030) and by the Texas A&M Center for Environmental and Rural Health via a grant from the National Institute of Environmental Health Sciences (P30-ES09106). A SAS macro is available from the Web site http://www.stat.sinica.edu.tw/yhchen/download.html.

APPENDIX A: BASIC LEMMAS

The following lemmas are required to derive the asymptotic distribution of the proposed estimator. Lemma A.1 below is in fact Lemma 3 of Chatterjee and Carroll (2005).

Lemma A.1

Under the case–control sampling design where the total sample size n = n1 + n0 tends to infinity but the sampling proportions for the cases and controls, that is, n1/n and n0/n, remain fixed at μ1 and μ0, we have for any measurable function M(D, Hdi, X) of data (D, Hdi, X),

equation M44

where E*(·|X) denotes the expectation with respect to the joint distribution of (D, Hdi) given X defined by

equation M45
(A.1)

and λ(x) = Σd Σh, S(d, h, x, B, γ1, γ0)

Lemma A.2 below provides an explicit expression for the estimating function LΦ(·).

Lemma A.2

Write

equation M46

Then

equation M47
(A.2)

Proof

By definition,

equation M48

and direct calculation yields

equation M49

which proves the result.

Lemma A.3 provides explicit forms for the information matrices.

Lemma A.3

Let x2110ΩΩ = Ecc{−LΩΩ(D, Hdi, X, Ω)}, where LΩΩ is the second derivative of L(·) with respective to Ω and x2110ΩΩ = Ecc{−[partial differential]LΩ(·)/[partial differential]ΩT}. Then

equation M50

Proof

The first identity has been given in Lemma 4 of Chatterjee and Carroll (2005). To show the second identity, applying the chain rule we have

equation M51
(A.3)

The first term of (A.3) equals

equation M52

By the definition of ω(h, Ω) given in (3.8), it easy to see that

equation M53

where the joint density p*(D, Hdi|X, Ω) is defined in (A.1). Recalling further that LΩ(D, h, X, Ω) = [partial differential] log p*(D, Hdi = h|X, Ω)/[partial differential]Ω and LΩ (D, G, X, Ω) = Σh ω(h, Ω)LΩ(D, h, X, Ω), we have the identity

equation M54

Hence, the second term of (A.3) leads to

equation M55

The desired result thus follows.

APPENDIX B: PROOF OF THEOREM

We will first obtain the asymptotic expansion of the proposed estimating equation (3.6), by which the large sample distribution theory for [Phi w/ hat] can be derived immediately from the central limit theorem.

For our estimator [Phi w/ hat] = (B, [gamma with circumflex]1, [theta w/ hat]) a standard Taylor series expansion of (3.6) yields

equation M56
(B.1)

In view of (2.5) and (3.5), let γ0 = An external file that holds a picture, illustration, etc.
Object name is nihms108120ig2.jpg (θ, γ1, F) and [gamma with circumflex]0 = An external file that holds a picture, illustration, etc.
Object name is nihms108120ig2.jpg(θ, γ1, Femp). Recall that [gamma with circumflex]0 is solved from (3.5):

equation M57

Making a further Taylor series expansion, we have

equation M58

Note that

equation M59

Hence,

equation M60

An explicit expression for Q is given in Appendix C. Consequently, the middle 2 terms in the expansion (B.1) reduce to

equation M61

Therefore, we have

equation M62

Note that in the last equality above, we have used the fact that

equation M63

and equation M64, which follow directly from Lemmas A.1 and A.2. This completes the proof.

APPENDIX C: EXPRESSIONS FOR THE DERIVATIVES OF An external file that holds a picture, illustration, etc.
Object name is nihms108120ig2.jpg (θ, γ1, F) AND An external file that holds a picture, illustration, etc.
Object name is nihms108120ig2.jpg(θ, γ1, Femp)

Recall that An external file that holds a picture, illustration, etc.
Object name is nihms108120ig2.jpg (θ, γ1, F)is the solution of γ0 to (2.5). Here, we define γ1 to be the vectorized form for the diplotype-effect parameters subject to certain mode of effects (e.g. additive effect) of haplotypes, and define η to be the vectorized form for the full set of diplotype effects {γ1j1j2}. Differentiating both sides of (2.5) with respect to γ1, we have

equation M65

Let

equation M66

and Q = ∫ Q(x)dF(x). Then

equation M67

Similarly, letting

equation M68

we have

equation M69

The derivatives of An external file that holds a picture, illustration, etc.
Object name is nihms108120ig2.jpg(θ, γ1, Femp) can be obtained by replacing dF in the above quantities with dFemp.

Footnotes

Conflict of Interest: None declared.

Contributor Information

YI-HAU CHEN, Institute of Statistical Science, Academia Sinica, Taipei 11529, Taiwan, People’s Republic of China.

NILANJAN CHATTERJEE, Biostatistics Branch, Division of Cancer Epidemiology and Genetics, National Cancer Institute, 6120 Executive Boulevard, EPS 8038, Rockville, MD 20852, USA, chattern/at/mail.nih.gov.

RAYMOND J. CARROLL, Department of Statistics, Texas A&M University, TAMU 3143, College Station, TX 77843-3143, USA.

References

  • Andersen JB. Asymptotic properties of conditional maximum likelihood estimators. Journal of the Royal Statistical Society, Series B. 1970;32:283–301.
  • Chatterjee N, Carroll RJ. Semiparametric maximum likelihood estimation in case-control studies of gene-environmental interactions. Biometrika. 2005;92:399–418.
  • Chatterjee N, Chen J, Spinka C, Carroll RJ. Comment on the paper likelihood based inference on haplotype effects in genetic association studies by D. J. Lin and D. Zeng. Journal of the American Statistical Association. 2006;101:108–110.
  • Epstein M, Satten G. Inference of haplotype effects in case-control studies using unphased genotype data. American Journal of Human Genetics. 2003;73:1316–1329. [PubMed]
  • Gohagan JK, Prorok PC, Hayes RB, Kramer BS. The Prostate, Lung, Colorectal and Ovarian (PLCO) Cancer Screening Trial of the National Cancer Institute: History, organization, and status. Controlled Clinical Trials. 2000;21:251S–272S. [PubMed]
  • Lake SL, Lyon H, Tantisira K, Silverman EK, Weiss ST, Laird NM, Schaid DJ. Estimation and tests of haplotype-environment interaction when linkage phase is ambiguous. Human Heredity. 2003;55:56–65. [PubMed]
  • Li SS, Khalid N, Carlson C, Zhao LP. Estimating haplotype frequencies and standard errors for multiple single nucleotide polymorphisms. Biostatistics. 2003;4:513–522. [PubMed]
  • Lin DY, Zeng D. Likelihood-based inference on haplotype effects in genetic association studies (with discussion) Journal of the American Statistical Association. 2006;101:89–118.
  • Moslehi R, Chatterjee N, Church TR, Chen J, Yeager M, Weissfield J, Hein DW, Hayes RB. Cigarette smoking, n-acetyltransferase genes and the risk of advanced colorectal adenoma. Pharmacogenomics. 2006;7:819–829. [PubMed]
  • Prentice RL, Pyke R. Logistic disease incidence models and case-control studies. Biometrika. 1979;66:403–411.
  • Roeder K, Carroll RJ, Lindsay BG. A nonparametric mixture approach to case-control studies with errors in covariables. Journal of the American Statistical Association. 1996;91:722–732.
  • Satten GA, Epstein MP. Comparison of prospective and retrospective methods for haplotype inference in case-control studies. Genetic Epidemiology. 2004;27:192–201. [PubMed]
  • Schaid DJ. Evaluating associations of haplotypes with traits. Genetic Epidemiology. 2004;27:348–364. [PubMed]
  • Spinka C, Carroll RJ, Chatterjee N. Analysis of case-control studies of genetic and environmental factors with missing genetic information and haplotype-phase ambiguity. Genetic Epidemiology. 2005;29:108–127. [PMC free article] [PubMed]
  • Wallenstein S, Hodge SE, Weston A. Logistic regression model for analyzing extended haplotype data. Genetic Epidemiology. 1998;15:173–181. [PubMed]
  • Zhao LP, Li SS, Khalid NA. Method for the assessment of disease associations with single-nucleotide polymorphism haplotypes and environmental variables in case-control studies. American Journal of Human Genetics. 2003;72:1231–1250. [PubMed]