An empirical Bayes method for genetic association analysis using case-control mother-child pair data

1.   Introduction

It is well known that most obstetric and early-life diseases or phenotypes have a multifactorial etiology involving genetic factors, environmental exposures, and interactions between them^[1–4]. For example, low maternal body mass index and genes of the mother and child were found to be associated with the risk of preterm delivery^{[5, 6]}. Moreover, maternal COL24A1 variants were shown to have significant genome-wide interactions with maternal pre-pregnancy overweight/ obesity on preterm birth risk using genetic data of 1733 African-American women from the Boston Birth Cohort^[7]. Identifying potential genetic risk factors is important for understanding their biological effects and for developing public health strategies for prevention. A popular design for the integrative investigation of the effects affecting obstetrical and early life phenotypes is to collect risk factor information from children and their mothers^{[1, 8, 9]}. In this study, the outcome of interest was a dichotomous phenotype of the mother or child, and cases and controls were ascertained based on this phenotype. This is the so-called mother-child case-control design^{[10, 11]}.

The standard prospective logistic regression method^[12] is adopted widely in case-control studies. However, this method is not efficient because it does not utilize information such as Mendelian inheritance between parental and child genotypes. Alternative methods have been developed to assess the genetic effects on both children and mothers. Shi et al.^[13] and Chen et al.^[8] proposed the log-linear method and retrospective method for fitting logistic regression models, respectively, and improved the statistical inference efficiency by incorporating constraints on the genotype distribution, such as Mendelian inheritance, random mating, and Hardy-Weinberg equilibrium (HWE) in the population of interest. Chen et al.^[10] extended the retrospective likelihood method to allow assessment of environmental effects and gene-environment interactions through a semiparametric maximum likelihood estimation (MLE) method, and improved statistical inference efficiency by incorporating an additional constraint that children’s genotype is conditionally independent of maternal environmental risk factors given maternal genotype.

Using the retrospective likelihood method to analyze case-control mother-child pair data can result in a noticeable efficiency gain in estimating odds ratio (OR) parameters under the model assumptions of HWE and maternal gene-environment independence in the population of interest; however, it would produce serious bias if the independence assumption or HWE is violated. Chen et al.^[10] alleviated the bias concern by allowing the conditional distribution of environmental risk factors for a given maternal genotype to be nonparametric, but this strategy is not sufficiently efficient. Alternatively, Zhang et al.^[14] proposed a modified profile likelihood estimation method by relating the maternal genotype to environmental risk factors through a novel “double-additive” logistic regression (daLOG) model, which is computationally robust but not sufficiently robust because estimation bias may be introduced when the daLOG model is seriously violated.

Two novel estimators are proposed in this study: the weighted averages of an efficient estimator and a robust estimator. The second estimator is robust in that it is nearly unbiased, and the corresponding significance test has Type I error controlled around the nominal level. Both estimators are closely related to the modified profile MLE of Zhang et al.^[14]. The efficient one imposes the assumptions of HWE and maternal gene-environment independence (referred to as “independent-model” estimator hereafter), which is biased if at least one of the assumptions is violated. In contrast, the robust one does not require the two assumptions, so it is relatively robust but might not be efficient (refer to as “dependent-model” estimator hereafter). Our weighting method was motivated by an empirical Bayes-type (EB-type) shrinkage estimator for estimating the effects of haplotypes and haplotype-environmental interactions with standard case-control data^[15]. Here, such shrinkage technique is extended to assess OR parameters associated with maternal genotype and environmental risk factors using case-control mother-child genotype data. Briefly, the robust “dependent-model” estimator is shrunk toward the more efficient but possibly biased “independent-model” estimator. The magnitude of “shrinkage” is data-dependent, which tends to be large if and only if HWE and maternal gene-environmental independence are satisfied.

The remainder of this paper is organized as follows. In Section 2, two MLEs are derived based on dependent and independent models. Subsequently, two EB-type estimators are constructed based on the two MLEs, and their asymptotic properties are established. In Sections 3 and 4, the desired finite-sample performance of the proposed EB-type estimators is demonstrated through extensive simulations and real-data applications. Concluding remarks are provided in Section 5.

2.   Materials and methods

2.1   Model, likelihood, and assumptions

Let $Y$ denote the dichotomous disease status of either mother (e.g., pre-eclampsia) or child (e.g., low birth weight), that is, $Y=0$ for a control pair and $Y=1$ for a case pair, $(G^{m},G^{c})$ denote the genotypes of a mother-child pair at a target genetic locus, and $X$ denotes a vector of $p$ maternal environmental risk factors. $Y$ and $(G^c, G^m, X)$ are related using the following logistic penetrance model:

where $\text{logit}(t)=\log(t/(1-t))$ is the logit function, $H(\cdot)$ is an arbitrarily specified function, and ${\boldsymbol{\beta}}$ is the $q$-vector of the regression parameters. $G^{c}$ and $G^{m}$ can be coded according to the mode of inheritance. The code set for a specific mode of inheritance is denoted as ${\cal{M}}$. For the additive mode of inheritance, the genotype is coded as the number of minor alleles (${\cal{M}}=\{0, 1, 2\}$); for the dominant mode of inheritance, the genotype is coded as 1 if at least one of the minor alleles is present and 0 otherwise (${\cal{M}}=\{0, 1\}$); for the recessive mode of inheritance, the genotype is coded as 1 if both alleles are minor alleles and 0 otherwise (${\cal{M}}=\{0, 1\}$).

Suppose that $(G^{c}, G^{m}, X)$ is collected from the $n_0$ control pairs $(Y=0)$ and $n_1$ case pairs $(Y=1)$. Let $n=n_0+n_1$. The retrospective likelihood for case-control mother-child data is

Assume that the phenotype prevalence $\text{pr}(Y=1)$ is known a priori to be $f$ or that it can be estimated using extra data. According to the Bayesian theorem, under constraint

maximizing the likelihood function above is equivalent to maximizing:

In this study, all the statistical inference procedures were based on the likelihood function (3). Under Mendelian inheritance and random mating, the joint probability $\text{pr}(G^c_{i},G^m_{i},X_{i})$ in (3) can be further decomposed according to the two distributional assumptions on $(G^c, G^m,X)$ (C1 and C2.1/C2.2).

C1. The child genotype is independent of maternal environmental risk factors conditional on the maternal genotype, i.e., $\text{pr}(G^c_{i}|G^m_{i},X_{i})=\text{pr}(G^c_{i}|G^m_{i})$.

C2.1. The joint distribution $(G^m_{i},X_{i})$ is completely unspecified (dependent model).

C2.2. $G^m_i$ is independent of $X_i$ (independent model).

Maternal environmental risk factors are generally maternal-related characteristics; therefore, they should be independent of the child genotype given the maternal genotype. Thus, C1 is a reasonable assumption in studies of gene-environment associations with obstetrical/maternal phenotypes using case control mother-child pair data. In contrast, maternal genotype and environmental factors can be correlated or uncorrelated, corresponding to Conditions C2.1 and C2.2, respectively.

In the following subsections, we modify the estimator proposed by Zhang et al.^[14] to obtain two estimators of the regression parameter vector ${\boldsymbol{\beta}}=(\beta_{0}, {\boldsymbol{\beta}}^{\mathsf{T}}_{1})^{\mathsf{T}}$ under Conditions C1 and C2.1/C2.2. The first estimator based on C2.1 is robust but inefficient, and the second estimator is efficient but can be biased if C2.2 is violated. We then extend the EB theory of Chen et al.^[15] to derive a shrinkage estimator that can adaptively balance bias and efficiency.

2.2   A dependent-model based estimator

Under Condition C1, the joint probability $\text{pr}(G^c, G^m, X)$ can be written as

If Mendelian inheritance and random mating hold in the population of interest, the conditional distribution $\text{pr}(G^c|G^m)$ is a function of MAF $\theta$ and fixation index $\rho$ (a parameter measuring deviation from HWE); refer to Ref. [14, Table 1] for the conditional distributions under various modes of inheritance. These distributional constraints are incorporated into the likelihood function (3): Under C2.1, the empirical likelihood method^[16] can be adopted by introducing the probability mass $\pi_{i}$ for $\text{pr}(G^m_{i}, X_{i})$ that satisfies the constraint

Denote ${\boldsymbol{\pi}}=(\pi_{1}, \dots, \pi_{n})$. The empirical log-likelihood function can be written as:

As shown in Supporting Information Appendix S1, under constraints (2) and (4), the profile likelihood function of ${\boldsymbol{\varTheta}}=(\theta,\rho, {\boldsymbol{\beta}}^{\mathsf{T}})^ {\mathsf{T}}$ can be obtained using the Lagrange multiplier method, which takes the form:

where

and $\lambda$ satisfies the equation:

which is exactly the “score” equation derived from profile likelihood $l^{ \text{J}}$ with respect to $\lambda$. Then, ${\boldsymbol{\varTheta}}$ and $\lambda$ can be estimated by solving Eq. (7) and ${\partial l^{\text{J}}({{\boldsymbol{\varTheta}}},\lambda)}/{\partial \boldsymbol{\boldsymbol{\varTheta}}}=0$,

As show in Ref. [17], it is computationally unstable to estimate ${\boldsymbol{\varTheta}}$ and $\lambda$ using the above method as the solution is a saddle point of $l^{\text{J}}$. The profile-likelihood function (6) can be modified to resolve this numerical problem. Specifically, $\lambda$ in Eq. (6) can be replaced by its limiting value:^[17]

The validity of this modification is as follows. Let ${\boldsymbol{\varTheta}}_{0}$ be the true value of ${\boldsymbol{\varTheta}}$. The “true” value $\lambda_0$ is the solution to the equations (refer to Supporting Information Appendix S2 for the detailed proof):

As a result, the so-called dependent-model estimator of ${\boldsymbol{\varTheta}}$ can be obtained by solving

Asymptotic properties of the dependent-model estimator can be established under some regularity conditions^[18]. First, with a probability tending to one as $n \to \infty$, there exists a solution $\hat{{\boldsymbol{\varTheta}}}_{\text{J}}$ to Eq. (10), which is consistent with ${\boldsymbol{\varTheta}}_{0}$. Second, $\hat{{\boldsymbol{\varTheta}}}_{\text{J}}$ is asymptotically normally distributed (refer to Supporting Information Appendix S3 for a detailed derivation).

where

Obviously, $I_{\text{J}}({\boldsymbol{\varTheta}}_{0})$ can be consistently estimated by the empirical Fisher information matrix, that is,

Among $n$ selected subjects, without loss of generality, assume that the first $n_0$ subjects are controls (i.e., $Y_i=0$ for $i=1,\cdots,n_0$) and the remaining $n_1 = n - n_0$ subjects are cases (i.e., $Y_i=1$ for $i=n_0+1,\cdots,n=n_0+n_1$). Then, $\varSigma_{\text{J}}({\boldsymbol{\varTheta}}_{0})$ can be estimated consistently using

where $a^{\otimes2}=aa^{\mathsf{T}}$ for any vector $a$. Therefore, the limiting variance-covariance matrix of $\sqrt{n}\hat{{\boldsymbol{\varTheta}}}_{\text{J}}$ can be consistently estimated using $\hat{I}^{-1}_{\text{J}}(\hat{{\boldsymbol{\varTheta}}}_{\text{J}}) \hat{\varSigma}_{\text{J}}(\hat{{\boldsymbol{\varTheta}}}_{\text{J}}) \hat{I}^{-1}_{\text{J}}(\hat{{\boldsymbol{\varTheta}}}_{\text{J}})$. Consequently, significance tests of the genetic effects can be constructed.

2.3   An independent-model based estimator

Now, we consider obtaining an independent-model estimator of ${\boldsymbol{\varTheta}}$ under Conditions C1 and C2.2. In this case, the joint probability $\text{pr}(G^c, G^m, X)$ can be written as $\text{pr}(G^c, G^m)\text{pr}(X)$. Under Mendelian inheritance and random mating, $\text{pr}(G^c, G^m)$ is a function of $\theta$ or ($\theta$, $\rho$) depending on whether HWE holds or not. The empirical likelihood method was adopted, as in the dependent model. Here, the distribution of $X$, instead of the joint distribution of $(G^m, X)$, is allowed to be nonparametric. Let $\mu_{i}$ denote the probability mass for $\text{pr}(X_{i})$ that satisfies

Denote ${\boldsymbol{\mu}}=(\mu_{1}, \cdots, \mu_{n})$. Then, the corresponding empirical log-likelihood function is

Similar to the derivation of $l^{\text{J}}({\boldsymbol{\varTheta}}, \lambda)$, under constraints (2) and (14), the profile likelihood function corresponding to Eq. (15) takes the form (refer to Supporting Information Appendix S1 for a detailed derivation):

where

and $\lambda^*$ satisfies the equation:

The MLEs of ${\boldsymbol{\varTheta}}$ and $\lambda^{*}$ can then be obtained by solving the “score” equations

Again, $\lambda^{*}$ can be replaced by $\lambda_{0}$ defined by Eq. (8). Similarly, $\lambda_{0}$ is the exact solution to the following equations:

As a result, the so-called independent-model estimator of ${\boldsymbol{\varTheta}}$ can be obtained by solving

As Supporting Information Appendix S3 shows, the independent model estimator, denoted by $\hat{{\boldsymbol{\varTheta}}}_{\text{I}}$, is asymptotically normal under some regularity conditions:

where $I_{\text{I}}({\boldsymbol{\varTheta}}_{0})$ and $\varSigma_{\text{I}}({\boldsymbol{\varTheta}}_{0})$ are the analogs of $I_{\text{J}}({\boldsymbol{\varTheta}}_{0})$ and $\varSigma_{\text{J}}({\boldsymbol{\varTheta}}_{0})$, respectively. As in Eqs. (12) and (13), $I_{\text{I}}({\boldsymbol{\varTheta}}_{0})$ and $\varSigma_{\text{I}}({\boldsymbol{\varTheta}}_{0})$ can be estimated consistently.

2.4   Two empirical Bayes-type shrinkage estimators

Compared with the dependent-model estimator, the independent-model estimator is more efficient, but serious bias could be produced if environmental risk factors are correlated with maternal genotype. Chen et al.^[15] were fully aware of the two possibly conflicting goals of improving efficiency and maintaining robustness when estimating the effects of haplotypes and haplotype-environmental interactions with standard case-control data and thus proposed an EB-type shrinkage estimator to balance efficiency and robustness in a data-adaptive fashion. Here, we adopt this “shrinkage” estimation technique to balance the estimation robustness and efficiency of the effects of maternal and child genetic variants and environmental risk factors in case-control mother-child studies.

Let ${\boldsymbol{\hat{\beta}}}_{\text{I}}$ and ${\boldsymbol{\hat{\beta}}}_{\text{J}}$ denote all genetic and environmental related effect estimates in the independent-model and dependent-model estimators, respectively. Intuitively, a weighted estimator has the form

where $I_{q}$ denotes an identity matrix of size $q = \text{dim}({\boldsymbol{\beta}})$, K is a weight factor that favors ${\boldsymbol{\hat{\beta}}}_{\text{I}}$ or ${\boldsymbol{\hat{\beta}}}_{\text{J}}$ depending on the bias of independent-model estimator. Obviously, K should have the property that it favors the robust dependent-model estimator if the independent-model estimator is biased, and vice versa. In the following, we derive the optimal weight matrix K based on the EB theory.

Throughout this study, we assume that the dependent model described in Section 2.2 is true. We denote the true value of ${\boldsymbol{\beta}}$ as ${\boldsymbol{\beta}}_{\text{J}}$. Let ${\boldsymbol{\beta}}_{\text{I}}$ be the value of ${\boldsymbol{\beta}}$ that minimizes the Kullback-Leibler discrepancy between dependent and independent models. Then, ${\boldsymbol{\hat{\beta}}}_{\text{I}}$ converges to ${\boldsymbol{\beta}}_{\text{I}}$ under certain regularity conditions. The prior distribution of ${\boldsymbol{\beta}}_{\text{J}}$ is assumed to be a multivariate normal distribution with expectation ${\boldsymbol{\beta}}_{\text{I}}$ and variance-covariance matrix ${\boldsymbol{A}}$, where ${\boldsymbol{A}}$ is a $q \times q$ positive definite matrix. In view of the asymptotic normality of ${\boldsymbol{\hat{\beta}}}_{\text{J}}$, conditional on the value of ${\boldsymbol{\beta}}_{\text{J}}$, ${\boldsymbol{\hat{\beta}}}_{\text{J}}$ approximately follows a multivariate normal distribution with mean vector ${\boldsymbol{\beta}}_{\text{J}}$ and variance-covariance matrix ${\boldsymbol{V}}$, where ${\boldsymbol{V}}$ is a $q \times q$ positive definite matrix. Consequently, the posterior expectation of ${\boldsymbol{\beta}}_{\text{J}}$ is given by

where ${\boldsymbol{\beta}}_{\text{I}}$ is estimated using $\hat{{\boldsymbol{\beta}}}_{\text{I}}$. The appropriate estimators of ${\boldsymbol{A}}$ and ${\boldsymbol{V}}$ are discussed as follows. Obviously, ${\boldsymbol{\psi}}:={\boldsymbol{\beta}}_{\text{J}}-{\boldsymbol{\beta}}_{\text{I}}$ follows a multivariate normal distribution with a zero mean vector and variance-covariance matrix ${\boldsymbol{A}}$. The prior hyperparameter ${\boldsymbol{A}}$ can be conservatively estimated using $\hat{{\boldsymbol{\psi}}}\hat{{\boldsymbol{\psi}}}^{\mathsf{T}}$, where $\hat{{\boldsymbol{\psi}}}=\hat{{\boldsymbol{\beta}}}_{\text{J}}-\hat{{\boldsymbol{\beta}}}_{\text{I}}$ ^{[15, 19]}. An estimator of ${\boldsymbol{V}}$, denoted by ${\boldsymbol{\hat{V}}}$, can be derived from the variance-covariance matrix of ${\bf{\hat{\varTheta}}}_{\text{J}}$ described in Section 2.2. The resulting EB-type shrinkage estimator can be expressed as

where the shrinkage factor, $K_1$, is

As in Ref. [15], another shrinkage factor $K_2$ alternative to $K_1$ is considered, which is a diagonal matrix with the $i{\text{th}}$ diagonal element $k_i={\hat\psi_i^2}/(\hat{v}_i+\hat\psi_i^2)$:

where $\hat{v}_i$ and $\hat{\psi_i}$ are the $i{\text{th}}$ diagonal elements of ${\boldsymbol{\hat{V}}}$ and ${\boldsymbol{\hat{\psi}}}$, respectively. This results in an alternative EB estimator:

If the observed data support the independence of environmental risk factors and maternal genotype, then ${\boldsymbol{\hat{\psi}}} \approx {\bf{0}}_{q}$ such that both ${\boldsymbol{\hat{\beta}}}_{\text{EB1}}$ and ${\boldsymbol{\hat{\beta}}}_{\text{EB2}}$ approximates ${\boldsymbol{\hat{\beta}}}_{\text{I}}$. If the observed data support the correlation between environmental risk factors and maternal genotype, then ${\boldsymbol{\hat{\psi}}} {\boldsymbol{\hat{\psi}}}^{\mathsf{T}}$ tends to be a non-zero $q$-dimensional matrix, which increases with the degree of deviation of the independent model; thus, ${\boldsymbol{\hat{\beta}}}_{\text{EB1}}$ and ${\boldsymbol{\hat{\beta}}}_{\text{EB2}}$ would tend to ${\boldsymbol{\hat{\beta}}}_{\text{J}}$.

In the following, the variance-covariance matrix of ${\boldsymbol{\hat{\beta}}}_{\text{EB1}}$ and ${\boldsymbol{\hat{\beta}}}_{\text{EB2}}$ is derived by applying the delta method. Note that they are functions of $({\boldsymbol{\hat{\beta}}}_{\text{I}}, {\boldsymbol{\hat{\beta}}}_{\text{J}})$. We first derive the asymptotic distribution of $({\boldsymbol{\hat{\beta}}}_{\text{I}}, {\boldsymbol{\hat{\beta}}}_{\text{J}})$. Based on the asymptotic expressions of ${\boldsymbol{\hat{\varTheta}}}_{\text{J}}$ and ${\boldsymbol{\hat{\varTheta}}}_{\text{I}}$ derived in Sections 2.2 and 2.3, we obtain the following asymptotic normality:

where the variance-covariance matrix

can be consistently estimated by

Here, ${\bf{0}}$ is a $q \times 2$ matrix of zeros, and $I_{\text{J}}({\boldsymbol{\varTheta}}_{0})$, $I_{\text{I}}({\boldsymbol{\varTheta}}_{0})$, $\hat{I}_{\text{J}}({\boldsymbol{\hat{\varTheta}}}_{\text{J}})$, and $\hat{I}_{\text{I}}({\boldsymbol{\hat{\varTheta}}}_{\text{I}})$ are given in Sections 2.2 and 2.3. Furthermore, $S^0_{\text{I}}({\boldsymbol{\varTheta}})$ and $S^0_{\text{J}}({\boldsymbol{\varTheta}})$ are the partial derivative vectors of $l^{\text{I}}({\boldsymbol{\varTheta}}, \lambda_0)$ and $l^{\text{J}}({\boldsymbol{\varTheta}}, \lambda_0)$, respectively, with respect to ${\boldsymbol{\varTheta}}$.

According to the delta method, the asymptotic variance-covariance matrix of ${\boldsymbol{\hat{\beta}}}_{\text{EB1}}$ and ${\boldsymbol{\hat{\beta}}}_{\text{EB2}}$ can be consistently estimated by (See Supporting Information Appendix S4)

respectively, where

and

Note that $\hat{G}_1$ and $\hat{G}_2$ are the partial derivative matrices of ${\boldsymbol{\hat{\beta}}}_{\text{EB1}}$ and ${\boldsymbol{\hat{\beta}}}_{\text{EB2}}$ with respect to $({\boldsymbol{\hat{\beta}}}_{\text{I}}^{\mathsf{T}}, {\boldsymbol{\hat{\beta}}}_{\text{J}}^{\mathsf{T}})^{\mathsf{T}}$.

3.   Simulation studies

The finite sample performance of the proposed method was evaluated through extensive simulations. First, we evaluated the estimation performance for the two EB-type estimators (23) and (24). The powers and Type I error rates of the corresponding significance tests were then examined.

3.1   Simulation setups

Several parameter combinations were considered, as described in the following subsections. For each parameter combination, genotype data $(G^m, G^c)$ were generated for a population of $1\times10^7$ mother-child pairs, and two covariates were considered: a continuous covariate $X_1$ and a discrete covariate $X_2$, where $X_1$ was generated according to a linear model and could be either independent of or correlated with maternal genotype, and $X_2$ was sampled from the binomial distribution $B(1, 0.5)$ and was independent of both maternal genotype and $X_1$. Given each vector $(G^m, G^c, X_1, X_2)$, a binary disease outcome $Y$ was generated according to the logistic regression model:

For each parameter combination, 1000 datasets, each consisting of $n_1$ case pairs and $n_0$ control pairs, were sampled from the population.

In all simulations, both $G^m$ and $G^c$ were coded as minor allele counts, so that the additive mode of inheritance was adopted in the penetrance model (1). Note that $\text{pr}(G^c|G^m)$ does not depend on the fixation index in this situation (refer to Ref. [14, Table 1]). Hereafter, let Dep represent the dependent-model-based method that does not depend on the HWE assumption. Let IndHWE and IndHWD represent the independent-model-based methods with and without the HWE assumption, respectively. Here, the HWD denotes Hardy-Weinberg disequilibrium ($\rho>0$). Let $\theta$ be fixed at 0.2, and let $\rho$ be fixed at 0 for HWE or 0.1 HWD. The covariate $X_1$ was generated according to the following linear model:

where $e$ was independent of $G^m$, and follows a standard normal distribution. Note that $\eta$ characterizes the correlation strength between $X_1$ and $G^m$, and a zero value of $\eta$ indicates independence between $X_1$ and $G^m$ and vice versa. In addition, the prevalence of $Y$ was fixed at 0.03 and assumed to be known in the relevant methods.

The value of $\eta$ was set to 0, $\log(1.5)$, or $\log(2.5)$ for various correlation strengths between $X_1$ and $G^m$. The odds ratio parameters were set to $\beta_{c} = \text{log}(1.8)$, $\beta_{m} = \text{log}(2)$, $\beta_{x_1} = \beta_{x_2} = -\text{log}(1.2)$, $\beta_{cx_1} = \beta_{cx_2} = -\text{log}(1.5)$ and $\beta_{mx_1} = \beta_{mx_2} = \text{log}(1.5)$. We considered various sample-size combinations with equal or unequal numbers of case and control pairs. In the equal situation, the common number was set to 150, 300, 600, or 1000; in the unequal situation, the ratio of case pairs to control pairs was set to be $1:2$, and the number of case pairs was set to be 100, 200, 400, or 700.

In significance tests, to illustrate the impact of the correlation between maternal genotype and environmental risk factors, the value of $\eta$ was set to range between $0$ and $\text{log}(2.5)$, and the values of $\beta_{c}$, $\beta_{x}$, and $\beta_{cx}$ were fixed as above, $\beta_{m}$ was set as $\log(1.3)$, $\beta_{mx_2}$ was set to $\log(1.2)$, and the value of $\beta_{mx_1}$ ranged between $-\log(1.5)$ and $\log(1.5)$. The significance level was fixed at $5\%$, and the sample size was $n_1=n_0 = 150$, 300, or 1000.

3.2   Estimation results

Five considered methods, namely IndHWE, IndHWD, Dep, EB1, and EB2 were applied to the generated data. To compare the estimation efficiencies of two estimators, we introduced the so-called efficiency gain of estimator A against estimator B as

where MSE(A) and MSE(B) are the mean squared errors of Methods A and B, respectively. The estimation results are summarized in Tables 1–3 and the Supporting Information Tables S1–S9. The efficiency gains with equal numbers of case and control pairs are presented in Supporting Information Tables S10–S15.

Table  1.  MSE and bias (in parenthesis) of various estimators when case and control pairs are equal and both the independence assumption ($\eta=0$) and HWE ($\rho=0$) hold.

n0 (n1)	Method	βc	βm	βx1	βx2	βcx1	βcx2	βmx1	βmx1
	IndHWE	0.056(0.004)	0.051(0.006)	0.024(–0.002)	0.110(–0.021)	0.024(–0.007)	0.101(–0.008)	0.023(0.009)	0.100(0.014)
	IndHWD	0.099(0.006)	0.088(–0.038)	0.027(–0.001)	0.124(–0.082)	0.037(–0.017)	0.166(0.007)	0.031(0.017)	0.155(0.039)
150	Dep	0.067(0.008)	0.104(0.033)	0.032(–0.012)	0.133(–0.024)	0.025(–0.010)	0.103(–0.009)	0.053(0.027)	0.193(0.019)
	EB1	0.062(0.008)	0.082(0.024)	0.029(–0.010)	0.123(–0.024)	0.024(–0.009)	0.102(–0.009)	0.040(0.022)	0.152(0.018)
	EB2	0.057(0.006)	0.068(0.021)	0.026(–0.007)	0.113(–0.023)	0.024(–0.007)	0.101(–0.008)	0.033(0.019)	0.127(0.017)
	IndHWE	0.025(0.006)	0.025(0.004)	0.011(–0.002)	0.054(0.005)	0.011(–0.009)	0.047(–0.006)	0.011(0.005)	0.045(0.006)
	IndHWD	0.041(0.004)	0.044(–0.044)	0.012(0.003)	0.064(–0.061)	0.015(–0.013)	0.068(0.016)	0.014 (0.004)	0.063(0.033)
300	Dep	0.032(0.009)	0.048(0.015)	0.015(–0.006)	0.066(0.004)	0.012(–0.011)	0.048(–0.007)	0.025(0.012)	0.091(0.007)
	EB1	0.029(0.010)	0.038(0.010)	0.013(–0.005)	0.061(0.005)	0.012(–0.011)	0.048(–0.007)	0.019(0.011)	0.070(0.007)
	EB2	0.026(0.008)	0.031(0.008)	0.012(–0.004)	0.056(0.005)	0.011(–0.010)	0.047(–0.006)	0.016(0.010)	0.058(0.005)
	IndHWE	0.013(0.004)	0.013(0.008)	0.006(–0.000)	0.026(–0.002)	0.006(–0.005)	0.025(0.004)	0.006(0.001)	0.021(–0.003)
	IndHWD	0.019(–0.014)	0.020(–0.024)	0.007(0.004)	0.035(–0.065)	0.007(–0.005)	0.032(0.031)	0.006(–0.001)	0.030(0.022)
600	Dep	0.017(0.004)	0.027(0.012)	0.008(–0.004)	0.030(–0.005)	0.007(–0.006)	0.025(0.005)	0.012(0.008)	0.048(0.003)
	EB1	0.015(0.005)	0.021(0.010)	0.007(–0.003)	0.028(–0.004)	0.007(–0.006)	0.025(0.004)	0.010(0.006)	0.036(0.001)
	EB2	0.014(0.004)	0.017(0.010)	0.006(–0.002)	0.026(–0.002)	0.006(–0.005)	0.025(0.004)	0.008(0.005)	0.029(–0.001)
	IndHWE	0.008(0.004)	0.008(–0.001)	0.003(0.001)	0.016(–0.002)	0.004(–0.005)	0.015(0.004)	0.004(0.001)	0.014(0.003)
	IndHWD	0.012(–0.013)	0.012(–0.021)	0.004(0.004)	0.023(–0.056)	0.004(–0.004)	0.020(0.026)	0.004(–0.001)	0.019(0.025)
1000	Dep	0.010(0.003)	0.015(0.007)	0.004(–0.001)	0.018(–0.002)	0.004(–0.005)	0.015(0.004)	0.008(0.004)	0.028(0.000)
	EB1	0.009(0.003)	0.012(0.005)	0.004(–0.000)	0.017(–0.001)	0.004(–0.005)	0.015(0.004)	0.006(0.004)	0.022(0.001)
	EB2	0.008(0.003)	0.010(0.004)	0.004(0.000)	0.016(–0.002)	0.004(–0.005)	0.015(0.004)	0.005(0.003)	0.018(0.003)

 | Show Table

DownLoad: CSV

Table  2.  MSE and bias (in parenthesis) of various estimators when case and control pairs are equal and the independence assumption is violated ($\eta={\rm log}(1.5)$) but HWE holds ($\rho=0$).

$n_0$ $(n_1)$	Method	$\beta_{c}$	$\beta_{m}$	$\beta_{x_1}$	$\beta_{x_2}$	$\beta_{cx_1}$	$\beta_{cx_2}$	$\beta_{mx_1}$	$\beta_{mx_2}$
	IndHWE	0.055(–0.003)	0.063(–0.071)	0.051(–0.167)	0.099(0.006)	0.022(0.007)	0.114(–0.005)	0.190(0.408)	0.100(–0.001)
	IndHWD	0.096(0.005)	0.100(–0.100)	0.046(–0.128)	0.112(–0.064)	0.033(0.006)	0.167(0.021)	0.162(0.341)	0.141(0.031)
150	Dep	0.066(0.011)	0.110(0.018)	0.031(–0.011)	0.120(–0.008)	0.023(–0.009)	0.115(–0.001)	0.046(0.012)	0.191(0.017)
	EB1	0.064(0.011)	0.099(0.006)	0.031(–0.028)	0.115(–0.006)	0.023(–0.007)	0.115(–0.002)	0.050(0.057)	0.171(0.014)
	EB2	0.056(0.003)	0.076(–0.023)	0.036(–0.086)	0.101(0.002)	0.022(0.005)	0.114(–0.005)	0.060(0.097)	0.127(0.007)
	IndHWE	0.027(–0.004)	0.031(–0.078)	0.038(–0.164)	0.049(–0.006)	0.011(0.009)	0.049(–0.006)	0.176(0.406)	0.048(0.012)
	IndHWD	0.041(–0.018)	0.049(–0.110)	0.037(–0.146)	0.063(–0.082)	0.013(0.013)	0.069(0.024)	0.162(0.376)	0.069(0.052)
300	Dep	0.032(0.002)	0.053(0.001)	0.015(0.001)	0.061(–0.017)	0.011(–0.004)	0.049(–0.002)	0.022(0.009)	0.098(0.029)
	EB1	0.032(0.002)	0.049(–0.005)	0.015(–0.016)	0.059(–0.016)	0.011(–0.003)	0.049(–0.002)	0.024(0.037)	0.091(0.028)
	EB2	0.027(–0.002)	0.036(–0.035)	0.019(–0.061)	0.050(–0.009)	0.011(0.008)	0.049(–0.006)	0.028(0.058)	0.063(0.020)
	IndHWE	0.013(–0.006)	0.019(–0.080)	0.032(–0.162)	0.026(–0.007)	0.005(0.011)	0.025(–0.000)	0.170(0.406)	0.024(0.004)
	IndHWD	0.017(–0.030)	0.028(–0.104)	0.033(–0.161)	0.038(–0.081)	0.006(0.018)	0.031(0.028)	0.168(0.398)	0.035(0.045)
600	Dep	0.015(–0.002)	0.025(–0.001)	0.007(–0.002)	0.032(–0.015)	0.006(–0.001)	0.024(0.004)	0.010(0.006)	0.046(0.014)
	EB1	0.015(–0.002)	0.024(–0.005)	0.007(–0.009)	0.031(–0.015)	0.006(–0.000)	0.024(0.004)	0.011(0.021)	0.045(0.014)
	EB2	0.013(–0.004)	0.020(–0.038)	0.009(–0.038)	0.027(–0.010)	0.005(0.009)	0.025(0.000)	0.012(0.031)	0.031(0.009)
	IndHWE	0.007(–0.001)	0.015(–0.082)	0.029(–0.161)	0.016(0.002)	0.003(0.009)	0.016(–0.006)	0.168(0.405)	0.014(0.006)
	IndHWD	0.010(–0.027)	0.021(–0.100)	0.031(–0.164)	0.025(–0.066)	0.004(0.015)	0.020(0.025)	0.169(0.405)	0.021(0.405)
1000	Dep	0.009(0.002)	0.017(0.001)	0.004(–0.000)	0.020(–0.001)	0.003(–0.002)	0.016(–0.002)	0.006(0.004)	0.029(0.006)
	EB1	0.009(0.002)	0.016(–0.001)	0.004(–0.004)	0.020(–0.001)	0.003(–0.001)	0.016(–0.002)	0.007(0.013)	0.028(0.006)
	EB2	0.008(0.000)	0.014(–0.035)	0.005(–0.024)	0.017(0.001)	0.003(0.008)	0.016(–0.006)	0.007(0.019)	0.018(0.005)

 | Show Table

DownLoad: CSV

Table  3.  MSE and bias (in parenthesis) of various estimators when case and control pairs are equal and the independence assumption is violated ($\eta={\rm log}(2.5)$) but HWE holds ($\rho=0$).

$n_0$ $(n_1)$	Method	$\beta_{c}$	$\beta_{m}$	$\beta_{x_1}$	$\beta_{x_2}$	$\beta_{cx_1}$	$\beta_{cx_2}$	$\beta_{mx_1}$	$\beta_{mx_2}$
	IndHWE	0.061(–0.053)	0.161(–0.326)	0.104(–0.295)	0.097(0.024)	0.017(0.039)	0.123(–0.009)	0.539(0.720)	0.110(–0.013)
	IndHWD	0.100(–0.057)	0.176(–0.270)	0.104(–0.256)	0.107(–0.044)	0.024(0.039)	0.171(0.035)	0.471(0.624)	0.163(–0.003)
150	Dep	0.072(–0.003)	0.128(0.033)	0.028(–0.012)	0.121(0.013)	0.018(–0.009)	0.118(0.006)	0.037(0.029)	0.205(–0.009)
	EB1	0.071(–0.004)	0.124(0.021)	0.028(–0.021)	0.119(0.014)	0.018(–0.008)	0.118(0.005)	0.040(0.051)	0.200(–0.009)
	EB2	0.063(–0.036)	0.118(–0.104)	0.037(–0.084)	0.100(0.021)	0.017(0.027)	0.123(–0.008)	0.046(0.076)	0.136(–0.012)
	IndHWE	0.028(–0.052)	0.144(–0.341)	0.091(–0.287)	0.050(–0.287)	0.010(0.042)	0.053(–0.007)	0.504(0.703)	0.047(–0.003)
	IndHWD	0.049(–0.092)	0.144(–0.323)	0.107(–0.296)	0.063(–0.079)	0.013(0.060)	0.071(0.040)	0.500(0.685)	0.067(0.009)
300	Dep	0.031(–0.012)	0.059(0.009)	0.013(–0.002)	0.064(–0.018)	0.009(–0.004)	0.049(0.010)	0.017(0.011)	0.093(0.006)
	EB1	0.031(–0.012)	0.058(0.003)	0.013(–0.007)	0.064(–0.018)	0.009(–0.003)	0.049(0.010)	0.018(0.022)	0.091(0.006)
	EB2	0.028(–0.038)	0.067(–0.094)	0.016(–0.043)	0.052(–0.008)	0.009(0.029)	0.052(–0.006)	0.019(0.034)	0.060(0.001)
	IndHWE	0.016(–0.044)	0.130(–0.339)	0.087(–0.287)	0.023(0.011)	0.006(0.043)	0.029(–0.013)	0.494(0.699)	0.026(–0.008)
	IndHWD	0.029(–0.104)	0.131(–0.331)	0.109(–0.318)	0.031(–0.072)	0.009(0.062)	0.039(0.045)	0.519(0.712)	0.033(0.005)
600	Dep	0.016(–0.005)	0.031(0.006)	0.006(–0.003)	0.029(–0.001)	0.005(–0.001)	0.027(0.004)	0.009(0.007)	0.044(0.000)
	EB1	0.016(–0.005)	0.031(0.003)	0.006(–0.006)	0.029(–0.001)	0.005(–0.001)	0.027(0.004)	0.009(0.013)	0.043(–0.000)
	EB2	0.015(–0.030)	0.037(–0.062)	0.007(–0.025)	0.024(0.008)	0.005(0.027)	0.029(–0.011)	0.009(0.019)	0.031(–0.005)
	IndHWE	0.010(–0.044)	0.121(–0.336)	0.084(–0.286)	0.015(0.007)	0.004(0.043)	0.017(–0.004)	0.485(0.694)	0.015(–0.013)
	IndHWD	0.022(–0.109)	0.125(–0.339)	0.107(–0.322)	0.025(–0.083)	0.007(0.065)	0.024(0.054)	0.519(0.717)	0.019(0.005)
1000	Dep	0.011(–0.007)	0.018(0.007)	0.004(–0.001)	0.020(–0.008)	0.003(–0.001)	0.017(0.012)	0.005(0.002)	0.026(–0.001)
	EB1	0.011(–0.007)	0.018(0.005)	0.004(–0.003)	0.020(–0.008)	0.003(–0.000)	0.017(0.012)	0.005(0.005)	0.026(–0.001)
	EB2	0.010(–0.030)	0.021(–0.039)	0.004(–0.014)	0.016(0.003)	0.003(0.023)	0.017(–0.003)	0.005(0.009)	0.018(–0.007)

 | Show Table

DownLoad: CSV

As shown in Tables 1 and S4, all five considered methods are essentially unbiased for all sample size combinations when both HWE and independence assumptions are satisfied. According to MSEs, IndHWE is more efficient than EB1 and EB2, and the latter two methods are more efficient than Dep. For example, when $n_1=n_0=150$, the efficiency gains of EB1 and EB2 against Dep were approximately $0.95\% - 24.29\%$ and $2.16\%-37.98\%$, respectively. IndHWE appears to be much more efficient than IndHWD with efficiency gains ranging between $9.04\%$ and $44.01\%$ for $n_1=n_0=150$, although IndHWD incorporates only one more parameter $\rho$ compared with IndHWE (Supporting Information Table S10). When the independence assumption holds but HWE is violated, IndHWE could produce considerable biases on maternal genetic effects but not on the other effects in our simulation situation, so that IndHWD outperforms IndHWE in estimating maternal genetic effects (Supporting Information Tables S1, S7, and S13). However, the estimation results of the remaining four methods appeared to be unaffected by the violation of HWE.

When HWE holds, but the $G^m$-$X_1$ independence assumption is violated ($\eta=\log(1.5)$), both IndHWE and IndHWD could be seriously biased, especially for main effect of X₁ and the interaction effect between and $G ^ m$ and $X_1$ (Tables 2 and S5). As expected, Dep, which does not require HWE or independence assumptions, is inherently unbiased. Correspondingly, the MSEs of the above two effects were generally much larger for IndHWE than for Dep. Furthermore, EB1 and EB2 can effectively balance the bias and efficiency. The efficiency gains of EB1 and EB2 against IndHWE are considerable for the above two effects. For example, when $n_1= n_0=150$, the efficiency gains of EB1 and EB2 against IndHWE for the two effects are approximately $29.67\%-40.04\%$ and $68.67\%-73.91\%$, respectively. For the other effects, EB1 and EB2 are less efficient than IndHWE but more efficient than Dep, with efficiency gains against Dep ranging between $0.22\%$ and $10.56\%$ for EB1 and 0.86% and 33.52% for EB2 when the sample size is 150 (Supporting Information Table S11). In addition, EB1 and EB2 are slightly more efficient than IndHWD in the simulation situations. When both HWE and the independence assumption are violated, both IndHWE and IndHWD are biased, and their performances are poorer than those of the other three methods (Supporting Information Tables S2, S8, and S14).

With a larger $\eta$ of $\log(2.5)$, IndHWE and IndHWD are more biased toward the maternal effect, $X_1$ environmental effect, and $G^m$-$X_1$ interaction (Tables 3 and S6), in contrast with the results for $\eta=\log(1.5)$ (Tables 2 and S5). Consequently, the efficiency gains of EB1 and EB2 against IndHWE and IndHWD increased. For example, with a sample size of 150, the efficiency gains of EB1 and EB2 against IndHWE for the three effects increased to $23.00\% - 26.30\%$, $64.75\% - 72.99\%$, and $91.51\% - 92.62\%$, respectively (Supporting Information Table S12). When HWE was violated, the efficiency gains of EB1, EB2, and Dep against IndHWE and IndHWD were even greater (Supporting Information Tables S3, S9, and S15).

In summary, IndHWE is the preferred method when the assumptions of HWE and maternal gene-environment independence hold. However, when any one of the two assumptions is violated, serious estimation biases can be produced for maternal genetic effects, environmental effects, and maternal gene-environment interactions. Dep is the most robust with respect to the two assumptions, but it is much less efficient than IndHWE, and its performance is much poorer than IndHWE in terms of MSE when independence is satisfied. In contrast, EB1 and EB2 can balance estimation bias and efficiency by appropriately weighting IndHWE and Dep. In the small-sized situation, EB1 is generally more unbiased than EB2, although the former is slightly less efficient. Consequently, EB1 may be preferable to EB2 with a small sample size; however, EB2 may be preferable to EB1 with a moderate or large sample size.

The five considered methods require specification of the prevalence $f$. Therefore, it would be interesting to study the impact of misspecifying $f$ on the estimation results. A sensitivity analysis was conducted for this purpose. The true parameters (including $f$) were the same as those of Table 1, but $f$ was misspecified to be 0.01 or 0.1 for all of the five methods. The simulation results shown in Supporting Information Tables S16–S21 suggest that all methods were not affected by misspecifying the prevalence, except that the MSEs were slightly larger, with the prevalence being overspecified. This finding is consistent with those of Refs. [14, 20].

3.3   Significance test results

In this subsection, the significance test results (type-I error rates and powers at a significance level of 0.05) of the five considered methods are summarized for maternal gene-environment interactions in the simulation settings described in Section 3.1. Fig. 1 and Supporting Information Table S22 show type-I error rates for the significance tests of $\beta_{mx_1}$ for various $\eta$ ($0-\log(2.5)$), sample sizes ($n_1=n_0=150$, $300$, and $1000$), and $\rho$ (0 for HWE and 0.1 for HWD). As expected, type-I error rates of IndHWE and IndHWD increase dramatically with $\eta$, because both methods rely on the independence assumption. In contrast, type-I error rates of both EB1 and EB2 are only slightly inflated (type I error rates: $4.4\% - 7.5\%$ and $3.9\%- 9.9\%$ for EB1 and EB2, respectively, $n_1= n_0=150$ and $\rho=0$), controlled around the nominal level even for large $\eta$, demonstrating the robustness of EB1 and EB2. Not surprisingly, Dep has well-controlled type-I error rates ($4.4\% - 6.2\%$, $n_1=n_0=150$, and $\rho=0$).

Figure  1.  Type-I error rates for the significance tests of $\beta_{mx_1}$ with $\rho=0$ (HWE) or $\rho=0.1$ (HWD), various sample sizes ($n_1=n_0=150$, $n_1=n_0=300$, and $n_1=n_0=1000$), and various $\eta$ values ($0$ through $\log(2.5)$). The other parameters were fixed: $\beta_{c}=\text{log}(1.8)$, $\beta_{m}=\text{log}(1.3)$, $\beta_{x_1}=\beta_{x_2}=-\text{log}(1.2)$, $\beta_{cx_1}=\beta_{cx_2}=-\text{log}(1.5)$, and $\beta_{mx_2}=\text{log}(1.2)$.

DownLoad: Full-Size Img PowerPoint

Figs. 2 and 3 show the powers for testing the $G^m$-$X_1$ interaction for various $\beta_{mx_1}$ ($\text{from} -\log(1.5) \; \text{to}\; \log(1.5)$), sample sizes ($n_1=n_0=150$, $300$, and $1\;000$), and $\rho$ (0 for HWE and 0.1 for HWD). Because Type I error rates of IndHWE and IndHWD cannot be controlled when $\eta>0$, powers are not reported for IndHWE and IndHWD in such a situation. As expected, IndHWE is the most powerful among the five methods when the independence assumption holds ($\eta=0$), followed by IndHWD ranks the second. Evidently, EB1 and EB2 were more powerful than Dep in this situation. When the independence assumption does not hold, EB1 and EB2 can be either more or less powerful than Dep, depending on the sign of $\beta_{mx_1}$.

Figure  2.  Powers for the significance tests of maternal gene-environment interaction ($\beta_{mx_1}$) with HWE ($\rho=0$), various sample sizes ($n_1=n_0=150$, $n_1=n_0=300$, and $n_1=n_0=1000$), and various $\eta$ values ($0$, $\log(1.2)$, and $\log(1.5)$). The other parameters were fixed: $\beta_{c}=\text{log}(1.8)$, $\beta_{m}=\text{log}(1.3)$, $\beta_{x_1}=\beta_{x_2}=-\text{log}(1.2)$, $\beta_{cx_1}=\beta_{cx_2}=-\text{log}(1.5)$, and $\beta_{mx_2}=\text{log}(1.2)$.

DownLoad: Full-Size Img PowerPoint

Figure  3.  Powers for the significance tests of maternal gene-environment interaction ($\beta_{mx_1}$) with HWD ($\rho=0.1$), various sample sizes ($n_1=n_0=150$, $n_1=n_0=300$, and $n_1=n_0=1000$), and various $\eta$ values ($0$, $\log(1.2)$, and $\log(1.5)$). The other parameters were fixed: $\beta_{c}=\text{log}(1.8)$, $\beta_{m}=\text{log}(1.3)$, $\beta_{x_1}=\beta_{x_2}=-\text{log}(1.2)$, $\beta_{cx_1}=\beta_{cx_2}=-\text{log}(1.5)$, and $\beta_{mx_2}=\text{log}(1.2)$.

DownLoad: Full-Size Img PowerPoint

The diagrams of Type I errors and powers for the significance tests of the combined effect of maternal gene-environment interactions are shown in the Supporting Information Figures S1–S3, which have patterns similar to those of the $G^m$-$X_1$ interaction effect.

4.   Application to the Danish National Birth Cohort study

Previous studies have shown that spontaneous preterm birth (SPTB) is the single largest cause of premature birth and is thought to be caused by a variety of factors, including infection or inflammation, stress, malnutrition, and genetic factors^[5]. Using case-control mother-child pair data nested in the Danish National Birth Cohort (DNBC), we examined the association between SPTB and candidate SNPs by adjusting for maternal pre-pregnancy BMI (pp-BMI).

The DNBC is a cohort study that began in 1996. Extensive biomaterials and epidemiological data were collected from more than 10 0000 mothers and their children in Denmark. It has been well established that both fetal and maternal genotypes contribute to the risk of preterm birth, and genetic loci associated with SPTB have been identified through genetic association studies^[21–24]. Furthermore, low maternal pp-BMI was believed to be associated with a high risk of SPTB^[25]. In this study, 720 mothers gave birth prematurely (gestational days $< 259$lt; 259$) and 906 term deliveries ($280 \leqslant$ gestational days $\leqslant 286$) among 1626 eligible mother-child pairs.

The SNPs rs9939609 in genes FTO and rs2684811 in gene IGF1R have been identified as associated with BMI and SPTB, respectively^{[26, 23]}. In the current study, we focused on the common SNPs located within 50kb of rs9939609 (eight SNPs in this region) and rs2684811 (24 SNPs in this region). All these 32 SNPs appear to be in HWE according to our analysis results by PLINK^[27]; refer to Supporting Information Table S24 for the estimated MAFs and the p-values for testing HWE. Moreover, maternal pp-BMI can be assumed to be independent of the child genotype conditional on the maternal genotype since it is a maternal phenotype. Independent tests were also performed between pp-BMI and maternal genotype using the Kruskal-Wallis test. We applied the five methods mentioned in Section 3 to this dataset and adopted an additive mode of inheritance for each method. As the five methods are robust with respect to misspecification of disease prevalence according to our simulation results described in Section 3, we fixed the preterm birth prevalence at 7% because the prevalence of preterm birth in Denmark is approximately 5%–9% ^{[28, 29]} according to a survey of birth data in developed countries.

All the considered methods are expected to be nearly unbiased if $G^{m}$ is independent of $X$, and EB1 and EB2 should be preferable to Dep in this situation. On the contrary, if $G^{m}$ is strongly correlated to $X$, then IndHWE is generally biased, and EB1 and EB2 should be preferable to IndHWE and IndHWD. To verify this assertion, we studied the relationship between the shrinkage strengths of EB estimators and the association strength between $X$ and $G^m$. First, a linear transformation was performed such that the transformed effects of Dep and IndHWE were 0 and 1, respectively. Then, the shrinkage strengths of the two EB estimators were measured using the regression coefficient between the transformed EB estimators and the p-value for testing the independence between $X$ and $G^m$. As shown in Supporting Information Fig. S4, the desired shrinkage strength of EB1 is significant (regression coefficient $= 0.398$; 95% confidence interval $= [0.130, 0.665]$), which is effect-independent (refer to Supporting Information Appendix S5 for a rigorous proof). In contrast, the shrinkage strength of EB2 is significant only for the $G^m$-$X$ interaction effect $\beta_{mx}$ (regression ${\rm coefficient} = 0.277$; 95% confidence interval $= [0.057, 0.497]$). In the following, we describe the analysis results for two selected SNPs (rs16945088 and rs9941349) in detail, which are also summarized in Table 4. The results of the analysis of the other 30 candidate SNPs are summarized in Supporting Information Tables S25–S28.

Table  4.  Genetic effect estimates of two common SNPs in genes FTO and IGF1R (Danish National Birth Cohort).

SNPa	log(OR)	Dep				IndHWE				EB1				EB2
		Estb	SEc	Pd		Estb	SEc	Pd		Estb	SEc	Pd		Estb	SEc	Pd
rs16945088	$\beta_c$	–0.151	0.140	0.283		–0.146	0.127	0.250		–0.146	0.127	0.249		–0.146	0.127	0.250
(0.636)	$\beta_m$	0.055	0.147	0.709		0.049	0.125	0.694		0.050	0.125	0.692		0.049	0.125	0.694
	$\beta_{cx}$	0.045	0.118	0.702		0.044	0.115	0.703		0.044	0.116	0.703		0.044	0.115	0.703
	$\beta_{mx}$	–0.081	0.160	0.613		–0.103	0.114	0.367		–0.102	0.116	0.382		–0.103	0.115	0.373
rs9941349	$\beta_c$	0.053	0.081	0.511		0.109	0.071	0.127		0.059	0.079	0.455		0.091	0.077	0.236
(0.039)	$\beta_m$	0.034	0.086	0.693		–0.031	0.074	0.676		0.027	0.083	0.747		–0.007	0.082	0.929
	$\beta_{cx}$	–0.070	0.073	0.341		–0.071	0.073	0.327		–0.070	0.073	0.340		–0.071	0.073	0.327
	$\beta_{mx}$	0.073	0.083	0.380		0.176	0.074	0.017		0.084	0.082	0.304		0.113	0.085	0.180
[Note] a Selected SNPs (the first one is BMI-related candidate SNP in gene FTO, and the second one is SPTB-related candidate SNP in gene IGF1R). Presented in parentheses are p-values for testing the independence between the SNPs and maternal pp-BMI using Kruskal-Wallis test. b Estimated effects. c Estimated standard errors. d Significance test p-values.

 | Show Table

DownLoad: CSV

For the first SNP rs16945088, the maternal genotype is shown to be independent of maternal pp-BMI ($\text{p-value} = 0.636$). Consequently, the effect estimates of Dep and IndHWE are close to each other, IndHWE is generally more efficient than Dep in terms of standard errors, and the estimates of EB1 and EB2 are much closer to those of IndHWE than to Dep. For example, the $\beta_{mx}$ estimates of Dep, IndHWE, EB1, and EB2 are $-0.081$, $-0.103$, $-0.102$, and $-0.103$, respectively, indicating that the estimates of EB1 and EB2 are much closer to that of IndHWE than those of Dep. In addition, the standard errors of the $\beta_{mx}$ estimates by IndHWE, EB1, and EB2 are 0.114, 0.116, and 0.115, respectively, compared to 0.160 by Dep, indicating that the former three methods are more efficient than Dep. For the second SNP rs9941349, maternal genotype is significantly associated with pp-BMI (p-value = 0.039). As a result, the estimates of EB1 are closer to those of Dep compared with IndHWE, whereas the estimates of EB2 are closer to those of IndHWE compared with Dep except that the $\beta_{mx}$ estimate by EB2 lies in the middle of the $\beta_{mx}$ estimates of Dep and IndHWE. For example, the $\beta_{mx}$ estimates of Dep, EB1, and EB2 are 0.073, 0.084, and 0.113, respectively, compared to 0.176 by IndHWE. Furthermore, the $G^m$-$X$ interaction is not significant for Dep, EB1, and EB2 ($\text{p-values} > 0.3$), whereas it is significant for IndHWE ($\text{p-value}= 0.017$). This suggests that the seemingly significant result of IndHWE might be a false positive finding.

5.   Conclusions

In this study, two EB-type estimators (EB1 and EB2) are proposed by appropriately weighting an independent model estimator and a dependent-model estimator, which can adaptively balance robustness and efficiency. The independent-model estimator is an efficient estimator that improves efficiency by imposing HWE and $G^m$-$X$ independence assumptions in the model, but it could produce a serious estimation bias if at least one of the assumptions is violated. In contrast, the dependent-model estimator Dep is a robust estimator that does not require the above two assumptions but might not be efficient. The EB-type estimators EB1 and EB2 developed in this study adaptively utilize the HWE assumption and the independence assumption between maternal genotype and maternal environmental risk factors, which appropriately balance robustness and efficiency. As in Dep, EB1 and EB2 fully utilize Mendelian inheritance and the conditional independence assumption between maternal environmental risk factors and children’s genotype given maternal genotype.

These properties were verified by our simulation results in various situations. Specifically, if both HWE and $G^m$-$X$ independence assumptions hold, EB1 and EB2 are more efficient than Dep, although they are less efficient than IndHWE. In contrast, IndHWE could be seriously biased if HWE or $G^m$-$X$ independence assumptions are violated, while EB1 and EB2 still perform well. If the sample size is small, EB1 is more robust than EB2, whereas EB2 appears more suitable for moderate or large samples in terms of estimation efficiency. The desired properties of EB1 and EB2 are also demonstrated using a real-data application. In practice, it is difficult to obtain prior information regarding the relationship between maternal genotype and environmental risk factors. In this situation, EB1 and EB2 serve as robust alternatives to IndHWE and Dep.

Our novel method has been designed to analyze common genetic variants (MAF$\geqslant 5\%$) and the low-dimensional covariate $X$. In practice, there could be many mother-related characteristics, such as age, pre-pregnancy BMI, and smoking. The dimension of the parameters could be at least three times that of $X$ when studying gene-environment interactions in case-control mother-child pair data. With high-dimensional parameters, we can adopt any regularization method to select associated covariates. However, increasing evidence shows that rare genetic variants (MAF$< 1\%$lt; 1\%$) play an important role in complex diseases and traits^[30–32]. Recently, some methods have been developed to analyze family based multilocus rare variants^[33–36]. However, these methods are not based on retrospective likelihood functions; therefore, they cannot fully utilize family information, such as Mendelian inheritance and gene-environment independence. Further investigation is thus needed to extend our method to analyze multilocus rare variants under a case-control mother-child pair design.

Data availability statement

The DNBC data are available at the dbGAP website (study accession: phs000103.v1.p1).

Supporting information

The supporting information for this article is available online at https://doi.org/10.52396/JUSTC-2022-0007. It includes Appendices S1–S5, Figures S1–S4, and Tables S1–S28.

Acknowledgements

This work was supported by the National Natural Science Foundation of China (11771096, 12171451, and 72091212). Funding support for the GWAS of Prematurity and Its Complications study was provided by the NIH Genes, Environment, and Health Initiative [GEI] (U01HG004423). The GWAS of Prematurity and its Complications study is one of the genome-wide association studies funded as part of the Gene Environment Association Studies (GENEVA) under GEI. Assistance with phenotype harmonization and genotype cleaning, as well as general study coordination, was provided by the GENEVA Coordinating Center (U01HG004446). The National Center for Biotechnology Information provided assistance with data cleaning. Funding support for genotyping, which was performed at the Johns Hopkins University Center for Inherited Disease Research, was provided by the NIH GEI (U01HG004438) and the NIH contract “High-throughput Genotyping for Studying the Genetic Contributions to Human Disease” (HHSN268200782096C).

Conflict of interest

The authors declare that they have no conflict of interest.