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

doi:10.1093/biostatistics/kxm011

Biostatistics. Author manuscript; available in PMC 2009 May 17.

Published in final edited form as:

Biostatistics. 2008 January; 9(1): 81–99.

Published online 2007 May 8. doi: 10.1093/biostatistics/kxm011

PMCID: PMC2683243

NIHMSID: NIHMS108120

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

YI-HAU CHEN, NILANJAN CHATTERJEE,^* and RAYMOND J. CARROLL

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

Contributor Information.

^* To whom correspondence should be addressed

The publisher's final edited version of this article is available free at Biostatistics

See other articles in PMC that cite the published article.

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 H^di = (H₁, H₂), where H₁ and H₂ denote the constituent haplotypes. We assume that there are J possible haplotypes indexed by h_j for j = 1, …, J. The diplotypes are then indexed by equation M1

, j₁ = 1, … j₁, j₂ = 1, …, j₂. 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

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 H^di and a set of environmental covariate X, we assume that the risk of the disease is given by the logistic regression model

(2.1)

for some known function m(·, β₁). Often one further imposes structural assumptions on the odds ratio parameters β₁ 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

(2.2)

where β_X is the main effect of X, β_{h_{j_k}} is the main effect of haplotype h_{j_k}, k = 1, 2, and β_{h_{j_k} X} is the interaction effect of X with haplotype h_{j_k}, 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 H^di and X, we assume a general polytomous logistic regression for the conditional distribution of H^di given X:

(2.3)

where

is a chosen reference diplotype. Observe that model (2.3) allows association between H^di and X through the regression parameters γ₁_j_₁_j_₂. Let γ₀ and γ₁ denote the vectorized forms for the parameters γ₀ _j_₁_j_₂ and γ₁_j_₁_j_₂. Let qhap(h^di|x, γ₀, γ₁) denote pr(H^di = h^di|X = x) as defined by model (2.3). We allow the marginal distribution of X, denoted by F(x), to remain completely unspecified. If H^di were directly observable, then, in principle, no further assumptions are necessary, and one can estimate γ₀ and γ₁ 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 γ₀ and γ₁ 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 γ₁ can be interpreted as measures of “diplotype effects” on the distribution of exposure. Thus, it is natural to specify the γ₁ parameters according to certain mode of effects of the underlying haplotypes. For example, assuming an additive effect for the haplotypes, one can write γ₁ _j_₁_j_₂ = γ₁, _j_₁ + γ₁, _j_₂, which allows the diplotype effects to be determined by a reduced set of “haplotype effect” parameters γ₁, j; in this case, γ₁ 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 γ₁ 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

This class of semiparametric models includes the parametric submodel where X|H^di = h^di 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 γ₀ 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 θ = (θ₂, …, θ_J) so that

(2.4)

where h₁ denotes the chosen reference haplotype and θ₁ = 0. Let

be the marginal frequency for the diplotype h^di. Recall that in the proposed model, γ₀ is defined as an implicit function of γ₁, θ, and F(x) through the relationship

(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 h^di.

Let

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

Recall that q_hap(h|X, γ₀, γ₁) = pr(H^di = h|X; γ₀, γ₁) is the conditional model of H^di given X that is specified as in (2.3).

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

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

(3.1)

and the profile log-likelihood

(3.2)

where

and

with

= (β₀, β₁, κ)^T and κ = β₀ + log(n₁/n₀) − 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 β₀ is separable from κ and hence is theoretically identifiable. In practice, however, there is usually little information about β₀ 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

and S(d, h, x,

, γ₀, γ₁). When pr(D = 1) is known, κ depends on β₀ only, hence here we define

= (β₀, β₁)^T. When the disease is rare so that

(3.3)

we have

Note that β₀ does not appear in this expression, and hence we define

= (κ, β₁)^T.

Our goal is to estimate the parameters (

, θ, γ₁) based on the profile log-likelihood (3.2), where γ₀ is defined as an implicit function of (θ, γ₁, F) through (2.5), and we write γ₀ = An external file that holds a picture, illustration, etc.
Object name is nihms108120ig2.jpg

An external file that holds a picture, illustration, etc.
Object name is nihms108120ig2.jpg

(θ, γ₁, F). Let Ω = (

, γ₁, γ₀), Ω^* = {

, γ₁,

(θ, γ₁, F)}, and Φ = (

, γ₁, θ). Let

_Ω(·) and

_Φ(·) be, respectively, the derivatives of

(·) with respect to Ω and Φ, and An external file that holds a picture, illustration, etc.
Object name is nihms108120ig2.jpg

_θ and

_γ_₁ the derivatives of An external file that holds a picture, illustration, etc.
Object name is nihms108120ig2.jpg

(·)with respect to and θ and γ₁ We then have

Explicit expressions for An external file that holds a picture, illustration, etc.
Object name is nihms108120ig2.jpg

_γ_₁ and

_θ; are given in Appendix C. Also, the information matrix is given by

where

_ΩΩ = E(−

_ΩΩ), with

_ΩΩ the second derivative of

with respect to Ω; note that the terms involving second derivatives of Ω^* do not appear in the information matrix because E(

_Ω) = 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

(3.4)

where we have substituted an estimate for F in An external file that holds a picture, illustration, etc.
Object name is nihms108120ig2.jpg

(·); that is, for each fixed value of (θ, γ₁), we solve γ₀ from

(3.5)

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

for the case where pr(D = 1) is known, where _emp,1(x) and _emp,0(x) are the empirical distributions of X in the case and control samples, and is given by _emp(x) = _emp,0(x) for the case where the rare-disease assumption can be made. An alternative choice of (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 _emp.

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 G_i denote the observed unphased genotype of subject i and An external file that holds a picture, illustration, etc.
Object name is nihms108120ig1.jpg

(G_i) the set of diplotypes that are consistent with G_i. When only G_i instead of equation M24

is observed for each subject, we propose to obtain the estimate [Phi w/ hat]

for Φ = (

, γ₁, θ) as the solution of the weighted version of (3.4):

(3.6)

where using the short-hand notation that ₀ = An external file that holds a picture, illustration, etc.
Object name is nihms108120ig2.jpg

(θ, γ₁, _emp), the weights are given by

(3.7)

The limiting version of the weights is given as

(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 γ₀, we

calculate the weights {_i} from (3.7) and
solve (3.6) to obtain an updated estimate of Φ using the weights {_i} given in (i); note that within this step we also need to solve (3.5) to obtain updated value of γ₀.

The algorithm is iterated between the 2 steps until convergence. Note that the weights {_i} are only used in solving Φ from (3.6) and are not required in solving γ₀ from (3.5).

3.3 Asymptotic theory

Make the following series of definitions. Expectations denoted as E_cc(·) are taken under the case–control sampling design, that is, for any random vector Y, equation M28

, where μ_d = plim n_d/n, d = 0, 1. Then, define

(3.9)

Note that the second derivative of Ω^* does not appear in

since E(

_Φ) = E(

_Φ) = 0, and the last identity in (3.9) is given by Lemma A.3 in Appendix A.

Define _emp(D_i) to be the mass of _emp(X_i), which is equal to equation M30

if π_d = pr (D = d) is known and is equal to I (D_i = 0)/n₀ when the rare-disease approximation is used. Let q_hap(X, γ₁, γ₀) = {q_hap(h|X, γ₁, γ₀)} be the vector collection over h of q_hap(h|X, γ₁, γ₀) for all diplotypes except the reference diplotype, and let q_HWE(θ) be defined similarly. Define

where

_Ω_γ_₀ is the obvious submatrix of x2110

_ΩΩ.

Theorem 3.1

Let

Suppose that E_cc{

(·)

^T (·)} exists and the matrix

is invertible. Then, n^1/2( [Phi w/ hat]

− Φ) is asymptotically normal with mean zero and covariance matrix

(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 _emp 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

⁻¹ −

⁻¹Ψ

⁻¹, where

. 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 γ₁_j_₁_j_₂= γ_1, _j_₁ + γ_1, _j_₂, where j₁ and j₂ 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 γ₀_j_₁_j_₂ ’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

where Z(5) denotes the number of the copies of the causal haplotype contained in H^di. The true value of the parameter vector (β₀, β_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 γ₁, 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 γ₁ 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 (more ...)

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 [H^di|X] is misspecified.

Following the argument of causality in Section 2, if [X|H^di] 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

where the diplotype data are again generated from the distribution in Table 1, ε is a t-distribution (d.f. = 3) truncated at ±5, λ₄ = λ₅ = −1.2, λ₁₂ = 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, [H^di|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(H^di|X) = pr(H^di) = q_HWE(H^di), 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 [H^di|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|H^di] have nonzero mean; for example, the estimate for the interaction parameter between h₅ 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 [H^di|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 (more ...)

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 β₁, the haplotype–environment odds ratio parameters γ₁, 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 γ₁, 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 γ₁ 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 H^di and environmental factors W = (X, S) is given by

(6.1)

where the stratum-specific intercept parameters γ₀_j_₁_j_₂ (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 = n₁ + n₀ tends to infinity but the sampling proportions for the cases and controls, that is, n₁/n and n₀/n, remain fixed at μ₁ and μ₀, we have for any measurable function M(D, H^di, X) of data (D, H^di, X),

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

(A.1)

and λ(x) = Σ_d Σ_h, S(d, h, x,

, γ₁, γ₀)

Lemma A.2 below provides an explicit expression for the estimating function

_Φ(·).

Lemma A.2

Write

Then

(A.2)

Proof

By definition,

and direct calculation yields

which proves the result.

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

Lemma A.3

Let

_ΩΩ = E_cc{−

_ΩΩ(D, H^di, X, Ω)}, where

_ΩΩ is the second derivative of

(·) with respective to Ω and x2110

_ΩΩ = E_cc{− [partial differential]

_Ω(·)/

Ω^T}. Then

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

(A.3)

The first term of (A.3) equals

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

where the joint density p^*(D, H^di|X, Ω) is defined in (A.1). Recalling further that

_Ω(D, h, X, Ω) = [partial differential]

log p^*(D, H^di = h|X, Ω)/ [partial differential]

Ω and

_Ω (D, G, X, Ω) = Σ_h ω(h, Ω)

_Ω(D, h, X, Ω), we have the identity

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

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]

= (

, ₁, ) a standard Taylor series expansion of (3.6) yields

(B.1)

In view of (2.5) and (3.5), let γ₀ = An external file that holds a picture, illustration, etc.
Object name is nihms108120ig2.jpg

(θ, γ₁, F) and ₀ = An external file that holds a picture, illustration, etc.
Object name is nihms108120ig2.jpg

(θ, γ₁,

_emp). Recall that ₀ is solved from (3.5):

Making a further Taylor series expansion, we have

Note that

Hence,

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

Therefore, we have

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

and

, 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

(θ, γ₁, F) AND An external file that holds a picture, illustration, etc.
Object name is nihms108120ig2.jpg

(θ, γ₁, _emp)

Recall that An external file that holds a picture, illustration, etc.
Object name is nihms108120ig2.jpg

(θ, γ₁, F)is the solution of γ₀ to (2.5). Here, we define γ₁ 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 {γ₁_j_₁_j_₂}. Differentiating both sides of (2.5) with respect to γ₁, we have

Let

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

Similarly, letting

we have

The derivatives of An external file that holds a picture, illustration, etc.
Object name is nihms108120ig2.jpg

(θ, γ₁, _emp) can be obtained by replacing dF in the above quantities with d_emp.

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, Email: 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]