An Efficient Test for Gene-Environment Interaction in Generalized Linear Mixed Models with Family Data

Gene-environment (GE) interaction has important implications in the etiology of complex diseases that are caused by a combination of genetic factors and environment variables. Several authors have developed GE analysis in the context of independent subjects or longitudinal data using a gene-set. In this paper, we propose to analyze GE interaction for discrete and continuous phenotypes in family studies by incorporating the relatedness among the relatives for each family into a generalized linear mixed model (GLMM) and by using a gene-based variance component test. In addition, we deal with collinearity problems arising from linkage disequilibrium among single nucleotide polymorphisms (SNPs) by considering their coefficients as random effects under the null model estimation. We show that the best linear unbiased predictor (BLUP) of such random effects in the GLMM is equivalent to the ridge regression estimator. This equivalence provides a simple method to estimate the ridge penalty parameter in comparison to other computationally-demanding estimation approaches based on cross-validation schemes. We evaluated the proposed test using simulation studies and applied it to real data from the Baependi Heart Study consisting of 76 families. Using our approach, we identified an interaction between BMI and the Peroxisome Proliferator Activated Receptor Gamma (PPARG) gene associated with diabetes.


Introduction
Linear mixed models (LMM) have been used to find associations between continuous phenotypes and genetic variants, genes, and gene-environment (GE) interactions in unrelated and related subjects in genome-wide association (GWA) analysis. For unrelated subjects, the analysis can be performed within the generalized linear model framework, however, for related subjects as in the case of family data, one has to include the kinship matrix to take into account the correlation among the relatives for each family. In this paper, we are interested in testing GE interaction for discrete phenotypes. Generalized linear mixed models (GLMM) proposed by Breslow and Clayton [1] is an ideal statistical approach to detect such an interaction with non-continuous phenotypes, because it can treat the familiar effect on the phenotype as a random effect.
Gene-based GE interaction tests have previously been proposed for independent subjects [2][3][4]. While each GE interaction can be tested individually using one single nucleotide polymorphism (SNP) at a time, it is known that single SNP association is not as powerful as the gene-based analysis [2] due to the linkage disequilibrium (LD) present among the SNPs in a gene. Lin et al. [2] proposed a variance component test (VCT) of the interactions by treating the interactions as a random effect.
This approach was extended to sequencing data with rare variants [3,5]. To overcome multicollinearity of the coefficients of the genetic markers, Lin et al. [2,3] applied ridge regression penalization of SNP coefficients and estimated the ridge penalty parameter with generalized cross-validation. However, this method is computationally demanding and their final test ignores the tuning of the ridge penalty parameter. Coombes [6] instead proposed treating the genetic coefficients as a random effect in a linear mixed model framework to perform the ridge penalization. This equivalence was initially proposed by Bishop and Tipping [7,8] for Bayesian ridge regression in linear models framework. While this approach was able to incorporate the ridge penalty into the test statistic, it was only developed for a quantitative phenotype [6].
Here, we propose a GLMM GE interaction framework for discrete and continuous phenotypes that treats the coefficients of genetic markers as random effects. Also, because the correlation among relatives cannot be ignored, this modeling framework incorporates the kinship matrix in the GLMM [9]. We test for GE interaction between a set of SNPs and an environment by treating interaction coefficients as random effects using a VCT. Our proposed model includes three random effects: the first for genetic variants, the second for gene-environment interaction, and the third for the inclusion of families. In the methods section, we develop the framework for this model. We present the VCT as proposed by Lin [10] in the presence of several random effects and adapt this VCT to accommodate the interactions [11]. We also prove that the corresponding best linear unbiased predictor (BLUP) in our GLMM model is equivalent to the ridge regression estimator, as proposed by Shen et al. [12]. In our simulations, we show that our model can be efficiently computed using the GMMAT package [13] in R and maintains appropriate type I error as well as sufficient power. Finally, we apply our methodology to test for genetic interactions with BMI associated with diabetes among the Baependi Heart Study [14], which consists of 76 families of varying pedigree size.

Generalized Linear Mixed Model
GLMMs have been widely applied in situations where the outcome is discrete and random components are involved in the linear predictor. GLMMs, as proposed by Breslow and Clayton [1], link a response variable y i , for i = 1, . . . , n, with vectors x i and z i of explanatory variables associated with the fixed and random effects. Given a r-dimensional vector d of random effects, the model is given is known as the link function and µ d i = E(y i |d) is the conditional mean. The conditional variance is given by Var(y i |d) = φω i ν µ d i , with ν(.) is a known function, φ is a scale parameter and ω i are known weights. Denoting the observation vector by y = (y 1 , . . . , y n ) T and the design matrices with rows x T i and z T i by X and Z, the general formulation of GLMM is given by with µ d = (µ d 1 , . . . , µ d n ) T , and where d is assumed multivariate normal distributed with mean 0 and covariance matrix D = D(π) depending on an unknown vector π of variance components.

Generalized Linear Mixed Model with Gene-Environment
To set up the GLMM to model GE interaction in families, assume we have a random sample of N independent families from a study population, with n i members in the ith family such the total number of individuals is n = ∑ N i=1 n i . For the jth member in the ith family, let Y ij be a discrete or continuous response variable for the phenotype of interest, X ij = (X 1 ij , . . . , X p ij ) T , with p equal to the number of non-environmental covariates, G ij = (G 1 ij , . . . , G q ij ) T with q equal to number of observed genotypes for SNPs in a gene, E ij the environmental variable of interest, and S ij = (E ij G 1 ij , . . . , E ij G q ij ) T the GE interaction. The GE interaction GLMM for families may be written as where g(.) is a monotone known function, Y ij |α ij follows a distribution in the exponential family, ν(.) is a known function, φ is a scale parameter that may be known or may need to be estimated, ω ij are known weights (commonly equal to 1), Φ i is the kinship matrix, and σ 2 is a parameter to be estimated.
Equation (2) can be rewritten using the Cholesky decomposition of the kinship matrix for the ith family 2Φ i = K i K T i and assuming that Then, the family model is given by To simplify the computational burden in the estimation process, we generalize Equation (3), by redefining its vectors and matrices as: Our goal is to test the null hypothesis H 0 : γ = 0 that the gene of interest has no GE interaction associated with response.

Proposed GE Interaction Test
As mentioned by Lin et al. [2], to treat γ as a fixed vector and proceed with a p degrees of freedom (DF) score test can result in loss of power to test for interaction. Another common strategy is to use a single SNP analysis of GE interaction, which assumes all SNPs are uncorrelated, however, this is usually not the case. In the majority of cases, the SNPs within a gene are highly correlated, thus, here we propose a test that accounts for correlation among SNPs and has uses less DF than the score test.
Using Equation (4), our proposed test assumes γ is a random vector following an arbitrary distribution with mean 0 and variance τ I q . Thus, a test of the null hypothesis H 0 : τ = 0 is equivalent to testing H 0 : γ = 0. For simplicity, we assume γ ∼ N(0, τI q ).
In order to account for LD among SNPs in a gene and avoid estimation issues related to multicollinearity, we use a ridge regression approach to impose a penalty on θ in the PQL proposed for GLMM [1]. However, the selection of a penalty parameter can be computationally demanding [2]. Thus, to expedite and incorporate the selection of a penalty parameter used in our proposed test, we specify d 1 = Gθ, as a random effect in Equation (4), with θ ∼ N(0, σ 2 θ I q ). By using this approach, we demonstrate later in Section 3.1, under the null model, the best linear unbiased predictor (BLUP) of d 1 is equivalent to the ridge regression estimator. Denoting d 2 = Kb and d 3 = Sγ, and assuming d 1 , d 2 and d 3 to be independent, we can write Equation (4) in the GLMM form as with Z = [I n I n I n ], where I n is the (n × n) identity matrix and Based on the penalized quasi-likelihood (PQL) [1], Lin [10] developed a VCT for independent subjects in the framework of a GLMM. To test the null hypothesis H 0 : τ = 0, the VCT uses the score statistic where β and π are the maximum likelihood (ML) estimators for β and π = (σ 2 θ , σ 2 , φ) T , respectively, under the null model as described in Section 3.1. In addition, is calculated under the null model, where g (.) denotes the derivate of function g(.). The covariance matrix for the working vectorỸ is given by Since SS T does not have a block diagonal structure, the score U τ β, π cannot be written as a sum of N independent random variables, corresponding to the families. Therefore, the asymptotic distribution of U τ β, π is not a normal distribution, in contrast to the VCT of Lin [10]. Instead, we follow the approach developed by Zhang and Lin [11], and propose, as score statistic, the first term in Equation (6), which corresponds to the quadratic form To correct for bias, we use the restricted maximum likelihood (REML) estimators [1] in the GLMM framework to obtain β and π under the null hypothesis.
Zhang and Lin [11] showed that under H 0 : τ = 0, U τ follows approximately a mixture of one degree of freedom, independent chi-square distributions. However, for computational ease, we use the Satterthwaite method [15] to approximate the distribution of U τ by a scaled chi-square distribution κχ 2 ξ , where the scale parameter κ and the degrees of freedom ξ can be calculated by equating the mean and variance of U τ to those of κχ 2 ξ . When REML estimates are used to calculate U τ , Zhang and Lin [11] showed that the mean and variance of U τ can be approximated, respectively, by Dashed lines in J and M represent the cases: (i) φ known, implying a (2 × 1) vector and (2 × 2) matrix, respectively and (ii) φ unknown, implying a (3 × 1) vector and (3 × 3) matrix, respectively.
As the mean and variance of κχ 2 ξ are given by κξ and 2κ 2 ξ, respectively, we obtain the equations tr( PSS T ) = κξ and I τ = 2κ 2 ξ, where P denotes the matrix P evaluated in π. By solving these equations, we demonstrate that κ = I τ /[2 tr( PSS T )] and ξ = 2 tr( PSS T ) 2 /I τ . Therefore, to test H 0 : τ = 0, we propose the statistic which follows approximately a chi-square distribution with ξ degrees of freedom.

Null Model Estimation
Our proposed score test requires that we first fit the null model. Under the null hypothesis H 0 : τ = 0, Equation (5) becomes To estimate the parameters in Equation (8), Breslow and Clayton [1] proposed a Fisher scoring solution that may be expressed as the iterative solution to the system This system is equivalent to the so called Henderson equations [16] for computing the best linear unbiased estimator (BLUE) of β and the best linear unbiased predictor (BLUP) of d 0 .
By re-expressing d 0 to (d T 1 , d T 2 ) T in Equation (9), we obtain the system   X This new system is equivalent to using a ridge regression penalization for parameter θ in Equation (4) (see Appendix A for details). Assuming that π = (σ θ , σ, φ) T is known, it can be shown that the solution of Equation (10) is given by the following equations: with Σ = σ 2 θ GG T + σ 2 KK T + φW −1 . Chen et al. [17] fitted their GLMM by defining, d s = (d 1 + d 2 ), and assuming d s ∼ N 0, σ 2 θ GG T + σ 2 KK T . However, the BLUP for d s is identical to the sum of BLUPs in Equation (11). Therefore, we use their PQL-based estimation algorithm implemented in the R package GMMAT (see Chen et al. [17] and Chen and Conomos [13] for details) to estimate our null model. The iterative process uses REML to estimate the variance parameters vector π used in the score statistic Equation (7).

Simulations
For our simulations, we used SimPed [18] to generate 100 SNPs (50 independent and 50 in LD) for 1000 independent families with identical pedigree structures of size 10 Figure 1. For each simulation, we randomly sampled without replacement 100 families to obtain a sample of 1000 individuals.
The environment was simulated to be correlated within a family and depend on age and sex of the subject using the following model with parameters chosen to mimic our real data example: where I(·) is the indicator function, ε i = (ε i1 , . . . , ε i10 ) T ∼ N(0, 4I 10 ) where I 10 is the (10 × 10) identity matrix, and γ i ∼ N(0, 4). In our simulations, we considered a SNP-set with 50 independent SNPs or in LD. The correlation structure for the SNPs in LD is shown in Figure 2. Using the R package SimCorMultRes [19], we simulated a correlated binary phenotype dependent on family using the following mean structure: where Φ i is the kinship matrix corresponding to the family pedigree in Figure 1. Given SimCorMultRes only allows for specification of the correlation matrix, σ 2 is set equal to 1. G1 ij and G2 ij are either independent SNPs (MAF = 0.3 and 0.1) or the first and fifth SNPs with MAF = 0.3 and 0.17 respectively from Figure 2. Note that only the SNPs with a main effect interact with the environment in our model. We generated 10,000 and 1000 datasets to estimate type I error and empirical power, respectively, at an α = 0.05 level. Using these datasets, we compared the performances of the score test, MinP test, and our proposed VCT. As previously mentioned, the score test treats γ as a fixed vector and results in a p DF test. The null model for this test was estimated as specified in Section 3.1. For the MinP test, which represents a single SNP analysis of GE interaction, we independently model the single SNP-by-environment interaction using Equation (4) where G and S are a vector, rather than a matrix, for a SNP and GE interaction, respectively. For each model, we calculated the p-value for the test of interaction and found the minimum p-value among all tests. We corrected for multiple testing by multiplying by the number of effective SNPs in the gene [20]. In the case of independence, the number of effective SNPs would be equivalent to the number of SNPs in the gene.  3  4  5  6  7  8  9  10  11  12  13  14  15  16  17  18  19  20  21  22  23  24  25  26  27  28  29  30  31  32  33  34  35  36  37  38  39  40  41  42  43  44  45  46  47  48  49  50  1  2  3  4  5  6  7  8  9  10  11  12  13  14  15  16  17  18  19  20  21  22  23

Type I Error
We first compared the empirical type I error of the different methods at 0.05 α-level. To evaluate type I error, we set γ 1 = γ 2 = 0 and varied the number of SNPs q in the gene. The SNPs were either independent or in LD. The empirical type I error rates are shown in Table 1 as well as the mean of the fittedσ 2 ,σ 2 θ , and λ = 1/ σ 2 G parameters from Equation (8) across all simulations. While the variance component for the random effect defined by the kinship matrix stays approximately constant, the penalization term λ increases as the number of SNPs in the model increases. Thus, like in ridge regression, increasing the number of parameters results in an increased penalization of the model. All of the methods were conservative in our simulations, but as q increased, the score test became useless due to the large DF of the test.

Empirical Power for Independent SNPs
We next compared the empirical power of the different methods with either five, 10, or 50 independent SNPs. We varied the amount of interaction for the two selected SNPs by varying γ 1 = γ 2 from 0 to 0.1 by 0.01. Figure 3 shows that as the number of SNPs in the model increases, each of the methods lose power due to the increase in DF of each test. While the VCT performs best with five or 10 SNPs, the MinP test performs best for 50 SNPs because only two out of 50 SNPs have interaction. The MinP test will always perform best if very few of the SNPs in the set have interaction.

Empirical Power for SNPs in LD
Finally, we compared the empirical power of the different methods with either five, 10, or 50 SNPs in LD. We varied the amount of interaction for the two selected SNPs by varying γ 1 = γ 2 from 0 to 0.1 by 0.01. Figure 4 shows that as before, each of the methods lose power as q increases. The VCT outperforms the MinP test in all scenarios because the SNPs are correlated which the MinP test fails to account for.

Application to Baependi Data
We use our proposed method to test for GE interactions in the Baependi Heart Study [14] between BMI and three different candidate genes that may be associated with type II diabetes (T2D). The first candidate gene we studied was the Peroxisome-Proliferator-Activated Receptors gamma (PPARG) gene, which is a key regulator of adipocyte differentiation and energy balance. Two of the mutations in the PPARG gene have been shown to be associated with obesity or diabetes-related phenotypes in different populations [21]. PPARG2, the predominantly isoform of PPARG, is expressed selectively and at a higher level in adipose tissue, where it modulates the expression of target genes implicated in adipocyte differentiation and glucose homeostasis [22]. Thus, the PPARG2 gene is a major candidate gene for T2D or obesity, both being complex phenotypes determined by the combination of multiple genetic and environmental factors [23,24]. The second candidate gene studied was the Fat Mass and Obesity associated protein (FTO), which confers risk for obesity and BMI. Since obesity is known to be a predisposing factor for the development of T2D, it is not surprising that variants in FTO have been also found in T2D GWAS [25]. The final candidate gene studied was the cyclin-dependent kinase 5 regulatory subunit associated protein 1-like 1 (CDKAL1) gene which confers risk for obesity and T2D [26]. In our study, the PPARG, FTO, and CDKAL1 genes had 16, 149, and 186 genetic variants genotyped, respectively. However these numbers of SNPs do not represent the number of effective SNPs discussed in Gao et al. [20], that is equivalent of the number of principal components to reach 99.5% of their total variation. Then, the effective number of SNPs associated with the PPARG, FTO, and CDKAL1 genes are 10, 93, and 92, respectively. We used Equation (5) to specify a logistic GLMM to test for GE interaction of the aforementioned genes with BMI associated with T2D status (case/control). Our model included age, sex, and the first two principal components of the entire genotype data of Baependi data as covariates. Due to some individuals missing genotype information for some SNPs, the tests of each gene had different sample sizes. Table 2 describes the sample sizes (number of subjects and number of families) for cases and controls included in analysis of each gene. Table 2. Summary of cases per subjects and families.

Gene
Subjects Families PPARG  845  83  928  43  42  85  FTO  712  71  783  47  38  85  CDKAL1  661  69  730  47  38  85 In Table 3, we report the p-values for each method. By comparing the p-values with respect to the corresponding α level, only the VCT identifies a significant GE interaction of BMI with PPARG. All other tests were non-significant for this gene as well as for other candidate genes. The variance estimates for families σ 2 and the ridge penalty σ 2 θ are reported in Table 3. In Section 4, we showed that the ridge penalty increases as the number of SNPs increases, however, our results for PPARG, FTO and CDKAL1 suggest that λ also depends on the number of subjects. Finally, Table 3 also shows the execution times using the R version 3.3.1 and a processor Intel(R) Core(TM) i5-6500 CPU @ 3.20 GHz with a RAM memory 8.00 GB and operating system 64-bits. Computation times for the VCT and the score test were considerably lower than those for the MinP test. The time to compute each test increased with the increase in number of SNPs in a gene and number of subjects in the analysis.

Conclusions
We have proposed a variance component score test for testing for interactions between a set of SNPs in a gene and an environmental variable with family data. We specified the interaction coefficients as random variables with common variance and evaluated the null hypothesis that the variance is equal to zero. Given the LD among some SNPs in a gene, we fit the null model assuming the SNP coefficients as random effects and showed that the corresponding BLUP was equivalent to ridge regression estimator. This approach gives a direct estimation for the ridge penalization parameter in comparison with other computationally demanding procedures based on cross validation [2]. We compared, via simulations and a real data application, our approach with the so called MinP test and also with the traditional q degrees of freedom score test. The results showed that the proposed test is robust and performs well, with considerable power. Simulations and application presented in this paper were done assuming a binary phenotype and a continuous environmental variable, however, the GLMM admits phenotypes with distribution belonging to the exponential family. These other distributions are currently unexplored in this paper. In addition, it is possible to have multiple environmental factors as well as environmental factors that are discrete. The proposed model can easily incorporate these cases. It is important to note that the proposed model does not account for a possible correlation of the environment among family members. A possible extension of this work is to include in the GLMM a shared household factor by adding a random effect that follows a normal distribution with mean vector zero and with variance the matrix that characterizes household sharing. Finally, using GLMMs can be computationally intensive and may experience convergence issues. In the future, we plan to explore using generalized estimating equations as an alternative approach to testing for interactions in families.

Appendix A. Ridge Regression and BLUP Equivalence in GLMM
Using the same notation of Section 2.2, where d 2 = Kb, the Equation (4) under the null hypothesis H 0 : γ = 0 becomes to g(µ d 2 ) =X β + Gθ + d 2 with µ d 2 = E(Y|d 2 ). Following Breslow and Clayton [1], the penalized quasi-likelihood (PQL) for this model is given by 2 withd 2 is choosen to maximize the sum of the last two terms, W R = diag ω ij / ν(µ d 2 ij )g (µ d 2 ij ) 2 and Subindex R denotes that the parameters are being estimated under Ridge regression method. We maximize the PQL with respect to β, θ and d 2 =d 2 jointly, obtaining the partial derivates (see Chen et al. [17] for details). Ridge regression estimator of θ is obtained by minimizing the function