Parametric and semi-parametric approaches modeling approaches

Many researchers have built upon the comforting regression framework in developing customized approaches to detect G-E interactions, including ordinary regression, penalized regression (Park and Hastie 2008) and logic regression (Schwender and Ruczinski 2010). In general, the joint effect of a genetic variant G and a given exposure E on a phenotype Y is often defined with the simple model:

where G is the number of allele (coded 0,1,2), E is continuous or categorical, Z represent a set of covariates one may adjust for, β are the linear effects of each component and g() the link function is the logit for dichotomous Y and the identity for quantitative Y. This model is a simplification, in that it ignores possible dominance effects. Still, just as the additive model has good power over a wide range of possible dominance models and has become the primary test statistic used in most GWAS (Lettre et al. 2007), the additive main and interaction effects will be detectably non-zero for a wide range of true dominance models, and the proportion of variance explained by the missing dominance effects will be quite small for most models.

Simplification is common in classical frequentist approaches, where adding degree of freedom can reduce statistical power. Or to quote the parallel from Kooperberg and Leblanc (Kooperberg and Leblanc 2008) with a cake: “if we want to divide the power over all possible interactions, nobody will get more than a crumb, and no-one will taste how good the cake is; we are better off dividing the cake among those people we believe to enjoy it.” For example a saturated linear model for a trichotomous E will have nine degree of freedom (df) compared to four df for the equation (1). In fact the same strategy has been used in most GWAs of marginal effect for the same reason.

It is important to note that even a simple model as equation (1) may encounter statistical issues. Especially, recent works from Tchetgen Tchetgen and Kraft (Tchetgen and Kraft 2011) have shown that when the main effect of continuous E, β_E in equation (1), is mis-specified, the likelihood ratio test, score test, and Wald test statistics of the main effect of G and the interaction effect can have incorrect type-1 error rates. This issue, which has been shown to be due to underestimation of the variance of β_GE, can be solved using different techniques (Cornelis et al. 2011): a) using a more flexible model for the environmental main effect (e.g., adding quadratic and cubic term for the exposure); b) using a robust “sandwich” estimator of the variance and c) modeling a continuous exposure by using general categorical variables.

A Bayesian framework gives the opportunity to make a step further in modeling the complexity of interaction effects. It provides a rational and quantitative way to consider a range of hypothesis in a single analysis. For example, Bayesian methods can be used to consider simultaneously multiple genetic models, some of them including diverse interaction effects, and to evaluate the posterior probability of each of these models (e.g., Crainiceanu et al. (2009) and Zhang and Liu (2007)). They also allow for multiple assumptions, which can be used to build composite estimators. If one wants to quantify the relevance of the G-E independence assumption (discussed in further sections), they offer solutions to trade off between bias and efficiency in a data adaptive way (Li and Conti 2009; Mukherjee et al. 2010). Finally, they allow incorporating biological information and knowledge accumulated in previous association studies, so that interaction effects can be weighted by their plausibility. However, despite their potential advantages, Bayesian approaches have been only sparsely used in genetic association studies and their advantages and limits from a modeling point of view need to be studied further. In particular, many hypothesized models are likely to be roughly equally consistent with the observed data for realistic sample sizes, making it difficult to infer which model provides the best fit: the cake will be split among so many people that nobody will get more than a crumb.

Screening for variants involved in interaction when interacting factor are unknown

Most genetic variants having effect through interactions with other risk factors are also likely to display marginal linear effect. For example, using random parameters for model (1) to simulate data —specifically, generating main effects and interaction effects independently of each other— will produce genetic variants with marginal effect almost 100% of time. This suggests one can simply test for marginal effect with power being almost only related to sample size, unless (as discussed below) the state of nature is such that most true models include interaction effects, but these are offset by the main effects so that the marginal genetic effects are quite small. This is especially useful if potential interacting factor are unmeasured or when interaction effects are expected to be difficult to assess.

Interaction models with small or no marginal genetic effects are theoretically possible(Culverhouse et al. 2002; Song et al. 2010). If such interactions are common then this will have significant consequences for how we go about searching for the genetic basis of complex phenotypes and will obviously limit the interest of screening for marginal effect. However such models have not yet been observed and confirmed in real data. This has led some to suggest that increasing sample size and testing for the marginal linear effect in agnostic GWAs scans might be the most powerful approach in most cases, while using more complex models might have only limited advantages (Clayton and McKeigue 2001; Hirschhorn and Daly 2005; Wang et al. 2005). The large success of this strategy in detecting genetic variants in GWAs has provided arguments in this direction, but the small amount of heritability explained by the “GWAs variants” is a potential rebuttal to the efficiency of this strategy.

When searching for quantitative trait loci (QTLs) an alternative for screening for the presence of interactions without using potential interacting factors is to test for homogeneity of variances across genotypic classes (Pare et al. 2010). The rationale is that, if the magnitude and the direction of the effect of a QTL differ depending on other genetic or non-genetic factors, the variability of the phenotypic outcome among individuals carrying the risk allele is likely to be larger than among the non-carrier. Hence, under the assumption that the main effect of the QTL affect neither the within genotype variance nor the between genotype variance, testing for heteroscedasticity will test for the presence of potential interactions. Note that heterogeneity of variances may be explained not only by the presence of interactions, but also by other biological mechanisms or other association patterns such as linkage disequilibrium with variants with large effect size (Takeuchi et al. 2010). Simulation studies have shown that the power of the test was limited when applying genome-wide significance threshold (Pare et al. 2010; Struchalin et al. 2010). It is also highly dependent on the main effect of the unknown interacting factors, having an optimal power for specific magnitude of main effect of E (Struchalin et al. 2010). Despite these limitations, testing for homogeneity of variance can be built in a two step approach. Since the power of the test, which highly depends on the main effect of the unknown interacting factors, is limited when applying genome-wide significance threshold, testing for homogeneity of variance has been proposed to preselect variants of interest that can be tested further for interactions. The potential of this approach has been recently demonstrated in a genome-wide association study of C-reactin and soluble ICAM-1 conducted in the Women’s Genome Health Study (Pare et al. 2010). Interestingly one of the identified GxE interaction was replicated in an independent study (Dehghan et al. 2011).

Leveraging interaction effect to improve detection of marginal effect

When a locus is expected to have residual marginal effects conditional on others factors tested for interaction, an efficient strategy is to use composite null hypothesis where both main effect and interaction effects are tested jointly (Kraft et al. 2007). Explicitly, testing the null hypothesis that the genetic variant has no effect in any strata or based on equation (1) H₀: β_G =0 or β_GE =0. This can be done by using a multivariate Wald test or a likelihood ratio test comparing a model including effect of E and Z only versus a model including effects of G, GE, E and Z. A simple alternative when exposure is binary or categorical is to test for marginal genetic effect in strata defined by exposure E. The joint test can then be computed as the sum of chi-squared for association derived from each stratum. Since the samples are independents, the sum follows a chi-square with the degree being equal to the of strata for E.

For case-control studies, the test for such joint effect can be performed using standard logistic regression, the more powerful retrospective likelihood approach (Chatterjee and Carroll 2005; Cornelis et al. 2011) can exploit an underlying gene-environment independence assumption or using the empirical Bayes approach (Chen et al. 2009a; Mukherjee and Chatterjee 2008) that can data adaptively relax the independence assumption. An extension from the family-based test for the joint test of gene main effect and G-E interaction (FBAT-J) has been recently proposed for dichotomous traits in trios and sibships (Hoffmann et al. 2009). The test assumes the genotype and the environment are independent conditional on the parental mating type. If the assumption does not hold, the test will have an inflated type I error rate (Weinberg and Umbach 2000).

By allowing for heterogeneous genetic effect among genetic or environmental strata one can maximize the statistical power to detect the locus while minimizing the loss of power when genetic effect is homogeneous. Simulation studies have shown that a joint test for a main genetic effect and interaction effect is likely to have higher statistical power than the marginal test or the standard one degree of freedom test in presence of moderate interaction effect or when interaction effect are in opposite direction to the main effect (Kraft et al. 2007). Conversely, in the presence of a small interaction effect, the marginal test may conserve the highest power.

Methods for meta-analysis of multiple parameters have been recently described so that estimates of effects from the joint test can also be combined across independent sample. In particular, Manning et al (2011) have described a general approach, while Aschard and colleagues (2011) have extended the aforementioned principle of analyzing sample stratified by environmental factors. The first approach should be used when analyzing quantitative exposures and in situations where the samples within each cohort have to be analyzed as a whole (e.g. in family data where one has to account for correlation among individuals). The second approach essentially offers practical advantage and it can be more flexible in situations where environmental categories may differ among the cohort analyzed. The first genome-wide application of the joint test has been published recently by Hazma et al (2011). They identified a new genetic variant associated with Parkinson’s disease and replicated the signal in independent samples.

As any test modeling interaction effect per se, the joint test is limited by the multiple testing issues in large scale data. Hence, it is only applicable in situations where there is a measured factor that might interact with the tested locus. Nevertheless some have shown that the joint test can be built in framework where multiple potential effect modifiers can be considered for a single locus. Strategy for testing can then be defined by averaging the effect of a given locus over other factors (Ferreira et al. 2007) or by testing the maximum joint test over a range of possible model (Chapman and Clayton 2007). It has been also suggested that degree-of-freedom for such joint tests can be reduced using Tukey style one-degrees-of-freedom model for interaction between groups of related genetic or/and environmental variables (Chapman and Clayton 2007; Chatterjee et al. 2006; Ciampa et al. 2011).

Testing for interaction per se

Besides TDT-like extension for G-E interaction as FBAT-I and its extension (Hoffmann et al. 2009; Lake and Laird 2004; Moerkerke et al. 2010) that are applicable to nuclear families data only, the traditional test for interaction consists in evaluating the term β_GE from equation (1). This test is relatively robust compared to many other approaches, although as described previously, misspecification of the main effect of a continuous E may increase type I error rate. The main concern when applying this simple test in GWAs data is its limited power (see “power” section above). Two types of strategies have been discussed to increase detection: a) to use multi-stage approaches to reduce multiple testing burden; and b) to leverage additional assumption on the data analyzed to improve efficiency.

Since the seminal paper from Marchini et al. (2005), multi-stage approaches using sequential test are considered as realistic approaches in GWAs. Even if not demonstrated, their work suggests that such strategy may improve the power of identifying interaction effects in GWAs. Since then diverse analysis strategies have been proposed, most of them focusing on the gene-gene interaction which face a strong multiple testing issues in GWAs. However these approaches can also be applied in the context of G-E. Examples include screening on genetic marginal effects (Kooperberg and Leblanc 2008; Macgregor and Khan 2006), or screening on a test that models the G-E association induced by an interaction in the combined case-control sample (Murcray et al. 2009). Simulation studies suggest that such approaches can be more powerful than traditional single-stage approach in which a huge penalty needs to be paid for multiple testing. Using a two-step strategy allows for less stringent thresholds of significance in the second step, since genetic markers have been prioritized in step one for their likely involvement in G-E interactions. While these methods became popular, questions have risen on how power and type 1 error are influenced by the correction among the two steps. While the two stages have been shown to be virtually independent in simulation study when screening on marginal effect (Kooperberg and Leblanc 2008; Marchini et al. 2005), recent work from Dai et al. (2010), provides proof of asymptotic independence of marginal association statistics and interaction statistics in linear regression, logistic regression, and Cox proportional hazard models when analyzing rare disease. Hence, in many situations the family-wise type I error rate might be controlled using classical Bonferroni correction for number of interaction tested at the second step only or by using permutation when markers considered at the second step are correlated.

Making assumption about the data analyzed to increase power of statistical test is a common principle. For binary trait such as disease status, the most popular one is the G-E independence assumption that allows testing for interaction in case-only data by testing for association between G and E among the cases using

Under the assumption of G-E independence in the whole population or G-E independence in controls for rare disease, testing for H₀: γ_E = 0 is equivalent to testing for H₀: β_GE = 0 from equation (1). When the assumption holds this method has the maximum power compared to most other approaches that leverage the G-E independence, except in the situation where the main effect of G or E is in opposite direction to the interaction effect (Mukherjee et al. 2011; Murcray et al. 2011). However it has also disadvantages: the main effect of G and E cannot be estimated and the type I error can be highly inflated when the assumption does not hold.

A range of other approaches have been proposed to leverage this assumption while providing a trade-off between increased power and controlled type I error rate (Chatterjee et al. 2005; Chen et al. 2009a; Cheng 2006; Mukherjee and Chatterjee 2008; Mukherjee et al. 2007). For example, when data on both cases and controls are available in a study, then one can be much more flexible than case-only analysis in studies of gene-environment interaction whether or not the independence assumption is valid. One can use a retrospective likelihood approach (Chatterjee et al. 2005) under the gene-environment independence assumption to obtain very efficient estimate all of the parameters of a general logistic regression model. On the other hand, if violation of the gene-environment independence assumption is suspected, one can perform data adaptive methods such as an empirical Bayes technique (Chen et al. 2009a; Mukherjee and Chatterjee 2008), which can be robust to violation of the independence assumption and yet can be more powerful than traditional case-control analysis when the independence assumption is valid. Other alternatives to the case-only test include multi-step approaches in a single sample (Gauderman et al. 2010; Murcray et al. 2009), multi-sample design (Chen et al. 2009b), and approaches that use Bayesian framework (Li and Conti 2009; Mukherjee et al. 2010). One should note that, based on recent reports, differences in performances between these methods only exist at the margin and they always depend on the type of model simulated (see Mukherjee and Chatterjee (2008) for a detailed comparison of several of these methods).

Exploratory or agnostic approaches

Traditional statistical methods such as multivariable linear or logistic regression are ill-equipped to incorporate all possible pairwise interactions among a large number of markers and exposures, let alone higher-order interactions. However, for complex diseases or traits the influence of non-linear or higher-order gene-gene and G-E interactions may be appreciable. Therefore, researchers are faced with difficult decisions to make their analysis practically feasible within computational and modeling restrictions (Maenner et al. 2009). The common alternative is to move away from the classical hypothesis testing framework and estimation of statistical significance level, and to use “model free” approaches or to adopt an agnostic approach to identify gene-environment interactions. Different analysis approaches from machine learning or data mining are needed to manage the high dimensionality of genome wide analysis studies and large scale data collections.

Interdisciplinary collaborations have led to the adoption of approaches from one community to another, especially in the field of gene-gene interactions. These include data segmentation methods (Tryon 1939), tree-based methods (Breiman et al. 1984), pattern recognition methods (Ripley 1996) and (non-)linear dimension reduction methods (Fodor 2002). A list of examples of these in the context of gene-gene interactions is given in Van Steen (Van Steen 2012). Unfortunately, the adoption of these methods in genome-wide based G-E interaction detection is not as “frequent” as it is genome-wide epistasis studies. In the following, we elaborate on two techniques that deserve more attention in the context of GWEI studies: Tree-based and Multifactor Dimensionality Reduction derived techniques.

Because the number of possible genetic model can be quite large, exploratory methods are often built on a trade-off and assume or favor some specific interaction models. Recursive partitioning approaches, such as Random Forests (Breiman et al. 1984; Schwarz et al. 2010) – a flexible and efficient data mining method based on regression or classification trees – also face such issues. Random Forests do not model interaction variables per se but they allow for interactions (or complex non-linear relationships) in the sense that they evaluate classification ability of particular combination of values taken by sets of predictor variables. Because of the independence assumption used during node splitting of “trees” these methods have been shown to have limited ability to detect pure interaction effects (McKinney et al. 2009). Notably, the recent SNPInterForest approach (Yoshida and Koike 2011) performed very well in successfully identifying pure epistatic interactions with high precision and was still more than capable of concurrently identifying multiple interactions under the existence of genetic heterogeneity. Hence, extensions that relax the independence assumption within a conditional inference framework (Hothorn et al. 2006) and improved procedures to extract interaction patterns from Random Forest (Yoshida and Koike 2011) make the Random Forest methodology particularly attractive for GWEI studies. Different variable importance measures have been proposed in the literature, including a joint importance measure which extends the idea of single importance to multiple importance and can be useful especially for interactions (Bureau et al. 2005). Note that correlated predictors and varying predictor categories or measurement scales are likely to exist in G-E studies and that care needs to be taken in the selection of the importance criterion. For instance, Strobl et al. (2008) identified the mechanisms causing the bias for permutation importance scores and developed a conditional variable importance which reflects the true impact of each predictor variable more reliable than the original permutation variable importance measure.

As an application example, Maenner et al. (Maenner et al. 2009) analyzed coronary heart disease cases from the Framingham Heart Study by firstly identifying influential SNPs for age of onset of early coronary heart disease using a random forest approach. Variable importance scores from a RF analysis provide measures to determine important SNPs and environmental exposures taking into account interactions without specifying a genetic model (Lunetta et al. 2004). Secondly, generalized estimating equations were used to evaluate the statistical significance of main effects and interactions of previously detected SNPs and smoking status (Maenner et al. 2009) (however note that such significance level should be take with caution since the selection at the “mining step” potentially overfits the data). The authors used a simple solution to handle family structure within their data by considering a binary family indicator as covariate for building the random forest. Similarly, Zhai et al. (2011) performed a 2-step approach with initial screening for SNPs associated to environmental measures by random forest and further analysis based on case-only logistic regression to obtain parameter estimates for the selected variables.

Tree-based methods might be an relevant alternative to logistic regression methods for identifying genes without strong marginal effects and of robustness to genetic heterogeneity where different subsets of genes can lead to a phenotype of interest (Lunetta et al. 2004). Random forests outperformed Fisher’s exact test when several risk SNPs interact (Lunetta et al. 2004) and behaved more robust when a high number of unassociated noise SNPs is present (Bureau et al. 2005). Another interesting approach combining regression models and tree-based methodology is a semi-parametric regression model, named partially linear tree-based regression model (PLTR) (Chen et al. 2007). The linear regression part of their model can control efficiently for confounders and provides the possibility to correct for linear main effects of variables so that a parsimonious summary of the joint effect of genetic and environmental variables is obtained.

Also non-parametric data mining methods such as Multifactor Dimensionality Reduction (MDR) (Ritchie et al. 2001), are the subjects of a trade-off. In contrast to logistic regression and random forests, MDR can be used to detect G-E interactions in the absence of any main effects. MDR can be applied to smaller sample sizes than logistic regression which needs enough observations to model all main and interactions effects. However, the “reduction” step consists in splitting the different combination of two variables (defined by E and G) in two groups of high risk versus low risk. This allows a range of model to be tested. But the interaction is summarized in a single binary parameter and is therefore unlikely to capture the full complexity of interactions (e.g., a gradient of effect across different combinations). Several extensions and variations of the MDR method have been proposed to address initial shortcomings of MDR (including the lack of correction for lower order effects, and the too stringent reduction into two risk groups). Model-Based MDR (MB-MDR) (Calle et al. 2010) and its extension to family data, family MDR (FAM-MDR) (Cattaert et al. 2010), enables adjustments for possible confounders and the handling of various phenotypes, e.g., continuous, categorical or censored. In particular, MB-MDR uses a reduction into a one-dimensional variable with three levels, i.e. high risk, no evidence, low risk, and potentially a continuum of risk groups (Calle et al. 2010; Cattaert et al. 2010). While comparing MB-MDR to MDR in the presence of noise, i.e. genotyping error, phenocopy and genetic heterogeneity, MB-MDR was found to have increased power in most situations, especially for genetic heterogeneity, phenocopies and minor allele frequencies. Previous to applying the MB-MDR method, FAM-MDR uses a preparation step where familial correlation free traits are obtained as residuals from a polygenic model (hence, hereby adjusting for potential population stratification). FAM-MDR outperformed Pedigree-based MDR (PGMDR) (Lou et al. 2008) in terms of handling multiple testing, empirical power and efficient use of available information from complex and extended pedigrees (Cattaert et al. 2010) and is therefore a promising alternative to the classical MDR derivatives to explore gene-environment interactions. One disadvantage of MDR is that its computational burden increases with the number of SNPs and the order of considered interactions. A parallel algorithm of MDR and MB-MDR has been implemented by Bush et al. (2006) and Van Lishout et al. (2011), respectively. Despite these efforts, filtering methods to preselect a subset of candidate factors and stochastic search algorithms (e.g., simulated annealing and evolutionary algorithms) are needed to assist researchers in the exhaustive search for interactions in genome-wide association studies. Knowledge about the pros and cons of these filtering approaches (as applied to genome-wide epistasis settings) will be most beneficial for GWEI studies and the availability of an entire exposeome.

Duell et al. (2008) compared MDR to focused interaction testing framework and logistic regression for identification of higher-order interaction effects in a case-control study using 26 polymorphisms and smoking as possible environmental risk factor. Little concordance existed between MDR and interaction testing framework with regard to the interaction factors. This finding may be caused by the different interaction modeling methodologies behind the approaches. The authors recommend using multiple approaches for data screening and analysis to detect potentially new risk factor combinations. More comparative studies are actually needed, examining differences between traditional (often regression-based) approaches with untraditional (often data-mining) methods in the context of GWEI studies. The study from Duell et al.(2008) also highlight the difficulties in computing a comprehensive significance level for exploratory methods. Overall, one should remind that there is no straightforward way to define a null hypothesis and to test it in these exploratory approaches. However strategies to statistically evaluate the significance of models obtained through data mining procedures are now discussed in the literature (e.g. (Pattin et al. 2009)) and more might be developed in the future.