Robust genetic interaction analysis

Robust genetic interaction analysis Abstract For the risk, progression, and response to treatment of many complex diseases, it has been increasingly recognized that genetic interactions (including gene–gene and gene–environment interactions) play important roles beyond the main genetic and environmental effects. In practical genetic interaction analyses, model mis-specification and outliers/contaminations in response variables and covariates are not uncommon, and demand robust analysis methods. Compared with their nonrobust counterparts, robust genetic interaction analysis methods are significantly less popular but are gaining attention fast. In this article, we provide a comprehensive review of robust genetic interaction analysis methods, on their methodologies and applications, for both marginal and joint analysis, and for addressing model mis-specification as well as outliers/contaminations in response variables and covariates. genetic interaction, robustness, model mis-specification, outlier/contamination Introduction Genetic interactions, including both gene–gene (G-G) interactions and gene–environment (G-E) interactions, have been shown to play crucial roles for the etiology, prognosis and response to treatment of many complex diseases beyond the main effects. In published studies, G factors include gene expressions, methylation, SNPs, as well as several other types of omics measurements. Beyond ‘real’ environmental variables, published studies have also examined demographic, clinical/medical, socioeconomic and other E factors. Many significant findings on genetic interactions have been made. For example, targeting the interaction between genes BIG3 and PHB2 has been suggested as a potential new treatment venue for ERα-positive breast cancer patients [1]. Researchers have found that the associations between ERCC2 polymorphisms and lung cancer risk are modified by the cumulative cigarette smoking [2]. A large number of statistical and machine learning methods have been developed for identifying important genetic interactions. Comprehensive discussions are available in [3–6], and many others. Take G-E interaction as an example. Most of the existing methods, although differ in technical details, share the following common strategy [7, 8]. Consider the model response∼f(G+E+G×E). In this model, the form of f is usually prespecified, for example, as logistic, linear and Cox for binary, continuous and survival responses, respectively; and linear covariate effects are assumed. Standard likelihood-based estimation is usually conducted, and conclusions on the importance of interactions are based on the estimates and significance levels of the G-E interaction terms. This strategy, despite many advantages, is nonrobust. In genetic interaction analysis, the demand for robustness arises in multiple aspects. Model mis-specification. There are at least two aspects. (a) Take binary responses, for example, disease/healthy status and response/resistance to treatment, as an example. Many of the published studies adopt the logistic regression model without sufficient justification. In ‘classic’ biomedical studies, it has been suggested that the logistic model assumption may fail, and alternative models, for example the probit and log models, may be needed [9, 10]. Genetic data and interactions are much more complex, and it is unlikely that they can be described using a single ‘default’ model. (b) The assumption of linear covariate effects may not always hold. For example, it has been shown that some E factors, such as age and body mass index (BMI), often have nonlinear effects [11]. A few publications have conducted exploratory analysis and suggested that some G effects may also deviate from linear [12]. Outliers/contaminations in predictors. E factors have been extensively studied in classic biomedical studies. It has been suggested that outliers/data contaminations are not uncommon, which can be caused by errors in data collection/recording as well as existence of influential observations [13]. For some types of G factors, for example gene expressions, outliers/contaminations can be caused by technical problems in sample preparation and profiling, human errors, as well as genetic abnormalities [14]. Outliers/contaminations in responses. Most if not all response variables studied in genetic studies have been investigated in classic biomedical studies, where the problem of outliers/contaminations in response has been well recognized [15]. For example, in prognosis studies, it has been recognized that human errors can be responsible for extremely long/short survivals as well as misclassification in the cause of death [16]. Specifically, consider the lung adenocarcinoma (LUAD) data collected by The Cancer Genome Atlas (TCGA) [17]. TCGA has generated comprehensive clinical, environmental and genetic data with high quality for multiple cancer types. The LUAD data have been analyzed in multiple published studies, including some interaction analyses [18, 19]. For the 504 subjects with complete mRNA expressions and survival time, there are six subjects with survival times 156.54, 162.98 163.99, 221.16, 232.00 and 238.11 months, respectively, while the rest 498 subjects have survival times ranging from 0.13 to 129.43. Some outliers are observed in response.A unique characteristic of genetic studies is that usually they cannot afford conducting strict subject selections, and the studied cohorts are often composed of multiple disease subtypes [20]. Different subtypes often behave differently in multiple aspects [21] and cannot be described using the same model. In this case, it can be viewed that data on the major subtype are ‘contaminated’ by the small subtypes. The possibility of model mis-specification and existence of outliers/contaminations in predictors and responses call for robust analysis methods. In low-dimensional biomedical studies, it has been well established that model mis-specification and failing to account for outliers and data contamination lead to biased estimation, false marker identification and inference and suboptimal models [22]. Robust analysis methods have been extensively developed and shown to have important scientific impact. With the high dimensionality and other unique features of genetic data, robust methods developed for low-dimensional biomedical data cannot be directly applied, and new development is needed. Our literature search and published reviews [23] suggest that, for the analysis of high-dimensional genetic data, the development of robust analysis methods is still much limited; however, it is gaining popularity fast. The majority of robust genetic data analysis methods focus on the main effect models, and the development on robust genetic interactions has been scattered. Considering the increasing importance of genetic interactions and demand for robust analysis methods in practical data analyses, our goal is to conduct a structured and systematic review of robust genetic interaction methods developed in recent literature. To the best of our knowledge, such a review is not available in the literature. This study targets to fill this knowledge gap. Methods In the analysis of genetic interactions (main effects as well), there are two generic analysis paradigm [3]. In marginal analysis, one or a small number of genes (or other omics measurement units) are analyzed at a time. In joint analysis, a large number of genes are analyzed under a single model. Although linked, the two analysis paradigms differ significantly in multiple aspects. Thus, below we review for each paradigm separately. In addition, the three robustness aspects, namely, model mis-specification, outliers/contaminations in predictors and outliers/contaminations in responses, are discussed separately. A summary of these methods is presented in Table 1. In addition, the following two considerations are taken into account. Table 1. Summary of the robust genetic analysis approaches (partial list) Approach  M/Ja  Robustb  G-E/G-G  Software  Applicationc  Conclusion  MDR [24, 25]  M  MisSpec  SNP-SNP  https://cran.r-project.org/ web/packages/MDR/ index.html  A sporadic breast cancer case-control data set, n=200,P=10  For the real data set, without any statistically significant main effects, MDR identifies one statistically significant high-order interaction among four SNPs from three different estrogen-metabolism genes, which is missed by the previous studies using logistic regression  W-test [26]  M  MisSpec  SNP-SNP  http://www2.ccrb.cuhk. edu.hk/wtest/download. html.  Two bipolar disorder GWAS data sets: (1) n=5000,P=414 682, (2) n=2420,P=729 304  For the two real data sets, W-test identifies significant interaction pairs that can be replicated. The genes in pairs play central roles in neurotransmission and synapse formation, and are undiscoverable by main effects with low frequency variants  RankCorr [27]  M  MisSpec  G-E    A lung cancer prognosis data set, n=336,P=2500 (22 283),q=5.  RankCorr identifies genes (SFRP1, IGFBP1, AHSG, AGR2 and others) and their interactions with important implications, which are different from those using alternatives. Most of the identified interactions are between genes and smoking status  3G-SPA [28]  M  MisSpec  G-G with SNPs  http://cran.r-project. org/package=SPA3G  (1) A birth weight data set, n = 1511 with 797 SNPs in 186 genes. (2) A yeast eQTL mapping data set, n = 112 with 6229 SNPs in 2956 genes  For the first data set, 23 G-G interactions are detected with the cutoff P-value 0.001, of which the genes include PDGFC, ANG, PTGER3, IL9 and others. For the second data set, 87 pairs are found to have significant interactions at an experiment-wise significant level of 0.05, where the corresponding genes include URA3, BAT2, LUE2 and others. For both data sets, the findings are shown to be biological sensible in the literature  R-GMDR [29]  M  Out-Res  SNP-SNP    A KARE data set with a HOMA-IR phenotype, n=8577,P=585 (327 872)  For the real data set, R-GMDR identifies a larger number of hub genes, which were reported in the literature than that from the ordinary GMDR method. For example, the KCNH1 gene has been demonstrated to have strong interaction effects with other genes on the function of insulin secretion  QR-based [30]  M  Out-Res  Gene/SNP-E    A lung cancer prognosis data set, n=246,P=2605 (22 283),q=5  For the real data set, some interesting findings different from the alternatives are observed. Specifically, with quantile = 0.5, the numbers of identified genes are 20 and corresponding interactions are 7 (with age), 7 (gender), 10 (smoking), 15 (chemotherapy treatment) and 22 (stage), respectively  SNPInter- Forest [31]  J  MisSpec  SNP-SNP  https://gwas. biosciencedbc. jp/SNPInterForest/ index.html  A GWAS data set of rheumatoid arthritis from the Wellcome Trust Case Control Consortium, n=3,499,P=10 000  For the real data set, SNPInterForest identifies important interactions with genes having some functions related to rheumatoid arthritis, such as PROK2 and PCDH15  BNN [32]  J  MisSpec  SNP-SNP  https://github.com/ beamandrew/BNN  A GWAS tuberculosis data set, n=104,P=16 925  The top five SNPs selected using BNN are shown located on the same chromosome and within 10–50 MB of loci previously reported as having a statistically significant association with pulmonary tuberculosis susceptibility  Nonpara- Scr [33]  J  MisSpec  SNP-SNP    A framingham heart study data set, n=977,P=349 985  Nonpara-Scr identifies 41 main effect SNPs and 36 SNP-SNP interactions in the real data set. The analysis indicates the importance of epistatic interactions in explaining BMI  SODA [34]  J  Out-Pre  G-G  https://cran.r-project.org/ web/packages/sodavis/ index.html  A Michigan lung cancer data set, n=86,P=5217  For the real data set, SODA selects two main effects and two interaction effects, which contribute substantially to classification accuracy compared with Lasso-Logistic  RankReg [35]  J  Out-Res  Gene/SNP-E  http://works.bepress. com/Shuangge/50/  A lung cancer prognosis data set, n=348,P=200 (22 283),q=2  In the analysis of a lung cancer study, RankReg identifies 11 genes with main effects, among which four and five genes have interactions with age and smoking status. A brief literature search shows that the identified genes may have important implications, including CYP2B6, CPS1, 1E2, HLA-DQA1 and others. Their better prediction performances are also observed  LAD-based [36]  J  Out-Res  Gene/SNP-E  https://github.com/cenwu.  A lung cancer prognosis data set, n=348,P=500 (22 283),q=5  In data analysis, the LAD-based method identifies 21 main gene effects and 10 G-E interactions, with seven G-age interactions, one G-smoke interaction and two G-stage interactions. The identified genes (PSPH, HMGA2, CHL1, PON3 and others) and their interactions with E factors are different from those of alternatives and are biologically sensible with better prediction performance  Approach  M/Ja  Robustb  G-E/G-G  Software  Applicationc  Conclusion  MDR [24, 25]  M  MisSpec  SNP-SNP  https://cran.r-project.org/ web/packages/MDR/ index.html  A sporadic breast cancer case-control data set, n=200,P=10  For the real data set, without any statistically significant main effects, MDR identifies one statistically significant high-order interaction among four SNPs from three different estrogen-metabolism genes, which is missed by the previous studies using logistic regression  W-test [26]  M  MisSpec  SNP-SNP  http://www2.ccrb.cuhk. edu.hk/wtest/download. html.  Two bipolar disorder GWAS data sets: (1) n=5000,P=414 682, (2) n=2420,P=729 304  For the two real data sets, W-test identifies significant interaction pairs that can be replicated. The genes in pairs play central roles in neurotransmission and synapse formation, and are undiscoverable by main effects with low frequency variants  RankCorr [27]  M  MisSpec  G-E    A lung cancer prognosis data set, n=336,P=2500 (22 283),q=5.  RankCorr identifies genes (SFRP1, IGFBP1, AHSG, AGR2 and others) and their interactions with important implications, which are different from those using alternatives. Most of the identified interactions are between genes and smoking status  3G-SPA [28]  M  MisSpec  G-G with SNPs  http://cran.r-project. org/package=SPA3G  (1) A birth weight data set, n = 1511 with 797 SNPs in 186 genes. (2) A yeast eQTL mapping data set, n = 112 with 6229 SNPs in 2956 genes  For the first data set, 23 G-G interactions are detected with the cutoff P-value 0.001, of which the genes include PDGFC, ANG, PTGER3, IL9 and others. For the second data set, 87 pairs are found to have significant interactions at an experiment-wise significant level of 0.05, where the corresponding genes include URA3, BAT2, LUE2 and others. For both data sets, the findings are shown to be biological sensible in the literature  R-GMDR [29]  M  Out-Res  SNP-SNP    A KARE data set with a HOMA-IR phenotype, n=8577,P=585 (327 872)  For the real data set, R-GMDR identifies a larger number of hub genes, which were reported in the literature than that from the ordinary GMDR method. For example, the KCNH1 gene has been demonstrated to have strong interaction effects with other genes on the function of insulin secretion  QR-based [30]  M  Out-Res  Gene/SNP-E    A lung cancer prognosis data set, n=246,P=2605 (22 283),q=5  For the real data set, some interesting findings different from the alternatives are observed. Specifically, with quantile = 0.5, the numbers of identified genes are 20 and corresponding interactions are 7 (with age), 7 (gender), 10 (smoking), 15 (chemotherapy treatment) and 22 (stage), respectively  SNPInter- Forest [31]  J  MisSpec  SNP-SNP  https://gwas. biosciencedbc. jp/SNPInterForest/ index.html  A GWAS data set of rheumatoid arthritis from the Wellcome Trust Case Control Consortium, n=3,499,P=10 000  For the real data set, SNPInterForest identifies important interactions with genes having some functions related to rheumatoid arthritis, such as PROK2 and PCDH15  BNN [32]  J  MisSpec  SNP-SNP  https://github.com/ beamandrew/BNN  A GWAS tuberculosis data set, n=104,P=16 925  The top five SNPs selected using BNN are shown located on the same chromosome and within 10–50 MB of loci previously reported as having a statistically significant association with pulmonary tuberculosis susceptibility  Nonpara- Scr [33]  J  MisSpec  SNP-SNP    A framingham heart study data set, n=977,P=349 985  Nonpara-Scr identifies 41 main effect SNPs and 36 SNP-SNP interactions in the real data set. The analysis indicates the importance of epistatic interactions in explaining BMI  SODA [34]  J  Out-Pre  G-G  https://cran.r-project.org/ web/packages/sodavis/ index.html  A Michigan lung cancer data set, n=86,P=5217  For the real data set, SODA selects two main effects and two interaction effects, which contribute substantially to classification accuracy compared with Lasso-Logistic  RankReg [35]  J  Out-Res  Gene/SNP-E  http://works.bepress. com/Shuangge/50/  A lung cancer prognosis data set, n=348,P=200 (22 283),q=2  In the analysis of a lung cancer study, RankReg identifies 11 genes with main effects, among which four and five genes have interactions with age and smoking status. A brief literature search shows that the identified genes may have important implications, including CYP2B6, CPS1, 1E2, HLA-DQA1 and others. Their better prediction performances are also observed  LAD-based [36]  J  Out-Res  Gene/SNP-E  https://github.com/cenwu.  A lung cancer prognosis data set, n=348,P=500 (22 283),q=5  In data analysis, the LAD-based method identifies 21 main gene effects and 10 G-E interactions, with seven G-age interactions, one G-smoke interaction and two G-stage interactions. The identified genes (PSPH, HMGA2, CHL1, PON3 and others) and their interactions with E factors are different from those of alternatives and are biologically sensible with better prediction performance  aM and J stand for marginal analysis and joint analysis. bMisSpec, Out-Pre and Out-Res stand for the robustness properties toward model mis-specification, outliers/contaminations in predictors and outliers/contaminations in responses, respectively. cn and p are the number of subjects and G factors in the application data set. For the G-E interaction analysis, q is the number of the E factors. For the value of P, the dimension after preprocessing is provided, together with the original dimension reported in parenthesis. Table 1. Summary of the robust genetic analysis approaches (partial list) Approach  M/Ja  Robustb  G-E/G-G  Software  Applicationc  Conclusion  MDR [24, 25]  M  MisSpec  SNP-SNP  https://cran.r-project.org/ web/packages/MDR/ index.html  A sporadic breast cancer case-control data set, n=200,P=10  For the real data set, without any statistically significant main effects, MDR identifies one statistically significant high-order interaction among four SNPs from three different estrogen-metabolism genes, which is missed by the previous studies using logistic regression  W-test [26]  M  MisSpec  SNP-SNP  http://www2.ccrb.cuhk. edu.hk/wtest/download. html.  Two bipolar disorder GWAS data sets: (1) n=5000,P=414 682, (2) n=2420,P=729 304  For the two real data sets, W-test identifies significant interaction pairs that can be replicated. The genes in pairs play central roles in neurotransmission and synapse formation, and are undiscoverable by main effects with low frequency variants  RankCorr [27]  M  MisSpec  G-E    A lung cancer prognosis data set, n=336,P=2500 (22 283),q=5.  RankCorr identifies genes (SFRP1, IGFBP1, AHSG, AGR2 and others) and their interactions with important implications, which are different from those using alternatives. Most of the identified interactions are between genes and smoking status  3G-SPA [28]  M  MisSpec  G-G with SNPs  http://cran.r-project. org/package=SPA3G  (1) A birth weight data set, n = 1511 with 797 SNPs in 186 genes. (2) A yeast eQTL mapping data set, n = 112 with 6229 SNPs in 2956 genes  For the first data set, 23 G-G interactions are detected with the cutoff P-value 0.001, of which the genes include PDGFC, ANG, PTGER3, IL9 and others. For the second data set, 87 pairs are found to have significant interactions at an experiment-wise significant level of 0.05, where the corresponding genes include URA3, BAT2, LUE2 and others. For both data sets, the findings are shown to be biological sensible in the literature  R-GMDR [29]  M  Out-Res  SNP-SNP    A KARE data set with a HOMA-IR phenotype, n=8577,P=585 (327 872)  For the real data set, R-GMDR identifies a larger number of hub genes, which were reported in the literature than that from the ordinary GMDR method. For example, the KCNH1 gene has been demonstrated to have strong interaction effects with other genes on the function of insulin secretion  QR-based [30]  M  Out-Res  Gene/SNP-E    A lung cancer prognosis data set, n=246,P=2605 (22 283),q=5  For the real data set, some interesting findings different from the alternatives are observed. Specifically, with quantile = 0.5, the numbers of identified genes are 20 and corresponding interactions are 7 (with age), 7 (gender), 10 (smoking), 15 (chemotherapy treatment) and 22 (stage), respectively  SNPInter- Forest [31]  J  MisSpec  SNP-SNP  https://gwas. biosciencedbc. jp/SNPInterForest/ index.html  A GWAS data set of rheumatoid arthritis from the Wellcome Trust Case Control Consortium, n=3,499,P=10 000  For the real data set, SNPInterForest identifies important interactions with genes having some functions related to rheumatoid arthritis, such as PROK2 and PCDH15  BNN [32]  J  MisSpec  SNP-SNP  https://github.com/ beamandrew/BNN  A GWAS tuberculosis data set, n=104,P=16 925  The top five SNPs selected using BNN are shown located on the same chromosome and within 10–50 MB of loci previously reported as having a statistically significant association with pulmonary tuberculosis susceptibility  Nonpara- Scr [33]  J  MisSpec  SNP-SNP    A framingham heart study data set, n=977,P=349 985  Nonpara-Scr identifies 41 main effect SNPs and 36 SNP-SNP interactions in the real data set. The analysis indicates the importance of epistatic interactions in explaining BMI  SODA [34]  J  Out-Pre  G-G  https://cran.r-project.org/ web/packages/sodavis/ index.html  A Michigan lung cancer data set, n=86,P=5217  For the real data set, SODA selects two main effects and two interaction effects, which contribute substantially to classification accuracy compared with Lasso-Logistic  RankReg [35]  J  Out-Res  Gene/SNP-E  http://works.bepress. com/Shuangge/50/  A lung cancer prognosis data set, n=348,P=200 (22 283),q=2  In the analysis of a lung cancer study, RankReg identifies 11 genes with main effects, among which four and five genes have interactions with age and smoking status. A brief literature search shows that the identified genes may have important implications, including CYP2B6, CPS1, 1E2, HLA-DQA1 and others. Their better prediction performances are also observed  LAD-based [36]  J  Out-Res  Gene/SNP-E  https://github.com/cenwu.  A lung cancer prognosis data set, n=348,P=500 (22 283),q=5  In data analysis, the LAD-based method identifies 21 main gene effects and 10 G-E interactions, with seven G-age interactions, one G-smoke interaction and two G-stage interactions. The identified genes (PSPH, HMGA2, CHL1, PON3 and others) and their interactions with E factors are different from those of alternatives and are biologically sensible with better prediction performance  Approach  M/Ja  Robustb  G-E/G-G  Software  Applicationc  Conclusion  MDR [24, 25]  M  MisSpec  SNP-SNP  https://cran.r-project.org/ web/packages/MDR/ index.html  A sporadic breast cancer case-control data set, n=200,P=10  For the real data set, without any statistically significant main effects, MDR identifies one statistically significant high-order interaction among four SNPs from three different estrogen-metabolism genes, which is missed by the previous studies using logistic regression  W-test [26]  M  MisSpec  SNP-SNP  http://www2.ccrb.cuhk. edu.hk/wtest/download. html.  Two bipolar disorder GWAS data sets: (1) n=5000,P=414 682, (2) n=2420,P=729 304  For the two real data sets, W-test identifies significant interaction pairs that can be replicated. The genes in pairs play central roles in neurotransmission and synapse formation, and are undiscoverable by main effects with low frequency variants  RankCorr [27]  M  MisSpec  G-E    A lung cancer prognosis data set, n=336,P=2500 (22 283),q=5.  RankCorr identifies genes (SFRP1, IGFBP1, AHSG, AGR2 and others) and their interactions with important implications, which are different from those using alternatives. Most of the identified interactions are between genes and smoking status  3G-SPA [28]  M  MisSpec  G-G with SNPs  http://cran.r-project. org/package=SPA3G  (1) A birth weight data set, n = 1511 with 797 SNPs in 186 genes. (2) A yeast eQTL mapping data set, n = 112 with 6229 SNPs in 2956 genes  For the first data set, 23 G-G interactions are detected with the cutoff P-value 0.001, of which the genes include PDGFC, ANG, PTGER3, IL9 and others. For the second data set, 87 pairs are found to have significant interactions at an experiment-wise significant level of 0.05, where the corresponding genes include URA3, BAT2, LUE2 and others. For both data sets, the findings are shown to be biological sensible in the literature  R-GMDR [29]  M  Out-Res  SNP-SNP    A KARE data set with a HOMA-IR phenotype, n=8577,P=585 (327 872)  For the real data set, R-GMDR identifies a larger number of hub genes, which were reported in the literature than that from the ordinary GMDR method. For example, the KCNH1 gene has been demonstrated to have strong interaction effects with other genes on the function of insulin secretion  QR-based [30]  M  Out-Res  Gene/SNP-E    A lung cancer prognosis data set, n=246,P=2605 (22 283),q=5  For the real data set, some interesting findings different from the alternatives are observed. Specifically, with quantile = 0.5, the numbers of identified genes are 20 and corresponding interactions are 7 (with age), 7 (gender), 10 (smoking), 15 (chemotherapy treatment) and 22 (stage), respectively  SNPInter- Forest [31]  J  MisSpec  SNP-SNP  https://gwas. biosciencedbc. jp/SNPInterForest/ index.html  A GWAS data set of rheumatoid arthritis from the Wellcome Trust Case Control Consortium, n=3,499,P=10 000  For the real data set, SNPInterForest identifies important interactions with genes having some functions related to rheumatoid arthritis, such as PROK2 and PCDH15  BNN [32]  J  MisSpec  SNP-SNP  https://github.com/ beamandrew/BNN  A GWAS tuberculosis data set, n=104,P=16 925  The top five SNPs selected using BNN are shown located on the same chromosome and within 10–50 MB of loci previously reported as having a statistically significant association with pulmonary tuberculosis susceptibility  Nonpara- Scr [33]  J  MisSpec  SNP-SNP    A framingham heart study data set, n=977,P=349 985  Nonpara-Scr identifies 41 main effect SNPs and 36 SNP-SNP interactions in the real data set. The analysis indicates the importance of epistatic interactions in explaining BMI  SODA [34]  J  Out-Pre  G-G  https://cran.r-project.org/ web/packages/sodavis/ index.html  A Michigan lung cancer data set, n=86,P=5217  For the real data set, SODA selects two main effects and two interaction effects, which contribute substantially to classification accuracy compared with Lasso-Logistic  RankReg [35]  J  Out-Res  Gene/SNP-E  http://works.bepress. com/Shuangge/50/  A lung cancer prognosis data set, n=348,P=200 (22 283),q=2  In the analysis of a lung cancer study, RankReg identifies 11 genes with main effects, among which four and five genes have interactions with age and smoking status. A brief literature search shows that the identified genes may have important implications, including CYP2B6, CPS1, 1E2, HLA-DQA1 and others. Their better prediction performances are also observed  LAD-based [36]  J  Out-Res  Gene/SNP-E  https://github.com/cenwu.  A lung cancer prognosis data set, n=348,P=500 (22 283),q=5  In data analysis, the LAD-based method identifies 21 main gene effects and 10 G-E interactions, with seven G-age interactions, one G-smoke interaction and two G-stage interactions. The identified genes (PSPH, HMGA2, CHL1, PON3 and others) and their interactions with E factors are different from those of alternatives and are biologically sensible with better prediction performance  aM and J stand for marginal analysis and joint analysis. bMisSpec, Out-Pre and Out-Res stand for the robustness properties toward model mis-specification, outliers/contaminations in predictors and outliers/contaminations in responses, respectively. cn and p are the number of subjects and G factors in the application data set. For the G-E interaction analysis, q is the number of the E factors. For the value of P, the dimension after preprocessing is provided, together with the original dimension reported in parenthesis. First, it is noted that different methods may be applicable to different types of G factors. This is mainly because of the distribution of G factor, which can be continuous or discrete. For example, SNPs are usually coded with three categories (0, 1 and 2) for genotypes (aa, Aa and AA), respectively. Therefore, some SNP interaction analysis methods are often based on contingency tables and statistical tests. On the other hand, gene expressions generally have continuous measurements, where the regression models are often adopted and heavy-tailed distributions are observed more frequently. To provide a deeper insight of each approach, we discuss their specific applicable types of G factors, which are also listed in Table 1. Second, despite similarities, there are still certain distinctions between models for G-E and G-G interaction analysis. E factors studied in most existing G-E interaction analysis are usually preselected with evidence from previous studies and have low dimensions. Thus, many methods always include E factors in the model without further selection. Compared with G-E interaction analysis, G-G interaction analysis often has much higher computational cost. It also has more ‘symmetrical’ form, as the interactions are between two G factors with similar distributions and other properties. These make many methods for identifying G-G interactions not directly applicable to analyze G-E interactions, and vice versa. The types of interactions studied in each model are demonstrated in the following section. Marginal analysis The most prominent feature of marginal analysis is that, each time, the number of variables (interactions + main effects) is considerably smaller than the sample size. Thus, standard estimation and inference approaches can be applied. In biomedical studies, marginal analysis is still highly popular because of its computational simplicity and stability. Robust to model mis-specification As mentioned in the first section and also discussed thoroughly in low-dimensional studies [37], with practical data, simply assuming a popular model (for example, logistic model for binary data) can be ‘dangerous’. With high-dimensional data, model diagnostics techniques, which may suggest model forms with low-dimensional data, have not been well developed. Thus, a strategy to achieve robustness to model mis-specification is to ‘avoid’ making specific model form assumptions. Three representative methods are reviewed below. Multifactor dimensionality reduction: Multifactor dimensionality reduction (MDR) is a model-free approach first proposed by [24] to identify genetic interactions for case-control studies. The most common application of MDR is in detecting SNP-SNP interactions and proceeds as follows. First, for each pair of genetic factors, their possible multifactor classes are represented in a two-dimensional space. Then, using the training data generated by V-fold cross-validation, each multifactor class is labeled as either ‘high-risk’ or ‘low-risk’ based on the ratio of the number of cases to the number of controls. Next, the prediction error of this model is estimated using the testing data. Finally, a permutation test is conducted to compare the average cross-validation consistency from the observed data to the distribution of average consistency under the null hypothesis of no association. The identification of important interactions can then be based on the P-values. The empirical studies demonstrate that MDR outperforms the parametric logistic regression. Software of MDR written in R is available at https://cran.r-project.org/web/packages/MDR/index.html [25]. Despite considerable successes, MDR is shown to be more computationally intensive than the generalized linear model, especially when the number of analyzed SNPs is >10 [24]. Many researchers have extended the framework of MDR to be more flexible and handle different types of outcomes. For example, [38] introduces a generalized MDR (GMDR) framework using the residual-based score of a generalized linear model, and this approach can accommodate dichotomous and continuous responses and has lead to some interesting findings [39]. The survival MDR (Surv-MDR) is proposed by [40] for survival data, which replaces the case-control ratios by the logrank test statistics. W-test and others: W-test is introduced by [26] to detect G-G interaction effects and measures the distributional difference between cases and controls through a combined log-odds ratio. It is a model free test in that it makes no assumption on the genotypic effect model and distribution. Similar to MDR, at the first stage, a 2 × 9 contingency table is generated to test the association for a pair of SNPs. For each cell of the table, the distribution is computed based on the number of subjects for the case and control groups, respectively. Then, a W-test statistic is developed to measure the discordance between the two distributions using the normalized log odds ratio adjusted by a data-adaptive factor. It follows a chi-squared distribution with the degree of freedom estimated from the covariance structure of the contingency table. Different from the logistic regression under which the odds ratio has a linear relationship with the interaction, the W-test approach can accommodate nonlinear interactions, as no specific genetic model is assumed. It is shown in simulation that it has robust power and reasonable type I error compared with logistic regression. The R code implementing this approach is available at http://www2.ccrb.cuhk.edu.hk/wtest/download.html. Using a laptop computer with 2.4 GHz CPU and 8 GB memory, W-test takes about 7.3 s to accomplish the analysis of one simulation with 1225 SNP-SNP interactions and 1000 subjects, compared with 7.7 s for chi-squared test and 77.7 s for MDR. It is observed that W-test has significant advantage over MDR in terms of computational cost. There are other robust tests for evaluating the associations of SNP-SNP interactions with a trait, such as the weighted rank test [41], mutual information-based test [42] and forward U-test [43]. They can be used to identify interactions without making stringent model assumption and are robust toward model mis-specification. Rank correlation-based methods [27]: The methods propose a robust rank correlation-based method (RankCorr) under a general semiparametric transformation model for the identification of G-E interactions. Strictly speaking, this method is not completely model-free. However, it makes much weaker assumptions on the model form and thus can be more robust than methods based on, for example, logistic or Cox models. In numerical studies, the continuous G factors are analyzed. However, this model can also accommodate discrete G factors. Denote y as the response variable, which can be a continuous marker, categorical disease status or survival time. Denote x=(x1,…,xp) as the p G factors and z=(z1,…,zq) as the q E factors. Assume n iid observations. For G factor j(=1,…,p), the following semiparametric transformation model is assumed:   E(gj(y))=∑k=1qzkαjk+xjβj+∑k=1qxjzkγjk=θj′wj, where gj(·) is a monotone increasing function with an unknown form, wj=(z,xj,zxj)′, αj=(αj1,…,αjk)′ and βj represent the main E and G effects, respectively, γj=(γj1,…,γjk)′ corresponds to G-E interactions, and θj=(αj′,βj,γj)′. A penalized loss function, O(θj)+ρ(θj), is proposed for parameter estimation, where O(θj) is the robust loss function, and ρ(θj)=∑k=12q+1ρ(θjk;λk) is the MCP penalty function with ρ(t;λ)=λ∫0|t|(1−xλξ)+dx. This method has two key differences from the existing alternatives. The first is that, without a model form assumption and hence likelihood function, a robust rank correlation-based loss function is constructed. For example, when y is continuous,   O(θj)=1n(n−1)∑i≠lI(yi≥yl)I(θj′wij≥θj′wlj), where I(·) is the indicator function. Robust loss functions for other types of responses can be constructed in a similar manner. The second difference is that a penalization approach is adopted for stabilizing estimates and selecting important effects. Simulation shows that this method outperforms the nonrobust model-based alternatives with more accurate identification, where it can achieve robustness to model mis-specification as well as outliers in response. A coordinate descent algorithm is developed for the optimization of RankCorr, where smoothed approximation and local linear approximation are adopted for the rank correlation objective function and MCP penalty, respectively. With fixed tunings, the analysis for one gene, five E factors and their interactions takes 0.46 s on a regular desktop computer. Kernel machine methods and others: [28] Kernel machine methods and others propose a model-based kernel machine method to identify significant G-G interactions for continuously valued quantitative traits. Different from the single SNP-based pairwise interaction analyses, this method takes genes as the testing units and identifies interactions at the gene level. The proposed gene-centric G-G Smoothing-sPline analysis of variance (ANOVA) model (3 G-SPA) can accommodate complex nonlinear SNP effects within a gene and imposes fewer assumptions on the underlying effects, making it less sensitive to model mis-specification. A two-step procedure is developed, which includes an exhaustive search for all pairwise interactions and assessment of significance of interactions for the identified gene pairs. In the first step, let yi and xi=(xi1,…,xiL) be the phenotype and genotype vector of the gene pair for subject i, where L=L1+L2 is the total number of SNPs in the two genes under investigation. A nonparametric regression model is adopted:   yi=m(xi)+εi,  i=1,2,…,n, where m(·) is an unknown function, and εi is a random error. This method is robust in that the form of m(·) does not need to be specified. The function m(·) is then decomposed into the main effects of the two genes and the interaction effect between them based on a functional ANOVA decomposition and reproducing kernels. In the second step, a score test statistic is developed to detect important interactions. Simulation shows that 3 G-SPA has better performance than a regression-based principle component (PC) analysis method with the linearity assumption. The R package SPA3G for implementing 3 G-SPA is available at http://cran.r-project.org/package=SPA3G. For 1000 simulation replicates, which include 500 subjects and two separate genes with 10 SNPs each, the computational time for 3 G-SPA and two PC-based approaches partial PCA (pPCA) and full PCA (fPCA) are 28, 24.5 and 24 min, respectively. 3 G-SPA has been extended in [44] to analyze case-control data based on a mixed effects logistic model. In addition, a multiple-kernel method is introduced in [45] to assess G-E interactions and can accommodate various trait types, such as continuous, binary and survival. Besides the kernel methods, other nonparametric techniques have also been developed for the marginal analysis of genetic interactions, including, for example, the multivariate adaptive regression splines [46] and regularized isotonic regression [47]. Robust to outliers/contaminations in predictors As described in ‘Introduction’ section and in the literature, outliers and contaminations in the predictors can happen to both E and G factors. However, our limited literature search suggests that there is still a lack of genetic interaction analysis methods tailored to accommodate outliers and contaminations in predictors. With the link between marginal and joint analyses, it is possible that methods developed for joint analysis can be ‘marginalized’ for marginal analysis. In addition, a few studies might have potential robustness properties toward outliers and contaminations in predictors but have not been well examined. Specifically, the L-estimator GMDR is proposed by [29] with the goal of accommodating outliers in response, and is developed based on least trimmed squared regression (LTS). LTS has been shown to be robust to outliers in both predictor and response spaces in the literature [48]. Thus, it is reasonable to conjecture that the L-estimator GMDR may also achieve robustness toward outliers in predictors. More theoretical and numerical analyses are needed to examine this robustness. The details of this approach are discussed in next subsection. Robust to outliers/contaminations in responses Robust GMDR [29]: It develops a robust GMDR (R-GMDR) approach for SNP-SNP interaction analysis based on the L-estimator and M-estimator estimation methods. The proposed method can reduce the effects caused by outlying responses. Specifically, for the L-estimator GMDR, the LTS regression, as opposed to the least squared regression, is conducted for generalized linear models. A studentized trimmed residual is proposed as the score for GMDR, where a threshold value is used to determine whether the ith sample should be regarded as an outlier and removed. For the M-estimator GMDR, an M-estimator type score based on the Tukey’s biweight function is developed, where a threshold T is used to shrink all residuals exceeding T. In extensive simulations where the response follows a skewed distribution or a normal distribution with a huge variance, both the L-estimator GMDR and M-estimator GMDR are shown to outperform the original GMDR. The original GMDR software is available at http://www.ssg.uab.edu/gmdr/. It is expected that the two robust methods can be implemented by modifying the original GMDR software. In addition, it is demonstrated that M-estimator GMDR has similar computational cost with the ordinary GMDR, and L-estimator GMDR takes slightly more computational time since an additional step to compute residuals is needed. Researchers have also proposed other improvements of MDR-related methods to accommodate long-tailed errors and outliers in responses. For example, for the accelerated failure time (AFT)-MDR designed for survival phenotypes, two improvements are introduced to reduce the effects of extreme values on the sum of standardized residuals in the AFT-MDR [49]. The first approach transforms the continuous standardized residuals into binary variables and then uses the original MDR instead of AFT-MDR. Under the second approach, the extreme values of the standardized residuals beyond the specified bounds are replaced by either the lower or upper bound. Both methods are demonstrated to have higher powers than the original AFT-MDR with the presence of outliers. Quantile regression and others [30]: Quantile regression and others propose a quantile regression-based (QR-based) method to identify G-E interactions for prognosis responses. As with other QR-based methods [50], it is robust to contamination in the responses. Both continuous and discrete G factors are analyzed in numerical studies, mimicking gene expression and SNP data, respectively. Here, we use notations similar to those for the rank correlation-based method reviewed above. Denote y as the survival time. The following check loss is adopted for G factor j(=1,…,p):   ρτ(y−∑k=1qzkαjk−xjβj−∑k=1qxjzkγjk), where the function ρτ=u{τ−I(u<0)} with 0<τ<1. Based on this loss function, two weighted methods are proposed to accommodate right censoring of the survival time. Important interactions can be identified based on significance level as well as using a penalization approach. It is observed in simulation that under data contamination and heavy-tailed errors, the robust methods can significantly outperform the nonrobust ones. Two techniques are involved in the computation process, including the majorize–minimization (MM) approach and local quadratic approximation. With more complex algorithm, the QR-based methods tend to be slower than the nonrobust methods, such as that based on least squared loss. In the literature [51, 52], the robust check loss function has also been used in the framework of regression tree for interaction analysis, and such an approach can be potentially adapted for genetic interaction analysis. The above ‘robust loss function + penalized identification’ strategy has also been adopted with other robust loss functions. For instance, [19] adopts an exponential squared loss for the G-E interaction analysis with prognosis data. This approach can down-weigh the influence of observations with larger residuals. Joint analysis Compared with marginal analysis, joint analysis may better reflect the fact that molecular mechanisms under complex diseases involves multiple genes (and their interactions). When a large number of interactions and main effects are jointly analyzed, because of the high dimensionality, complex statistical methods are usually needed. In recent statistical literature, joint analysis has been popular [53, 54]. In biomedical studies, its popularity has been increasing fast. Recently, the ‘main effects, interactions’ hierarchical constraint is often used in joint analysis, where an interaction is allowed into the model only if the corresponding main effects are also identified [53–56]. Although this constraint is not strictly validated in biomedical studies, it is suggested to achieve better estimation and interpretation [53]. More detailed discussions are available in the literature [3, 57]. For the following approaches, their properties of respecting hierarchy are discussed briefly. Robust to model mis-specification Random forest: Random forest is an ensemble of regression trees and can predict an outcome based on a large number of predictors. A novel method SNPInterForest is proposed by [31] for detecting G-G interactions by modifying the construction of the random forest. It is a model-free approach and does not require (semi)parametric model specification. Specifically, in the first step, a random forest is constructed using case-control data with SNPs as categorical variables. Different from the original random forest under which variables are selected based on individual effects, a combination of multiple SNPs as well as a single SNP is considered when choosing a split variable at each node. Then, the interactions of variables are represented within branches of trees, which are the paths from the root nodes to the leaf nodes. Next, for each pair of SNPs, its interaction strength is measured by the number of simultaneous appearance in the same branch over all trees. That is, an SNP-SNP interaction is relatively more important in affecting the disease outcome if it appears more frequently in the same branch. Finally, the above interaction score is normalized by using its respective baseline level, which is based on the hypothesis that the SNP combinations do not interact. The random forest methods usually do not consider the ‘main effects, interactions’ hierarchy and are even more focusing on interactions without main effects. In extensive simulations under a wide spectrum of disease models, SNPInterForest is shown to achieve better performance than the model-based approach BOOST. Software for SNPInterForest is available at https://gwas.biosciencedbc.jp/SNPInterForest/index.html. For the analysis of one real GWAS data with 3499 individuals and 10 000 SNPs, the running time of SNPInterForest is about 98 h using a computer with 2.67 GHz CPU and 6 GB memory, compared with 11 min for BOOST and 5 h for SNPHarvester. There are several other genetic interaction studies that are also under the framework of random forest [58, 59]. Different approaches may have minor difference in the construction of forest and definition of interaction score. Neural network: Researchers have developed a series of neural network methods to model and detect the functional relationships between genetic interactions and response variables. Neural network is a model-free approach and is capable of dealing with a large amount of data. It has been widely applied to identify SNP-SNP and SNP-E interactions. Similar to random forest, the neural network approaches also usually focus on identifying interactions with weak main effects or without main effects and do not take hierarchical constraints into account. In general, the underlying structure of the neural network for genetic interaction analysis is a weighted and directed graph, where the nodes represent genetic factors and the edges represent genetic interactions. The graph is divided into multiple layers. The input layer contains all genetic factors under consideration, and each node is connected to one or more nodes (genetic factors) in a hidden layer. Nodes in the hidden layer are further connected to either additional hidden layers or to an output node representing the response. Signals pass from the input layer to the output layer through the interactions with weights and activation functions. Different algorithms are proposed to learn the weights and activation functions in genetic interaction studies, including for example the genetic programming neural network [60], grammatical evolution neural networks [61], feed-forward multilayer perceptron [62] and Bayesian neural network (BNN) [32]. Simulations in [62] and [63] show that the neural network is more accurate in identifying interactions than the nonrobust logistic regression model. There are several publicly available softwares for the above approaches. For example, the code implementing BNN is available at https://github.com/beamandrew/BNN. Compared with MDR, BNN is shown to be able to analyze a larger number of SNPs in a relatively short time. Specifically, the analysis of a data set with 104 subjects and 500 000 SNPs can be accomplished within 30 min. Nonparametric screening [33]: It proposes a nonparametric joint modeling method for detecting SNP-SNP interactions based on variable screening and variable selection (Nonpara-Scr). The associations between the phenotype and genotypes are characterized by an unknown function. Robustness is achieved by not assuming the specific model form. Use notations similar to those above. This approach considers the following generic nonparametric model:   yi=μ+∑j=1pfj(xij)+∑j=k+1p∑k=1phjk(xijxik)+εi, where μ is the population mean, fj(·)’s and hjk(·)’s are the unknown functions representing the main effects and interactions, respectively, and εi is the random error. The basis expansion technique is used to approximate the nonparametric functions, that is:   fj(u)=∑lcjlϕl(u),  hjk(u)=∑ldjklϕl(u), (1) where ϕl(·) is the lth basis function, and cjl and djkl are the corresponding coefficients. A hybrid of variable screening and variable selection is proposed for interaction detection. First, consider the model with only one main effect:   yi=μ+∑lcjlϕl(xij)+εi,  for j=1,…,p. The subset of important main effects is then defined as M^1={1≤j≤p:∑lc^jl2>C} with a threshold C. Next, the following model is considered to screen interactions:   yi=μ+∑j∈M^1p∑lcjlϕl(xij)+∑j=1p∑k∈M^1p∑ldjklϕl(xijxik)+εi. The subset of selected interactions is estimated by M^2={1≤j≤p,k∈M^1:∑ld^jkl2>C′} for some C′>0. A variable selection method based on LASSO is further used for the identification of final interactions following the above screening procedure. In the variable screening step, the weak hierarchy is taken into account, where the interactions can be included in the model if at least one of the corresponding main effects is selected. However, Lasso is then adopted directly in the variable selection step and cannot guarantee the hierarchical structure of ‘main effects, interactions’. In simulation studies, it is shown that the proposed method has an outstanding performance in selecting important interactions compared with parametric models, such as the forward LASSO. In addition, the proposed procedure is demonstrated to be computationally efficient. Robust to outliers/contaminations in predictors Stepwise cOnditional likelihood variable selection for Discriminant Analysis [34]: Stepwise cOnditional likelihood variable selection for Discriminant Analysis (SODA) proposes SODA to detect both main and quadratic interaction effects. The continuous mRNA expression measurements are analyzed in this study. This method does not require the normality (or sub-Gaussian) assumption on predictors, which is usually made in other studies, leading to much enhanced robustness towards long-tailed distributed predictors and outliers. The SODA approach consists of three stages: a preliminary forward main effect selection, a forward variable selection (for both main and interaction effects) and a backward elimination. A key point of the above procedure is that the extended Bayesian information criterion (EBIC) is used as the variable selection/elimination criterion at each iteration. Specifically, for the K-class classification problem, let y∈{1,…,K} denote the class label and x be the p-dimensional predictor vector. The following logistic model is adopted:   P(y=k|x,θ)= exp ⁡[δk(x|θ)]1+∑l=1K−1 exp ⁡[δl(x|θ)], (2) where δk(x|θ) is the discriminant function for class k, and θ denotes the parameter vector. The function δk(x|θ) is further defined as:   δk(x|θ)={αk+βk′x+x′Akx,k=1,…,K−1,0,k=K. (3) This method does not need to model the distribution of x and includes the quadratic discriminant analysis as a special case, under which the predictors are required to follow a multivariate normal distribution. At each iteration with the considered predictor set S, the EBIC is estimated based on the log-likelihood function derived from equation (2) with the predictors in S. The proposed method does not respect the hierarchical constraints. In simulations with multiple predictor distributions, the proposed method is observed to achieve higher identification and classification accuracy, compared with the nonrobust methods including the backward procedure and forward–backward method. The R package sodavis for SODA is available at https://cran.r-project.org/web/packages/sodavis/index.html. To analyze one simulated data set with 1000 subjects and 1000 genes, SODA takes about 4 min, the backward procedure takes 22 min and Hier [53] takes >24 h. Robust to outliers/contaminations in responses A rank-based penalization approach: In [35] for the identification of G-E interactions, a rank-based regression (RankReg) is adopted, which does not make stringent distributional assumptions on the random error and can accommodate long-tailed errors and hence outliers/contaminations in the response. With the regression framework, this method can be applied to analyze both continuous and discrete G factors. For subject i, let yi be the response, xi=(xi0,xi1,…,xip) be the (p+1)-dimensional vector of G factors with xi0=1 and Ui and Ei be the continuous and discrete E factors, respectively. The following partially linear varying coefficient (VC) model is adopted for joint analysis:   yi=∑j=0pαk(Ui)xij+∑j=0pβjxijEi+εi, where αk(·) is the smooth varying-coefficient (VC) function and εi is the random error. The VC function is approximated using a similar basis expansion technique as in equation (1). Beyond being robust to outliers/contaminations, the proposed method can also accommodate unspecified nonlinear effects of E factors, and thus, has a second robustness property. The rescaled rank-based loss function   1n∑i<l|εi−εl| is considered for estimating the unknown functions (coefficients). For regularized estimation and more importantly selection of important effects, a penalization approach is adopted. The proposed method does not respect the ‘main effects, interactions’ hierarchy because of computation consideration. However, it is demonstrated that one can first apply the proposed method and then add back the main effects if the hierarchy is needed. Simulation shows that the proposed method has satisfactory performance and outperforms several competing alternatives including the nonrobust approach based on the ordinary least squared loss function. The computer code written in R is available at http://works.bepress.com/Shuangge/50/. The proposed objective function is optimized using the group coordinate descent algorithm. The analysis of one simulated replicate with n = 200, P = 200 and two E factors takes about 1.6 min using a regular laptop. It is also shown that the computational cost of proposed method grows linearly with n and p. A least absolute deviation based method: The least absolute deviation (LAD) loss is a special case of the check loss and well known for its robustness to long-tailed errors and outliers in responses. This method [36] proposes a penalized robust approach to dissect G-E interactions based on the LAD loss for prognosis data. It has been applied to both continuous and categorical predictors, corresponding to gene expression and SNP data. Use notations similar to those above. Let T be the logarithm of survival time. The AFT model is assumed under joint modeling, that is:   T=∑k=1qαkzk+∑j=1pβjxj+∑j=1p∑k=1qxjzkγjk+ε, where ε is the random error with an unknown distribution function. It is especially noted that the random error can have a long-tailed distribution or be contaminated. Based on the Kaplan–Meier weights wi(i=1,2,…,n) for accommodating the right censoring of survival time, the following weighted LAD loss function is proposed:   ∑i=1nwi|yi−∑k=1qαkzik−∑j=1pβjxij−∑j=1p∑k=1qxijzikγjk|, where yi=min⁡(Ti,Ci), and Ci is the log censoring time. A two-stage penalty proposed in [64] is further added for the identification of interactions respecting the strong ‘main effects, interactions’ hierarchy. Simulation demonstrates the advantage of the proposed approach over the alternatives without the robustness property, which are based on least squared loss. Research code written in R is publicly available at https://github.com/cenwu. A group coordinate descent algorithm is developed for the proposed method. For one simulated replicate with 200 subjects, 200 G factors and 10 E factors, the proposed LAD-based method takes 12.8 s to accomplish the analysis given fixed tunings, compared with 56.5, 57.4 and 57.3 s taken by the other three alternatives. Some insights into the numerical performance of different methods The implementation of all the reviewed methods is a prohibitive task and beyond the scope of this review. As there are not publicly available softwares for many approaches, strictly numerical comparisons of different methods cannot be conducted. Following the literature [23], we summarize the simulation settings and results reported in some of the reviewed papers in Table 2, to provide some insights into their numerical performance and partly show their advantages and shortcomings. It is observed that the performance of these methods depends on the underlying genetic models, the distributions of predictors and response (error), data dimensionality, sample size and others. Among them, no method can perform consistently better than others under all scenarios. In general, the robust methods have relatively poor or competitive performance compared with the nonrobust ones when the models are not mis-specified and the data do not have contaminations and outliers. This result is as expected, as the robustness can be achieved usually at the cost of some loss in efficiency [65]. When there are model mis-specification and data contamination, the robust interaction analysis methods mostly have significant advantages. Table 2. Summary of comparisons of the robust genetic analysis approaches against alternatives (partial list) References  Proposed/alternative  Main settings  Major conclusions  Wang et al. [26]  W-test/logistic regression; chi-squared test; MDR  n=1000,P=50; Discrete predictors: six different genetic architectures with minor allele frequency (MAF) in the common range and in the low frequency range, and LD in the high range, medium range and low range; binary response: (1) d=(4,2)a, linear interaction model with main effect, (2) d=(0,2), nonlinear interaction model without main effect  (1) The model-free W-test and MDR outperform the logistic regression even under the linear model, with higher power and lower type I error rate  (2) Their advantages become more prominent when the underlying model is nonlinear. For one genetic architecture, the power is 5.9% (logistic), 72.6% (Chi-squared), 95.5% (MDR) and 88% (W-test). W-test behaves much better than MDR when 1%<MAF<5%, where the average power of W-test and MDR are 87.6 and 73.2%, respectively  Shi et al. [27]  RankCorr/P-value-based significance analysis (Sig); (1) Penalized least squares (PeLS), (2) Penalized Cox model (PeCOX)  n=200,P=500 1000,q=5, d=(3,4,20)b; Continuous predictors: follow multivariate normal distribution with two correlation structures and two correlations; (1) Continuous response under linear model, (2) Censored response under transformation model: g(t)=t3,t5,t and t1/3; Error distributions: N(0, 1), 0.7N(0,1)+0.3Cauchy(0,1) and t(1)  (1) When the error has a normal distribution, Sig has the best performance. However, when the error is contaminated or has a long tail, RankCorr can significantly outperform the alternatives (2) RankCorr has the largest AUC under all simulation scenarios, including those with linear transformation function and normal distribution error, followed by the nonrobust alternatives PeCOX and Sig  Li and Cui [28]  3G-SPA/Partial PCA (pPCA); Full PCA (fPCA)  n = 500, 1000, two separate genes with 10 SNPs each; discrete predictors: generated via MS program; Continuous response under linear model: (1) without genetic effects, (2) with main effects but no interaction, (3) with interaction effect smaller than main effects, (4) with interaction effect larger than main effects  3G-SPA has the highest power with the type I error reasonably controlled. The PCA-based methods are more conservative. The power of 3G-SPA is not depending on whether the genetic effect is because of its main effect or interaction, whereas the PCA-based methods prefer the scenarios with larger interactions  Kim and Park [29]  L-estimator GMDR; M-estimator GMDR/Ordinary GMDR  (1) n=1000,P=100, (2) n=3000,P=1000, d=(2,1); Discrete predictors: follow multinomial distribution with two different probabilities; Continuous response under linear model; Error distributions: (1−τ)N(0,1)+τN(0,400) and (1−τ)N(0,1)+τGamma(α,β) with τ=0,0.033 and 0.067  When there are no outliers (τ = 0), the ordinary GMDR and M-estimator GMDR outperform L- estimator GMDR. For all the designs with outliers, the performance of L-estimator and M-estimator GMDR are better than that of ordinary GMDR, with higher power. In some designs, L-estimator GMDR outperforms M-estimator GMDR  Yoshida and Koike [31]  SNPInterForest/SNPHarvester; BOOST  (1) n=4000,P=1000; Discrete predictors: the genetic model with different values of LD and MAFs; Binary response: the additive, multiplicative, and threshold models of interactions with weak marginal effects (2) n=400,P=1000; Discrete predictors: same as (1); Binary response: the pure interaction model without marginal effects  (1) The performance of SNPInterForest is superior to two alternatives in terms of their power SNPInterForest and SNPHarvester do not identify any spurious associations, whereas BOOST falsely identifies five interactions (2) The SNPInterForest and BOOST outperform SNPHarvester, with higher power and better identification accuracy The false discovery rates of SNPInterForest and SNPHarvester are much lower than that of BOOST  Beam et al. [32]  (1) BNN/BEAM; χ2 test (2)BNN/BEAM; χ2 test; GBM; MDR  n=2000,P=1000; Discrete predictors: the genetic model with different values of MAF; (1) Binary response: three biallelic models with main effects, (2) Binary response: pure interaction models without main effects  (1) BNN is uniformly more powerful than the χ2 test under all the simulated genetic models, and it outperforms BEAM in most instances (2) BNN performs better than BEAM, χ2 test and GBM under most of the simulations, and has similar performance with MDR. But MDR is much more time-consuming  Li et al. [33]  Nonpara-Scr/Screen and Clean (SC); SIS+LASSO; Forward LASSO; MDR  n=1000, P=3512, d=(12,8); Discrete predictors: dichotomizing the continuous variables generated from normal distribution with different correlations; Continuous response: a genetic model with both additive effects and dominant effects, considering three different noise levels  Nonpara-Scrbehaves better than the other four alternatives, with higher power and much smaller false-positive rates. Their average power are 97% (Nonpara-Scr), 90% (SIS+LASSO), 80% (SC), 48% (Forward LASSO) and 64% (MDR)  Li and Liu [34]  SODA/ Backward procedure (Back); Forward–backward method (ForBack); hierNet; IIS-SQDA  Ten different sample sizes with n=100, 100×101/9,…,1000; (1) P=50,d=(1,4); Continuous predictors: generated from normal, non-normal and heteroskedastic distributions, (2) P=1000,d=(1,4); Continuous predictors: non-normal distribution (3) P=50,d=(0,4), Continuous predictors: same as (2), (4) P=50,d=(2,4), Continuous predictors: same as (2); Binary response: Bayes classification rule  With Gaussian predictors, SODA, Back and ForBack can detect all important predictors as n increases, with small false positives, whereas IIS-SQDA and hierNet select too many false positives, resulting in high test error rates. Under other predictor settings, SODA performs much better than the nonrobust alternatives and achieves near-oracle classification accuracy when sample size is large enough  Wu et al. [35]  RankReg/LSregression-based method; Robust marginal analysis (RobMA)  n=200, P=200, q=2, d=(0,0,20); (1) Discrete predictors: MS program and Hardy–Weinberg equilibrium-based method, (2) Continuous predictors: follow multivariate normal distribution with an auto-regressive correlation structure; Continuous response under the linear model; Error distributions: N(0,1),t(3),0.85N(0,1)+0.15 log Normal(0,1) and 0.85N(0,1)+0.15Cauchy(0,1)  With N(0, 1) error distribution, RankReg is slightly inferior to the unrobust alternative LS, but still outperforms RobMA. With other three error types, RankReg has significant advantages, with more true positives, much fewer false positives, and higher estimation accuracy  Wu et al. [36]  LAD-hier/LAD-nhier; LAD-group; LS-hier; LS-nhier; LS-group  n=200, P=200, q=10, d=(10, 10, 15); (1) Continuous predictors: follow multivariate normal distribution with the banded and autoregressive correlation structure, (2) Discrete predictors: dichotomizing the above continuous values at the first and third quartiles; Censored response under the AFT model; Error distributions: N(0,1),0.80N(0,1)+0.20 log Normal(0,1),t(2) and 0.80N(0,1)+0.20Cauchy(0,1)  Without data contamination, the robust LAD-hier and the nonrobust counterpart LS-hier have comparable performances. With heavy-tailed errors, LAD-hier is significantly superior to the alternatives, in terms of both identification and prediction performance  References  Proposed/alternative  Main settings  Major conclusions  Wang et al. [26]  W-test/logistic regression; chi-squared test; MDR  n=1000,P=50; Discrete predictors: six different genetic architectures with minor allele frequency (MAF) in the common range and in the low frequency range, and LD in the high range, medium range and low range; binary response: (1) d=(4,2)a, linear interaction model with main effect, (2) d=(0,2), nonlinear interaction model without main effect  (1) The model-free W-test and MDR outperform the logistic regression even under the linear model, with higher power and lower type I error rate  (2) Their advantages become more prominent when the underlying model is nonlinear. For one genetic architecture, the power is 5.9% (logistic), 72.6% (Chi-squared), 95.5% (MDR) and 88% (W-test). W-test behaves much better than MDR when 1%<MAF<5%, where the average power of W-test and MDR are 87.6 and 73.2%, respectively  Shi et al. [27]  RankCorr/P-value-based significance analysis (Sig); (1) Penalized least squares (PeLS), (2) Penalized Cox model (PeCOX)  n=200,P=500 1000,q=5, d=(3,4,20)b; Continuous predictors: follow multivariate normal distribution with two correlation structures and two correlations; (1) Continuous response under linear model, (2) Censored response under transformation model: g(t)=t3,t5,t and t1/3; Error distributions: N(0, 1), 0.7N(0,1)+0.3Cauchy(0,1) and t(1)  (1) When the error has a normal distribution, Sig has the best performance. However, when the error is contaminated or has a long tail, RankCorr can significantly outperform the alternatives (2) RankCorr has the largest AUC under all simulation scenarios, including those with linear transformation function and normal distribution error, followed by the nonrobust alternatives PeCOX and Sig  Li and Cui [28]  3G-SPA/Partial PCA (pPCA); Full PCA (fPCA)  n = 500, 1000, two separate genes with 10 SNPs each; discrete predictors: generated via MS program; Continuous response under linear model: (1) without genetic effects, (2) with main effects but no interaction, (3) with interaction effect smaller than main effects, (4) with interaction effect larger than main effects  3G-SPA has the highest power with the type I error reasonably controlled. The PCA-based methods are more conservative. The power of 3G-SPA is not depending on whether the genetic effect is because of its main effect or interaction, whereas the PCA-based methods prefer the scenarios with larger interactions  Kim and Park [29]  L-estimator GMDR; M-estimator GMDR/Ordinary GMDR  (1) n=1000,P=100, (2) n=3000,P=1000, d=(2,1); Discrete predictors: follow multinomial distribution with two different probabilities; Continuous response under linear model; Error distributions: (1−τ)N(0,1)+τN(0,400) and (1−τ)N(0,1)+τGamma(α,β) with τ=0,0.033 and 0.067  When there are no outliers (τ = 0), the ordinary GMDR and M-estimator GMDR outperform L- estimator GMDR. For all the designs with outliers, the performance of L-estimator and M-estimator GMDR are better than that of ordinary GMDR, with higher power. In some designs, L-estimator GMDR outperforms M-estimator GMDR  Yoshida and Koike [31]  SNPInterForest/SNPHarvester; BOOST  (1) n=4000,P=1000; Discrete predictors: the genetic model with different values of LD and MAFs; Binary response: the additive, multiplicative, and threshold models of interactions with weak marginal effects (2) n=400,P=1000; Discrete predictors: same as (1); Binary response: the pure interaction model without marginal effects  (1) The performance of SNPInterForest is superior to two alternatives in terms of their power SNPInterForest and SNPHarvester do not identify any spurious associations, whereas BOOST falsely identifies five interactions (2) The SNPInterForest and BOOST outperform SNPHarvester, with higher power and better identification accuracy The false discovery rates of SNPInterForest and SNPHarvester are much lower than that of BOOST  Beam et al. [32]  (1) BNN/BEAM; χ2 test (2)BNN/BEAM; χ2 test; GBM; MDR  n=2000,P=1000; Discrete predictors: the genetic model with different values of MAF; (1) Binary response: three biallelic models with main effects, (2) Binary response: pure interaction models without main effects  (1) BNN is uniformly more powerful than the χ2 test under all the simulated genetic models, and it outperforms BEAM in most instances (2) BNN performs better than BEAM, χ2 test and GBM under most of the simulations, and has similar performance with MDR. But MDR is much more time-consuming  Li et al. [33]  Nonpara-Scr/Screen and Clean (SC); SIS+LASSO; Forward LASSO; MDR  n=1000, P=3512, d=(12,8); Discrete predictors: dichotomizing the continuous variables generated from normal distribution with different correlations; Continuous response: a genetic model with both additive effects and dominant effects, considering three different noise levels  Nonpara-Scrbehaves better than the other four alternatives, with higher power and much smaller false-positive rates. Their average power are 97% (Nonpara-Scr), 90% (SIS+LASSO), 80% (SC), 48% (Forward LASSO) and 64% (MDR)  Li and Liu [34]  SODA/ Backward procedure (Back); Forward–backward method (ForBack); hierNet; IIS-SQDA  Ten different sample sizes with n=100, 100×101/9,…,1000; (1) P=50,d=(1,4); Continuous predictors: generated from normal, non-normal and heteroskedastic distributions, (2) P=1000,d=(1,4); Continuous predictors: non-normal distribution (3) P=50,d=(0,4), Continuous predictors: same as (2), (4) P=50,d=(2,4), Continuous predictors: same as (2); Binary response: Bayes classification rule  With Gaussian predictors, SODA, Back and ForBack can detect all important predictors as n increases, with small false positives, whereas IIS-SQDA and hierNet select too many false positives, resulting in high test error rates. Under other predictor settings, SODA performs much better than the nonrobust alternatives and achieves near-oracle classification accuracy when sample size is large enough  Wu et al. [35]  RankReg/LSregression-based method; Robust marginal analysis (RobMA)  n=200, P=200, q=2, d=(0,0,20); (1) Discrete predictors: MS program and Hardy–Weinberg equilibrium-based method, (2) Continuous predictors: follow multivariate normal distribution with an auto-regressive correlation structure; Continuous response under the linear model; Error distributions: N(0,1),t(3),0.85N(0,1)+0.15 log Normal(0,1) and 0.85N(0,1)+0.15Cauchy(0,1)  With N(0, 1) error distribution, RankReg is slightly inferior to the unrobust alternative LS, but still outperforms RobMA. With other three error types, RankReg has significant advantages, with more true positives, much fewer false positives, and higher estimation accuracy  Wu et al. [36]  LAD-hier/LAD-nhier; LAD-group; LS-hier; LS-nhier; LS-group  n=200, P=200, q=10, d=(10, 10, 15); (1) Continuous predictors: follow multivariate normal distribution with the banded and autoregressive correlation structure, (2) Discrete predictors: dichotomizing the above continuous values at the first and third quartiles; Censored response under the AFT model; Error distributions: N(0,1),0.80N(0,1)+0.20 log Normal(0,1),t(2) and 0.80N(0,1)+0.20Cauchy(0,1)  Without data contamination, the robust LAD-hier and the nonrobust counterpart LS-hier have comparable performances. With heavy-tailed errors, LAD-hier is significantly superior to the alternatives, in terms of both identification and prediction performance  a d=(b,c): b and c represent the number of true non-zero G effects and G-G interactions. b d=(a,b,c): a, b and c represent the number of true non-zero E effects, G effects and G-E interactions. Table 2. Summary of comparisons of the robust genetic analysis approaches against alternatives (partial list) References  Proposed/alternative  Main settings  Major conclusions  Wang et al. [26]  W-test/logistic regression; chi-squared test; MDR  n=1000,P=50; Discrete predictors: six different genetic architectures with minor allele frequency (MAF) in the common range and in the low frequency range, and LD in the high range, medium range and low range; binary response: (1) d=(4,2)a, linear interaction model with main effect, (2) d=(0,2), nonlinear interaction model without main effect  (1) The model-free W-test and MDR outperform the logistic regression even under the linear model, with higher power and lower type I error rate  (2) Their advantages become more prominent when the underlying model is nonlinear. For one genetic architecture, the power is 5.9% (logistic), 72.6% (Chi-squared), 95.5% (MDR) and 88% (W-test). W-test behaves much better than MDR when 1%<MAF<5%, where the average power of W-test and MDR are 87.6 and 73.2%, respectively  Shi et al. [27]  RankCorr/P-value-based significance analysis (Sig); (1) Penalized least squares (PeLS), (2) Penalized Cox model (PeCOX)  n=200,P=500 1000,q=5, d=(3,4,20)b; Continuous predictors: follow multivariate normal distribution with two correlation structures and two correlations; (1) Continuous response under linear model, (2) Censored response under transformation model: g(t)=t3,t5,t and t1/3; Error distributions: N(0, 1), 0.7N(0,1)+0.3Cauchy(0,1) and t(1)  (1) When the error has a normal distribution, Sig has the best performance. However, when the error is contaminated or has a long tail, RankCorr can significantly outperform the alternatives (2) RankCorr has the largest AUC under all simulation scenarios, including those with linear transformation function and normal distribution error, followed by the nonrobust alternatives PeCOX and Sig  Li and Cui [28]  3G-SPA/Partial PCA (pPCA); Full PCA (fPCA)  n = 500, 1000, two separate genes with 10 SNPs each; discrete predictors: generated via MS program; Continuous response under linear model: (1) without genetic effects, (2) with main effects but no interaction, (3) with interaction effect smaller than main effects, (4) with interaction effect larger than main effects  3G-SPA has the highest power with the type I error reasonably controlled. The PCA-based methods are more conservative. The power of 3G-SPA is not depending on whether the genetic effect is because of its main effect or interaction, whereas the PCA-based methods prefer the scenarios with larger interactions  Kim and Park [29]  L-estimator GMDR; M-estimator GMDR/Ordinary GMDR  (1) n=1000,P=100, (2) n=3000,P=1000, d=(2,1); Discrete predictors: follow multinomial distribution with two different probabilities; Continuous response under linear model; Error distributions: (1−τ)N(0,1)+τN(0,400) and (1−τ)N(0,1)+τGamma(α,β) with τ=0,0.033 and 0.067  When there are no outliers (τ = 0), the ordinary GMDR and M-estimator GMDR outperform L- estimator GMDR. For all the designs with outliers, the performance of L-estimator and M-estimator GMDR are better than that of ordinary GMDR, with higher power. In some designs, L-estimator GMDR outperforms M-estimator GMDR  Yoshida and Koike [31]  SNPInterForest/SNPHarvester; BOOST  (1) n=4000,P=1000; Discrete predictors: the genetic model with different values of LD and MAFs; Binary response: the additive, multiplicative, and threshold models of interactions with weak marginal effects (2) n=400,P=1000; Discrete predictors: same as (1); Binary response: the pure interaction model without marginal effects  (1) The performance of SNPInterForest is superior to two alternatives in terms of their power SNPInterForest and SNPHarvester do not identify any spurious associations, whereas BOOST falsely identifies five interactions (2) The SNPInterForest and BOOST outperform SNPHarvester, with higher power and better identification accuracy The false discovery rates of SNPInterForest and SNPHarvester are much lower than that of BOOST  Beam et al. [32]  (1) BNN/BEAM; χ2 test (2)BNN/BEAM; χ2 test; GBM; MDR  n=2000,P=1000; Discrete predictors: the genetic model with different values of MAF; (1) Binary response: three biallelic models with main effects, (2) Binary response: pure interaction models without main effects  (1) BNN is uniformly more powerful than the χ2 test under all the simulated genetic models, and it outperforms BEAM in most instances (2) BNN performs better than BEAM, χ2 test and GBM under most of the simulations, and has similar performance with MDR. But MDR is much more time-consuming  Li et al. [33]  Nonpara-Scr/Screen and Clean (SC); SIS+LASSO; Forward LASSO; MDR  n=1000, P=3512, d=(12,8); Discrete predictors: dichotomizing the continuous variables generated from normal distribution with different correlations; Continuous response: a genetic model with both additive effects and dominant effects, considering three different noise levels  Nonpara-Scrbehaves better than the other four alternatives, with higher power and much smaller false-positive rates. Their average power are 97% (Nonpara-Scr), 90% (SIS+LASSO), 80% (SC), 48% (Forward LASSO) and 64% (MDR)  Li and Liu [34]  SODA/ Backward procedure (Back); Forward–backward method (ForBack); hierNet; IIS-SQDA  Ten different sample sizes with n=100, 100×101/9,…,1000; (1) P=50,d=(1,4); Continuous predictors: generated from normal, non-normal and heteroskedastic distributions, (2) P=1000,d=(1,4); Continuous predictors: non-normal distribution (3) P=50,d=(0,4), Continuous predictors: same as (2), (4) P=50,d=(2,4), Continuous predictors: same as (2); Binary response: Bayes classification rule  With Gaussian predictors, SODA, Back and ForBack can detect all important predictors as n increases, with small false positives, whereas IIS-SQDA and hierNet select too many false positives, resulting in high test error rates. Under other predictor settings, SODA performs much better than the nonrobust alternatives and achieves near-oracle classification accuracy when sample size is large enough  Wu et al. [35]  RankReg/LSregression-based method; Robust marginal analysis (RobMA)  n=200, P=200, q=2, d=(0,0,20); (1) Discrete predictors: MS program and Hardy–Weinberg equilibrium-based method, (2) Continuous predictors: follow multivariate normal distribution with an auto-regressive correlation structure; Continuous response under the linear model; Error distributions: N(0,1),t(3),0.85N(0,1)+0.15 log Normal(0,1) and 0.85N(0,1)+0.15Cauchy(0,1)  With N(0, 1) error distribution, RankReg is slightly inferior to the unrobust alternative LS, but still outperforms RobMA. With other three error types, RankReg has significant advantages, with more true positives, much fewer false positives, and higher estimation accuracy  Wu et al. [36]  LAD-hier/LAD-nhier; LAD-group; LS-hier; LS-nhier; LS-group  n=200, P=200, q=10, d=(10, 10, 15); (1) Continuous predictors: follow multivariate normal distribution with the banded and autoregressive correlation structure, (2) Discrete predictors: dichotomizing the above continuous values at the first and third quartiles; Censored response under the AFT model; Error distributions: N(0,1),0.80N(0,1)+0.20 log Normal(0,1),t(2) and 0.80N(0,1)+0.20Cauchy(0,1)  Without data contamination, the robust LAD-hier and the nonrobust counterpart LS-hier have comparable performances. With heavy-tailed errors, LAD-hier is significantly superior to the alternatives, in terms of both identification and prediction performance  References  Proposed/alternative  Main settings  Major conclusions  Wang et al. [26]  W-test/logistic regression; chi-squared test; MDR  n=1000,P=50; Discrete predictors: six different genetic architectures with minor allele frequency (MAF) in the common range and in the low frequency range, and LD in the high range, medium range and low range; binary response: (1) d=(4,2)a, linear interaction model with main effect, (2) d=(0,2), nonlinear interaction model without main effect  (1) The model-free W-test and MDR outperform the logistic regression even under the linear model, with higher power and lower type I error rate  (2) Their advantages become more prominent when the underlying model is nonlinear. For one genetic architecture, the power is 5.9% (logistic), 72.6% (Chi-squared), 95.5% (MDR) and 88% (W-test). W-test behaves much better than MDR when 1%<MAF<5%, where the average power of W-test and MDR are 87.6 and 73.2%, respectively  Shi et al. [27]  RankCorr/P-value-based significance analysis (Sig); (1) Penalized least squares (PeLS), (2) Penalized Cox model (PeCOX)  n=200,P=500 1000,q=5, d=(3,4,20)b; Continuous predictors: follow multivariate normal distribution with two correlation structures and two correlations; (1) Continuous response under linear model, (2) Censored response under transformation model: g(t)=t3,t5,t and t1/3; Error distributions: N(0, 1), 0.7N(0,1)+0.3Cauchy(0,1) and t(1)  (1) When the error has a normal distribution, Sig has the best performance. However, when the error is contaminated or has a long tail, RankCorr can significantly outperform the alternatives (2) RankCorr has the largest AUC under all simulation scenarios, including those with linear transformation function and normal distribution error, followed by the nonrobust alternatives PeCOX and Sig  Li and Cui [28]  3G-SPA/Partial PCA (pPCA); Full PCA (fPCA)  n = 500, 1000, two separate genes with 10 SNPs each; discrete predictors: generated via MS program; Continuous response under linear model: (1) without genetic effects, (2) with main effects but no interaction, (3) with interaction effect smaller than main effects, (4) with interaction effect larger than main effects  3G-SPA has the highest power with the type I error reasonably controlled. The PCA-based methods are more conservative. The power of 3G-SPA is not depending on whether the genetic effect is because of its main effect or interaction, whereas the PCA-based methods prefer the scenarios with larger interactions  Kim and Park [29]  L-estimator GMDR; M-estimator GMDR/Ordinary GMDR  (1) n=1000,P=100, (2) n=3000,P=1000, d=(2,1); Discrete predictors: follow multinomial distribution with two different probabilities; Continuous response under linear model; Error distributions: (1−τ)N(0,1)+τN(0,400) and (1−τ)N(0,1)+τGamma(α,β) with τ=0,0.033 and 0.067  When there are no outliers (τ = 0), the ordinary GMDR and M-estimator GMDR outperform L- estimator GMDR. For all the designs with outliers, the performance of L-estimator and M-estimator GMDR are better than that of ordinary GMDR, with higher power. In some designs, L-estimator GMDR outperforms M-estimator GMDR  Yoshida and Koike [31]  SNPInterForest/SNPHarvester; BOOST  (1) n=4000,P=1000; Discrete predictors: the genetic model with different values of LD and MAFs; Binary response: the additive, multiplicative, and threshold models of interactions with weak marginal effects (2) n=400,P=1000; Discrete predictors: same as (1); Binary response: the pure interaction model without marginal effects  (1) The performance of SNPInterForest is superior to two alternatives in terms of their power SNPInterForest and SNPHarvester do not identify any spurious associations, whereas BOOST falsely identifies five interactions (2) The SNPInterForest and BOOST outperform SNPHarvester, with higher power and better identification accuracy The false discovery rates of SNPInterForest and SNPHarvester are much lower than that of BOOST  Beam et al. [32]  (1) BNN/BEAM; χ2 test (2)BNN/BEAM; χ2 test; GBM; MDR  n=2000,P=1000; Discrete predictors: the genetic model with different values of MAF; (1) Binary response: three biallelic models with main effects, (2) Binary response: pure interaction models without main effects  (1) BNN is uniformly more powerful than the χ2 test under all the simulated genetic models, and it outperforms BEAM in most instances (2) BNN performs better than BEAM, χ2 test and GBM under most of the simulations, and has similar performance with MDR. But MDR is much more time-consuming  Li et al. [33]  Nonpara-Scr/Screen and Clean (SC); SIS+LASSO; Forward LASSO; MDR  n=1000, P=3512, d=(12,8); Discrete predictors: dichotomizing the continuous variables generated from normal distribution with different correlations; Continuous response: a genetic model with both additive effects and dominant effects, considering three different noise levels  Nonpara-Scrbehaves better than the other four alternatives, with higher power and much smaller false-positive rates. Their average power are 97% (Nonpara-Scr), 90% (SIS+LASSO), 80% (SC), 48% (Forward LASSO) and 64% (MDR)  Li and Liu [34]  SODA/ Backward procedure (Back); Forward–backward method (ForBack); hierNet; IIS-SQDA  Ten different sample sizes with n=100, 100×101/9,…,1000; (1) P=50,d=(1,4); Continuous predictors: generated from normal, non-normal and heteroskedastic distributions, (2) P=1000,d=(1,4); Continuous predictors: non-normal distribution (3) P=50,d=(0,4), Continuous predictors: same as (2), (4) P=50,d=(2,4), Continuous predictors: same as (2); Binary response: Bayes classification rule  With Gaussian predictors, SODA, Back and ForBack can detect all important predictors as n increases, with small false positives, whereas IIS-SQDA and hierNet select too many false positives, resulting in high test error rates. Under other predictor settings, SODA performs much better than the nonrobust alternatives and achieves near-oracle classification accuracy when sample size is large enough  Wu et al. [35]  RankReg/LSregression-based method; Robust marginal analysis (RobMA)  n=200, P=200, q=2, d=(0,0,20); (1) Discrete predictors: MS program and Hardy–Weinberg equilibrium-based method, (2) Continuous predictors: follow multivariate normal distribution with an auto-regressive correlation structure; Continuous response under the linear model; Error distributions: N(0,1),t(3),0.85N(0,1)+0.15 log Normal(0,1) and 0.85N(0,1)+0.15Cauchy(0,1)  With N(0, 1) error distribution, RankReg is slightly inferior to the unrobust alternative LS, but still outperforms RobMA. With other three error types, RankReg has significant advantages, with more true positives, much fewer false positives, and higher estimation accuracy  Wu et al. [36]  LAD-hier/LAD-nhier; LAD-group; LS-hier; LS-nhier; LS-group  n=200, P=200, q=10, d=(10, 10, 15); (1) Continuous predictors: follow multivariate normal distribution with the banded and autoregressive correlation structure, (2) Discrete predictors: dichotomizing the above continuous values at the first and third quartiles; Censored response under the AFT model; Error distributions: N(0,1),0.80N(0,1)+0.20 log Normal(0,1),t(2) and 0.80N(0,1)+0.20Cauchy(0,1)  Without data contamination, the robust LAD-hier and the nonrobust counterpart LS-hier have comparable performances. With heavy-tailed errors, LAD-hier is significantly superior to the alternatives, in terms of both identification and prediction performance  a d=(b,c): b and c represent the number of true non-zero G effects and G-G interactions. b d=(a,b,c): a, b and c represent the number of true non-zero E effects, G effects and G-E interactions. Remarks Our literature review suggests that there have been significant developments in robust genetic interaction analysis methods, especially in the recent literature. However, compared with low-dimensional studies and high-dimensional studies on main effects, research tailored to genetic interactions is still limited. Based on strategies developed in published studies, much robust genetic interaction research can be conducted. For example, the ‘robust loss function + penalization’ joint analysis can be extended to other robust loss functions, such as the Huber’s loss [66] and Tukey’s bi-square loss [67]. Some other robust estimations such as S-estimation and MM-estimation, which have been shown to be powerful in low-dimensional settings, can be potentially adapted to detect genetic interactions. The high-dimensional robust screening techniques may be potentially coupled with the nonparametric feature screening [68], rank correlation screening [69] and others. Most of the existing methods, including those reviewed above, are limited to one aspect of robustness. They often make no or weaker assumptions on one of the following aspects: specific model form, conditional distribution of the response (distribution of the error) or distribution of the predictors, resulting in different types of robustness. In practical data analysis, for example, model mis-specification and outliers may happen together and demand robustness in multiple aspects. A few studies that can accommodate both model mis-specification and outliers in response include the rank correlation-based method [27] and R-GMDR [29]. For high-dimensional main effect analysis, some methods have been proposed for variable selection and can accommodate outliers/contaminations in both responses and predictors, such as the sparse least trimmed squares regression [48], the robust L2 boosting method [70] and others. Similar development is also needed for interaction analysis. Our review also suggests that for some aspects of robustness, there are multiple methods available. However, there is still a lack of direct and comprehensive comparisons. In addition, most of the existing methods focus on methodological development. Theoretical studies are limited. It has been noted that theoretical investigation with high-dimensional data is challenging. Significantly more challenges can be brought by interactions. The computational algorithms are critically dependent on the methodology frameworks. However, some similar techniques are involved. For example, to address the problem that the robust objective functions are often not differentiable, the smoothed approximations are adopted. In addition, the (group) coordinate descent algorithms are usually developed for penalization approaches. The computational cost differs widely across these approaches because of their methodologies, adopted algorithms, sizes of the analyzed data sets and others. The robust methods are observed to be in general more computationally intensive than nonrobust methods, as they often have more complex objective functions and optimization algorithms. However, the results suggest that these robust interaction learning methods are still computational affordable. For some of the reviewed methods, software is publicly available. However, the development of user-friendly software is still badly needed. Applications The methods reviewed above and other robust methods in the literature have been applied to the analysis of multiple diseases and made sensible discoveries different from those detected by nonrobust methods. A few representative findings are briefly reviewed below. A sporadic breast cancer case-control data set is analyzed using MDR [24]. This case study aims at detecting SNP-SNP interactions associated with sporadic breast cancer using SNPs from selected genes. The data set contains data on 200 White women with sporadic primary invasive breast cancer treated at the Vanderbilt University Medical Center during 1982–96. Ten SNPs in five selected genes (COMT, CYP1A1, CYP1B1, GSTM1 and GSTT1) are analyzed. A four-locus interaction among four SNPs (COMT, CYP1A1m1, CYP1B1 codon 48 and CYP1B1 codon 432) is detected to be significantly associated with the risk of sporadic breast cancer. The protein products of the corresponding genes of these four SNPs have been shown to interact as enzymes in the metabolism of estrogens in breast tissues. Estrogens have been recognized as one of the principal factors related to breast cancer. Thus, it is expected that the interactions between the SNPs in these four genes have independent contributions to the risk of breast cancer. However, using the parametric logistic regression analysis, previous studies fail to find any association between these genes or interactions and an increased risk of breast cancer. Similar promising findings have also been made in other MDR studies. The nonparametric W-test proposed by [26] is applied to two bipolar disorder genome-wide association data sets. The first data set is from the Wellcome Trust Case Control Consortium (WTCCC), including 2000 bipolar cases and 3000 controls of European Caucasians with 414 682 SNPs. The second data set is collected by the Genetic Association Information Network (GAIN) project and consists of 653 cases and 1767 controls with 729 304 SNPs. For these two real GWAS data sets, a three-step procedure is used to search for SNP-SNP interactions and their corresponding G-G interactions. First, the main effects of the whole SNP set are evaluated using W-test and then a candidate SNP set is identified based on the P-value. Second, the W-test for detecting important SNP-SNP interactions is conducted on the candidate SNP set. Third, the G-G interactions corresponding to the identified SNP-SNP interactions are obtained using Web-based databases. It is interesting that the W-test approach identifies a number of genes that can be replicated in the two independent GWAS data sets, and that the interactions among 18 identified genes naturally form two networks. The identification results are different from the previous GWAS studies. Most of the findings play important roles in brain and neuronal functions and are highly relevant to neuron disorders, including RTN4R, DPP10, CMSD1, SLIT3, PARK2, TMEM132D, HNT, CNTNAP2, ACCN1, A2BP1 and others. For instance, RTN4R has been shown to encode a Nogo receptor that mediates axonal growth inhibition and plays a role in regulating axonal regeneration and plasticity in the central nervous system. Kim and Park [29] apply the R-GMDR approach to the Korea Association Resource (KARE) data, which is collected by the large-scale genome-wide association analyses in the Gyeonggi Province of South Korea. The outcome of interest is the (log-transformed) Homeostasis Model Assessment of Insulin Resistance (HOMA-IR) level, which is continuous and widely used to estimate insulin resistance. After removing samples with the missing BMI or HOMA-IR score, this study includes of 8577 individuals and 327 872 SNPs. To reduce the computational burden on the R-GMDR, a pre-screening procedure based on the P-value of a marginal SNP analysis and the linkage disequilibrium (LD) coefficient is conducted, resulting in 585 SNPs for the downstream analysis. Then, the R-GMDR is applied to identify important SNP-SNP interactions between these 585 SNPs. The identified SNPs are matched to their corresponding genes to facilitate the analysis of possible biological implications. The two R-GMDR methods identify sensible genes, especially hub genes, that are missed by the original GMDR. Specifically, the KCNH1 gene is identified to have strong interaction effects with other genes on the function of insulin secretion. The other hub gene RYR2 has been demonstrated to be related to diabetes. Outliers are observed in the tail part of the response, leading to an adverse influence on the original GMDR. Wu et al. [35] conduct a G-E interaction analysis for lung cancer prognosis. Data from four independent studies, conducted by the Dana-Farber Cancer Institute, HLM (Moffitt Cancer Center), UM (University of Michigan Cancer Center) and MSKCC, are pooled and analyzed. The combined sample size is 351. Different from the former three applications, which focus on SNP-SNP interactions, the continuous gene expression measurements on 22 283 probes are analyzed as G factors. Two clinical factors, age and smoking status, are analyzed as E factors. The response variable of interest is right censored overall survival. The robust method identifies 11 genes as having significant associations with prognosis, among which four and five genes are detected to having interactions with age and smoking status, respectively. The nonrobust approach based on the ordinary least squared loss function identifies 19 genes with main effects (6 overlap with the robust method) and three interactions with smoking status. A brief literature search is conducted and suggests that most of the genes identified by the robust approach have important implications for lung cancer prognosis. Findings missed by the nonrobust approach include genes CYP2B6, WIF1, HLA-DQA1 and others. Take gene WIF1 as an example. This gene has been implicated in lung cancer pathogenesis, and the methylation silencing of WIF1 is a usual and potentially central mechanism of the aberrant activation of the Wnt signaling pathway. Remarks In the aforementioned studies, robust methods have led to promising new discoveries. However, as robust methods are relatively new and limited, their practical impact may still be limited. In particular, their applications in large-scale cutting-edge studies and biomedical publications are still rare. Significant effort in applications is needed. Discussion Extensive research on genetic interactions has been conducted, and tailored analysis methods have been developed in a large number of, especially recent, studies. In this review, we have focused on robust methods, which have important practical implications but have not received enough attention in the literature. To the best of our knowledge, this article is the first to comprehensively review robust genetic interaction analysis methods. With our limited knowledge and literature search, it is inevitable the selection of reviewed methods is biased and incomplete. With the limited scope of this review, it is not our goal to provide the ‘full picture’. The review has been conducted in a rather ‘statistical’ manner because of the high complexity of the methods. In addition, we have tried to review the generically applicable (as opposed to case specific) methods. Even with these limitations, being the first of its kind, this review can still be much useful for researchers in the field of genetic interaction analysis. Funding This work was partly supported by the National Institutes of Health (grant numbers CA216017, CA191383, CA204120); and the National Natural Science Foundation of China (grant numbers 61402276 and 91546202). Mengyun Wu is an Assistant Professor in School of Statistics and Management, Shanghai University of Finance and Economics. She is also a Postdoctoral Fellow in the Department of Biostatistics at Yale University. Shuangge Ma is a Professor in the Department of Biostatistics at Yale University. References 1 Yoshimaru T, Komatsu M, Matsuo T, et al.   Targeting BIG3-PHB2 interaction to overcome tamoxifen resistance in breast cancer cells. Nat Commun  2013; 4: 2443. Google Scholar CrossRef Search ADS PubMed  2 Zhou W, Liu G, Miller DP, et al.   Gene-environment interaction for the ERCC2 polymorphisms and cumulative cigarette smoking exposure in lung cancer. Cancer Res  2002; 62: 1377– 81. Google Scholar PubMed  3 Cordell HJ. Detecting gene-gene interactions that underlie human diseases. Nat Rev Genet  2009; 10( 6): 392– 404. Google Scholar CrossRef Search ADS PubMed  4 Van Steen K. Travelling the world of gene-gene interactions. Brief Bioinform  2012; 13( 1): 1– 19. Google Scholar CrossRef Search ADS PubMed  5 Thomas D. Gene-environment-wide association studies: emerging approaches. Nat Rev Genet  2010; 11( 4): 259– 72. Google Scholar CrossRef Search ADS PubMed  6 Simonds NI, Ghazarian AA, Pimentel CB, et al.   Review of the gene-environment interaction literature in cancer: what do we know? Genet Epidemiol  2016; 40( 5): 356– 65. Google Scholar CrossRef Search ADS PubMed  7 Purcell S, Neale B, Todd-Brown K, et al.   PLINK: a tool set for whole-genome association and population-based linkage analyses. Am J Hum Genet  2007; 81( 3): 559– 75. Google Scholar CrossRef Search ADS PubMed  8 Zhao J, Zhu Y, Xiong M. Genome-wide gene-gene interaction analysis for next-generation sequencing. Eur J Hum Genet  2016; 24( 3): 421– 8. Google Scholar CrossRef Search ADS PubMed  9 Bild AH, Yao G, Chang JT, et al.   Oncogenic pathway signatures in human cancers as a guide to targeted therapies. Nature  2006; 439( 7074): 353– 7. Google Scholar CrossRef Search ADS PubMed  10 Weissbrod O, Lippert C, Geiger D, et al.   Accurate liability estimation improves power in ascertained case-control studies. Nat Methods  2015; 12( 4): 332. Google Scholar CrossRef Search ADS PubMed  11 Stark A, Stahl MS, Kirchner HL, et al.   Body mass index at the time of diagnosis and the risk of advanced stages and poorly differentiated cancers of the breast: findings from a case-series study. Int J Obes  2010; 34( 9): 1381– 6. Google Scholar CrossRef Search ADS   12 Stephan J, Stegle O, Beyer A. A random forest approach to capture genetic effects in the presence of population structure. Nat Commun  2015; 6( 1): 7432. Google Scholar CrossRef Search ADS PubMed  13 Osborne JW, Overbay A. The power of outliers (and why researchers should always check for them). Pract Assess Res Eval  2010; 9: 1– 12. 14 Li C, Wong WH. Model-based analysis of oligonucleotide arrays: expression index computation and outlier detection. Proc Natl Acad Sci USA  2001; 98( 1): 31– 6. Google Scholar CrossRef Search ADS PubMed  15 Shieh AD, Hung YS. Detecting outlier samples in microarray data. Stat Appl Genet Mol  2009; 8( 1): 1– 24. Google Scholar CrossRef Search ADS   16 Rampatige R, Gamage S, Peiris S, et al.   Assessing the reliability of causes of death reported by the vital registration system in Sri Lanka: medical records review in Colombo. Health Inf Manag  2013; 42( 3): 20– 8. Google Scholar PubMed  17 Cancer Genome Atlas Research Network. Comprehensive molecular profiling of lung adenocarcinoma. Nature  2014; 511: 543– 50. CrossRef Search ADS PubMed  18 Wu M, Zang Y, Zhang S, et al.   Accommodating missingness in environmental measurements in gene-environment interaction analysis. Genet Epidemiol  2017; 41( 6): 523– 54. Google Scholar CrossRef Search ADS PubMed  19 Chai H, Zhang Q, Jiang Y, et al.   Identifying gene-environment interactions for prognosis using a robust approach. Econ Stat  2017; 4: 105– 20. 20 Burgess DJ. Cancer genetics: initially complex, always heterogeneous. Nat Rev Cancer  2011; 11( 3): 153. Google Scholar CrossRef Search ADS PubMed  21 Haibe-Kains B, Desmedt C, Loi S, et al.   A three-gene model to robustly identify breast cancer molecular subtypes. J Natl Cancer Inst  2012; 104( 4): 311– 25. Google Scholar CrossRef Search ADS PubMed  22 Huber PJ, Ronchetti EM. Robust Statistics. Wiley Series in Probability and Statistics , 2nd ed. Hoboken, NJ: Wiley, 2009. 23 Wu C, Ma S. A selective review of robust variable selection with applications in bioinformatics. Brief Bioinform  2015; 16( 5): 873– 83. Google Scholar CrossRef Search ADS PubMed  24 Ritchie MD, Hahn LW, Roodi N, et al.   Multifactor-dimensionality reduction reveals high-order interactions among estrogen-metabolism genes in sporadic breast cancer. Am J Hum Genet  2001; 69( 1): 138– 47. Google Scholar CrossRef Search ADS PubMed  25 Winham SJ, Motsinger-Reif AA. An R package implementation of multifactor dimensionality reduction. Biodata Min  2011; 4( 1): 24. Google Scholar CrossRef Search ADS PubMed  26 Wang MH, Sun R, Guo J, et al.   A fast and powerful W-test for pairwise epistasis testing. Nucleic Acids Res  2016; 44( 12): e115. Google Scholar CrossRef Search ADS PubMed  27 Shi X, Liu J, Huang J, et al.   A penalized robust method for identifying gene-environment interactions. Genet Epidemiol  2014; 38( 3): 220– 30. Google Scholar CrossRef Search ADS PubMed  28 Li S, Cui Y. Gene-centric gene-gene interaction: a model-based kernel machine method. Ann Appl Stat  2012; 6( 3): 1134– 61. Google Scholar CrossRef Search ADS   29 Kim Y, Park T. Robust gene-gene interaction analysis in genome wide association studies. PLoS One  2015; 10( 8): e0135016. Google Scholar CrossRef Search ADS PubMed  30 Wang G, Zhao Y, Zhang Q, et al.   Identifying gene-environment interactions associated with prognosis using penalized quantile regression. In: Big and Complex Data Analysis . Springer International Publishing, 2017, 347– 67. Google Scholar CrossRef Search ADS   31 Yoshida M, Koike A. SNPInterForest: a new method for detecting epistatic interactions. BMC Bioinform  2011; 12: 469. Google Scholar CrossRef Search ADS   32 Beam AL, Motsinger-Reif A, Doyle J. Bayesian neural networks for detecting epistasis in genetic association studies. BMC Bioinform  2014; 15( 1): 368. Google Scholar CrossRef Search ADS   33 Li J, Dan J, Li C, et al.   A model-free approach for detecting interactions in genetic association studies. Brief Bioinform  2014; 15( 6): 1057– 8. Google Scholar CrossRef Search ADS PubMed  34 Li Y, Liu JS. Robust variable and interaction selection for logistic regression and general index models. J Am Stat Assoc  2017, in press. 35 Wu C, Shi X, Cui Y, et al.   A penalized robust semiparametric approach for gene-environment interactions. Stat Med  2015; 34( 30): 4016– 30. Google Scholar CrossRef Search ADS PubMed  36 Wu C, Jiang Y, Ren J, et al.   Dissecting gene-environment interactions: a penalized robust approach accounting for hierarchical structures. Stat Med  2018; 37( 3): 437– 56. Google Scholar CrossRef Search ADS PubMed  37 Heagerty PJ, Kurland BF. Misspecified maximum likelihood estimates and generalised linear mixed models. Biometrika  2001; 88( 4): 973– 85. Google Scholar CrossRef Search ADS   38 Lou XY, Chen GB, Yan L, et al.   A generalized combinatorial approach for detecting gene-by-gene and gene-by-environment interactions with application to nicotine dependence. Am J Hum Genet  2007; 80( 6): 1125– 37. Google Scholar CrossRef Search ADS PubMed  39 Li MD, Burmeister M. New insights into the genetics of addiction. Nat Rev Cancer  2009; 10( 4): 225. Google Scholar CrossRef Search ADS   40 Gui J, Moore JH, Kelsey KT, et al.   A novel survival multifactor dimensionality reduction method for detecting gene-gene interactions with application to bladder cancer prognosis. Hum Genet  2011; 129( 1): 101– 10. Google Scholar CrossRef Search ADS PubMed  41 Gao X, Alvo M. A unified nonparametric approach for unbalanced factorial designs. J Am Stat Assoc  2005; 100( 471): 926– 41. Google Scholar CrossRef Search ADS   42 Wu X, Jin L, Xiong M. Mutual information for testing gene-environment interaction. PLoS One  2009; 4( 2): e4578. Google Scholar CrossRef Search ADS PubMed  43 Li M, Ye C, Fu W, et al.   Detecting genetic interactions for quantitative traits with U-statistics. Genet Epidemiol  2011; 35( 6): 457– 68. Google Scholar PubMed  44 Larson NB, Schaid DJ. A kernel regression approach to gene-gene interaction detection for case-control studies. Genet Epidemiol  2013; 37( 7): 695– 703. Google Scholar CrossRef Search ADS PubMed  45 Marceau R, Lu W, Holloway S, et al.   A fast multiple-kernel method with applications to detect gene-environment interaction. Genet Epidemiol  2015; 39( 6): 456– 68. Google Scholar CrossRef Search ADS PubMed  46 Lin HY, Wang W, Liu YH, et al.   Comparison of multivariate adaptive regression splines and logistic regression in detecting SNP-SNP interactions and their application in prostate cancer. J Hum Genet  2008; 53( 9): 802– 11. Google Scholar CrossRef Search ADS PubMed  47 Luss R, Rosset S, Shahar M. Efficient regularized isotonic regression with application to gene-gene interaction search. Ann Appl Stat  2012; 6( 1): 253– 83. Google Scholar CrossRef Search ADS   48 Alfons A, Croux C, Gelper S. Sparse least trimmed squares regression for analyzing high-dimensional large data sets. Ann Appl Stat  2013; 7( 1): 226– 48. Google Scholar CrossRef Search ADS   49 Lee S, Kim Y, Kwon MS, et al.   A comparative study on multifactor dimensionality reduction methods for detecting gene-gene interactions with the survival phenotype. BioMed Res Int  2015; 2015: 671859. Google Scholar PubMed  50 Fan J, Xue L, Zou H. Multitask quantile regression under the transnormal model. J Am Stat Assoc  2016; 111( 516): 1726– 35. Google Scholar CrossRef Search ADS PubMed  51 Chaudhuri P, Loh WY. Nonparametric estimation of conditional quantiles using quantile regression trees. Bernoulli  2002; 8: 561– 76. 52 Zhu Q, Hu Y, Tian M. Identifying interaction effects via additive quantile regression models. Stat Its Interface  2017; 10( 2): 255– 65. Google Scholar CrossRef Search ADS   53 Bien J, Taylor J, Tibshirani R. A lasso for hierarchical interactions. Ann Stat  2013; 41( 3): 1111– 41. Google Scholar CrossRef Search ADS PubMed  54 Lim M, Hastie T. Learning interactions via hierarchical group-lasso regularization. J Comput Graph Stat  2015; 24( 3): 627– 54. Google Scholar CrossRef Search ADS PubMed  55 Zhu R, Zhao H, Ma S. Identifying gene-environment and gene-gene interactions using a progressive penalization approach. Genet Epidemiol  2014; 38( 4): 353– 68. Google Scholar CrossRef Search ADS PubMed  56 Wu M, Huang J, Ma S. Identifying gene-gene interactions using penalized tensor regression. Stat Med  2018; 37( 4): 598– 610. Google Scholar CrossRef Search ADS PubMed  57 Hao N, Zhang HH. A note on high-dimensional linear regression with interactions. Am Stat  2017; 71( 4): 291– 7. Google Scholar CrossRef Search ADS   58 Winham SJ, Colby CL, Freimuth RR, et al.   SNP interaction detection with random forests in high-dimensional genetic data. BMC Bioinform  2012; 13( 1): 164. Google Scholar CrossRef Search ADS   59 Li J, Malley JD, Andrew AS, et al.   Detecting gene-gene interactions using a permutation-based random forest method. BioData Min  2016; 9( 1): 14. Google Scholar CrossRef Search ADS PubMed  60 Ritchie MD, White BC, Parker JS, et al.   Optimizationof neural network architecture using genetic programming improvesdetection and modeling of gene-gene interactions in studies of human diseases. BMC Bioinform  2003; 4: 28. Google Scholar CrossRef Search ADS   61 Motsinger-Reif AA, Dudek SM, Hahn LW, et al.   Comparison of approaches for machine-learning optimization of neural networks for detecting gene-gene interactions in genetic epidemiology. Genet Epidemiol  2008; 32( 4): 325– 40. Google Scholar CrossRef Search ADS PubMed  62 Günther F, Wawro N, Bammann K. Neural networks for modeling gene-gene interactions in association studies. BMC Genet  2009; 10( 1): 87. Google Scholar CrossRef Search ADS PubMed  63 Ritchie MD, Motsinger AA, Bush WS, et al.   Genetic programming neural networks: a powerful bioinformatics tool for human genetics. Appl Soft Comput  2007; 7( 1): 471– 9. Google Scholar CrossRef Search ADS PubMed  64 Liu J, Huang J, Zhang Y, et al.   Identification of gene-environment interactions in cancer studies using penalization. Genomics  2013; 102( 4): 189– 94. Google Scholar CrossRef Search ADS PubMed  65 Maronna RA, Yohai VJ. High finite-sample efficiency and robustness based on distance-constrained maximum likelihood. Comput Stat Data Anal  2015; 83: 262– 74. Google Scholar CrossRef Search ADS   66 Huber PJ. Robust estimation of a location parameter. Ann Math Stat  1964; 35( 1): 73– 101. Google Scholar CrossRef Search ADS   67 Beaton A, Tukey J. The fitting of power series, meaning polynomials, illustrated on band-spectroscopic data. Technometrics  1974; 16( 2): 147– 85. Google Scholar CrossRef Search ADS   68 Lin L, Sun J, Zhu L. Nonparametric feature screening. Comput Stat Data Anal  2013; 67: 162– 74. Google Scholar CrossRef Search ADS   69 Li G, Peng H, Zhang J, et al.   Robust rank correlation based screening. Ann Stat  2012; 40( 3): 1846– 77. Google Scholar CrossRef Search ADS   70 Lutz RW, Kalisch M, Bühlmann P. Robustified L2 boosting. Comput Stat Data Anal  2008; 52( 7): 3331– 41. Google Scholar CrossRef Search ADS   © The Author(s) 2018. Published by Oxford University Press. All rights reserved. For permissions, please email: journals.permissions@oup.com This article is published and distributed under the terms of the Oxford University Press, Standard Journals Publication Model (https://academic.oup.com/journals/pages/about_us/legal/notices) http://www.deepdyve.com/assets/images/DeepDyve-Logo-lg.png Briefings in Bioinformatics Oxford University Press

Robust genetic interaction analysis

Loading next page...
 
/lp/ou_press/robust-genetic-interaction-analysis-98UB3NyMEm
Publisher
Oxford University Press
Copyright
© The Author(s) 2018. Published by Oxford University Press. All rights reserved. For permissions, please email: journals.permissions@oup.com
ISSN
1467-5463
eISSN
1477-4054
D.O.I.
10.1093/bib/bby033
Publisher site
See Article on Publisher Site

Abstract

Abstract For the risk, progression, and response to treatment of many complex diseases, it has been increasingly recognized that genetic interactions (including gene–gene and gene–environment interactions) play important roles beyond the main genetic and environmental effects. In practical genetic interaction analyses, model mis-specification and outliers/contaminations in response variables and covariates are not uncommon, and demand robust analysis methods. Compared with their nonrobust counterparts, robust genetic interaction analysis methods are significantly less popular but are gaining attention fast. In this article, we provide a comprehensive review of robust genetic interaction analysis methods, on their methodologies and applications, for both marginal and joint analysis, and for addressing model mis-specification as well as outliers/contaminations in response variables and covariates. genetic interaction, robustness, model mis-specification, outlier/contamination Introduction Genetic interactions, including both gene–gene (G-G) interactions and gene–environment (G-E) interactions, have been shown to play crucial roles for the etiology, prognosis and response to treatment of many complex diseases beyond the main effects. In published studies, G factors include gene expressions, methylation, SNPs, as well as several other types of omics measurements. Beyond ‘real’ environmental variables, published studies have also examined demographic, clinical/medical, socioeconomic and other E factors. Many significant findings on genetic interactions have been made. For example, targeting the interaction between genes BIG3 and PHB2 has been suggested as a potential new treatment venue for ERα-positive breast cancer patients [1]. Researchers have found that the associations between ERCC2 polymorphisms and lung cancer risk are modified by the cumulative cigarette smoking [2]. A large number of statistical and machine learning methods have been developed for identifying important genetic interactions. Comprehensive discussions are available in [3–6], and many others. Take G-E interaction as an example. Most of the existing methods, although differ in technical details, share the following common strategy [7, 8]. Consider the model response∼f(G+E+G×E). In this model, the form of f is usually prespecified, for example, as logistic, linear and Cox for binary, continuous and survival responses, respectively; and linear covariate effects are assumed. Standard likelihood-based estimation is usually conducted, and conclusions on the importance of interactions are based on the estimates and significance levels of the G-E interaction terms. This strategy, despite many advantages, is nonrobust. In genetic interaction analysis, the demand for robustness arises in multiple aspects. Model mis-specification. There are at least two aspects. (a) Take binary responses, for example, disease/healthy status and response/resistance to treatment, as an example. Many of the published studies adopt the logistic regression model without sufficient justification. In ‘classic’ biomedical studies, it has been suggested that the logistic model assumption may fail, and alternative models, for example the probit and log models, may be needed [9, 10]. Genetic data and interactions are much more complex, and it is unlikely that they can be described using a single ‘default’ model. (b) The assumption of linear covariate effects may not always hold. For example, it has been shown that some E factors, such as age and body mass index (BMI), often have nonlinear effects [11]. A few publications have conducted exploratory analysis and suggested that some G effects may also deviate from linear [12]. Outliers/contaminations in predictors. E factors have been extensively studied in classic biomedical studies. It has been suggested that outliers/data contaminations are not uncommon, which can be caused by errors in data collection/recording as well as existence of influential observations [13]. For some types of G factors, for example gene expressions, outliers/contaminations can be caused by technical problems in sample preparation and profiling, human errors, as well as genetic abnormalities [14]. Outliers/contaminations in responses. Most if not all response variables studied in genetic studies have been investigated in classic biomedical studies, where the problem of outliers/contaminations in response has been well recognized [15]. For example, in prognosis studies, it has been recognized that human errors can be responsible for extremely long/short survivals as well as misclassification in the cause of death [16]. Specifically, consider the lung adenocarcinoma (LUAD) data collected by The Cancer Genome Atlas (TCGA) [17]. TCGA has generated comprehensive clinical, environmental and genetic data with high quality for multiple cancer types. The LUAD data have been analyzed in multiple published studies, including some interaction analyses [18, 19]. For the 504 subjects with complete mRNA expressions and survival time, there are six subjects with survival times 156.54, 162.98 163.99, 221.16, 232.00 and 238.11 months, respectively, while the rest 498 subjects have survival times ranging from 0.13 to 129.43. Some outliers are observed in response.A unique characteristic of genetic studies is that usually they cannot afford conducting strict subject selections, and the studied cohorts are often composed of multiple disease subtypes [20]. Different subtypes often behave differently in multiple aspects [21] and cannot be described using the same model. In this case, it can be viewed that data on the major subtype are ‘contaminated’ by the small subtypes. The possibility of model mis-specification and existence of outliers/contaminations in predictors and responses call for robust analysis methods. In low-dimensional biomedical studies, it has been well established that model mis-specification and failing to account for outliers and data contamination lead to biased estimation, false marker identification and inference and suboptimal models [22]. Robust analysis methods have been extensively developed and shown to have important scientific impact. With the high dimensionality and other unique features of genetic data, robust methods developed for low-dimensional biomedical data cannot be directly applied, and new development is needed. Our literature search and published reviews [23] suggest that, for the analysis of high-dimensional genetic data, the development of robust analysis methods is still much limited; however, it is gaining popularity fast. The majority of robust genetic data analysis methods focus on the main effect models, and the development on robust genetic interactions has been scattered. Considering the increasing importance of genetic interactions and demand for robust analysis methods in practical data analyses, our goal is to conduct a structured and systematic review of robust genetic interaction methods developed in recent literature. To the best of our knowledge, such a review is not available in the literature. This study targets to fill this knowledge gap. Methods In the analysis of genetic interactions (main effects as well), there are two generic analysis paradigm [3]. In marginal analysis, one or a small number of genes (or other omics measurement units) are analyzed at a time. In joint analysis, a large number of genes are analyzed under a single model. Although linked, the two analysis paradigms differ significantly in multiple aspects. Thus, below we review for each paradigm separately. In addition, the three robustness aspects, namely, model mis-specification, outliers/contaminations in predictors and outliers/contaminations in responses, are discussed separately. A summary of these methods is presented in Table 1. In addition, the following two considerations are taken into account. Table 1. Summary of the robust genetic analysis approaches (partial list) Approach  M/Ja  Robustb  G-E/G-G  Software  Applicationc  Conclusion  MDR [24, 25]  M  MisSpec  SNP-SNP  https://cran.r-project.org/ web/packages/MDR/ index.html  A sporadic breast cancer case-control data set, n=200,P=10  For the real data set, without any statistically significant main effects, MDR identifies one statistically significant high-order interaction among four SNPs from three different estrogen-metabolism genes, which is missed by the previous studies using logistic regression  W-test [26]  M  MisSpec  SNP-SNP  http://www2.ccrb.cuhk. edu.hk/wtest/download. html.  Two bipolar disorder GWAS data sets: (1) n=5000,P=414 682, (2) n=2420,P=729 304  For the two real data sets, W-test identifies significant interaction pairs that can be replicated. The genes in pairs play central roles in neurotransmission and synapse formation, and are undiscoverable by main effects with low frequency variants  RankCorr [27]  M  MisSpec  G-E    A lung cancer prognosis data set, n=336,P=2500 (22 283),q=5.  RankCorr identifies genes (SFRP1, IGFBP1, AHSG, AGR2 and others) and their interactions with important implications, which are different from those using alternatives. Most of the identified interactions are between genes and smoking status  3G-SPA [28]  M  MisSpec  G-G with SNPs  http://cran.r-project. org/package=SPA3G  (1) A birth weight data set, n = 1511 with 797 SNPs in 186 genes. (2) A yeast eQTL mapping data set, n = 112 with 6229 SNPs in 2956 genes  For the first data set, 23 G-G interactions are detected with the cutoff P-value 0.001, of which the genes include PDGFC, ANG, PTGER3, IL9 and others. For the second data set, 87 pairs are found to have significant interactions at an experiment-wise significant level of 0.05, where the corresponding genes include URA3, BAT2, LUE2 and others. For both data sets, the findings are shown to be biological sensible in the literature  R-GMDR [29]  M  Out-Res  SNP-SNP    A KARE data set with a HOMA-IR phenotype, n=8577,P=585 (327 872)  For the real data set, R-GMDR identifies a larger number of hub genes, which were reported in the literature than that from the ordinary GMDR method. For example, the KCNH1 gene has been demonstrated to have strong interaction effects with other genes on the function of insulin secretion  QR-based [30]  M  Out-Res  Gene/SNP-E    A lung cancer prognosis data set, n=246,P=2605 (22 283),q=5  For the real data set, some interesting findings different from the alternatives are observed. Specifically, with quantile = 0.5, the numbers of identified genes are 20 and corresponding interactions are 7 (with age), 7 (gender), 10 (smoking), 15 (chemotherapy treatment) and 22 (stage), respectively  SNPInter- Forest [31]  J  MisSpec  SNP-SNP  https://gwas. biosciencedbc. jp/SNPInterForest/ index.html  A GWAS data set of rheumatoid arthritis from the Wellcome Trust Case Control Consortium, n=3,499,P=10 000  For the real data set, SNPInterForest identifies important interactions with genes having some functions related to rheumatoid arthritis, such as PROK2 and PCDH15  BNN [32]  J  MisSpec  SNP-SNP  https://github.com/ beamandrew/BNN  A GWAS tuberculosis data set, n=104,P=16 925  The top five SNPs selected using BNN are shown located on the same chromosome and within 10–50 MB of loci previously reported as having a statistically significant association with pulmonary tuberculosis susceptibility  Nonpara- Scr [33]  J  MisSpec  SNP-SNP    A framingham heart study data set, n=977,P=349 985  Nonpara-Scr identifies 41 main effect SNPs and 36 SNP-SNP interactions in the real data set. The analysis indicates the importance of epistatic interactions in explaining BMI  SODA [34]  J  Out-Pre  G-G  https://cran.r-project.org/ web/packages/sodavis/ index.html  A Michigan lung cancer data set, n=86,P=5217  For the real data set, SODA selects two main effects and two interaction effects, which contribute substantially to classification accuracy compared with Lasso-Logistic  RankReg [35]  J  Out-Res  Gene/SNP-E  http://works.bepress. com/Shuangge/50/  A lung cancer prognosis data set, n=348,P=200 (22 283),q=2  In the analysis of a lung cancer study, RankReg identifies 11 genes with main effects, among which four and five genes have interactions with age and smoking status. A brief literature search shows that the identified genes may have important implications, including CYP2B6, CPS1, 1E2, HLA-DQA1 and others. Their better prediction performances are also observed  LAD-based [36]  J  Out-Res  Gene/SNP-E  https://github.com/cenwu.  A lung cancer prognosis data set, n=348,P=500 (22 283),q=5  In data analysis, the LAD-based method identifies 21 main gene effects and 10 G-E interactions, with seven G-age interactions, one G-smoke interaction and two G-stage interactions. The identified genes (PSPH, HMGA2, CHL1, PON3 and others) and their interactions with E factors are different from those of alternatives and are biologically sensible with better prediction performance  Approach  M/Ja  Robustb  G-E/G-G  Software  Applicationc  Conclusion  MDR [24, 25]  M  MisSpec  SNP-SNP  https://cran.r-project.org/ web/packages/MDR/ index.html  A sporadic breast cancer case-control data set, n=200,P=10  For the real data set, without any statistically significant main effects, MDR identifies one statistically significant high-order interaction among four SNPs from three different estrogen-metabolism genes, which is missed by the previous studies using logistic regression  W-test [26]  M  MisSpec  SNP-SNP  http://www2.ccrb.cuhk. edu.hk/wtest/download. html.  Two bipolar disorder GWAS data sets: (1) n=5000,P=414 682, (2) n=2420,P=729 304  For the two real data sets, W-test identifies significant interaction pairs that can be replicated. The genes in pairs play central roles in neurotransmission and synapse formation, and are undiscoverable by main effects with low frequency variants  RankCorr [27]  M  MisSpec  G-E    A lung cancer prognosis data set, n=336,P=2500 (22 283),q=5.  RankCorr identifies genes (SFRP1, IGFBP1, AHSG, AGR2 and others) and their interactions with important implications, which are different from those using alternatives. Most of the identified interactions are between genes and smoking status  3G-SPA [28]  M  MisSpec  G-G with SNPs  http://cran.r-project. org/package=SPA3G  (1) A birth weight data set, n = 1511 with 797 SNPs in 186 genes. (2) A yeast eQTL mapping data set, n = 112 with 6229 SNPs in 2956 genes  For the first data set, 23 G-G interactions are detected with the cutoff P-value 0.001, of which the genes include PDGFC, ANG, PTGER3, IL9 and others. For the second data set, 87 pairs are found to have significant interactions at an experiment-wise significant level of 0.05, where the corresponding genes include URA3, BAT2, LUE2 and others. For both data sets, the findings are shown to be biological sensible in the literature  R-GMDR [29]  M  Out-Res  SNP-SNP    A KARE data set with a HOMA-IR phenotype, n=8577,P=585 (327 872)  For the real data set, R-GMDR identifies a larger number of hub genes, which were reported in the literature than that from the ordinary GMDR method. For example, the KCNH1 gene has been demonstrated to have strong interaction effects with other genes on the function of insulin secretion  QR-based [30]  M  Out-Res  Gene/SNP-E    A lung cancer prognosis data set, n=246,P=2605 (22 283),q=5  For the real data set, some interesting findings different from the alternatives are observed. Specifically, with quantile = 0.5, the numbers of identified genes are 20 and corresponding interactions are 7 (with age), 7 (gender), 10 (smoking), 15 (chemotherapy treatment) and 22 (stage), respectively  SNPInter- Forest [31]  J  MisSpec  SNP-SNP  https://gwas. biosciencedbc. jp/SNPInterForest/ index.html  A GWAS data set of rheumatoid arthritis from the Wellcome Trust Case Control Consortium, n=3,499,P=10 000  For the real data set, SNPInterForest identifies important interactions with genes having some functions related to rheumatoid arthritis, such as PROK2 and PCDH15  BNN [32]  J  MisSpec  SNP-SNP  https://github.com/ beamandrew/BNN  A GWAS tuberculosis data set, n=104,P=16 925  The top five SNPs selected using BNN are shown located on the same chromosome and within 10–50 MB of loci previously reported as having a statistically significant association with pulmonary tuberculosis susceptibility  Nonpara- Scr [33]  J  MisSpec  SNP-SNP    A framingham heart study data set, n=977,P=349 985  Nonpara-Scr identifies 41 main effect SNPs and 36 SNP-SNP interactions in the real data set. The analysis indicates the importance of epistatic interactions in explaining BMI  SODA [34]  J  Out-Pre  G-G  https://cran.r-project.org/ web/packages/sodavis/ index.html  A Michigan lung cancer data set, n=86,P=5217  For the real data set, SODA selects two main effects and two interaction effects, which contribute substantially to classification accuracy compared with Lasso-Logistic  RankReg [35]  J  Out-Res  Gene/SNP-E  http://works.bepress. com/Shuangge/50/  A lung cancer prognosis data set, n=348,P=200 (22 283),q=2  In the analysis of a lung cancer study, RankReg identifies 11 genes with main effects, among which four and five genes have interactions with age and smoking status. A brief literature search shows that the identified genes may have important implications, including CYP2B6, CPS1, 1E2, HLA-DQA1 and others. Their better prediction performances are also observed  LAD-based [36]  J  Out-Res  Gene/SNP-E  https://github.com/cenwu.  A lung cancer prognosis data set, n=348,P=500 (22 283),q=5  In data analysis, the LAD-based method identifies 21 main gene effects and 10 G-E interactions, with seven G-age interactions, one G-smoke interaction and two G-stage interactions. The identified genes (PSPH, HMGA2, CHL1, PON3 and others) and their interactions with E factors are different from those of alternatives and are biologically sensible with better prediction performance  aM and J stand for marginal analysis and joint analysis. bMisSpec, Out-Pre and Out-Res stand for the robustness properties toward model mis-specification, outliers/contaminations in predictors and outliers/contaminations in responses, respectively. cn and p are the number of subjects and G factors in the application data set. For the G-E interaction analysis, q is the number of the E factors. For the value of P, the dimension after preprocessing is provided, together with the original dimension reported in parenthesis. Table 1. Summary of the robust genetic analysis approaches (partial list) Approach  M/Ja  Robustb  G-E/G-G  Software  Applicationc  Conclusion  MDR [24, 25]  M  MisSpec  SNP-SNP  https://cran.r-project.org/ web/packages/MDR/ index.html  A sporadic breast cancer case-control data set, n=200,P=10  For the real data set, without any statistically significant main effects, MDR identifies one statistically significant high-order interaction among four SNPs from three different estrogen-metabolism genes, which is missed by the previous studies using logistic regression  W-test [26]  M  MisSpec  SNP-SNP  http://www2.ccrb.cuhk. edu.hk/wtest/download. html.  Two bipolar disorder GWAS data sets: (1) n=5000,P=414 682, (2) n=2420,P=729 304  For the two real data sets, W-test identifies significant interaction pairs that can be replicated. The genes in pairs play central roles in neurotransmission and synapse formation, and are undiscoverable by main effects with low frequency variants  RankCorr [27]  M  MisSpec  G-E    A lung cancer prognosis data set, n=336,P=2500 (22 283),q=5.  RankCorr identifies genes (SFRP1, IGFBP1, AHSG, AGR2 and others) and their interactions with important implications, which are different from those using alternatives. Most of the identified interactions are between genes and smoking status  3G-SPA [28]  M  MisSpec  G-G with SNPs  http://cran.r-project. org/package=SPA3G  (1) A birth weight data set, n = 1511 with 797 SNPs in 186 genes. (2) A yeast eQTL mapping data set, n = 112 with 6229 SNPs in 2956 genes  For the first data set, 23 G-G interactions are detected with the cutoff P-value 0.001, of which the genes include PDGFC, ANG, PTGER3, IL9 and others. For the second data set, 87 pairs are found to have significant interactions at an experiment-wise significant level of 0.05, where the corresponding genes include URA3, BAT2, LUE2 and others. For both data sets, the findings are shown to be biological sensible in the literature  R-GMDR [29]  M  Out-Res  SNP-SNP    A KARE data set with a HOMA-IR phenotype, n=8577,P=585 (327 872)  For the real data set, R-GMDR identifies a larger number of hub genes, which were reported in the literature than that from the ordinary GMDR method. For example, the KCNH1 gene has been demonstrated to have strong interaction effects with other genes on the function of insulin secretion  QR-based [30]  M  Out-Res  Gene/SNP-E    A lung cancer prognosis data set, n=246,P=2605 (22 283),q=5  For the real data set, some interesting findings different from the alternatives are observed. Specifically, with quantile = 0.5, the numbers of identified genes are 20 and corresponding interactions are 7 (with age), 7 (gender), 10 (smoking), 15 (chemotherapy treatment) and 22 (stage), respectively  SNPInter- Forest [31]  J  MisSpec  SNP-SNP  https://gwas. biosciencedbc. jp/SNPInterForest/ index.html  A GWAS data set of rheumatoid arthritis from the Wellcome Trust Case Control Consortium, n=3,499,P=10 000  For the real data set, SNPInterForest identifies important interactions with genes having some functions related to rheumatoid arthritis, such as PROK2 and PCDH15  BNN [32]  J  MisSpec  SNP-SNP  https://github.com/ beamandrew/BNN  A GWAS tuberculosis data set, n=104,P=16 925  The top five SNPs selected using BNN are shown located on the same chromosome and within 10–50 MB of loci previously reported as having a statistically significant association with pulmonary tuberculosis susceptibility  Nonpara- Scr [33]  J  MisSpec  SNP-SNP    A framingham heart study data set, n=977,P=349 985  Nonpara-Scr identifies 41 main effect SNPs and 36 SNP-SNP interactions in the real data set. The analysis indicates the importance of epistatic interactions in explaining BMI  SODA [34]  J  Out-Pre  G-G  https://cran.r-project.org/ web/packages/sodavis/ index.html  A Michigan lung cancer data set, n=86,P=5217  For the real data set, SODA selects two main effects and two interaction effects, which contribute substantially to classification accuracy compared with Lasso-Logistic  RankReg [35]  J  Out-Res  Gene/SNP-E  http://works.bepress. com/Shuangge/50/  A lung cancer prognosis data set, n=348,P=200 (22 283),q=2  In the analysis of a lung cancer study, RankReg identifies 11 genes with main effects, among which four and five genes have interactions with age and smoking status. A brief literature search shows that the identified genes may have important implications, including CYP2B6, CPS1, 1E2, HLA-DQA1 and others. Their better prediction performances are also observed  LAD-based [36]  J  Out-Res  Gene/SNP-E  https://github.com/cenwu.  A lung cancer prognosis data set, n=348,P=500 (22 283),q=5  In data analysis, the LAD-based method identifies 21 main gene effects and 10 G-E interactions, with seven G-age interactions, one G-smoke interaction and two G-stage interactions. The identified genes (PSPH, HMGA2, CHL1, PON3 and others) and their interactions with E factors are different from those of alternatives and are biologically sensible with better prediction performance  Approach  M/Ja  Robustb  G-E/G-G  Software  Applicationc  Conclusion  MDR [24, 25]  M  MisSpec  SNP-SNP  https://cran.r-project.org/ web/packages/MDR/ index.html  A sporadic breast cancer case-control data set, n=200,P=10  For the real data set, without any statistically significant main effects, MDR identifies one statistically significant high-order interaction among four SNPs from three different estrogen-metabolism genes, which is missed by the previous studies using logistic regression  W-test [26]  M  MisSpec  SNP-SNP  http://www2.ccrb.cuhk. edu.hk/wtest/download. html.  Two bipolar disorder GWAS data sets: (1) n=5000,P=414 682, (2) n=2420,P=729 304  For the two real data sets, W-test identifies significant interaction pairs that can be replicated. The genes in pairs play central roles in neurotransmission and synapse formation, and are undiscoverable by main effects with low frequency variants  RankCorr [27]  M  MisSpec  G-E    A lung cancer prognosis data set, n=336,P=2500 (22 283),q=5.  RankCorr identifies genes (SFRP1, IGFBP1, AHSG, AGR2 and others) and their interactions with important implications, which are different from those using alternatives. Most of the identified interactions are between genes and smoking status  3G-SPA [28]  M  MisSpec  G-G with SNPs  http://cran.r-project. org/package=SPA3G  (1) A birth weight data set, n = 1511 with 797 SNPs in 186 genes. (2) A yeast eQTL mapping data set, n = 112 with 6229 SNPs in 2956 genes  For the first data set, 23 G-G interactions are detected with the cutoff P-value 0.001, of which the genes include PDGFC, ANG, PTGER3, IL9 and others. For the second data set, 87 pairs are found to have significant interactions at an experiment-wise significant level of 0.05, where the corresponding genes include URA3, BAT2, LUE2 and others. For both data sets, the findings are shown to be biological sensible in the literature  R-GMDR [29]  M  Out-Res  SNP-SNP    A KARE data set with a HOMA-IR phenotype, n=8577,P=585 (327 872)  For the real data set, R-GMDR identifies a larger number of hub genes, which were reported in the literature than that from the ordinary GMDR method. For example, the KCNH1 gene has been demonstrated to have strong interaction effects with other genes on the function of insulin secretion  QR-based [30]  M  Out-Res  Gene/SNP-E    A lung cancer prognosis data set, n=246,P=2605 (22 283),q=5  For the real data set, some interesting findings different from the alternatives are observed. Specifically, with quantile = 0.5, the numbers of identified genes are 20 and corresponding interactions are 7 (with age), 7 (gender), 10 (smoking), 15 (chemotherapy treatment) and 22 (stage), respectively  SNPInter- Forest [31]  J  MisSpec  SNP-SNP  https://gwas. biosciencedbc. jp/SNPInterForest/ index.html  A GWAS data set of rheumatoid arthritis from the Wellcome Trust Case Control Consortium, n=3,499,P=10 000  For the real data set, SNPInterForest identifies important interactions with genes having some functions related to rheumatoid arthritis, such as PROK2 and PCDH15  BNN [32]  J  MisSpec  SNP-SNP  https://github.com/ beamandrew/BNN  A GWAS tuberculosis data set, n=104,P=16 925  The top five SNPs selected using BNN are shown located on the same chromosome and within 10–50 MB of loci previously reported as having a statistically significant association with pulmonary tuberculosis susceptibility  Nonpara- Scr [33]  J  MisSpec  SNP-SNP    A framingham heart study data set, n=977,P=349 985  Nonpara-Scr identifies 41 main effect SNPs and 36 SNP-SNP interactions in the real data set. The analysis indicates the importance of epistatic interactions in explaining BMI  SODA [34]  J  Out-Pre  G-G  https://cran.r-project.org/ web/packages/sodavis/ index.html  A Michigan lung cancer data set, n=86,P=5217  For the real data set, SODA selects two main effects and two interaction effects, which contribute substantially to classification accuracy compared with Lasso-Logistic  RankReg [35]  J  Out-Res  Gene/SNP-E  http://works.bepress. com/Shuangge/50/  A lung cancer prognosis data set, n=348,P=200 (22 283),q=2  In the analysis of a lung cancer study, RankReg identifies 11 genes with main effects, among which four and five genes have interactions with age and smoking status. A brief literature search shows that the identified genes may have important implications, including CYP2B6, CPS1, 1E2, HLA-DQA1 and others. Their better prediction performances are also observed  LAD-based [36]  J  Out-Res  Gene/SNP-E  https://github.com/cenwu.  A lung cancer prognosis data set, n=348,P=500 (22 283),q=5  In data analysis, the LAD-based method identifies 21 main gene effects and 10 G-E interactions, with seven G-age interactions, one G-smoke interaction and two G-stage interactions. The identified genes (PSPH, HMGA2, CHL1, PON3 and others) and their interactions with E factors are different from those of alternatives and are biologically sensible with better prediction performance  aM and J stand for marginal analysis and joint analysis. bMisSpec, Out-Pre and Out-Res stand for the robustness properties toward model mis-specification, outliers/contaminations in predictors and outliers/contaminations in responses, respectively. cn and p are the number of subjects and G factors in the application data set. For the G-E interaction analysis, q is the number of the E factors. For the value of P, the dimension after preprocessing is provided, together with the original dimension reported in parenthesis. First, it is noted that different methods may be applicable to different types of G factors. This is mainly because of the distribution of G factor, which can be continuous or discrete. For example, SNPs are usually coded with three categories (0, 1 and 2) for genotypes (aa, Aa and AA), respectively. Therefore, some SNP interaction analysis methods are often based on contingency tables and statistical tests. On the other hand, gene expressions generally have continuous measurements, where the regression models are often adopted and heavy-tailed distributions are observed more frequently. To provide a deeper insight of each approach, we discuss their specific applicable types of G factors, which are also listed in Table 1. Second, despite similarities, there are still certain distinctions between models for G-E and G-G interaction analysis. E factors studied in most existing G-E interaction analysis are usually preselected with evidence from previous studies and have low dimensions. Thus, many methods always include E factors in the model without further selection. Compared with G-E interaction analysis, G-G interaction analysis often has much higher computational cost. It also has more ‘symmetrical’ form, as the interactions are between two G factors with similar distributions and other properties. These make many methods for identifying G-G interactions not directly applicable to analyze G-E interactions, and vice versa. The types of interactions studied in each model are demonstrated in the following section. Marginal analysis The most prominent feature of marginal analysis is that, each time, the number of variables (interactions + main effects) is considerably smaller than the sample size. Thus, standard estimation and inference approaches can be applied. In biomedical studies, marginal analysis is still highly popular because of its computational simplicity and stability. Robust to model mis-specification As mentioned in the first section and also discussed thoroughly in low-dimensional studies [37], with practical data, simply assuming a popular model (for example, logistic model for binary data) can be ‘dangerous’. With high-dimensional data, model diagnostics techniques, which may suggest model forms with low-dimensional data, have not been well developed. Thus, a strategy to achieve robustness to model mis-specification is to ‘avoid’ making specific model form assumptions. Three representative methods are reviewed below. Multifactor dimensionality reduction: Multifactor dimensionality reduction (MDR) is a model-free approach first proposed by [24] to identify genetic interactions for case-control studies. The most common application of MDR is in detecting SNP-SNP interactions and proceeds as follows. First, for each pair of genetic factors, their possible multifactor classes are represented in a two-dimensional space. Then, using the training data generated by V-fold cross-validation, each multifactor class is labeled as either ‘high-risk’ or ‘low-risk’ based on the ratio of the number of cases to the number of controls. Next, the prediction error of this model is estimated using the testing data. Finally, a permutation test is conducted to compare the average cross-validation consistency from the observed data to the distribution of average consistency under the null hypothesis of no association. The identification of important interactions can then be based on the P-values. The empirical studies demonstrate that MDR outperforms the parametric logistic regression. Software of MDR written in R is available at https://cran.r-project.org/web/packages/MDR/index.html [25]. Despite considerable successes, MDR is shown to be more computationally intensive than the generalized linear model, especially when the number of analyzed SNPs is >10 [24]. Many researchers have extended the framework of MDR to be more flexible and handle different types of outcomes. For example, [38] introduces a generalized MDR (GMDR) framework using the residual-based score of a generalized linear model, and this approach can accommodate dichotomous and continuous responses and has lead to some interesting findings [39]. The survival MDR (Surv-MDR) is proposed by [40] for survival data, which replaces the case-control ratios by the logrank test statistics. W-test and others: W-test is introduced by [26] to detect G-G interaction effects and measures the distributional difference between cases and controls through a combined log-odds ratio. It is a model free test in that it makes no assumption on the genotypic effect model and distribution. Similar to MDR, at the first stage, a 2 × 9 contingency table is generated to test the association for a pair of SNPs. For each cell of the table, the distribution is computed based on the number of subjects for the case and control groups, respectively. Then, a W-test statistic is developed to measure the discordance between the two distributions using the normalized log odds ratio adjusted by a data-adaptive factor. It follows a chi-squared distribution with the degree of freedom estimated from the covariance structure of the contingency table. Different from the logistic regression under which the odds ratio has a linear relationship with the interaction, the W-test approach can accommodate nonlinear interactions, as no specific genetic model is assumed. It is shown in simulation that it has robust power and reasonable type I error compared with logistic regression. The R code implementing this approach is available at http://www2.ccrb.cuhk.edu.hk/wtest/download.html. Using a laptop computer with 2.4 GHz CPU and 8 GB memory, W-test takes about 7.3 s to accomplish the analysis of one simulation with 1225 SNP-SNP interactions and 1000 subjects, compared with 7.7 s for chi-squared test and 77.7 s for MDR. It is observed that W-test has significant advantage over MDR in terms of computational cost. There are other robust tests for evaluating the associations of SNP-SNP interactions with a trait, such as the weighted rank test [41], mutual information-based test [42] and forward U-test [43]. They can be used to identify interactions without making stringent model assumption and are robust toward model mis-specification. Rank correlation-based methods [27]: The methods propose a robust rank correlation-based method (RankCorr) under a general semiparametric transformation model for the identification of G-E interactions. Strictly speaking, this method is not completely model-free. However, it makes much weaker assumptions on the model form and thus can be more robust than methods based on, for example, logistic or Cox models. In numerical studies, the continuous G factors are analyzed. However, this model can also accommodate discrete G factors. Denote y as the response variable, which can be a continuous marker, categorical disease status or survival time. Denote x=(x1,…,xp) as the p G factors and z=(z1,…,zq) as the q E factors. Assume n iid observations. For G factor j(=1,…,p), the following semiparametric transformation model is assumed:   E(gj(y))=∑k=1qzkαjk+xjβj+∑k=1qxjzkγjk=θj′wj, where gj(·) is a monotone increasing function with an unknown form, wj=(z,xj,zxj)′, αj=(αj1,…,αjk)′ and βj represent the main E and G effects, respectively, γj=(γj1,…,γjk)′ corresponds to G-E interactions, and θj=(αj′,βj,γj)′. A penalized loss function, O(θj)+ρ(θj), is proposed for parameter estimation, where O(θj) is the robust loss function, and ρ(θj)=∑k=12q+1ρ(θjk;λk) is the MCP penalty function with ρ(t;λ)=λ∫0|t|(1−xλξ)+dx. This method has two key differences from the existing alternatives. The first is that, without a model form assumption and hence likelihood function, a robust rank correlation-based loss function is constructed. For example, when y is continuous,   O(θj)=1n(n−1)∑i≠lI(yi≥yl)I(θj′wij≥θj′wlj), where I(·) is the indicator function. Robust loss functions for other types of responses can be constructed in a similar manner. The second difference is that a penalization approach is adopted for stabilizing estimates and selecting important effects. Simulation shows that this method outperforms the nonrobust model-based alternatives with more accurate identification, where it can achieve robustness to model mis-specification as well as outliers in response. A coordinate descent algorithm is developed for the optimization of RankCorr, where smoothed approximation and local linear approximation are adopted for the rank correlation objective function and MCP penalty, respectively. With fixed tunings, the analysis for one gene, five E factors and their interactions takes 0.46 s on a regular desktop computer. Kernel machine methods and others: [28] Kernel machine methods and others propose a model-based kernel machine method to identify significant G-G interactions for continuously valued quantitative traits. Different from the single SNP-based pairwise interaction analyses, this method takes genes as the testing units and identifies interactions at the gene level. The proposed gene-centric G-G Smoothing-sPline analysis of variance (ANOVA) model (3 G-SPA) can accommodate complex nonlinear SNP effects within a gene and imposes fewer assumptions on the underlying effects, making it less sensitive to model mis-specification. A two-step procedure is developed, which includes an exhaustive search for all pairwise interactions and assessment of significance of interactions for the identified gene pairs. In the first step, let yi and xi=(xi1,…,xiL) be the phenotype and genotype vector of the gene pair for subject i, where L=L1+L2 is the total number of SNPs in the two genes under investigation. A nonparametric regression model is adopted:   yi=m(xi)+εi,  i=1,2,…,n, where m(·) is an unknown function, and εi is a random error. This method is robust in that the form of m(·) does not need to be specified. The function m(·) is then decomposed into the main effects of the two genes and the interaction effect between them based on a functional ANOVA decomposition and reproducing kernels. In the second step, a score test statistic is developed to detect important interactions. Simulation shows that 3 G-SPA has better performance than a regression-based principle component (PC) analysis method with the linearity assumption. The R package SPA3G for implementing 3 G-SPA is available at http://cran.r-project.org/package=SPA3G. For 1000 simulation replicates, which include 500 subjects and two separate genes with 10 SNPs each, the computational time for 3 G-SPA and two PC-based approaches partial PCA (pPCA) and full PCA (fPCA) are 28, 24.5 and 24 min, respectively. 3 G-SPA has been extended in [44] to analyze case-control data based on a mixed effects logistic model. In addition, a multiple-kernel method is introduced in [45] to assess G-E interactions and can accommodate various trait types, such as continuous, binary and survival. Besides the kernel methods, other nonparametric techniques have also been developed for the marginal analysis of genetic interactions, including, for example, the multivariate adaptive regression splines [46] and regularized isotonic regression [47]. Robust to outliers/contaminations in predictors As described in ‘Introduction’ section and in the literature, outliers and contaminations in the predictors can happen to both E and G factors. However, our limited literature search suggests that there is still a lack of genetic interaction analysis methods tailored to accommodate outliers and contaminations in predictors. With the link between marginal and joint analyses, it is possible that methods developed for joint analysis can be ‘marginalized’ for marginal analysis. In addition, a few studies might have potential robustness properties toward outliers and contaminations in predictors but have not been well examined. Specifically, the L-estimator GMDR is proposed by [29] with the goal of accommodating outliers in response, and is developed based on least trimmed squared regression (LTS). LTS has been shown to be robust to outliers in both predictor and response spaces in the literature [48]. Thus, it is reasonable to conjecture that the L-estimator GMDR may also achieve robustness toward outliers in predictors. More theoretical and numerical analyses are needed to examine this robustness. The details of this approach are discussed in next subsection. Robust to outliers/contaminations in responses Robust GMDR [29]: It develops a robust GMDR (R-GMDR) approach for SNP-SNP interaction analysis based on the L-estimator and M-estimator estimation methods. The proposed method can reduce the effects caused by outlying responses. Specifically, for the L-estimator GMDR, the LTS regression, as opposed to the least squared regression, is conducted for generalized linear models. A studentized trimmed residual is proposed as the score for GMDR, where a threshold value is used to determine whether the ith sample should be regarded as an outlier and removed. For the M-estimator GMDR, an M-estimator type score based on the Tukey’s biweight function is developed, where a threshold T is used to shrink all residuals exceeding T. In extensive simulations where the response follows a skewed distribution or a normal distribution with a huge variance, both the L-estimator GMDR and M-estimator GMDR are shown to outperform the original GMDR. The original GMDR software is available at http://www.ssg.uab.edu/gmdr/. It is expected that the two robust methods can be implemented by modifying the original GMDR software. In addition, it is demonstrated that M-estimator GMDR has similar computational cost with the ordinary GMDR, and L-estimator GMDR takes slightly more computational time since an additional step to compute residuals is needed. Researchers have also proposed other improvements of MDR-related methods to accommodate long-tailed errors and outliers in responses. For example, for the accelerated failure time (AFT)-MDR designed for survival phenotypes, two improvements are introduced to reduce the effects of extreme values on the sum of standardized residuals in the AFT-MDR [49]. The first approach transforms the continuous standardized residuals into binary variables and then uses the original MDR instead of AFT-MDR. Under the second approach, the extreme values of the standardized residuals beyond the specified bounds are replaced by either the lower or upper bound. Both methods are demonstrated to have higher powers than the original AFT-MDR with the presence of outliers. Quantile regression and others [30]: Quantile regression and others propose a quantile regression-based (QR-based) method to identify G-E interactions for prognosis responses. As with other QR-based methods [50], it is robust to contamination in the responses. Both continuous and discrete G factors are analyzed in numerical studies, mimicking gene expression and SNP data, respectively. Here, we use notations similar to those for the rank correlation-based method reviewed above. Denote y as the survival time. The following check loss is adopted for G factor j(=1,…,p):   ρτ(y−∑k=1qzkαjk−xjβj−∑k=1qxjzkγjk), where the function ρτ=u{τ−I(u<0)} with 0<τ<1. Based on this loss function, two weighted methods are proposed to accommodate right censoring of the survival time. Important interactions can be identified based on significance level as well as using a penalization approach. It is observed in simulation that under data contamination and heavy-tailed errors, the robust methods can significantly outperform the nonrobust ones. Two techniques are involved in the computation process, including the majorize–minimization (MM) approach and local quadratic approximation. With more complex algorithm, the QR-based methods tend to be slower than the nonrobust methods, such as that based on least squared loss. In the literature [51, 52], the robust check loss function has also been used in the framework of regression tree for interaction analysis, and such an approach can be potentially adapted for genetic interaction analysis. The above ‘robust loss function + penalized identification’ strategy has also been adopted with other robust loss functions. For instance, [19] adopts an exponential squared loss for the G-E interaction analysis with prognosis data. This approach can down-weigh the influence of observations with larger residuals. Joint analysis Compared with marginal analysis, joint analysis may better reflect the fact that molecular mechanisms under complex diseases involves multiple genes (and their interactions). When a large number of interactions and main effects are jointly analyzed, because of the high dimensionality, complex statistical methods are usually needed. In recent statistical literature, joint analysis has been popular [53, 54]. In biomedical studies, its popularity has been increasing fast. Recently, the ‘main effects, interactions’ hierarchical constraint is often used in joint analysis, where an interaction is allowed into the model only if the corresponding main effects are also identified [53–56]. Although this constraint is not strictly validated in biomedical studies, it is suggested to achieve better estimation and interpretation [53]. More detailed discussions are available in the literature [3, 57]. For the following approaches, their properties of respecting hierarchy are discussed briefly. Robust to model mis-specification Random forest: Random forest is an ensemble of regression trees and can predict an outcome based on a large number of predictors. A novel method SNPInterForest is proposed by [31] for detecting G-G interactions by modifying the construction of the random forest. It is a model-free approach and does not require (semi)parametric model specification. Specifically, in the first step, a random forest is constructed using case-control data with SNPs as categorical variables. Different from the original random forest under which variables are selected based on individual effects, a combination of multiple SNPs as well as a single SNP is considered when choosing a split variable at each node. Then, the interactions of variables are represented within branches of trees, which are the paths from the root nodes to the leaf nodes. Next, for each pair of SNPs, its interaction strength is measured by the number of simultaneous appearance in the same branch over all trees. That is, an SNP-SNP interaction is relatively more important in affecting the disease outcome if it appears more frequently in the same branch. Finally, the above interaction score is normalized by using its respective baseline level, which is based on the hypothesis that the SNP combinations do not interact. The random forest methods usually do not consider the ‘main effects, interactions’ hierarchy and are even more focusing on interactions without main effects. In extensive simulations under a wide spectrum of disease models, SNPInterForest is shown to achieve better performance than the model-based approach BOOST. Software for SNPInterForest is available at https://gwas.biosciencedbc.jp/SNPInterForest/index.html. For the analysis of one real GWAS data with 3499 individuals and 10 000 SNPs, the running time of SNPInterForest is about 98 h using a computer with 2.67 GHz CPU and 6 GB memory, compared with 11 min for BOOST and 5 h for SNPHarvester. There are several other genetic interaction studies that are also under the framework of random forest [58, 59]. Different approaches may have minor difference in the construction of forest and definition of interaction score. Neural network: Researchers have developed a series of neural network methods to model and detect the functional relationships between genetic interactions and response variables. Neural network is a model-free approach and is capable of dealing with a large amount of data. It has been widely applied to identify SNP-SNP and SNP-E interactions. Similar to random forest, the neural network approaches also usually focus on identifying interactions with weak main effects or without main effects and do not take hierarchical constraints into account. In general, the underlying structure of the neural network for genetic interaction analysis is a weighted and directed graph, where the nodes represent genetic factors and the edges represent genetic interactions. The graph is divided into multiple layers. The input layer contains all genetic factors under consideration, and each node is connected to one or more nodes (genetic factors) in a hidden layer. Nodes in the hidden layer are further connected to either additional hidden layers or to an output node representing the response. Signals pass from the input layer to the output layer through the interactions with weights and activation functions. Different algorithms are proposed to learn the weights and activation functions in genetic interaction studies, including for example the genetic programming neural network [60], grammatical evolution neural networks [61], feed-forward multilayer perceptron [62] and Bayesian neural network (BNN) [32]. Simulations in [62] and [63] show that the neural network is more accurate in identifying interactions than the nonrobust logistic regression model. There are several publicly available softwares for the above approaches. For example, the code implementing BNN is available at https://github.com/beamandrew/BNN. Compared with MDR, BNN is shown to be able to analyze a larger number of SNPs in a relatively short time. Specifically, the analysis of a data set with 104 subjects and 500 000 SNPs can be accomplished within 30 min. Nonparametric screening [33]: It proposes a nonparametric joint modeling method for detecting SNP-SNP interactions based on variable screening and variable selection (Nonpara-Scr). The associations between the phenotype and genotypes are characterized by an unknown function. Robustness is achieved by not assuming the specific model form. Use notations similar to those above. This approach considers the following generic nonparametric model:   yi=μ+∑j=1pfj(xij)+∑j=k+1p∑k=1phjk(xijxik)+εi, where μ is the population mean, fj(·)’s and hjk(·)’s are the unknown functions representing the main effects and interactions, respectively, and εi is the random error. The basis expansion technique is used to approximate the nonparametric functions, that is:   fj(u)=∑lcjlϕl(u),  hjk(u)=∑ldjklϕl(u), (1) where ϕl(·) is the lth basis function, and cjl and djkl are the corresponding coefficients. A hybrid of variable screening and variable selection is proposed for interaction detection. First, consider the model with only one main effect:   yi=μ+∑lcjlϕl(xij)+εi,  for j=1,…,p. The subset of important main effects is then defined as M^1={1≤j≤p:∑lc^jl2>C} with a threshold C. Next, the following model is considered to screen interactions:   yi=μ+∑j∈M^1p∑lcjlϕl(xij)+∑j=1p∑k∈M^1p∑ldjklϕl(xijxik)+εi. The subset of selected interactions is estimated by M^2={1≤j≤p,k∈M^1:∑ld^jkl2>C′} for some C′>0. A variable selection method based on LASSO is further used for the identification of final interactions following the above screening procedure. In the variable screening step, the weak hierarchy is taken into account, where the interactions can be included in the model if at least one of the corresponding main effects is selected. However, Lasso is then adopted directly in the variable selection step and cannot guarantee the hierarchical structure of ‘main effects, interactions’. In simulation studies, it is shown that the proposed method has an outstanding performance in selecting important interactions compared with parametric models, such as the forward LASSO. In addition, the proposed procedure is demonstrated to be computationally efficient. Robust to outliers/contaminations in predictors Stepwise cOnditional likelihood variable selection for Discriminant Analysis [34]: Stepwise cOnditional likelihood variable selection for Discriminant Analysis (SODA) proposes SODA to detect both main and quadratic interaction effects. The continuous mRNA expression measurements are analyzed in this study. This method does not require the normality (or sub-Gaussian) assumption on predictors, which is usually made in other studies, leading to much enhanced robustness towards long-tailed distributed predictors and outliers. The SODA approach consists of three stages: a preliminary forward main effect selection, a forward variable selection (for both main and interaction effects) and a backward elimination. A key point of the above procedure is that the extended Bayesian information criterion (EBIC) is used as the variable selection/elimination criterion at each iteration. Specifically, for the K-class classification problem, let y∈{1,…,K} denote the class label and x be the p-dimensional predictor vector. The following logistic model is adopted:   P(y=k|x,θ)= exp ⁡[δk(x|θ)]1+∑l=1K−1 exp ⁡[δl(x|θ)], (2) where δk(x|θ) is the discriminant function for class k, and θ denotes the parameter vector. The function δk(x|θ) is further defined as:   δk(x|θ)={αk+βk′x+x′Akx,k=1,…,K−1,0,k=K. (3) This method does not need to model the distribution of x and includes the quadratic discriminant analysis as a special case, under which the predictors are required to follow a multivariate normal distribution. At each iteration with the considered predictor set S, the EBIC is estimated based on the log-likelihood function derived from equation (2) with the predictors in S. The proposed method does not respect the hierarchical constraints. In simulations with multiple predictor distributions, the proposed method is observed to achieve higher identification and classification accuracy, compared with the nonrobust methods including the backward procedure and forward–backward method. The R package sodavis for SODA is available at https://cran.r-project.org/web/packages/sodavis/index.html. To analyze one simulated data set with 1000 subjects and 1000 genes, SODA takes about 4 min, the backward procedure takes 22 min and Hier [53] takes >24 h. Robust to outliers/contaminations in responses A rank-based penalization approach: In [35] for the identification of G-E interactions, a rank-based regression (RankReg) is adopted, which does not make stringent distributional assumptions on the random error and can accommodate long-tailed errors and hence outliers/contaminations in the response. With the regression framework, this method can be applied to analyze both continuous and discrete G factors. For subject i, let yi be the response, xi=(xi0,xi1,…,xip) be the (p+1)-dimensional vector of G factors with xi0=1 and Ui and Ei be the continuous and discrete E factors, respectively. The following partially linear varying coefficient (VC) model is adopted for joint analysis:   yi=∑j=0pαk(Ui)xij+∑j=0pβjxijEi+εi, where αk(·) is the smooth varying-coefficient (VC) function and εi is the random error. The VC function is approximated using a similar basis expansion technique as in equation (1). Beyond being robust to outliers/contaminations, the proposed method can also accommodate unspecified nonlinear effects of E factors, and thus, has a second robustness property. The rescaled rank-based loss function   1n∑i<l|εi−εl| is considered for estimating the unknown functions (coefficients). For regularized estimation and more importantly selection of important effects, a penalization approach is adopted. The proposed method does not respect the ‘main effects, interactions’ hierarchy because of computation consideration. However, it is demonstrated that one can first apply the proposed method and then add back the main effects if the hierarchy is needed. Simulation shows that the proposed method has satisfactory performance and outperforms several competing alternatives including the nonrobust approach based on the ordinary least squared loss function. The computer code written in R is available at http://works.bepress.com/Shuangge/50/. The proposed objective function is optimized using the group coordinate descent algorithm. The analysis of one simulated replicate with n = 200, P = 200 and two E factors takes about 1.6 min using a regular laptop. It is also shown that the computational cost of proposed method grows linearly with n and p. A least absolute deviation based method: The least absolute deviation (LAD) loss is a special case of the check loss and well known for its robustness to long-tailed errors and outliers in responses. This method [36] proposes a penalized robust approach to dissect G-E interactions based on the LAD loss for prognosis data. It has been applied to both continuous and categorical predictors, corresponding to gene expression and SNP data. Use notations similar to those above. Let T be the logarithm of survival time. The AFT model is assumed under joint modeling, that is:   T=∑k=1qαkzk+∑j=1pβjxj+∑j=1p∑k=1qxjzkγjk+ε, where ε is the random error with an unknown distribution function. It is especially noted that the random error can have a long-tailed distribution or be contaminated. Based on the Kaplan–Meier weights wi(i=1,2,…,n) for accommodating the right censoring of survival time, the following weighted LAD loss function is proposed:   ∑i=1nwi|yi−∑k=1qαkzik−∑j=1pβjxij−∑j=1p∑k=1qxijzikγjk|, where yi=min⁡(Ti,Ci), and Ci is the log censoring time. A two-stage penalty proposed in [64] is further added for the identification of interactions respecting the strong ‘main effects, interactions’ hierarchy. Simulation demonstrates the advantage of the proposed approach over the alternatives without the robustness property, which are based on least squared loss. Research code written in R is publicly available at https://github.com/cenwu. A group coordinate descent algorithm is developed for the proposed method. For one simulated replicate with 200 subjects, 200 G factors and 10 E factors, the proposed LAD-based method takes 12.8 s to accomplish the analysis given fixed tunings, compared with 56.5, 57.4 and 57.3 s taken by the other three alternatives. Some insights into the numerical performance of different methods The implementation of all the reviewed methods is a prohibitive task and beyond the scope of this review. As there are not publicly available softwares for many approaches, strictly numerical comparisons of different methods cannot be conducted. Following the literature [23], we summarize the simulation settings and results reported in some of the reviewed papers in Table 2, to provide some insights into their numerical performance and partly show their advantages and shortcomings. It is observed that the performance of these methods depends on the underlying genetic models, the distributions of predictors and response (error), data dimensionality, sample size and others. Among them, no method can perform consistently better than others under all scenarios. In general, the robust methods have relatively poor or competitive performance compared with the nonrobust ones when the models are not mis-specified and the data do not have contaminations and outliers. This result is as expected, as the robustness can be achieved usually at the cost of some loss in efficiency [65]. When there are model mis-specification and data contamination, the robust interaction analysis methods mostly have significant advantages. Table 2. Summary of comparisons of the robust genetic analysis approaches against alternatives (partial list) References  Proposed/alternative  Main settings  Major conclusions  Wang et al. [26]  W-test/logistic regression; chi-squared test; MDR  n=1000,P=50; Discrete predictors: six different genetic architectures with minor allele frequency (MAF) in the common range and in the low frequency range, and LD in the high range, medium range and low range; binary response: (1) d=(4,2)a, linear interaction model with main effect, (2) d=(0,2), nonlinear interaction model without main effect  (1) The model-free W-test and MDR outperform the logistic regression even under the linear model, with higher power and lower type I error rate  (2) Their advantages become more prominent when the underlying model is nonlinear. For one genetic architecture, the power is 5.9% (logistic), 72.6% (Chi-squared), 95.5% (MDR) and 88% (W-test). W-test behaves much better than MDR when 1%<MAF<5%, where the average power of W-test and MDR are 87.6 and 73.2%, respectively  Shi et al. [27]  RankCorr/P-value-based significance analysis (Sig); (1) Penalized least squares (PeLS), (2) Penalized Cox model (PeCOX)  n=200,P=500 1000,q=5, d=(3,4,20)b; Continuous predictors: follow multivariate normal distribution with two correlation structures and two correlations; (1) Continuous response under linear model, (2) Censored response under transformation model: g(t)=t3,t5,t and t1/3; Error distributions: N(0, 1), 0.7N(0,1)+0.3Cauchy(0,1) and t(1)  (1) When the error has a normal distribution, Sig has the best performance. However, when the error is contaminated or has a long tail, RankCorr can significantly outperform the alternatives (2) RankCorr has the largest AUC under all simulation scenarios, including those with linear transformation function and normal distribution error, followed by the nonrobust alternatives PeCOX and Sig  Li and Cui [28]  3G-SPA/Partial PCA (pPCA); Full PCA (fPCA)  n = 500, 1000, two separate genes with 10 SNPs each; discrete predictors: generated via MS program; Continuous response under linear model: (1) without genetic effects, (2) with main effects but no interaction, (3) with interaction effect smaller than main effects, (4) with interaction effect larger than main effects  3G-SPA has the highest power with the type I error reasonably controlled. The PCA-based methods are more conservative. The power of 3G-SPA is not depending on whether the genetic effect is because of its main effect or interaction, whereas the PCA-based methods prefer the scenarios with larger interactions  Kim and Park [29]  L-estimator GMDR; M-estimator GMDR/Ordinary GMDR  (1) n=1000,P=100, (2) n=3000,P=1000, d=(2,1); Discrete predictors: follow multinomial distribution with two different probabilities; Continuous response under linear model; Error distributions: (1−τ)N(0,1)+τN(0,400) and (1−τ)N(0,1)+τGamma(α,β) with τ=0,0.033 and 0.067  When there are no outliers (τ = 0), the ordinary GMDR and M-estimator GMDR outperform L- estimator GMDR. For all the designs with outliers, the performance of L-estimator and M-estimator GMDR are better than that of ordinary GMDR, with higher power. In some designs, L-estimator GMDR outperforms M-estimator GMDR  Yoshida and Koike [31]  SNPInterForest/SNPHarvester; BOOST  (1) n=4000,P=1000; Discrete predictors: the genetic model with different values of LD and MAFs; Binary response: the additive, multiplicative, and threshold models of interactions with weak marginal effects (2) n=400,P=1000; Discrete predictors: same as (1); Binary response: the pure interaction model without marginal effects  (1) The performance of SNPInterForest is superior to two alternatives in terms of their power SNPInterForest and SNPHarvester do not identify any spurious associations, whereas BOOST falsely identifies five interactions (2) The SNPInterForest and BOOST outperform SNPHarvester, with higher power and better identification accuracy The false discovery rates of SNPInterForest and SNPHarvester are much lower than that of BOOST  Beam et al. [32]  (1) BNN/BEAM; χ2 test (2)BNN/BEAM; χ2 test; GBM; MDR  n=2000,P=1000; Discrete predictors: the genetic model with different values of MAF; (1) Binary response: three biallelic models with main effects, (2) Binary response: pure interaction models without main effects  (1) BNN is uniformly more powerful than the χ2 test under all the simulated genetic models, and it outperforms BEAM in most instances (2) BNN performs better than BEAM, χ2 test and GBM under most of the simulations, and has similar performance with MDR. But MDR is much more time-consuming  Li et al. [33]  Nonpara-Scr/Screen and Clean (SC); SIS+LASSO; Forward LASSO; MDR  n=1000, P=3512, d=(12,8); Discrete predictors: dichotomizing the continuous variables generated from normal distribution with different correlations; Continuous response: a genetic model with both additive effects and dominant effects, considering three different noise levels  Nonpara-Scrbehaves better than the other four alternatives, with higher power and much smaller false-positive rates. Their average power are 97% (Nonpara-Scr), 90% (SIS+LASSO), 80% (SC), 48% (Forward LASSO) and 64% (MDR)  Li and Liu [34]  SODA/ Backward procedure (Back); Forward–backward method (ForBack); hierNet; IIS-SQDA  Ten different sample sizes with n=100, 100×101/9,…,1000; (1) P=50,d=(1,4); Continuous predictors: generated from normal, non-normal and heteroskedastic distributions, (2) P=1000,d=(1,4); Continuous predictors: non-normal distribution (3) P=50,d=(0,4), Continuous predictors: same as (2), (4) P=50,d=(2,4), Continuous predictors: same as (2); Binary response: Bayes classification rule  With Gaussian predictors, SODA, Back and ForBack can detect all important predictors as n increases, with small false positives, whereas IIS-SQDA and hierNet select too many false positives, resulting in high test error rates. Under other predictor settings, SODA performs much better than the nonrobust alternatives and achieves near-oracle classification accuracy when sample size is large enough  Wu et al. [35]  RankReg/LSregression-based method; Robust marginal analysis (RobMA)  n=200, P=200, q=2, d=(0,0,20); (1) Discrete predictors: MS program and Hardy–Weinberg equilibrium-based method, (2) Continuous predictors: follow multivariate normal distribution with an auto-regressive correlation structure; Continuous response under the linear model; Error distributions: N(0,1),t(3),0.85N(0,1)+0.15 log Normal(0,1) and 0.85N(0,1)+0.15Cauchy(0,1)  With N(0, 1) error distribution, RankReg is slightly inferior to the unrobust alternative LS, but still outperforms RobMA. With other three error types, RankReg has significant advantages, with more true positives, much fewer false positives, and higher estimation accuracy  Wu et al. [36]  LAD-hier/LAD-nhier; LAD-group; LS-hier; LS-nhier; LS-group  n=200, P=200, q=10, d=(10, 10, 15); (1) Continuous predictors: follow multivariate normal distribution with the banded and autoregressive correlation structure, (2) Discrete predictors: dichotomizing the above continuous values at the first and third quartiles; Censored response under the AFT model; Error distributions: N(0,1),0.80N(0,1)+0.20 log Normal(0,1),t(2) and 0.80N(0,1)+0.20Cauchy(0,1)  Without data contamination, the robust LAD-hier and the nonrobust counterpart LS-hier have comparable performances. With heavy-tailed errors, LAD-hier is significantly superior to the alternatives, in terms of both identification and prediction performance  References  Proposed/alternative  Main settings  Major conclusions  Wang et al. [26]  W-test/logistic regression; chi-squared test; MDR  n=1000,P=50; Discrete predictors: six different genetic architectures with minor allele frequency (MAF) in the common range and in the low frequency range, and LD in the high range, medium range and low range; binary response: (1) d=(4,2)a, linear interaction model with main effect, (2) d=(0,2), nonlinear interaction model without main effect  (1) The model-free W-test and MDR outperform the logistic regression even under the linear model, with higher power and lower type I error rate  (2) Their advantages become more prominent when the underlying model is nonlinear. For one genetic architecture, the power is 5.9% (logistic), 72.6% (Chi-squared), 95.5% (MDR) and 88% (W-test). W-test behaves much better than MDR when 1%<MAF<5%, where the average power of W-test and MDR are 87.6 and 73.2%, respectively  Shi et al. [27]  RankCorr/P-value-based significance analysis (Sig); (1) Penalized least squares (PeLS), (2) Penalized Cox model (PeCOX)  n=200,P=500 1000,q=5, d=(3,4,20)b; Continuous predictors: follow multivariate normal distribution with two correlation structures and two correlations; (1) Continuous response under linear model, (2) Censored response under transformation model: g(t)=t3,t5,t and t1/3; Error distributions: N(0, 1), 0.7N(0,1)+0.3Cauchy(0,1) and t(1)  (1) When the error has a normal distribution, Sig has the best performance. However, when the error is contaminated or has a long tail, RankCorr can significantly outperform the alternatives (2) RankCorr has the largest AUC under all simulation scenarios, including those with linear transformation function and normal distribution error, followed by the nonrobust alternatives PeCOX and Sig  Li and Cui [28]  3G-SPA/Partial PCA (pPCA); Full PCA (fPCA)  n = 500, 1000, two separate genes with 10 SNPs each; discrete predictors: generated via MS program; Continuous response under linear model: (1) without genetic effects, (2) with main effects but no interaction, (3) with interaction effect smaller than main effects, (4) with interaction effect larger than main effects  3G-SPA has the highest power with the type I error reasonably controlled. The PCA-based methods are more conservative. The power of 3G-SPA is not depending on whether the genetic effect is because of its main effect or interaction, whereas the PCA-based methods prefer the scenarios with larger interactions  Kim and Park [29]  L-estimator GMDR; M-estimator GMDR/Ordinary GMDR  (1) n=1000,P=100, (2) n=3000,P=1000, d=(2,1); Discrete predictors: follow multinomial distribution with two different probabilities; Continuous response under linear model; Error distributions: (1−τ)N(0,1)+τN(0,400) and (1−τ)N(0,1)+τGamma(α,β) with τ=0,0.033 and 0.067  When there are no outliers (τ = 0), the ordinary GMDR and M-estimator GMDR outperform L- estimator GMDR. For all the designs with outliers, the performance of L-estimator and M-estimator GMDR are better than that of ordinary GMDR, with higher power. In some designs, L-estimator GMDR outperforms M-estimator GMDR  Yoshida and Koike [31]  SNPInterForest/SNPHarvester; BOOST  (1) n=4000,P=1000; Discrete predictors: the genetic model with different values of LD and MAFs; Binary response: the additive, multiplicative, and threshold models of interactions with weak marginal effects (2) n=400,P=1000; Discrete predictors: same as (1); Binary response: the pure interaction model without marginal effects  (1) The performance of SNPInterForest is superior to two alternatives in terms of their power SNPInterForest and SNPHarvester do not identify any spurious associations, whereas BOOST falsely identifies five interactions (2) The SNPInterForest and BOOST outperform SNPHarvester, with higher power and better identification accuracy The false discovery rates of SNPInterForest and SNPHarvester are much lower than that of BOOST  Beam et al. [32]  (1) BNN/BEAM; χ2 test (2)BNN/BEAM; χ2 test; GBM; MDR  n=2000,P=1000; Discrete predictors: the genetic model with different values of MAF; (1) Binary response: three biallelic models with main effects, (2) Binary response: pure interaction models without main effects  (1) BNN is uniformly more powerful than the χ2 test under all the simulated genetic models, and it outperforms BEAM in most instances (2) BNN performs better than BEAM, χ2 test and GBM under most of the simulations, and has similar performance with MDR. But MDR is much more time-consuming  Li et al. [33]  Nonpara-Scr/Screen and Clean (SC); SIS+LASSO; Forward LASSO; MDR  n=1000, P=3512, d=(12,8); Discrete predictors: dichotomizing the continuous variables generated from normal distribution with different correlations; Continuous response: a genetic model with both additive effects and dominant effects, considering three different noise levels  Nonpara-Scrbehaves better than the other four alternatives, with higher power and much smaller false-positive rates. Their average power are 97% (Nonpara-Scr), 90% (SIS+LASSO), 80% (SC), 48% (Forward LASSO) and 64% (MDR)  Li and Liu [34]  SODA/ Backward procedure (Back); Forward–backward method (ForBack); hierNet; IIS-SQDA  Ten different sample sizes with n=100, 100×101/9,…,1000; (1) P=50,d=(1,4); Continuous predictors: generated from normal, non-normal and heteroskedastic distributions, (2) P=1000,d=(1,4); Continuous predictors: non-normal distribution (3) P=50,d=(0,4), Continuous predictors: same as (2), (4) P=50,d=(2,4), Continuous predictors: same as (2); Binary response: Bayes classification rule  With Gaussian predictors, SODA, Back and ForBack can detect all important predictors as n increases, with small false positives, whereas IIS-SQDA and hierNet select too many false positives, resulting in high test error rates. Under other predictor settings, SODA performs much better than the nonrobust alternatives and achieves near-oracle classification accuracy when sample size is large enough  Wu et al. [35]  RankReg/LSregression-based method; Robust marginal analysis (RobMA)  n=200, P=200, q=2, d=(0,0,20); (1) Discrete predictors: MS program and Hardy–Weinberg equilibrium-based method, (2) Continuous predictors: follow multivariate normal distribution with an auto-regressive correlation structure; Continuous response under the linear model; Error distributions: N(0,1),t(3),0.85N(0,1)+0.15 log Normal(0,1) and 0.85N(0,1)+0.15Cauchy(0,1)  With N(0, 1) error distribution, RankReg is slightly inferior to the unrobust alternative LS, but still outperforms RobMA. With other three error types, RankReg has significant advantages, with more true positives, much fewer false positives, and higher estimation accuracy  Wu et al. [36]  LAD-hier/LAD-nhier; LAD-group; LS-hier; LS-nhier; LS-group  n=200, P=200, q=10, d=(10, 10, 15); (1) Continuous predictors: follow multivariate normal distribution with the banded and autoregressive correlation structure, (2) Discrete predictors: dichotomizing the above continuous values at the first and third quartiles; Censored response under the AFT model; Error distributions: N(0,1),0.80N(0,1)+0.20 log Normal(0,1),t(2) and 0.80N(0,1)+0.20Cauchy(0,1)  Without data contamination, the robust LAD-hier and the nonrobust counterpart LS-hier have comparable performances. With heavy-tailed errors, LAD-hier is significantly superior to the alternatives, in terms of both identification and prediction performance  a d=(b,c): b and c represent the number of true non-zero G effects and G-G interactions. b d=(a,b,c): a, b and c represent the number of true non-zero E effects, G effects and G-E interactions. Table 2. Summary of comparisons of the robust genetic analysis approaches against alternatives (partial list) References  Proposed/alternative  Main settings  Major conclusions  Wang et al. [26]  W-test/logistic regression; chi-squared test; MDR  n=1000,P=50; Discrete predictors: six different genetic architectures with minor allele frequency (MAF) in the common range and in the low frequency range, and LD in the high range, medium range and low range; binary response: (1) d=(4,2)a, linear interaction model with main effect, (2) d=(0,2), nonlinear interaction model without main effect  (1) The model-free W-test and MDR outperform the logistic regression even under the linear model, with higher power and lower type I error rate  (2) Their advantages become more prominent when the underlying model is nonlinear. For one genetic architecture, the power is 5.9% (logistic), 72.6% (Chi-squared), 95.5% (MDR) and 88% (W-test). W-test behaves much better than MDR when 1%<MAF<5%, where the average power of W-test and MDR are 87.6 and 73.2%, respectively  Shi et al. [27]  RankCorr/P-value-based significance analysis (Sig); (1) Penalized least squares (PeLS), (2) Penalized Cox model (PeCOX)  n=200,P=500 1000,q=5, d=(3,4,20)b; Continuous predictors: follow multivariate normal distribution with two correlation structures and two correlations; (1) Continuous response under linear model, (2) Censored response under transformation model: g(t)=t3,t5,t and t1/3; Error distributions: N(0, 1), 0.7N(0,1)+0.3Cauchy(0,1) and t(1)  (1) When the error has a normal distribution, Sig has the best performance. However, when the error is contaminated or has a long tail, RankCorr can significantly outperform the alternatives (2) RankCorr has the largest AUC under all simulation scenarios, including those with linear transformation function and normal distribution error, followed by the nonrobust alternatives PeCOX and Sig  Li and Cui [28]  3G-SPA/Partial PCA (pPCA); Full PCA (fPCA)  n = 500, 1000, two separate genes with 10 SNPs each; discrete predictors: generated via MS program; Continuous response under linear model: (1) without genetic effects, (2) with main effects but no interaction, (3) with interaction effect smaller than main effects, (4) with interaction effect larger than main effects  3G-SPA has the highest power with the type I error reasonably controlled. The PCA-based methods are more conservative. The power of 3G-SPA is not depending on whether the genetic effect is because of its main effect or interaction, whereas the PCA-based methods prefer the scenarios with larger interactions  Kim and Park [29]  L-estimator GMDR; M-estimator GMDR/Ordinary GMDR  (1) n=1000,P=100, (2) n=3000,P=1000, d=(2,1); Discrete predictors: follow multinomial distribution with two different probabilities; Continuous response under linear model; Error distributions: (1−τ)N(0,1)+τN(0,400) and (1−τ)N(0,1)+τGamma(α,β) with τ=0,0.033 and 0.067  When there are no outliers (τ = 0), the ordinary GMDR and M-estimator GMDR outperform L- estimator GMDR. For all the designs with outliers, the performance of L-estimator and M-estimator GMDR are better than that of ordinary GMDR, with higher power. In some designs, L-estimator GMDR outperforms M-estimator GMDR  Yoshida and Koike [31]  SNPInterForest/SNPHarvester; BOOST  (1) n=4000,P=1000; Discrete predictors: the genetic model with different values of LD and MAFs; Binary response: the additive, multiplicative, and threshold models of interactions with weak marginal effects (2) n=400,P=1000; Discrete predictors: same as (1); Binary response: the pure interaction model without marginal effects  (1) The performance of SNPInterForest is superior to two alternatives in terms of their power SNPInterForest and SNPHarvester do not identify any spurious associations, whereas BOOST falsely identifies five interactions (2) The SNPInterForest and BOOST outperform SNPHarvester, with higher power and better identification accuracy The false discovery rates of SNPInterForest and SNPHarvester are much lower than that of BOOST  Beam et al. [32]  (1) BNN/BEAM; χ2 test (2)BNN/BEAM; χ2 test; GBM; MDR  n=2000,P=1000; Discrete predictors: the genetic model with different values of MAF; (1) Binary response: three biallelic models with main effects, (2) Binary response: pure interaction models without main effects  (1) BNN is uniformly more powerful than the χ2 test under all the simulated genetic models, and it outperforms BEAM in most instances (2) BNN performs better than BEAM, χ2 test and GBM under most of the simulations, and has similar performance with MDR. But MDR is much more time-consuming  Li et al. [33]  Nonpara-Scr/Screen and Clean (SC); SIS+LASSO; Forward LASSO; MDR  n=1000, P=3512, d=(12,8); Discrete predictors: dichotomizing the continuous variables generated from normal distribution with different correlations; Continuous response: a genetic model with both additive effects and dominant effects, considering three different noise levels  Nonpara-Scrbehaves better than the other four alternatives, with higher power and much smaller false-positive rates. Their average power are 97% (Nonpara-Scr), 90% (SIS+LASSO), 80% (SC), 48% (Forward LASSO) and 64% (MDR)  Li and Liu [34]  SODA/ Backward procedure (Back); Forward–backward method (ForBack); hierNet; IIS-SQDA  Ten different sample sizes with n=100, 100×101/9,…,1000; (1) P=50,d=(1,4); Continuous predictors: generated from normal, non-normal and heteroskedastic distributions, (2) P=1000,d=(1,4); Continuous predictors: non-normal distribution (3) P=50,d=(0,4), Continuous predictors: same as (2), (4) P=50,d=(2,4), Continuous predictors: same as (2); Binary response: Bayes classification rule  With Gaussian predictors, SODA, Back and ForBack can detect all important predictors as n increases, with small false positives, whereas IIS-SQDA and hierNet select too many false positives, resulting in high test error rates. Under other predictor settings, SODA performs much better than the nonrobust alternatives and achieves near-oracle classification accuracy when sample size is large enough  Wu et al. [35]  RankReg/LSregression-based method; Robust marginal analysis (RobMA)  n=200, P=200, q=2, d=(0,0,20); (1) Discrete predictors: MS program and Hardy–Weinberg equilibrium-based method, (2) Continuous predictors: follow multivariate normal distribution with an auto-regressive correlation structure; Continuous response under the linear model; Error distributions: N(0,1),t(3),0.85N(0,1)+0.15 log Normal(0,1) and 0.85N(0,1)+0.15Cauchy(0,1)  With N(0, 1) error distribution, RankReg is slightly inferior to the unrobust alternative LS, but still outperforms RobMA. With other three error types, RankReg has significant advantages, with more true positives, much fewer false positives, and higher estimation accuracy  Wu et al. [36]  LAD-hier/LAD-nhier; LAD-group; LS-hier; LS-nhier; LS-group  n=200, P=200, q=10, d=(10, 10, 15); (1) Continuous predictors: follow multivariate normal distribution with the banded and autoregressive correlation structure, (2) Discrete predictors: dichotomizing the above continuous values at the first and third quartiles; Censored response under the AFT model; Error distributions: N(0,1),0.80N(0,1)+0.20 log Normal(0,1),t(2) and 0.80N(0,1)+0.20Cauchy(0,1)  Without data contamination, the robust LAD-hier and the nonrobust counterpart LS-hier have comparable performances. With heavy-tailed errors, LAD-hier is significantly superior to the alternatives, in terms of both identification and prediction performance  References  Proposed/alternative  Main settings  Major conclusions  Wang et al. [26]  W-test/logistic regression; chi-squared test; MDR  n=1000,P=50; Discrete predictors: six different genetic architectures with minor allele frequency (MAF) in the common range and in the low frequency range, and LD in the high range, medium range and low range; binary response: (1) d=(4,2)a, linear interaction model with main effect, (2) d=(0,2), nonlinear interaction model without main effect  (1) The model-free W-test and MDR outperform the logistic regression even under the linear model, with higher power and lower type I error rate  (2) Their advantages become more prominent when the underlying model is nonlinear. For one genetic architecture, the power is 5.9% (logistic), 72.6% (Chi-squared), 95.5% (MDR) and 88% (W-test). W-test behaves much better than MDR when 1%<MAF<5%, where the average power of W-test and MDR are 87.6 and 73.2%, respectively  Shi et al. [27]  RankCorr/P-value-based significance analysis (Sig); (1) Penalized least squares (PeLS), (2) Penalized Cox model (PeCOX)  n=200,P=500 1000,q=5, d=(3,4,20)b; Continuous predictors: follow multivariate normal distribution with two correlation structures and two correlations; (1) Continuous response under linear model, (2) Censored response under transformation model: g(t)=t3,t5,t and t1/3; Error distributions: N(0, 1), 0.7N(0,1)+0.3Cauchy(0,1) and t(1)  (1) When the error has a normal distribution, Sig has the best performance. However, when the error is contaminated or has a long tail, RankCorr can significantly outperform the alternatives (2) RankCorr has the largest AUC under all simulation scenarios, including those with linear transformation function and normal distribution error, followed by the nonrobust alternatives PeCOX and Sig  Li and Cui [28]  3G-SPA/Partial PCA (pPCA); Full PCA (fPCA)  n = 500, 1000, two separate genes with 10 SNPs each; discrete predictors: generated via MS program; Continuous response under linear model: (1) without genetic effects, (2) with main effects but no interaction, (3) with interaction effect smaller than main effects, (4) with interaction effect larger than main effects  3G-SPA has the highest power with the type I error reasonably controlled. The PCA-based methods are more conservative. The power of 3G-SPA is not depending on whether the genetic effect is because of its main effect or interaction, whereas the PCA-based methods prefer the scenarios with larger interactions  Kim and Park [29]  L-estimator GMDR; M-estimator GMDR/Ordinary GMDR  (1) n=1000,P=100, (2) n=3000,P=1000, d=(2,1); Discrete predictors: follow multinomial distribution with two different probabilities; Continuous response under linear model; Error distributions: (1−τ)N(0,1)+τN(0,400) and (1−τ)N(0,1)+τGamma(α,β) with τ=0,0.033 and 0.067  When there are no outliers (τ = 0), the ordinary GMDR and M-estimator GMDR outperform L- estimator GMDR. For all the designs with outliers, the performance of L-estimator and M-estimator GMDR are better than that of ordinary GMDR, with higher power. In some designs, L-estimator GMDR outperforms M-estimator GMDR  Yoshida and Koike [31]  SNPInterForest/SNPHarvester; BOOST  (1) n=4000,P=1000; Discrete predictors: the genetic model with different values of LD and MAFs; Binary response: the additive, multiplicative, and threshold models of interactions with weak marginal effects (2) n=400,P=1000; Discrete predictors: same as (1); Binary response: the pure interaction model without marginal effects  (1) The performance of SNPInterForest is superior to two alternatives in terms of their power SNPInterForest and SNPHarvester do not identify any spurious associations, whereas BOOST falsely identifies five interactions (2) The SNPInterForest and BOOST outperform SNPHarvester, with higher power and better identification accuracy The false discovery rates of SNPInterForest and SNPHarvester are much lower than that of BOOST  Beam et al. [32]  (1) BNN/BEAM; χ2 test (2)BNN/BEAM; χ2 test; GBM; MDR  n=2000,P=1000; Discrete predictors: the genetic model with different values of MAF; (1) Binary response: three biallelic models with main effects, (2) Binary response: pure interaction models without main effects  (1) BNN is uniformly more powerful than the χ2 test under all the simulated genetic models, and it outperforms BEAM in most instances (2) BNN performs better than BEAM, χ2 test and GBM under most of the simulations, and has similar performance with MDR. But MDR is much more time-consuming  Li et al. [33]  Nonpara-Scr/Screen and Clean (SC); SIS+LASSO; Forward LASSO; MDR  n=1000, P=3512, d=(12,8); Discrete predictors: dichotomizing the continuous variables generated from normal distribution with different correlations; Continuous response: a genetic model with both additive effects and dominant effects, considering three different noise levels  Nonpara-Scrbehaves better than the other four alternatives, with higher power and much smaller false-positive rates. Their average power are 97% (Nonpara-Scr), 90% (SIS+LASSO), 80% (SC), 48% (Forward LASSO) and 64% (MDR)  Li and Liu [34]  SODA/ Backward procedure (Back); Forward–backward method (ForBack); hierNet; IIS-SQDA  Ten different sample sizes with n=100, 100×101/9,…,1000; (1) P=50,d=(1,4); Continuous predictors: generated from normal, non-normal and heteroskedastic distributions, (2) P=1000,d=(1,4); Continuous predictors: non-normal distribution (3) P=50,d=(0,4), Continuous predictors: same as (2), (4) P=50,d=(2,4), Continuous predictors: same as (2); Binary response: Bayes classification rule  With Gaussian predictors, SODA, Back and ForBack can detect all important predictors as n increases, with small false positives, whereas IIS-SQDA and hierNet select too many false positives, resulting in high test error rates. Under other predictor settings, SODA performs much better than the nonrobust alternatives and achieves near-oracle classification accuracy when sample size is large enough  Wu et al. [35]  RankReg/LSregression-based method; Robust marginal analysis (RobMA)  n=200, P=200, q=2, d=(0,0,20); (1) Discrete predictors: MS program and Hardy–Weinberg equilibrium-based method, (2) Continuous predictors: follow multivariate normal distribution with an auto-regressive correlation structure; Continuous response under the linear model; Error distributions: N(0,1),t(3),0.85N(0,1)+0.15 log Normal(0,1) and 0.85N(0,1)+0.15Cauchy(0,1)  With N(0, 1) error distribution, RankReg is slightly inferior to the unrobust alternative LS, but still outperforms RobMA. With other three error types, RankReg has significant advantages, with more true positives, much fewer false positives, and higher estimation accuracy  Wu et al. [36]  LAD-hier/LAD-nhier; LAD-group; LS-hier; LS-nhier; LS-group  n=200, P=200, q=10, d=(10, 10, 15); (1) Continuous predictors: follow multivariate normal distribution with the banded and autoregressive correlation structure, (2) Discrete predictors: dichotomizing the above continuous values at the first and third quartiles; Censored response under the AFT model; Error distributions: N(0,1),0.80N(0,1)+0.20 log Normal(0,1),t(2) and 0.80N(0,1)+0.20Cauchy(0,1)  Without data contamination, the robust LAD-hier and the nonrobust counterpart LS-hier have comparable performances. With heavy-tailed errors, LAD-hier is significantly superior to the alternatives, in terms of both identification and prediction performance  a d=(b,c): b and c represent the number of true non-zero G effects and G-G interactions. b d=(a,b,c): a, b and c represent the number of true non-zero E effects, G effects and G-E interactions. Remarks Our literature review suggests that there have been significant developments in robust genetic interaction analysis methods, especially in the recent literature. However, compared with low-dimensional studies and high-dimensional studies on main effects, research tailored to genetic interactions is still limited. Based on strategies developed in published studies, much robust genetic interaction research can be conducted. For example, the ‘robust loss function + penalization’ joint analysis can be extended to other robust loss functions, such as the Huber’s loss [66] and Tukey’s bi-square loss [67]. Some other robust estimations such as S-estimation and MM-estimation, which have been shown to be powerful in low-dimensional settings, can be potentially adapted to detect genetic interactions. The high-dimensional robust screening techniques may be potentially coupled with the nonparametric feature screening [68], rank correlation screening [69] and others. Most of the existing methods, including those reviewed above, are limited to one aspect of robustness. They often make no or weaker assumptions on one of the following aspects: specific model form, conditional distribution of the response (distribution of the error) or distribution of the predictors, resulting in different types of robustness. In practical data analysis, for example, model mis-specification and outliers may happen together and demand robustness in multiple aspects. A few studies that can accommodate both model mis-specification and outliers in response include the rank correlation-based method [27] and R-GMDR [29]. For high-dimensional main effect analysis, some methods have been proposed for variable selection and can accommodate outliers/contaminations in both responses and predictors, such as the sparse least trimmed squares regression [48], the robust L2 boosting method [70] and others. Similar development is also needed for interaction analysis. Our review also suggests that for some aspects of robustness, there are multiple methods available. However, there is still a lack of direct and comprehensive comparisons. In addition, most of the existing methods focus on methodological development. Theoretical studies are limited. It has been noted that theoretical investigation with high-dimensional data is challenging. Significantly more challenges can be brought by interactions. The computational algorithms are critically dependent on the methodology frameworks. However, some similar techniques are involved. For example, to address the problem that the robust objective functions are often not differentiable, the smoothed approximations are adopted. In addition, the (group) coordinate descent algorithms are usually developed for penalization approaches. The computational cost differs widely across these approaches because of their methodologies, adopted algorithms, sizes of the analyzed data sets and others. The robust methods are observed to be in general more computationally intensive than nonrobust methods, as they often have more complex objective functions and optimization algorithms. However, the results suggest that these robust interaction learning methods are still computational affordable. For some of the reviewed methods, software is publicly available. However, the development of user-friendly software is still badly needed. Applications The methods reviewed above and other robust methods in the literature have been applied to the analysis of multiple diseases and made sensible discoveries different from those detected by nonrobust methods. A few representative findings are briefly reviewed below. A sporadic breast cancer case-control data set is analyzed using MDR [24]. This case study aims at detecting SNP-SNP interactions associated with sporadic breast cancer using SNPs from selected genes. The data set contains data on 200 White women with sporadic primary invasive breast cancer treated at the Vanderbilt University Medical Center during 1982–96. Ten SNPs in five selected genes (COMT, CYP1A1, CYP1B1, GSTM1 and GSTT1) are analyzed. A four-locus interaction among four SNPs (COMT, CYP1A1m1, CYP1B1 codon 48 and CYP1B1 codon 432) is detected to be significantly associated with the risk of sporadic breast cancer. The protein products of the corresponding genes of these four SNPs have been shown to interact as enzymes in the metabolism of estrogens in breast tissues. Estrogens have been recognized as one of the principal factors related to breast cancer. Thus, it is expected that the interactions between the SNPs in these four genes have independent contributions to the risk of breast cancer. However, using the parametric logistic regression analysis, previous studies fail to find any association between these genes or interactions and an increased risk of breast cancer. Similar promising findings have also been made in other MDR studies. The nonparametric W-test proposed by [26] is applied to two bipolar disorder genome-wide association data sets. The first data set is from the Wellcome Trust Case Control Consortium (WTCCC), including 2000 bipolar cases and 3000 controls of European Caucasians with 414 682 SNPs. The second data set is collected by the Genetic Association Information Network (GAIN) project and consists of 653 cases and 1767 controls with 729 304 SNPs. For these two real GWAS data sets, a three-step procedure is used to search for SNP-SNP interactions and their corresponding G-G interactions. First, the main effects of the whole SNP set are evaluated using W-test and then a candidate SNP set is identified based on the P-value. Second, the W-test for detecting important SNP-SNP interactions is conducted on the candidate SNP set. Third, the G-G interactions corresponding to the identified SNP-SNP interactions are obtained using Web-based databases. It is interesting that the W-test approach identifies a number of genes that can be replicated in the two independent GWAS data sets, and that the interactions among 18 identified genes naturally form two networks. The identification results are different from the previous GWAS studies. Most of the findings play important roles in brain and neuronal functions and are highly relevant to neuron disorders, including RTN4R, DPP10, CMSD1, SLIT3, PARK2, TMEM132D, HNT, CNTNAP2, ACCN1, A2BP1 and others. For instance, RTN4R has been shown to encode a Nogo receptor that mediates axonal growth inhibition and plays a role in regulating axonal regeneration and plasticity in the central nervous system. Kim and Park [29] apply the R-GMDR approach to the Korea Association Resource (KARE) data, which is collected by the large-scale genome-wide association analyses in the Gyeonggi Province of South Korea. The outcome of interest is the (log-transformed) Homeostasis Model Assessment of Insulin Resistance (HOMA-IR) level, which is continuous and widely used to estimate insulin resistance. After removing samples with the missing BMI or HOMA-IR score, this study includes of 8577 individuals and 327 872 SNPs. To reduce the computational burden on the R-GMDR, a pre-screening procedure based on the P-value of a marginal SNP analysis and the linkage disequilibrium (LD) coefficient is conducted, resulting in 585 SNPs for the downstream analysis. Then, the R-GMDR is applied to identify important SNP-SNP interactions between these 585 SNPs. The identified SNPs are matched to their corresponding genes to facilitate the analysis of possible biological implications. The two R-GMDR methods identify sensible genes, especially hub genes, that are missed by the original GMDR. Specifically, the KCNH1 gene is identified to have strong interaction effects with other genes on the function of insulin secretion. The other hub gene RYR2 has been demonstrated to be related to diabetes. Outliers are observed in the tail part of the response, leading to an adverse influence on the original GMDR. Wu et al. [35] conduct a G-E interaction analysis for lung cancer prognosis. Data from four independent studies, conducted by the Dana-Farber Cancer Institute, HLM (Moffitt Cancer Center), UM (University of Michigan Cancer Center) and MSKCC, are pooled and analyzed. The combined sample size is 351. Different from the former three applications, which focus on SNP-SNP interactions, the continuous gene expression measurements on 22 283 probes are analyzed as G factors. Two clinical factors, age and smoking status, are analyzed as E factors. The response variable of interest is right censored overall survival. The robust method identifies 11 genes as having significant associations with prognosis, among which four and five genes are detected to having interactions with age and smoking status, respectively. The nonrobust approach based on the ordinary least squared loss function identifies 19 genes with main effects (6 overlap with the robust method) and three interactions with smoking status. A brief literature search is conducted and suggests that most of the genes identified by the robust approach have important implications for lung cancer prognosis. Findings missed by the nonrobust approach include genes CYP2B6, WIF1, HLA-DQA1 and others. Take gene WIF1 as an example. This gene has been implicated in lung cancer pathogenesis, and the methylation silencing of WIF1 is a usual and potentially central mechanism of the aberrant activation of the Wnt signaling pathway. Remarks In the aforementioned studies, robust methods have led to promising new discoveries. However, as robust methods are relatively new and limited, their practical impact may still be limited. In particular, their applications in large-scale cutting-edge studies and biomedical publications are still rare. Significant effort in applications is needed. Discussion Extensive research on genetic interactions has been conducted, and tailored analysis methods have been developed in a large number of, especially recent, studies. In this review, we have focused on robust methods, which have important practical implications but have not received enough attention in the literature. To the best of our knowledge, this article is the first to comprehensively review robust genetic interaction analysis methods. With our limited knowledge and literature search, it is inevitable the selection of reviewed methods is biased and incomplete. With the limited scope of this review, it is not our goal to provide the ‘full picture’. The review has been conducted in a rather ‘statistical’ manner because of the high complexity of the methods. In addition, we have tried to review the generically applicable (as opposed to case specific) methods. Even with these limitations, being the first of its kind, this review can still be much useful for researchers in the field of genetic interaction analysis. Funding This work was partly supported by the National Institutes of Health (grant numbers CA216017, CA191383, CA204120); and the National Natural Science Foundation of China (grant numbers 61402276 and 91546202). Mengyun Wu is an Assistant Professor in School of Statistics and Management, Shanghai University of Finance and Economics. She is also a Postdoctoral Fellow in the Department of Biostatistics at Yale University. Shuangge Ma is a Professor in the Department of Biostatistics at Yale University. References 1 Yoshimaru T, Komatsu M, Matsuo T, et al.   Targeting BIG3-PHB2 interaction to overcome tamoxifen resistance in breast cancer cells. Nat Commun  2013; 4: 2443. Google Scholar CrossRef Search ADS PubMed  2 Zhou W, Liu G, Miller DP, et al.   Gene-environment interaction for the ERCC2 polymorphisms and cumulative cigarette smoking exposure in lung cancer. Cancer Res  2002; 62: 1377– 81. Google Scholar PubMed  3 Cordell HJ. Detecting gene-gene interactions that underlie human diseases. Nat Rev Genet  2009; 10( 6): 392– 404. Google Scholar CrossRef Search ADS PubMed  4 Van Steen K. Travelling the world of gene-gene interactions. Brief Bioinform  2012; 13( 1): 1– 19. Google Scholar CrossRef Search ADS PubMed  5 Thomas D. Gene-environment-wide association studies: emerging approaches. Nat Rev Genet  2010; 11( 4): 259– 72. Google Scholar CrossRef Search ADS PubMed  6 Simonds NI, Ghazarian AA, Pimentel CB, et al.   Review of the gene-environment interaction literature in cancer: what do we know? Genet Epidemiol  2016; 40( 5): 356– 65. Google Scholar CrossRef Search ADS PubMed  7 Purcell S, Neale B, Todd-Brown K, et al.   PLINK: a tool set for whole-genome association and population-based linkage analyses. Am J Hum Genet  2007; 81( 3): 559– 75. Google Scholar CrossRef Search ADS PubMed  8 Zhao J, Zhu Y, Xiong M. Genome-wide gene-gene interaction analysis for next-generation sequencing. Eur J Hum Genet  2016; 24( 3): 421– 8. Google Scholar CrossRef Search ADS PubMed  9 Bild AH, Yao G, Chang JT, et al.   Oncogenic pathway signatures in human cancers as a guide to targeted therapies. Nature  2006; 439( 7074): 353– 7. Google Scholar CrossRef Search ADS PubMed  10 Weissbrod O, Lippert C, Geiger D, et al.   Accurate liability estimation improves power in ascertained case-control studies. Nat Methods  2015; 12( 4): 332. Google Scholar CrossRef Search ADS PubMed  11 Stark A, Stahl MS, Kirchner HL, et al.   Body mass index at the time of diagnosis and the risk of advanced stages and poorly differentiated cancers of the breast: findings from a case-series study. Int J Obes  2010; 34( 9): 1381– 6. Google Scholar CrossRef Search ADS   12 Stephan J, Stegle O, Beyer A. A random forest approach to capture genetic effects in the presence of population structure. Nat Commun  2015; 6( 1): 7432. Google Scholar CrossRef Search ADS PubMed  13 Osborne JW, Overbay A. The power of outliers (and why researchers should always check for them). Pract Assess Res Eval  2010; 9: 1– 12. 14 Li C, Wong WH. Model-based analysis of oligonucleotide arrays: expression index computation and outlier detection. Proc Natl Acad Sci USA  2001; 98( 1): 31– 6. Google Scholar CrossRef Search ADS PubMed  15 Shieh AD, Hung YS. Detecting outlier samples in microarray data. Stat Appl Genet Mol  2009; 8( 1): 1– 24. Google Scholar CrossRef Search ADS   16 Rampatige R, Gamage S, Peiris S, et al.   Assessing the reliability of causes of death reported by the vital registration system in Sri Lanka: medical records review in Colombo. Health Inf Manag  2013; 42( 3): 20– 8. Google Scholar PubMed  17 Cancer Genome Atlas Research Network. Comprehensive molecular profiling of lung adenocarcinoma. Nature  2014; 511: 543– 50. CrossRef Search ADS PubMed  18 Wu M, Zang Y, Zhang S, et al.   Accommodating missingness in environmental measurements in gene-environment interaction analysis. Genet Epidemiol  2017; 41( 6): 523– 54. Google Scholar CrossRef Search ADS PubMed  19 Chai H, Zhang Q, Jiang Y, et al.   Identifying gene-environment interactions for prognosis using a robust approach. Econ Stat  2017; 4: 105– 20. 20 Burgess DJ. Cancer genetics: initially complex, always heterogeneous. Nat Rev Cancer  2011; 11( 3): 153. Google Scholar CrossRef Search ADS PubMed  21 Haibe-Kains B, Desmedt C, Loi S, et al.   A three-gene model to robustly identify breast cancer molecular subtypes. J Natl Cancer Inst  2012; 104( 4): 311– 25. Google Scholar CrossRef Search ADS PubMed  22 Huber PJ, Ronchetti EM. Robust Statistics. Wiley Series in Probability and Statistics , 2nd ed. Hoboken, NJ: Wiley, 2009. 23 Wu C, Ma S. A selective review of robust variable selection with applications in bioinformatics. Brief Bioinform  2015; 16( 5): 873– 83. Google Scholar CrossRef Search ADS PubMed  24 Ritchie MD, Hahn LW, Roodi N, et al.   Multifactor-dimensionality reduction reveals high-order interactions among estrogen-metabolism genes in sporadic breast cancer. Am J Hum Genet  2001; 69( 1): 138– 47. Google Scholar CrossRef Search ADS PubMed  25 Winham SJ, Motsinger-Reif AA. An R package implementation of multifactor dimensionality reduction. Biodata Min  2011; 4( 1): 24. Google Scholar CrossRef Search ADS PubMed  26 Wang MH, Sun R, Guo J, et al.   A fast and powerful W-test for pairwise epistasis testing. Nucleic Acids Res  2016; 44( 12): e115. Google Scholar CrossRef Search ADS PubMed  27 Shi X, Liu J, Huang J, et al.   A penalized robust method for identifying gene-environment interactions. Genet Epidemiol  2014; 38( 3): 220– 30. Google Scholar CrossRef Search ADS PubMed  28 Li S, Cui Y. Gene-centric gene-gene interaction: a model-based kernel machine method. Ann Appl Stat  2012; 6( 3): 1134– 61. Google Scholar CrossRef Search ADS   29 Kim Y, Park T. Robust gene-gene interaction analysis in genome wide association studies. PLoS One  2015; 10( 8): e0135016. Google Scholar CrossRef Search ADS PubMed  30 Wang G, Zhao Y, Zhang Q, et al.   Identifying gene-environment interactions associated with prognosis using penalized quantile regression. In: Big and Complex Data Analysis . Springer International Publishing, 2017, 347– 67. Google Scholar CrossRef Search ADS   31 Yoshida M, Koike A. SNPInterForest: a new method for detecting epistatic interactions. BMC Bioinform  2011; 12: 469. Google Scholar CrossRef Search ADS   32 Beam AL, Motsinger-Reif A, Doyle J. Bayesian neural networks for detecting epistasis in genetic association studies. BMC Bioinform  2014; 15( 1): 368. Google Scholar CrossRef Search ADS   33 Li J, Dan J, Li C, et al.   A model-free approach for detecting interactions in genetic association studies. Brief Bioinform  2014; 15( 6): 1057– 8. Google Scholar CrossRef Search ADS PubMed  34 Li Y, Liu JS. Robust variable and interaction selection for logistic regression and general index models. J Am Stat Assoc  2017, in press. 35 Wu C, Shi X, Cui Y, et al.   A penalized robust semiparametric approach for gene-environment interactions. Stat Med  2015; 34( 30): 4016– 30. Google Scholar CrossRef Search ADS PubMed  36 Wu C, Jiang Y, Ren J, et al.   Dissecting gene-environment interactions: a penalized robust approach accounting for hierarchical structures. Stat Med  2018; 37( 3): 437– 56. Google Scholar CrossRef Search ADS PubMed  37 Heagerty PJ, Kurland BF. Misspecified maximum likelihood estimates and generalised linear mixed models. Biometrika  2001; 88( 4): 973– 85. Google Scholar CrossRef Search ADS   38 Lou XY, Chen GB, Yan L, et al.   A generalized combinatorial approach for detecting gene-by-gene and gene-by-environment interactions with application to nicotine dependence. Am J Hum Genet  2007; 80( 6): 1125– 37. Google Scholar CrossRef Search ADS PubMed  39 Li MD, Burmeister M. New insights into the genetics of addiction. Nat Rev Cancer  2009; 10( 4): 225. Google Scholar CrossRef Search ADS   40 Gui J, Moore JH, Kelsey KT, et al.   A novel survival multifactor dimensionality reduction method for detecting gene-gene interactions with application to bladder cancer prognosis. Hum Genet  2011; 129( 1): 101– 10. Google Scholar CrossRef Search ADS PubMed  41 Gao X, Alvo M. A unified nonparametric approach for unbalanced factorial designs. J Am Stat Assoc  2005; 100( 471): 926– 41. Google Scholar CrossRef Search ADS   42 Wu X, Jin L, Xiong M. Mutual information for testing gene-environment interaction. PLoS One  2009; 4( 2): e4578. Google Scholar CrossRef Search ADS PubMed  43 Li M, Ye C, Fu W, et al.   Detecting genetic interactions for quantitative traits with U-statistics. Genet Epidemiol  2011; 35( 6): 457– 68. Google Scholar PubMed  44 Larson NB, Schaid DJ. A kernel regression approach to gene-gene interaction detection for case-control studies. Genet Epidemiol  2013; 37( 7): 695– 703. Google Scholar CrossRef Search ADS PubMed  45 Marceau R, Lu W, Holloway S, et al.   A fast multiple-kernel method with applications to detect gene-environment interaction. Genet Epidemiol  2015; 39( 6): 456– 68. Google Scholar CrossRef Search ADS PubMed  46 Lin HY, Wang W, Liu YH, et al.   Comparison of multivariate adaptive regression splines and logistic regression in detecting SNP-SNP interactions and their application in prostate cancer. J Hum Genet  2008; 53( 9): 802– 11. Google Scholar CrossRef Search ADS PubMed  47 Luss R, Rosset S, Shahar M. Efficient regularized isotonic regression with application to gene-gene interaction search. Ann Appl Stat  2012; 6( 1): 253– 83. Google Scholar CrossRef Search ADS   48 Alfons A, Croux C, Gelper S. Sparse least trimmed squares regression for analyzing high-dimensional large data sets. Ann Appl Stat  2013; 7( 1): 226– 48. Google Scholar CrossRef Search ADS   49 Lee S, Kim Y, Kwon MS, et al.   A comparative study on multifactor dimensionality reduction methods for detecting gene-gene interactions with the survival phenotype. BioMed Res Int  2015; 2015: 671859. Google Scholar PubMed  50 Fan J, Xue L, Zou H. Multitask quantile regression under the transnormal model. J Am Stat Assoc  2016; 111( 516): 1726– 35. Google Scholar CrossRef Search ADS PubMed  51 Chaudhuri P, Loh WY. Nonparametric estimation of conditional quantiles using quantile regression trees. Bernoulli  2002; 8: 561– 76. 52 Zhu Q, Hu Y, Tian M. Identifying interaction effects via additive quantile regression models. Stat Its Interface  2017; 10( 2): 255– 65. Google Scholar CrossRef Search ADS   53 Bien J, Taylor J, Tibshirani R. A lasso for hierarchical interactions. Ann Stat  2013; 41( 3): 1111– 41. Google Scholar CrossRef Search ADS PubMed  54 Lim M, Hastie T. Learning interactions via hierarchical group-lasso regularization. J Comput Graph Stat  2015; 24( 3): 627– 54. Google Scholar CrossRef Search ADS PubMed  55 Zhu R, Zhao H, Ma S. Identifying gene-environment and gene-gene interactions using a progressive penalization approach. Genet Epidemiol  2014; 38( 4): 353– 68. Google Scholar CrossRef Search ADS PubMed  56 Wu M, Huang J, Ma S. Identifying gene-gene interactions using penalized tensor regression. Stat Med  2018; 37( 4): 598– 610. Google Scholar CrossRef Search ADS PubMed  57 Hao N, Zhang HH. A note on high-dimensional linear regression with interactions. Am Stat  2017; 71( 4): 291– 7. Google Scholar CrossRef Search ADS   58 Winham SJ, Colby CL, Freimuth RR, et al.   SNP interaction detection with random forests in high-dimensional genetic data. BMC Bioinform  2012; 13( 1): 164. Google Scholar CrossRef Search ADS   59 Li J, Malley JD, Andrew AS, et al.   Detecting gene-gene interactions using a permutation-based random forest method. BioData Min  2016; 9( 1): 14. Google Scholar CrossRef Search ADS PubMed  60 Ritchie MD, White BC, Parker JS, et al.   Optimizationof neural network architecture using genetic programming improvesdetection and modeling of gene-gene interactions in studies of human diseases. BMC Bioinform  2003; 4: 28. Google Scholar CrossRef Search ADS   61 Motsinger-Reif AA, Dudek SM, Hahn LW, et al.   Comparison of approaches for machine-learning optimization of neural networks for detecting gene-gene interactions in genetic epidemiology. Genet Epidemiol  2008; 32( 4): 325– 40. Google Scholar CrossRef Search ADS PubMed  62 Günther F, Wawro N, Bammann K. Neural networks for modeling gene-gene interactions in association studies. BMC Genet  2009; 10( 1): 87. Google Scholar CrossRef Search ADS PubMed  63 Ritchie MD, Motsinger AA, Bush WS, et al.   Genetic programming neural networks: a powerful bioinformatics tool for human genetics. Appl Soft Comput  2007; 7( 1): 471– 9. Google Scholar CrossRef Search ADS PubMed  64 Liu J, Huang J, Zhang Y, et al.   Identification of gene-environment interactions in cancer studies using penalization. Genomics  2013; 102( 4): 189– 94. Google Scholar CrossRef Search ADS PubMed  65 Maronna RA, Yohai VJ. High finite-sample efficiency and robustness based on distance-constrained maximum likelihood. Comput Stat Data Anal  2015; 83: 262– 74. Google Scholar CrossRef Search ADS   66 Huber PJ. Robust estimation of a location parameter. Ann Math Stat  1964; 35( 1): 73– 101. Google Scholar CrossRef Search ADS   67 Beaton A, Tukey J. The fitting of power series, meaning polynomials, illustrated on band-spectroscopic data. Technometrics  1974; 16( 2): 147– 85. Google Scholar CrossRef Search ADS   68 Lin L, Sun J, Zhu L. Nonparametric feature screening. Comput Stat Data Anal  2013; 67: 162– 74. Google Scholar CrossRef Search ADS   69 Li G, Peng H, Zhang J, et al.   Robust rank correlation based screening. Ann Stat  2012; 40( 3): 1846– 77. Google Scholar CrossRef Search ADS   70 Lutz RW, Kalisch M, Bühlmann P. Robustified L2 boosting. Comput Stat Data Anal  2008; 52( 7): 3331– 41. Google Scholar CrossRef Search ADS   © The Author(s) 2018. Published by Oxford University Press. All rights reserved. For permissions, please email: journals.permissions@oup.com This article is published and distributed under the terms of the Oxford University Press, Standard Journals Publication Model (https://academic.oup.com/journals/pages/about_us/legal/notices)

Journal

Briefings in BioinformaticsOxford University Press

Published: Apr 19, 2018

There are no references for this article.

You’re reading a free preview. Subscribe to read the entire article.


DeepDyve is your
personal research library

It’s your single place to instantly
discover and read the research
that matters to you.

Enjoy affordable access to
over 18 million articles from more than
15,000 peer-reviewed journals.

All for just $49/month

Explore the DeepDyve Library

Search

Query the DeepDyve database, plus search all of PubMed and Google Scholar seamlessly

Organize

Save any article or search result from DeepDyve, PubMed, and Google Scholar... all in one place.

Access

Get unlimited, online access to over 18 million full-text articles from more than 15,000 scientific journals.

Your journals are on DeepDyve

Read from thousands of the leading scholarly journals from SpringerNature, Elsevier, Wiley-Blackwell, Oxford University Press and more.

All the latest content is available, no embargo periods.

See the journals in your area

DeepDyve

Freelancer

DeepDyve

Pro

Price

FREE

$49/month
$360/year

Save searches from
Google Scholar,
PubMed

Create lists to
organize your research

Export lists, citations

Read DeepDyve articles

Abstract access only

Unlimited access to over
18 million full-text articles

Print

20 pages / month

PDF Discount

20% off