Home | About | Journals | Submit | Contact Us | Français |

**|**HHS Author Manuscripts**|**PMC2683243

Formats

Article sections

- Summary
- 1. Introduction
- 2. Notations and proposed model
- 3. Semiparametric estimating equation inference
- 4. Simulations
- 5. Case–control study of colorectal adenoma study, NAT2 haplotype, and smoking
- 6. Concluding remarks
- References

Authors

Related links

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

Published in final edited form as:

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

PMCID: PMC2683243

NIHMSID: NIHMS108120

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

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

See other articles in PMC that cite the published article.

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.

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.

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*_{1}, *H*_{2}), where *H*_{1} and *H*_{2} denote the constituent haplotypes. We assume that there are *J* possible haplotypes indexed by *h _{j}* for

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

$$\text{logit}\{\text{pr}(D=1\mid {H}^{\text{di}},X)\}={\beta}_{0}+m({H}^{\text{di}},X,{\beta}_{1}),$$

(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

$$m\{{H}^{\text{di}}=({h}_{{j}_{1}},{h}_{{j}_{2}}),X,{\beta}_{1}\}={\beta}_{X}X+{\beta}_{{h}_{{j}_{1}}}+{\beta}_{{h}_{{j}_{2}}}+{\beta}_{{h}_{{j}_{1}}:X}X+{\beta}_{{h}_{{j}_{2}}:X}X,$$

(2.2)

where *β _{X}* is the main effect of

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*:

$$log\left[\frac{\text{pr}\{{H}^{\text{di}}=({h}_{{j}_{1}},{h}_{{j}_{2}})\mid X\}}{\text{pr}\left\{{H}^{\text{di}}=\left({h}_{{j}_{1}^{\ast}},{h}_{{j}_{2}^{\ast}}\right)\mid X\right\}}\right]={\gamma}_{0{j}_{1}{j}_{2}}+{\gamma}_{1{j}_{1}{j}_{2}}X,$$

(2.3)

where
${h}_{0}^{\text{di}}=({h}_{{j}_{1}^{\ast}},{h}_{{j}_{2}^{\ast}})$ is a chosen reference diplotype. Observe that model (2.3) allows association between *H*^{di} and *X* through the regression parameters *γ*_{1}_{j}_{1}_{j}_{2}. Let *γ*_{0} and *γ*_{1} denote the vectorized forms for the parameters *γ*_{0} _{j}_{1}_{j}_{2} and *γ*_{1}_{j}_{1}_{j}_{2}. Let *q*hap(*h*^{di}|*x*, *γ*_{0}, *γ*_{1}) 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 *γ*_{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} _{j}_{1}_{j}_{2} = *γ*_{1}, _{j}_{1} + *γ*_{1}, _{j}_{2}, 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

$$\text{pr}(x\mid h)=\frac{{q}_{\text{hap}}({h}^{\text{di}},x,{\gamma}_{0},{\gamma}_{1})\text{d}F(x)}{\int {q}_{\text{hap}}({h}^{\text{di}},v,{\gamma}_{0},{\gamma}_{1})\text{d}F(v)}.$$

This class of semiparametric models includes the parametric submodel where *X*|*H*^{di} = *h*^{di} follows a multivariate normal distribution with mean *μ _{h}*di and common variance–covariance matrix Σ. In this case, it is easy to see that
${\gamma}_{1{h}^{\text{di}}}={({\mu}_{{h}^{\text{di}}}-{\mu}_{{h}_{0}^{\text{di}}})}^{\text{T}}{\mathrm{\sum}}^{-1}$, which is a measure of the shift in the mean of the distribution of

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

$$log\left[\frac{\text{pr}\{{H}^{\text{di}}={h}^{\text{di}}=({h}_{{j}_{1}},{h}_{{j}_{2}})\}}{\text{pr}\{{H}^{\text{di}}=({h}_{1},{h}_{1})\}}\right]={\theta}_{{j}_{1}}+{\theta}_{{j}_{2}}\equiv {\theta}_{{h}^{\text{di}}},$$

(2.4)

where *h*_{1} denotes the chosen reference haplotype and *θ*_{1} = 0. Let

$${q}_{\text{HWE}}({h}^{\text{di}},\theta )=\text{pr}({H}^{\text{di}}={h}^{\text{di}},\theta )=\frac{exp({\theta}_{{h}^{\text{di}}})}{1+{\sum}_{{h}_{\ast}^{\text{di}}}exp({\theta}_{{h}_{\ast}^{\text{di}}})}$$

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

$${q}_{\text{HWE}}({h}^{\text{di}},\theta )=\int {q}_{\text{hap}}({H}^{\text{di}}={h}^{\text{di}}\mid X,{\gamma}_{0},{\gamma}_{1})\text{d}F(X).$$

(2.5)

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

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

$${p}_{\text{risk}}({H}^{\text{di}},D,X,{\beta}_{0},{\beta}_{1})={[\mathcal{H}\{{\beta}_{0}+m({H}^{\text{di}},X,{\beta}_{1})\}]}^{D}{[1-\mathcal{H}\{{\beta}_{0}+m({H}^{\text{di}},X,{\beta}_{1})\}]}^{1-D}.$$

Recall that *q*_{hap}(*h*|*X*, *γ*_{0}, *γ*_{1}) = pr(*H*^{di} = *h*|*X; γ*_{0}, *γ*_{1}) 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

$$\begin{array}{l}l=\sum _{d=0}^{1}\sum _{k=1}^{K}\sum _{h}{n}_{\mathit{dkh}}log\{{p}_{\text{risk}}(h,d,{x}_{k},{\beta}_{0},{\beta}_{1}){q}_{\text{hap}}(h\mid {x}_{k},{\gamma}_{0},{\gamma}_{1}){\delta}_{k}\}\\ -\sum _{d=0}^{1}{n}_{d}log\left\{\sum _{k=1}^{K}\sum _{h}{p}_{\text{risk}}(h,d,{x}_{k},{\beta}_{0},{\beta}_{1}){q}_{\text{hap}}(h\mid {x}_{k},{\gamma}_{0},{\gamma}_{1}){\delta}_{k}\right\}.\end{array}$$

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

$${\delta}_{k}=\frac{{\sum}_{i=1}^{n}I({X}_{i}={x}_{k})}{{\sum}_{d}{\sum}_{h}({n}_{d}/{\pi}_{d}){p}_{\text{risk}}(h,d,{x}_{k},{\beta}_{0},{\beta}_{1}){q}_{\text{hap}}(h\mid {x}_{k},{\gamma}_{1},{\gamma}_{0})}$$

(3.1)

and the profile log-likelihood

$$l(\omega ,\widehat{\delta}(\omega ))=\sum _{i=1}^{n}\mathcal{L}({D}_{i},{H}_{i}^{\text{di}},{X}_{i},\mathcal{B},{\gamma}_{1},{\gamma}_{0})=\sum _{i=1}^{n}log\left\{\frac{S({D}_{i},{H}_{i}^{\text{di}},{X}_{i},\mathcal{B},{\gamma}_{1},{\gamma}_{0})}{{\sum}_{{d}_{\ast}=0}^{1}{\sum}_{{h}_{\ast}}S({d}_{\ast},{h}_{\ast},{X}_{i},\mathcal{B},{\gamma}_{1},{\gamma}_{0})}\right\},$$

(3.2)

where

$${\pi}_{d}=\text{pr}(D=d)=\sum _{h,k}{p}_{\text{risk}}(h,d,{x}_{k},{\beta}_{0},{\beta}_{1}){q}_{\text{hap}}(h\mid {x}_{k},{\gamma}_{1},{\gamma}_{0}){\delta}_{k}$$

and

$$S(d,h,x,\mathcal{B},{\gamma}_{1},{\gamma}_{0})=exp\{d(\kappa -{\beta}_{0})\}{p}_{\text{risk}}(h,d,x,{\beta}_{0},{\beta}_{1}){q}_{\text{hap}}(h\mid x,{\gamma}_{0},{\gamma}_{1}),$$

with = (*β*_{0}, *β*_{1}, *κ*)^{T} and *κ* = *β*_{0} + log(*n*_{1}/*n*_{0}) − 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 and *S*(*d*, *h*, *x*, , *γ*_{0}, *γ*_{1}). When pr(*D* = 1) is known, *κ* depends on *β*_{0} only, hence here we define = (*β*_{0}, *β*_{1})^{T}. When the disease is rare so that

$${p}_{\text{risk}}(h,d,x,{\beta}_{0},{\beta}_{1})\approx exp[d\{{\beta}_{0}+m(h,x,{\beta}_{1})\}],$$

(3.3)

we have

$$S(d,h,x,\mathcal{B},{\gamma}_{1},{\gamma}_{0})=exp[d\{\kappa +m(h,x,{\beta}_{1})\}]{q}_{\text{hap}}(h\mid x,{\gamma}_{0},{\gamma}_{1}).$$

Note that *β*_{0} does not appear in this expression, and hence we define = (*κ*, *β*_{1})^{T}.

Our goal is to estimate the parameters (, *θ*, *γ*_{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} =
(*θ*, *γ*_{1}, *F*). Let Ω = (, *γ*_{1}, *γ*_{0}), Ω^{*} = {, *γ*_{1},
(*θ*, *γ*_{1}, *F*)}, and Φ = (, *γ*_{1}, *θ*). Let _{Ω}(·) and _{Φ}(·) be, respectively, the derivatives of (·) with respect to Ω and Φ, and
* _{θ}* and

$${\mathcal{L}}_{\mathrm{\Phi}}(\xb7)=\frac{\partial {\mathrm{\Omega}}^{\ast}}{\partial \mathrm{\Phi}}{\mathcal{L}}_{\mathrm{\Omega}}(\xb7)={\{{\mathcal{L}}_{\mathcal{B}}^{\text{T}}(\xb7),{({\mathcal{L}}_{{\gamma}_{1}}+{\mathcal{G}}_{{\gamma}_{1}}{\mathcal{L}}_{{\gamma}_{0}})}^{\text{T}}(\xb7),{({\mathcal{G}}_{\theta}{\mathcal{L}}_{{\gamma}_{0}})}^{\text{T}}(\xb7)\}}^{\text{T}}.$$

Explicit expressions for
_{γ}_{1} and
* _{θ;}* are given in Appendix C. Also, the information matrix is given by

$$\mathcal{I}=\frac{\partial {\mathrm{\Omega}}^{\ast}}{\partial \mathrm{\Phi}}{\mathcal{I}}_{\mathrm{\Omega}\mathrm{\Omega}}{\left(\frac{\partial {\mathrm{\Omega}}^{\ast}}{\partial \mathrm{\Phi}}\right)}^{\text{T}},$$

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

$$\begin{array}{l}0={n}^{-1/2}\sum _{i=1}^{n}{\mathcal{L}}_{\mathrm{\Phi}}\{{D}_{i},{H}_{i}^{\text{di}},{X}_{i},\mathcal{B},{\gamma}_{1},\mathcal{G}(\theta ,{\gamma}_{1},\widehat{F})\}\\ \equiv {n}^{-1/2}\sum _{i=1}^{n}{\mathcal{L}}_{\mathrm{\Phi}}({D}_{i},{H}_{i}^{\text{di}},{X}_{i},\mathrm{\Phi}),\end{array}$$

(3.4)

where we have substituted an estimate for *F* in (·); that is, for each fixed value of (*θ*, *γ*_{1}), we solve *γ*_{0} from

$${q}_{\text{HWE}}(h,\theta )=\int {q}_{\text{hap}}(h\mid x,{\gamma}_{0},{\gamma}_{1})\text{d}\widehat{F}(x).$$

(3.5)

One convenient choice of is the empirical estimate _{emp}, which is given by

$${\widehat{F}}_{\text{emp}}(x)={\widehat{F}}_{\text{emp},1}(x)\text{pr}(D=1)+{\widehat{F}}_{\text{emp},0}(x)\text{pr}(D=0)$$

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}.

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

$$\begin{array}{l}0={n}^{-1/2}\sum _{i=1}^{n}\sum _{h}{\widehat{\omega}}_{i}\{h,\mathcal{B},{\gamma}_{1},\mathcal{G}(\theta ,{\gamma}_{1},{\widehat{F}}_{\text{emp}})\}{\mathcal{L}}_{\mathrm{\Phi}}\{{D}_{i},h,{X}_{i},\mathcal{B},{\gamma}_{1},\mathcal{G}(\theta ,{\gamma}_{1},{\widehat{F}}_{\text{emp}})\}\\ \equiv {n}^{-1/2}\sum _{i=1}^{n}{\overline{\mathcal{L}}}_{\mathrm{\Phi}}\{{D}_{i},{G}_{i},{X}_{i},\mathrm{\Phi},\mathcal{G}(\theta ,{\gamma}_{1},{\widehat{F}}_{\text{emp}})\},\end{array}$$

(3.6)

where using the short-hand notation that _{0} = (*θ*, *γ*_{1}, _{emp}), the weights are given by

$$\begin{array}{l}{\widehat{\omega}}_{i}\{h,\mathcal{B},{\gamma}_{1},\mathcal{G}(\theta ,{\gamma}_{1},{\widehat{F}}_{\text{emp}})\}=\text{pr}({H}^{\text{di}}=h\mid {D}_{i},{X}_{i},{G}_{i},\mathcal{B},{\gamma}_{1},{\widehat{\gamma}}_{0})\\ =I\{h\in \mathcal{C}({G}_{i})\}\frac{{p}_{\text{risk}}(h,{D}_{i},{X}_{i},\mathcal{B}){q}_{\text{hap}}(h\mid {X}_{i},{\gamma}_{1},{\widehat{\gamma}}_{0})}{{\sum}_{{h}_{\ast}\in \mathcal{C}({G}_{i})}{p}_{\text{risk}}({h}_{\ast},{D}_{i},{X}_{i},\mathcal{B}){q}_{\text{hap}}({h}_{\ast}\mid {X}_{i},{\gamma}_{1},{\widehat{\gamma}}_{0})}.\end{array}$$

(3.7)

The limiting version of the weights is given as

$$\omega (h,\mathrm{\Omega})=\frac{I\{h\in C(G)\}S(D,h,X,\mathcal{B},{\gamma}_{1},{\gamma}_{0})}{{\sum}_{{h}_{\ast}\in C(G)}S(D,{h}_{\ast},X,\mathcal{B},{\gamma}_{1},{\gamma}_{0})}=\text{pr}({H}^{\text{di}}=h\mid D,G,X,\mathrm{\Omega}).$$

(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

- calculate the weights {
} from (3.7) and_{i}

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

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*,
${E}_{\text{cc}}(Y)={\sum}_{d=0}^{1}{\mu}_{d}E(Y\mid D=d)$, where *μ _{d}* = plim

$$\begin{array}{l}\overline{\mathcal{I}}={E}_{\text{cc}}\left\{-\frac{\partial}{\partial {\mathrm{\Phi}}^{\text{T}}}{\overline{\mathcal{L}}}_{\mathrm{\Phi}}(D,G,X,\mathrm{\Phi})\right\}=\frac{\partial {\mathrm{\Omega}}^{\ast}}{\partial \mathrm{\Phi}}{\overline{\mathcal{I}}}_{\mathrm{\Omega}\mathrm{\Omega}}{\left(\frac{\partial {\mathrm{\Omega}}^{\ast}}{\partial \mathrm{\Phi}}\right)}^{\text{T}},\\ {\overline{\mathcal{L}}}_{\mathrm{\Omega}}(D,G,X,\mathrm{\Omega})=\sum _{h}\omega (h,\mathrm{\Omega}){\mathcal{L}}_{\mathrm{\Omega}}(D,h,X,\mathrm{\Omega}),\\ {\overline{\mathcal{I}}}_{\mathrm{\Omega}\mathrm{\Omega}}={E}_{\text{cc}}\{-\partial {\overline{\mathcal{L}}}_{\mathrm{\Omega}}(\xb7)/\partial {\mathrm{\Omega}}^{T}\}={E}_{\text{cc}}\{{\overline{\mathcal{L}}}_{\mathrm{\Omega}}(D,G,X,\mathrm{\Omega}){\overline{\mathcal{L}}}_{\mathrm{\Omega}}^{\text{T}}(D,G,X,\mathrm{\Omega})\}.\end{array}$$

(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

$$\begin{array}{l}{Q}_{n}=\sum _{i=1}^{n}\frac{\partial}{\partial {\gamma}_{0}}{\mathbf{q}}_{\text{hap}}({X}_{i},{\gamma}_{0},{\gamma}_{1}){\widehat{p}}_{\text{emp}}({D}_{i}),\\ Q=\int \frac{\partial}{\partial {\gamma}_{0}}{\mathbf{q}}_{\text{hap}}(x,{\gamma}_{0},{\gamma}_{1})\text{d}F(x)=E\left\{\frac{\partial}{\partial {\gamma}_{0}}{\mathbf{q}}_{\text{hap}}(X,{\gamma}_{0},{\gamma}_{1})\right\},\\ k(X,D)=n{Q}^{-1}\{{\mathbf{q}}_{\text{hap}}(X,{\gamma}_{1},{\gamma}_{0})-{\mathbf{q}}_{\text{HWE}}(\theta )\}{\widehat{p}}_{\text{emp}}({D}_{i}),\\ {\overline{\mathcal{I}}}_{\mathrm{\Phi}{\gamma}_{0}}=(\partial {\mathrm{\Omega}}^{\ast}/\partial \mathrm{\Phi}){\overline{\mathcal{I}}}_{\mathrm{\Phi}{\gamma}_{0}},\\ {\overline{\mathcal{I}}}_{\mathrm{\Phi}{\gamma}_{0}}=E(-\partial {\overline{\mathcal{L}}}_{\mathrm{\Omega}}/\partial {{\gamma}_{0}}^{\text{T}}),\\ \mathcal{K}(X,D)={\overline{\mathcal{I}}}_{\mathrm{\Phi}{\gamma}_{0}}k(X,D),\end{array}$$

where _{Ω}_{γ}_{0} is the obvious submatrix of _{ΩΩ}.

Let

$$\overline{\mathcal{E}}(D,G,X,\mathrm{\Phi})={\overline{\mathcal{L}}}_{\mathrm{\Phi}}(D,G,X,\mathrm{\Phi})-E\{{\overline{\mathcal{L}}}_{\mathrm{\Phi}}(D,G,X,\mathrm{\Phi})\mid D\}-\{\mathcal{K}(X,D)-E(\mathcal{K}(X,D)\mid D)\}.$$

Suppose that *E*_{cc}{(·) ^{T} (·)} exists and the matrix is invertible. Then, *n*^{1/2}( − Φ) is asymptotically normal with mean zero and covariance matrix

$$\overline{\mathrm{\Gamma}}={\overline{\mathcal{I}}}^{-1}{E}_{\text{cc}}(\overline{\mathcal{E}}{\overline{\mathcal{E}}}^{\text{T}}){\overline{\mathcal{I}}}^{-1}.$$

(3.10)

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

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 (·) 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 ^{−1} − ^{−1}Ψ^{−1}, where
$\mathrm{\Psi}={\sum}_{d=0}^{1}({n}_{d}/n){[E\{{\overline{\mathcal{L}}}_{\mathrm{\Phi}}(D,G,X,\mathrm{\Phi})\mid D=d\}]}^{\otimes 2}$. Whether this naive estimate performs well in general is unknown, and we suggest using the estimate based on (3.10).

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 *γ*_{1}_{j}_{1}_{j}_{2}= *γ*_{1,} _{j}_{1} + *γ*_{1,} _{j}_{2}, where *j*_{1} and *j*_{2} 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

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

$$\text{logit}\{\text{pr}(D=1\mid {H}^{\text{di}},X)\}={\beta}_{0}+{\beta}_{H}Z(5)+{\beta}_{X}X+{\beta}_{HX}Z(5)X,$$

where *Z*(5) denotes the number of the copies of the causal haplotype contained in *H*^{di}. The true value of the parameter vector (*β*_{0}, *β _{H}*,

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}*

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.

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

$$\left[X\mid {H}^{\text{di}}=({h}_{{j}_{1}},{h}_{{j}_{2}})\right]\sim {\lambda}_{{j}_{1}}+{\lambda}_{{j}_{2}}+\epsilon ,$$

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, [

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*_{5} 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.

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.

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.

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

$$log\left[\frac{\text{pr}\{{H}^{\text{di}}=({h}_{{j}_{1}},{h}_{{j}_{2}})\mid X,S\}}{\text{pr}\left\{{H}^{\text{di}}=\left({h}_{{j}_{1}^{\ast}},{h}_{{j}_{2}^{\ast}}\right)\mid X,S\right\}}\right]={\gamma}_{0{j}_{1}{j}_{2}}(S)+{\gamma}_{1{j}_{1}{j}_{2}}X,$$

(6.1)

where the stratum-specific intercept parameters *γ*_{0}_{j}_{1}_{j}_{2} (*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.

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.

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).

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

$$\begin{array}{l}{E}_{\text{cc}}\{M(D,{H}^{\text{di}},X)\}=\sum _{d=0}^{1}{\mu}_{d}E\{M(D,{H}^{\text{di}},X)\mid D=d\}\\ =\int {E}^{\ast}\{M(D,{H}^{\text{di}},X)\mid X=x\}\lambda (x)\text{d}F(x),\end{array}$$

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

$${p}^{\ast}({H}^{\text{di}},D\mid X,\mathcal{B},{\gamma}_{1},{\gamma}_{0})=\frac{S(D,{H}^{\text{di}},X,\mathcal{B},{\gamma}_{1},{\gamma}_{0})}{{\sum}_{d}{\sum}_{h}S(d,h,X,\mathcal{B},{\gamma}_{1},{\gamma}_{0})}$$

(A.1)

and *λ*(*x*) = Σ* _{d}* Σ

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

Write

$$\begin{array}{l}S(D,{H}^{\text{di}},X,\mathrm{\Phi})=S\{D,{H}^{\text{di}},X,\mathcal{B},{\gamma}_{1},\mathcal{G}(\theta ,{\gamma}_{1},F)\},\\ {p}^{\ast}(D,{H}^{\text{di}}\mid X,\mathrm{\Phi})=\frac{S(D,{H}^{\text{di}},X,\mathrm{\Phi})}{{\sum}_{{d}_{\ast}}{\sum}_{{h}_{\ast}}S({d}_{\ast},{h}_{\ast},X,\mathrm{\Phi})},\\ \mathcal{R}(D,{H}^{\text{di}},X,\mathrm{\Phi})=\frac{\partial}{\partial \mathrm{\Phi}}logS(D,{H}^{\text{di}},X,\mathrm{\Phi}).\end{array}$$

Then

$${\mathcal{L}}_{\mathrm{\Phi}}\{D,{H}^{\text{di}},X,\mathcal{B},{\gamma}_{1},\mathcal{G}(\theta ,{\gamma}_{1},F)\}=\mathcal{R}(D,{H}^{\text{di}},X,\mathrm{\Phi})-{E}^{\ast}\{\mathcal{R}(\xb7)\mid X\}.$$

(A.2)

By definition,

$${\mathcal{L}}_{\mathrm{\Phi}}\{D,{H}^{\text{di}},X,\mathcal{B},{\gamma}_{1},\mathcal{G}(\theta ,{\gamma}_{1},F)\}=\frac{\partial}{\partial \mathrm{\Phi}}log{p}^{\ast}(D,{H}^{\text{di}}\mid X,\mathrm{\Phi}),$$

and direct calculation yields

$$\frac{\partial}{\partial \mathrm{\Phi}}log{p}^{\ast}(D,{H}^{\text{di}}\mid X,\mathrm{\Phi})=\mathcal{R}(D,{H}^{\text{di}},X,\mathrm{\Phi})-\sum _{d}\sum _{h}\frac{\mathcal{R}(d,h,X,\mathrm{\Phi})S(d,h,X,\mathrm{\Phi})}{{\sum}_{{d}_{\ast}}{\sum}_{{h}_{\ast}}S({d}_{\ast},{h}_{\ast},X,\mathrm{\Phi})},$$

which proves the result.

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

Let _{ΩΩ} = *E*_{cc}{−_{ΩΩ}(*D, H*^{di}, *X*, Ω)}, where _{ΩΩ} is the second derivative of (·) with respective to Ω and _{ΩΩ} = E_{cc}{−_{Ω}(·)/Ω^{T}}. Then

$$\begin{array}{l}{\mathcal{I}}_{\mathrm{\Omega}\mathrm{\Omega}}={E}_{\text{cc}}\{{\mathcal{L}}_{\mathrm{\Omega}}(D,{H}^{\text{di}},X,\mathrm{\Omega}){\mathcal{L}}_{\mathrm{\Omega}}^{\text{T}}(D,{H}^{\text{di}},X,\mathrm{\Omega})\},\\ {\overline{\mathcal{I}}}_{\mathrm{\Omega}\mathrm{\Omega}}={E}_{\text{cc}}\{{\overline{\mathcal{L}}}_{\mathrm{\Omega}}(D,G,X,\mathrm{\Omega}){\overline{\mathcal{L}}}_{\mathrm{\Omega}}^{\text{T}}(D,G,X,\mathrm{\Omega})\}.\end{array}$$

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

$$\begin{array}{l}{E}_{\text{cc}}\{-{\overline{\mathcal{L}}}_{\mathrm{\Omega}\mathrm{\Omega}}(D,G,X,\mathrm{\Omega})={E}_{\text{cc}}[E\{-{\mathcal{L}}_{\mathrm{\Omega}\mathrm{\Omega}}(D,{H}^{\text{di}},X,\mathrm{\Omega})\mid D,G,X\}]\\ -{E}_{\text{cc}}\left\{\sum _{h\in \mathcal{C}(G)}{\mathcal{L}}_{\mathrm{\Omega}}(D,h,X,\mathrm{\Omega})\frac{\partial \omega (h,\mathrm{\Omega})}{\partial {\mathrm{\Omega}}^{\text{T}}}\right\}.\end{array}$$

(A.3)

The first term of (A.3) equals

$${E}_{\text{cc}}\{-{\mathcal{L}}_{\mathrm{\Omega}\mathrm{\Omega}}(D,{H}^{\text{di}},X,\mathrm{\Omega})\}={E}_{\text{cc}}\{{\mathcal{L}}_{\mathrm{\Omega}}(D,{H}^{\text{di}},X){\mathcal{L}}_{\mathrm{\Omega}}^{\text{T}}(D,{H}^{\text{di}},X)\}.$$

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

$$\omega (h,\mathrm{\Omega})={p}^{\ast}({H}^{\text{di}}=h\mid D,X,\mathrm{\Omega}),$$

where the joint density *p*^{*}(*D*, *H*^{di}|*X*, Ω) is defined in (A.1). Recalling further that _{Ω}(*D, h, X*, Ω) = log *p*^{*}(*D*, *H*^{di} = *h*|*X*, Ω)/Ω and _{Ω} (*D, G, X*, Ω) = Σ* _{h} ω*(

$$\frac{\partial \omega (h,\mathrm{\Omega})}{\partial \mathrm{\Omega}}={\mathcal{L}}_{\mathrm{\Omega}}(D,h,X,\mathrm{\Omega})\omega (h,\mathrm{\Omega})-{\overline{\mathcal{L}}}_{\mathrm{\Omega}}(D,G,X,\mathrm{\Omega})\omega (h,\mathrm{\Omega}).$$

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

$${E}_{\text{cc}}\{{\mathcal{L}}_{\mathrm{\Omega}}(D,{H}^{\text{di}},X,\mathrm{\Omega}){\mathcal{L}}_{\mathrm{\Omega}}^{\text{T}}(D,{H}^{\text{di}},X,\mathrm{\Omega})\}-{E}_{\text{cc}}\{{\overline{\mathcal{L}}}_{\mathrm{\Omega}}(D,G,X,\mathrm{\Omega}){\overline{\mathcal{L}}}_{\mathrm{\Omega}}^{\text{T}}(D,G,X,\mathrm{\Omega})\}.$$

The desired result thus follows.

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

For our estimator = (, _{1}, ) a standard Taylor series expansion of (3.6) yields

$$\begin{array}{l}0={n}^{-1/2}\sum _{i=1}^{n}{\overline{\mathcal{L}}}_{\mathrm{\Phi}}\{{D}_{i},{G}_{i},{X}_{i},\mathcal{B},{\gamma}_{1},\mathcal{G}(\theta ,{\gamma}_{1},F)\}\\ +{n}^{-1/2}\sum _{i=1}^{n}[{\overline{\mathcal{L}}}_{\mathrm{\Phi}}\{{D}_{i},{G}_{i},{X}_{i},\mathcal{B},{\gamma}_{1},\mathcal{G}(\theta ,{\gamma}_{1},{\widehat{F}}_{\text{emp}})\}\\ -{\overline{\mathcal{L}}}_{\mathrm{\Phi}}\{{D}_{i},{G}_{i},{X}_{i},\mathcal{B},{\gamma}_{1},\mathcal{G}(\theta ,{\gamma}_{1},F)\}]\\ -\overline{\mathcal{I}}{n}^{1/2}(\widehat{\mathrm{\Phi}}-\mathrm{\Phi})+{o}_{p}(1).\end{array}$$

(B.1)

In view of (2.5) and (3.5), let *γ*_{0} =
(*θ, γ*_{1}, *F*) and _{0} = (*θ, γ*_{1}, _{emp}). Recall that _{0} is solved from (3.5):

$$0=\sum _{i=1}^{n}{\mathbf{q}}_{\text{hap}}({X}_{i},{\gamma}_{1},{\widehat{\gamma}}_{0}){\widehat{p}}_{\text{emp}}({D}_{i})-{\mathbf{q}}_{\text{HWE}}(\mathbf{h},\theta ).$$

Making a further Taylor series expansion, we have

$$0={n}^{1/2}\left\{\sum _{i=1}^{n}{\mathbf{q}}_{\text{hap}}({X}_{i},{\gamma}_{1},{\widehat{\gamma}}_{0}){\widehat{p}}_{\text{emp}}({D}_{i})-{\mathbf{q}}_{\text{HWE}}(\theta )\right\}+{Q}_{n}{n}^{1/2}({\widehat{\gamma}}_{0}-{\gamma}_{0})+{o}_{p}(1).$$

Note that

$$\begin{array}{l}{Q}_{n}=\sum _{i=1}^{n}\frac{\partial}{\partial {\gamma}_{0}}{\mathbf{q}}_{\text{hap}}({X}_{i},{\gamma}_{0},{\gamma}_{1}){\widehat{p}}_{\text{emp}}({D}_{i})=\int \frac{\partial}{\partial {\gamma}_{0}}{\mathbf{q}}_{\text{hap}}(X,{\gamma}_{0},{\gamma}_{1})\text{d}{\widehat{F}}_{\text{emp}}(X)\\ =Q+{o}_{p}(1).\end{array}$$

Hence,

$$\begin{array}{l}{n}^{1/2}({\widehat{\gamma}}_{0}-{\gamma}_{0})={n}^{1/2}{Q}^{-1}\sum _{i=1}^{n}\{{\mathbf{q}}_{\text{hap}}({X}_{i},{\gamma}_{1},{\gamma}_{0})-{\mathbf{q}}_{\text{HWE}}(\theta )\}{\widehat{p}}_{\text{emp}}({D}_{i})+{o}_{p}(1)\\ \equiv {n}^{-1/2}\sum _{i=1}^{n}k({X}_{i},{D}_{i})+{o}_{p}(1).\end{array}$$

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

$$\begin{array}{l}{n}^{-1/2}\sum _{i=1}^{n}[{\overline{\mathcal{L}}}_{\mathrm{\Phi}}\{{D}_{i},{G}_{i},{X}_{i},\mathcal{B},{\gamma}_{1},\mathcal{G}(\theta ,{\gamma}_{1},{\widehat{F}}_{\text{emp}})\}-{\overline{\mathcal{L}}}_{\mathrm{\Phi}}\{{D}_{i},{G}_{i},{X}_{i},\mathcal{B},{\gamma}_{1},\mathcal{G}(\theta ,{\gamma}_{1},F)\}]\\ =-{\overline{\mathcal{I}}}_{\mathrm{\Phi}{\gamma}_{0}}{n}^{1/2}({\widehat{\gamma}}_{0}-{\gamma}_{0})+{o}_{p}(1)\\ =-{n}^{-1/2}\sum _{i=1}^{n}\mathcal{K}({X}_{i},{D}_{i})+{o}_{p}(1).\end{array}$$

Therefore, we have

$$\begin{array}{l}{n}^{1/2}(\widehat{\mathrm{\Phi}}-\mathrm{\Phi})={n}^{-1/2}{\overline{\mathcal{I}}}^{-1}\sum _{i=1}^{n}[{\overline{\mathcal{L}}}_{\mathrm{\Phi}}\{{D}_{i},{G}_{i},{X}_{i},\mathcal{B},{\gamma}_{1},\mathcal{G}(\theta ,{\gamma}_{1},F)\}-\mathcal{K}({X}_{i},{D}_{i})]+{o}_{p}(1)\\ ={n}^{-1/2}{\overline{\mathcal{I}}}^{-1}\sum _{i=1}^{n}{\overline{\mathcal{L}}}_{\mathrm{\Phi}}\{{D}_{i},{G}_{i},{X}_{i},\mathcal{B},{\gamma}_{1},\mathcal{G}(\theta ,{\gamma}_{1},F)\}\\ -E[{\overline{\mathcal{L}}}_{\mathrm{\Phi}}\{{D}_{i},G,X,\mathcal{B},{\gamma}_{1},\mathcal{G}(\theta ,{\gamma}_{1},F)\}\mid {D}_{i}]\\ -[\mathcal{K}({X}_{i},{D}_{i})-E\{\mathcal{K}(X,{D}_{i})\mid {D}_{i}\}]+{o}_{p}(1).\end{array}$$

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

$$\sum _{i=1}^{n}E[{\overline{\mathcal{L}}}_{\mathrm{\Phi}}\{{D}_{i},G,X,\mathcal{B},{\gamma}_{1},\mathcal{G}(\theta ,{\gamma}_{1},F)\}\mid {D}_{i}]=0$$

and ${\sum}_{i=1}^{n}E\{\mathcal{K}(X,{D}_{i})\mid {D}_{i}\}=0$, which follow directly from Lemmas A.1 and A.2. This completes the proof.

Recall that
(*θ, γ*_{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 {*γ*_{1}_{j}_{1}_{j}_{2}}. Differentiating both sides of (2.5) with respect to *γ*_{1}, we have

$$0=\left\{\int \frac{\partial {\mathbf{q}}_{\text{hap}}(x,{\gamma}_{0},{\gamma}_{1})}{\partial {\gamma}_{0}}\text{d}F(x)\right\}\frac{\partial {\gamma}_{0}}{\partial {\gamma}_{1}}+\int \frac{\partial {\mathbf{q}}_{\text{hap}}(x,{\gamma}_{0},{\gamma}_{1})}{\partial \eta}\left(\frac{\partial \eta}{\partial {{\gamma}_{1}}^{\text{T}}}\otimes {x}^{\text{T}}\right)\text{d}F(x).$$

Let

$$Q(x)=\text{diag}\{{\mathbf{q}}_{\text{hap}}(x,{\gamma}_{0},{\gamma}_{1})\}-{\mathbf{q}}_{\text{hap}}(x,{\gamma}_{0},{\gamma}_{1}){\mathbf{q}}_{\text{hap}}{(x,{\gamma}_{0},{\gamma}_{1})}^{\text{T}}$$

and *Q* = ∫ *Q*(*x*)d*F*(*x*). Then

$${\mathcal{G}}_{{\gamma}_{1}}=\frac{\partial {\gamma}_{0}}{\partial {\gamma}_{1}}=-{Q}^{-1}\int Q(x)\left(\frac{\partial \eta}{\partial {{\gamma}_{1}}^{\text{T}}}\otimes {x}^{\text{T}}\right)\text{d}F(x).$$

Similarly, letting

$$R=\text{diag}\{{\mathbf{q}}_{\text{HWE}}(\theta )\}-{\mathbf{q}}_{\text{HWE}}(\theta ){\mathbf{q}}_{\text{HWE}}{(\theta )}^{\text{T}},$$

we have

$${\mathcal{G}}_{\theta}=\frac{\partial {\gamma}_{0}}{\partial \theta}={Q}^{-1}R\frac{\partial {\theta}_{{h}^{\text{di}}}}{\partial {\theta}^{\text{T}}}.$$

The derivatives of (*θ, γ*_{1}, _{emp}) can be obtained by replacing d*F* in the above quantities with d_{emp}.

*Conflict of Interest:* None declared.

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: vog.hin.liam@nrettahc.

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

- 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]

PubMed Central Canada is a service of the Canadian Institutes of Health Research (CIHR) working in partnership with the National Research Council's national science library in cooperation with the National Center for Biotechnology Information at the U.S. National Library of Medicine(NCBI/NLM). It includes content provided to the PubMed Central International archive by participating publishers. |