In this section, we first introduce the common and robust marginal modeling methods for G-E interactions. Consider a study with N independent subjects. For a subject i, suppose that there are q environmental/clinical variables ${𝒆}_i=(e_{i1},e_{i2},\mathrm{\dots },e_{iq})^{{}^{\prime }}{}_{1\times q}$, and p genes ${𝒙}_i=(x_{i1},x_{i2},\mathrm{\dots },x_{ip})^{{}^{\prime }}{}_{1\times p}$. Assume the survival outcome $T_i$ is related to the environmental/clinical variables ${𝒆}_i$, gene expression covariates ${𝒙}_i$, and their component-wise interactions

(1)$$\begin{array}{c}\hfill {𝒘}_{i..}={(w_{i11},\mathrm{\dots },w_{i1p},w_{i21},\mathrm{\dots },w_{iqp})}_{1\times qp}^{\prime }=\hfill \\ \hfill {(e_{i1}{x}_{i1},\mathrm{\dots },e_{i1}{x}_{ip},e_{i2}{x}_{i1},\mathrm{\dots },e_{iq}{x}_{ip})}_{1\times qp}^{\prime }.\hfill \end{array}$$

Since the survival time may be right-censored and incompletely observed. Define observe survival time $V=min(T,C)$, C is censoring time, and we use $\delta =I(T\le C)$ is the indicator of whether the survival time of subject is censored.

In a marginal Cox’s regression framework for G-E selection, the hazard function at time t for subject i ‘s survival given all environmental factors, one gene factor k and their corresponding interactions in a single model is modeled as

(2)$$\begin{array}{c}\hfill \lambda (\mathrm{t}\mid {𝐞}_{\mathrm{i}},{\mathrm{x}}_{\mathrm{ik}},{𝐰}_{\mathrm{i}\cdot \mathrm{k}})={\lambda }_0(\mathrm{t})\exp ({𝐞}_{\mathrm{i}}^{\prime }𝜶+{\mathrm{x}}_{\mathrm{ik}}{\beta }_{\mathrm{k}}+{𝐰}_{\mathrm{i}\cdot \mathrm{k}}^{\prime }{𝜸}_{\cdot \mathrm{k}}),\mathrm{i}=\hfill \\ \hfill 1,\mathrm{\dots },\mathrm{N};\mathrm{k}=1,\mathrm{\dots },\mathrm{p}.\hfill \end{array}$$

where ${\lambda }_0\left(t\right)$ is a non-negative deterministic baseline hazard function, ${𝒘}_{i\cdot k}=(e_{i1}{x}_{ik},\mathrm{\dots },e_{iq}{x}_{ik})^{\prime }{}_{1\times q}$ and $(𝜶,{\beta }_k,{𝜸}_{\cdot k})$ are corresponding parameters for these considering biomarkers. We define $𝜽={({\alpha }_1,\mathrm{\dots },{\alpha }_q,{\beta }_1,\mathrm{\dots },{\beta }_p,{\gamma }_{11},\mathrm{\dots },{\gamma }_{qp})}_{1\times (q+p+qp)}^{{}^{\prime }}$ as full parameters of the full model.

In a marginal accelerated failure time model framework for G-E selection, the log of survival time for subject i ‘s given all environmental factors, one gene factor k and their corresponding interactions in a single model is modeled as

(3)$$\begin{array}{c}\hfill \log \left({\mathrm{T}}_{\mathrm{i}}\right)={\alpha }_0+{𝐞}_{\mathrm{i}}^{\prime }𝜶+{\mathrm{x}}_{\mathrm{ik}}{\beta }_{\mathrm{k}}+{𝐰}_{\mathrm{i}\cdot \mathrm{k}}^{\prime }{𝜸}_{\cdot \mathrm{k}}+{ϵ}_{\mathrm{i}},\mathrm{i}=1,\mathrm{\dots },\mathrm{N};\hfill \\ \hfill \mathrm{k}=1,\mathrm{\dots },\mathrm{p}.\hfill \end{array}$$

Where ${\alpha }_0$ is an intercept term and $ϵ$ is the error term. Note that the significant G-E interaction can be selected based on the correspondence of the marginal p-value.

Although inverse probability-of-censoring weighted (IPCW) Kendall’s partial correlation [8] was originally developed for G-G interaction selection, this method can naturally be applied to the G-E interaction selection issue, as the concept is the same. Kendall [10] defined the partial rank correlation in the context of Kendall’s correlation, and showed that Pearson’s partial correlation formula still applies to Kendall’s correlation. For example, for four random variables K1; K2; K3; K4, Kendall’s partial correlation is calculated by the following formula

(4)$${\tau }_{12\cdot 34}=\frac{{\tau }_{12\cdot 3}-{\tau }_{14\cdot 3}{\tau }_{24\cdot 3}}{\sqrt{1-{\tau }_{14\cdot 3}^2}\sqrt{1-{\tau }_{24\cdot 3}^2}},$$

gives the Kendall’s partial correlation between K1 and K2 conditional on K3 and K4.

Therefore, the Kendall’s partial correlation of the survival trait with the G-E interaction terms can be obtained as follows,

(5)$$\begin{array}{c}\hfill {\tau }_{T,E_j{G}_k\cdot E_j,G_k}=\frac{{\tau }_{T,E_j{G}_K\cdot E_j}-{\tau }_{T,G_k\cdot E_j}{\tau }_{E_j{G}_k,G_k\cdot E_j}}{\sqrt{1-{\tau }_{T,G_k\cdot E_j}^2}\sqrt{1-{\tau }_{E_j{G}_k,G_k\cdot E_j}^2}},\hfill \\ \hfill j=1,\mathrm{\dots },q;k=1,\mathrm{\dots },p.\hfill \end{array}$$

To accommodate right-censored survival time data, we utilize the IPCW Kendall’s tau statistic proposed by Wang and Chen [8] and consider the resulting partial correlation statistics:

(6)$$\begin{array}{c}\hfill {\tilde{\tau }}_{T,E_j{G}_k\cdot E_j,G_k}=\frac{{\tilde{\tau }}_{T,E_j{G}_k\cdot E_j}-{\tilde{\tau }}_{T,G_k\cdot E_j}{\tau }_{E_j{G}_k,G_k\cdot E_j}}{\sqrt{1-{\tilde{\tau }}_{T,G_k\cdot E_j}^2}\sqrt{1-{\tau }_{E_j{G}_k,G_k\cdot E_j}^2}},\hfill \\ \hfill j=1,\mathrm{\dots },q;k=1,\mathrm{\dots },p.\hfill \end{array}$$

The full computation details of (${\tilde{\tau }}_{T,E_j{G}_k\cdot E_j},{\tilde{\tau }}_{T,G_k\cdot E_j},{\tau }_{E_j{G}_k,G_k\cdot E_j}$) can be seen as follows,

(7)$$\begin{array}{c}\hfill {\tilde{\tau }}_{\mathrm{T},{\mathrm{E}}_{\mathrm{j}}{\mathrm{G}}_{\mathrm{k}}\cdot {\mathrm{E}}_{\mathrm{j}}}=\frac{{\tilde{\tau }}_{\mathrm{T},{\mathrm{E}}_{\mathrm{j}}{\mathrm{G}}_{\mathrm{k}}}-{\tilde{\tau }}_{\mathrm{T},{\mathrm{E}}_{\mathrm{j}}}{\tau }_{{\mathrm{E}}_{\mathrm{j}}{\mathrm{G}}_{\mathrm{k}},{\mathrm{E}}_{\mathrm{j}}}}{\sqrt{1-{\tilde{\tau }}_{\mathrm{T},{\mathrm{E}}_{\mathrm{j}}}^2}\sqrt{1-{\tau }_{{\mathrm{E}}_{\mathrm{j}}{\mathrm{G}}_{\mathrm{k}},{\mathrm{E}}_{\mathrm{j}}}^2}},\mathrm{j}=1,\mathrm{\dots },\mathrm{q};\hfill \\ \hfill \mathrm{k}=1,\mathrm{\dots },\mathrm{p};\hfill \\ \hfill {\tilde{\tau }}_{\mathrm{T},{\mathrm{G}}_{\mathrm{k}}\cdot {\mathrm{E}}_{\mathrm{j}}}=\frac{{\tilde{\tau }}_{\mathrm{T},{\mathrm{G}}_{\mathrm{k}}}-{\tilde{\tau }}_{\mathrm{T},{\mathrm{E}}_{\mathrm{j}}}{\tau }_{{\mathrm{G}}_{\mathrm{k}},{\mathrm{E}}_{\mathrm{j}}}}{\sqrt{1-{\tilde{\tau }}_{\mathrm{T},{\mathrm{E}}_{\mathrm{j}}}^2}\sqrt{1-{\tau }_{{\mathrm{G}}_{\mathrm{k}},{\mathrm{E}}_{\mathrm{j}}}^2}},\mathrm{j}=1,\mathrm{\dots },\mathrm{q};\hfill \\ \hfill \mathrm{k}=1,\mathrm{\dots },\mathrm{p};\hfill \\ \hfill {\tau }_{{\mathrm{E}}_{\mathrm{j}}{\mathrm{G}}_{\mathrm{k}},{\mathrm{G}}_{\mathrm{k}}\cdot {\mathrm{E}}_{\mathrm{j}}}=\frac{{\tau }_{{\mathrm{E}}_{\mathrm{j}}{\mathrm{G}}_{\mathrm{k}},{\mathrm{G}}_{\mathrm{k}}}-{\tau }_{{\mathrm{E}}_{\mathrm{j}}{\mathrm{G}}_{\mathrm{k}},{\mathrm{E}}_{\mathrm{j}}}{\tau }_{{\mathrm{G}}_{\mathrm{k}},{\mathrm{E}}_{\mathrm{j}}}}{\sqrt{1-{\left({\tau }_{{\mathrm{E}}_{\mathrm{j}}{\mathrm{G}}_{\mathrm{k}},{\mathrm{E}}_{\mathrm{j}}}\right)}^2}\sqrt{1-{\left({\tau }_{{\mathrm{G}}_{\mathrm{k}},{\mathrm{E}}_{\mathrm{j}}}\right)}^2}},\mathrm{j}=1,\mathrm{\dots },\mathrm{q};\hfill \\ \hfill \mathrm{k}=1,\mathrm{\dots },\mathrm{p}.\hfill \end{array}$$

Note that the IPCW Kendall’s tau statistic [8] can be computed as follows

(8)

The CQPCorr method [5] consists of three steps, where quantile regression is used to accommodate long-tailed or contaminated responses, partial correlation is used to determine important interactions biomarkers, and the main G and E effects are appropriately controlled. The full detail can be seen in Xu, et al. [5]. Note that there is a tuning parameter in the CQPCorr approach, which is a quantile used in the first and third steps, with quantile $\tau =0.5$ being the most popular choice [5]. The CQPCorr approach can be performed by “QPCorr.matrix” R function of the R package “GEInter” (https://cran.r-project.org/web/packages/GEInter/) [11].

To evaluate the performance of G-E interaction selection, we considered seven measures, such as accuracy rate (ACC), true positive rate (TPR), precision rate (PRE), true negative rate (TNR), false positive rate (FPR), false negative rate (FNR) and minimum of model size (MMS), where the definitions of the first six metrics are displayed in Supplementary Table 1 and MMS is the minimum of model size of a selected set of associated interaction variables, including all potential valid interaction biomarkers. MMS measures the complexity of the selected model and reflects the accuracy of the screening process; larger ACC, TPR, PRE, TNR/smaller FPR, FNR, MMS values indicate higher accuracy of feature screening. We adopt the hard thresholding rule proposed by Fan and Lv [12] to select the candidate set of G-E interactions; that is, after ranking the G-E interaction predictors using a certain correlation measure, we select a prefixed number ($\frac{N}{2*log(N)}$) of top-ranked G-E interaction predictors as our candidate model. We then report the average number of these seven measures for each method in 200 replications.

In order to respect the “main effect, interaction” hierarchical constraint, if we choose a ${𝒘}_{\cdot jk}$ interaction biomarker in the model, then we assume the main factors ${𝒆}_{\cdot j}$ are ${𝒙}_{\cdot k}$ are related to the response, and all environmental factors are considered into the model, we then estimate the corresponding parameters using a maximum-likelihood estimation approach based on the AFT regression model. We report the root mean square error (RMSE) to measure the accuracy of the estimation, which is defined as

(9)$$\text{RMSE}=\sqrt{\frac{1}{S}{\sum }_{j=1}^{\mathrm{S}}{\left({\theta }_j-{\widehat{\theta }}_j\right)}^2},$$

where S is the full model size including all main and interaction covariates. In order to evaluate the estimation of the performance of the selected biomarkers, we report RMSE.M, the average of the root mean square error of 200 replications.

To evaluate the performance of survival prediction, let $\widehat{𝜽}$ be an estimator of the AFT regression parameter in a prediction model obtained from the training dataset and ($V_i^{*}$,${\delta }_i^{\mathrm{*}}$ ,${𝐞}_i^{{\mathrm{*}}^{\mathrm{\prime }}}$ , ${𝐱}_i^{{\mathrm{*}}^{\mathrm{\prime }}}$, ${𝐰}_i^{{\mathrm{*}}^{\mathrm{\prime }}}$) the survival and covariate data of subject i in the test data. Define (${𝒆}_i^{{*}^{\prime }}$, ${𝒙}_i^{{*}^{\prime }}$, ${𝒘}_i^{{*}^{\prime }}$)$\widehat{𝜽}$ as the prognosis index (PI) value for subject i. The AFT test is defined as the p-value of PI, where PI is used as the covariate in the univariate AFT model of survival outcomes in the test data. Similarly, LR test is the p-value of the log-rank test of the null hypothesis of equal survival between the “poor” and “good” prognostic groups in the test data, depending on whether the PI is higher or lower than the median PI value. Also, the c-index metric is considered to investigate survival prediction accuracy. Smaller AFT test and LR test values/larger c-index correspond to better predictive power.

In the following simulations, a series of simulation studies were conducted to compare the existing marginal modeling methods with IPCW Kendall’s partial correlation approach (IPCW-pcorr) in the selection of G-E interactions, the estimation of the selected biomarker effects, and prediction of the final survival prediction model. We also considered a simple measure, IPCW Kendall’s tau correlation (IPCW-tau). IPCW Kendall’s tau correlation [8] measures the association between survival characteristics and G-E interaction biomarkers without conditional main effects.

In order to compare the CQPCorr and IPCW Kendall’s partial correlation approaches clearly, we follow the simulation settings of Xu, et al. [5] by generating a cohort of 200 subjects and 100 subjects for training and testing data respectively. Each subject’s survival time follows the accelerated failure time model,

(10)$$\begin{array}{c}\hfill \log \left({\mathrm{T}}_i\right)=\sum _{j=1}^5{\alpha }_j{e}_{ij}+\sum _{k=1}^{1000}{\beta }_k{x}_{ik}+\sum _{j=1}^5\sum _{k=1}^{1000}{\gamma }_{jk}{w}_{ijk}+{\epsilon }_i,\hfill \\ \hfill i=1,\mathrm{\dots },200,\hfill \end{array}$$

where the covariates e jointly follow a 5-dimensional multivariate standard normal distribution with the first-order autoregressive (AR(1)) structure that is $corr({𝒆}_{.j},{𝒆}_{.k})={0.3}^{|j-k|}$, and the covariates jointly follow a 1000-dimensional multivariate standard normal distribution with the AR(1) structure that is $corr({𝒙}_{.j},{𝒙}_{.k})={0.5}^{|j-k|}$. Moreover, we assume that gene expression may be contaminated by outliers generated from a t-distribution with two degrees of freedom with a probability of 0.1. The outlier generation setting is the same as that of Wang and Chen [8]. Consider three error distributions: (Error 1) N(0, 1), (Error 2) 90% N(0, 1)+10% N($\pm$50, 1) and (Error 3) 80% N(0, 1)+20% N(0, 50). The last two error distributions lead to long-tailed distributions/contamination. The censoring time distribution follows a uniform distribution U(a,b) , which is chosen to control the censoring rate at about 30% (light censoring) and 60% (heavy censoring) respectively.

Five parameters scenarios are considered:

C1 setting has ${\gamma }_{jk}=2$, ${\alpha }_j=1$, ${\beta }_k=1$ for $j=1,2$ and $\text{k}=1,2,\mathrm{\dots },5.$, and ${\gamma }_{jk}=1$ for $j=3,4,5$ and $\text{k}=6,7$. All other coefficients are 0. Under this scenario, the first type interactions are stronger than the corresponding main effects.

C2 setting is the same as C1 setting except that the first type interactions and the corresponding main effects are at the same level. Specifically, ${\gamma }_{jk}={\alpha }_j={\beta }_k=1.5$ for $j=1,2$ and $\text{k}=1,2,\mathrm{\dots },5$.

C3 setting is the same as C1 setting except that the magnitudes of the main effects are larger. Specifically, ${\alpha }_1={\alpha }_2={\beta }_1=\mathrm{\dots }={\beta }_5=3$.

C4 setting is the same as C1 setting except that the magnitudes of the interactions are smaller. Specifically, ${\gamma }_{jk}=0.5$ for $j=1,2$ and $\text{k}=1,2,\mathrm{\dots },5$, and $j=3,4,5$ and $\text{k}=6,7$.

C5 setting is the same as C1 setting except that the first type interactions have negative effects. Specifically, ${\gamma }_{jk}=-2$ for $j=1,2$ and $\text{k}=1,2,\mathrm{\dots },5$.

Each scenario has sixteen important G-E interactions together with two main E effects and five main G effects. There are two types of important interactions. The first type includes ten interactions ${\gamma }_{jk},\text{j}=1,2\text{ and k}=1,2,\mathrm{\dots },5$ with both main $E({\alpha }_1,{\alpha }_2)$ and $G({\beta }_1,\mathrm{\dots },{\beta }_5)$ effects. The second type includes six interactions ${\gamma }_{jk},\text{j}=3,4,5\text{ and k}=6,7$ without main effects, which violates the “main effects, interactions” hierarchy. These simulated settings are close to the actual data. Note that the above simulated survival clinical genomic data can be generated by “simulated_data” function of the R package “GEInter” [11].

The numerical results are summarized in Tables 1,2,3,4,5 for scenarios C1 to C5 with a censoring rate of 30%. In each cell, mean (standard deviation, SD) values are based on 200 replicates. We observe that the performances of the IPCW Kendall’s partial correlation approach are always better than the CQPCorr, marginal Cox, marginal AFT approaches in all evaluation metrics, but the marginal Cox and marginal AFT approaches have smaller variability compared to the IPCW Kendall’s partial correlation approach in C2, C3, and C4 scenarios.

The numerical results are summarized in Supplementary Tables 2–6 for scenarios C1 to C5 with a censoring rate of 60%. In each cell, mean (SD) values are based on 200 replicates. We observe that the performances of the CQPCorr approach are always better than the alternative approaches in all evaluation metrics for scenarios C2, C3 and C4, and the IPCW Kendall’s partial or IPCW-tau correlation approach performs better than others for scenarios C1 and C5. In addition, we observe the marginal Cox and marginal AFT approaches have smaller variability compared to the IPCW Kendall’s partial correlation and CQPCorr approaches.

Table 6 presents simulation findings in tabular form based on censoring rates (30% and 60%) and five parameter scenarios to give readers a better understanding of how these approaches stack up.

After excluding patients with missing survival time data, our analysis is focused on the subset of the TCGA ESCA data with 368 patients and 20,501 gene expression variables. The censoring rate of the survival time in the data is about 58%.

The top 1000 genes with the smallest p-values based on the marginal (univariate) COX model are selected for downstream analysis, since the number of cancer-related genes is expected to be limited. The seven clinical variables whose E effects are analyzed include age, gender, esophageal tumor central location, person neoplasm cancer status, race, BMI and AJCC pathologic stage, and their summary information is reported in the Supplementary Table 7, with its source derived from Wang and Chen [13]. Some of the clinical variables contain missing values, and we use the sparse boosting method [14] in the R package “GEInter” to perform multiple imputation for the missing values in the clinical variables.

We take fifty random splits of the whole data into 258:110 training/test sets of the data to evaluate the performance of all methods for survival prediction via the significant proportion of the LR-test in the TCGA ESCA data. Higher significant proportion of LR-test corresponds to better prediction accuracy. From Table 7 we see that the performance of the IPCW Kendall’s partial correlation approach is better than the marginal Cox, marginal AFT and CQPCorr methods.

In addition, we apply the proposed IPCW Kendall’s partial correlation approach for whole data to identify several important G-E interaction biomarkers and estimate the corresponding parameters by AFT regression model. Table 8 provides the list of selected associated predictors with their correspondence weights, with the “*” notation meaning statistical significance (p-value 0.05). One candidate gene (GADD45B) has been shown to be related to ESCA [15] while Weygant, et al. [16] indicated that the gender factor can be used as prognosis biomarker for ESCA. In addition, we found the GADD45B-gender interaction biomarker to be significant from Table 8, and therefore we consider the GADD45B-gender biomarker to be a potential prognostic biomarker for ESCA, as the major gene GADD45B and gender factor have been documented to be associated with ESCA.

After excluding patients with missing survival time data, our analysis is focused on the subset of the TCGA PAAD data with 170 patients and 20,501 gene expression variables. The censoring rate of the survival time in the data is about 58%.

The top 1000 genes with the smallest p-values based on the marginal (univariate) COX model are selected for downstream analysis, since the number of cancer-related genes is expected to be limited. The seven clinical variables whose E effects are ethnicity, race, lymph node examined count, maximum tumor dimension, anatomic neoplasm subdivision, gender and age, and their summary information is reported in Supplementary Table 8, with its source being from Wang, et al. [13]. Some of the clinical variables contain missing values, and we use the sparse boosting method [14] in the R package “GEInter” to perform multiple imputation for the missing values in the clinical variables.

We take fifty random splits of the whole data into 119:51 training/test sets of the data to evaluate the performance of all methods for survival prediction via the significant proportion of the LR-test in the TCGA PAAD data. Higher significant proportion of LR-test corresponds to better prediction accuracy. From Table 7, we see that the performance of the IPCW Kendall’s partial correlation approach is better than the marginal Cox, marginal AFT and CQPCorr methods.

In addition, we apply the proposed IPCW Kendall’s partial correlation approach for whole data to identify several important G-E interaction biomarkers and estimate the corresponding parameters by AFT regression model. Table 9 provides the list of selected associated predictors with their correspondence weights, with the “*” notation meaning statistical significance (p-value 0.05). Two candidate genes (NRSN2 and TRIM59) have been shown to be related to PAAD [17, 18], while Chakladar, et al. [19] indicated that the gender factor can be used as prognosis biomarker for PAAD. In addition, we found the TRIM59-gender interaction biomarker to be significant from Table 9, and therefore we consider the TRIM59-gender biomarker to be a potential prognostic biomarker for PAAD, as the major gene TRIM59 and gender factor have been documented to be associated with PAAD.

After excluding patients with missing survival time data, our analysis is focused on the subset of the TCGA LUAD data with 505 patients and 20,501 gene expression variables. The censoring rate of the survival time in the data is about 64%.

The top 1000 genes with the smallest p-values based on the marginal (univariate) COX model are selected for downstream analysis, since the number of cancer-related genes is expected to be limited. The eight clinical variables whose E effects are analyzed include age at initial pathologic diagnosis, number pack years smoked, AJCC pathologic metastasis, AJCC pathologic nodes, AJCC pathologic stage, race, gender and AJCC pathologic tumor, with their summary information reported in Supplementary Table 9, sourced from Wang, et al. [13]. Some of the clinical variables contain missing values, and we use the sparse boosting method [14] in the R package “GEInter” to perform multiple imputation for the missing values in the clinical variables.

We take fifty random splits of the whole data into 354:151 training/test sets of the data to evaluate the performance of all methods for survival prediction via the significant proportion of the LR-test in the TCGA LUAD data. Higher significant proportion of LR-test corresponds to better prediction accuracy. From Table 7, we see that the performance of the IPCW Kendall’s partial correlation approach is better than the marginal Cox, marginal AFT and CQPCorr methods.

In addition, we apply the proposed IPCW Kendall’s partial correlation approach for whole data to identify several important G-E interaction biomarkers and estimate the corresponding parameters by AFT regression model. Table 10 provides the list of selected associated predictors with their correspondence weights, with the “*” notation meaning statistical significance (p-value 0.05). There are six candidate genes (ATP13A4, GPR116, HABP2, MBIP, SFTA3, and ZSCAN16) have been shown to be related to LUAD [20, 21, 22, 23, 24].