\section*{Introduction}

For many complex diseases such as cancer, diabetes, and cardiovascular diseases, prognosis is of essential interest. In the omics era, profiling studies have been extensively conducted, searching for genetic markers associated with prognosis. It has been suggested that, beyond the main effects of genetic (G) and environmental (E) risk factors, the gene-environment (G$\times$E) interactions also have important implications. Multiple statistical methods have been developed for G$\times$E interaction analysis. For reviews, refer to~\cite{Caspi2006gene},~\cite{Cordell2009detecting},~\cite{Thomas2010methods}, and others.

Denote $T$ as the prognosis time of interest, $X = (X_1, \ldots, X_q)$ as the $q$ environmental/clinical variables, and $Z=(Z_1, \ldots, Z_p)$ as the $p$ genetic variables. Assume $n$ independent subjects. Regression-based interaction analysis, with its broad applicability, has been extensively adopted and proceeds as follows. (a) For gene $j (=1, \ldots, p)$, consider the model $T \sim \phi( \alpha_{j,0}+\sum_{k=1}^qX_k\alpha_{j,k} + Z_j\beta_j + Z_j\sum_{k = 1}^qX_k\gamma_{j,k})$, where $\phi(\cdot)$ is the known link function (for example, the Cox or exponential model), and $\alpha_{j,0},\alpha_{j,k}, \beta_{j}, \gamma_{j,k}$ are the unknown regression coefficients. As usually $q<<n$, this is a low-dimensional model and can be fitted using standard, usually likelihood-based, techniques and software. Denote $p_{j,k}$ as the p-value for $\gamma_{j,k}$. (b) With $\{p_{j,k}: j=1, \ldots, p; k=1, \ldots, q \}$, conduct multiple comparison adjustment using the Bonferroini or FDR (false discovery rate) approach, and identify important interactions.

A limitation of the above approach is a lack of robustness properties. Usually it is assumed that all subjects satisfy the same prognosis models. In practice, most genetic studies cannot afford to conduct strict subject selection. Seemingly homogeneous subjects can have different disease subtypes~\citep{Haibe-Kains2012three} and different survival patterns~\citep{burgess2011cancer}. Cause of death can be misclassified, leading to contamination in disease-specific survival~\citep{fall2008reliability}. The survival times extracted from medical records are not always reliable (\cite{bowman2011doctors},~\cite{rampatige2013assessing}). With unrobust for example likelihood-based estimation, even a single contaminated observation can result in severely biased estimates~\citep{huber2009robust}, which can lead to false marker identification. Another limitation is that, significance level-based identification, although {\it asymptotically} valid, may generate unreliable results when the sample sizes are small to moderate, as in typical profiling studies. A recent study suggests that regularized estimation can lead to more reliable estimation and hence more accurate marker identification~\citep{shi2014gene}.

With low-dimensional biomedical data, robust statistical methods have demonstrated great power. As suggested in a recent review by~\cite{wu2014selective}, with genetic data, the development is limited and unsystematic and has been mostly on the analysis of main effects but not interactions. In the literature, relevant studies include~\cite{gui2011novel}, which identifies important interactions using the multifactor dimensionality reduction (MDR) technique. However, this method is limited to categorial data (such as SNPs) and not broadly applicable.~\cite{shi2014gene} developed a rank-based method, which is robust to model mis-specification but not data contamination. The most relevant study is perhaps~\cite{wang2015identifying}, which developed a quantile-regression based method. With that method, the quantile needs to be specified, which is not a trivial task in practical data analysis. The objective function is not differentiable, causing difficulty in estimation. In addition, as to be shown in this study, its numerical performance can be less satisfactory under many data settings. With low-dimensional data, it has been shown that no robust method dominates the others. It is thus prudent to develop alternative robust methods.

Consider prognosis data with both G and E measurements. The goal is to develop a new method for identifying important G$\times$E interactions. Significantly advancing from the existing studies, the proposed method is robust to contamination in the prognosis data. In addition, our simulation suggests that, under certain scenarios, its numerical performance is better than the quantile regression analysis, which is perhaps the most popular robust method for genetic data~\citep{wu2014selective}. A penalization approach is adopted for marker identification, which differs from the significance level-based approach and can have better numerical performance when the sample size is small to moderate.


\section*{Robust identification of interactions}

\subsection*{Data and model settings}

In the literature, there are multiple prognosis models. In this study, we consider the accelerated failure time (AFT) model, which has been adopted in~\cite{shi2014gene}, \cite{liu2013identification}, and many others. Compared to alternatives such as the Cox model, the AFT model can be preferred because of its more intuitive interpretations and lower computational cost, which are especially desirable with high-dimensional genetic data. With a slight abuse of notation, still $T$ to denote the logarithm of prognosis time. For gene $j$, the AFT model assumes that
$$T = \alpha_{j,0}+\sum_{k=1}^qX_k\alpha_{j,k} + Z_j\beta_j + Z_j\sum_{k = 1}^qX_k\gamma_{j,k}+\epsilon,$$
where $\epsilon$ is the random error.

Assume $n$ independent observations. Consider the scenario where a small subset of the random errors are contaminated, leading to contamination in the prognosis times. In a typical profiling study, $q<<n$, while $p$ can be comparable to or much larger than $n$. For subject $i$, denote $C_i$ as the logarithm of censoring time and $\mathbf{x}_i$ and  $\mathbf{z}_i$ as the observed $X$ and $Z$ values, respectively. Under right censoring, we observe $(y_i=min(T_i, C_i), \delta_i=I(T_i \leq C_i), \mathbf{x}_i, \mathbf{z}_i)$. Further denote $\mathbf{u}_{i,j}=(1,\mathbf{x}^{\prime}_i, z_{i,j},z_{i,j}x_{i,1},\ldots,z_{i,j}x_{i,q})'$,  $\bm\zeta_j=(\alpha_{j,0},\alpha_{j,1},\ldots,\alpha_{j,q},\beta_j,\gamma_{j,1},\ldots, \gamma_{j,q})'$, and $\mathbf{U}_j=(\mathbf{u}_{1,j}, \ldots, \mathbf{u}_{n,j})'$ which is a $n\times (2q+2)$ matrix. For gene $j$ and subject $i$, the AFT model can now be written as
$T_i = \mathbf{u}_{i,j}^\top\bm{\zeta}_j + \epsilon_{i,j}.$ Without loss of generality, assume that the data $\{(y_i,\delta_i,\mathbf{x}_i,\mathbf{z}_i), i=1,\ldots, n\}$ have been sorted according to $y_i$'s from the smallest to the largest.


\subsection*{Penalized robust identification}

A penalized marker identification method is defined by its loss function and penalty, which are defined separately as follows.

\noindent{\bf Loss function} Before defining the robust loss function, we take one step back and consider the scenario with no data contamination. When the distribution of $\epsilon$ is not specified, likelihood-based estimation cannot be adopted. A popular approach is the weighted least squared estimation proposed in~\cite{Stute1993consistent} and proceeds as follows. First compute the Kaplan-Meier weights as $$\omega_{1}=\frac{\delta_{1}}{n},~\omega_{i}=
\frac{\delta_{i}}{n-i+1}\prod_{j=1}^{i-1}
(\frac{n-j}{n-j+1})^{\delta_{j}},~i=2,\cdots,n.$$
For gene $j(=1,\ldots, p)$, the weighted least squared objective function is defined as
\begin{eqnarray}
\sum_{i=1}^n \omega_{i} (y_i - \mathbf{u}_{i,j}^\top\bm\zeta_j)^2,
\end{eqnarray}
This loss function has the following properties. With a simple linear regression model, the most commonly adopted loss function is the least squared one. Here to accommodate censoring, a weight function is imposed to re-weigh different observations according to their observed times and event status. When there is no censoring, $\omega_i=1/n$. With the quadratic form, the loss function is not robust to data contamination. If subject $i$ is not censored, then $\omega_i\neq 0$, and an arbitrarily large $y_i$ results in arbitrarily large estimates. Biased estimation can happen with data contamination, which may lead to false marker identification.

To accommodate data contamination, for gene $j(=1,\ldots, p)$, we propose the exponential squared loss function
\begin{eqnarray}
Q_{\theta}(\bm\zeta_j|\mathbf{y}, \mathbf{U}_j,\bm\omega)=
\sum_{i=1}^n\omega_i\exp(-(y_i - \mathbf{u}_{i,j}^\top\bm\zeta_j)^2 / \theta).
\end{eqnarray}
$\mathbf{y}$ and $\bm\omega$ denote the vectors composed of $y_i$'s and $\omega_i$'s, respectively. $\theta>0$ is a tuning parameter. The rationale of this approach is as follows. For contaminated subjects with $y_i$s deviate from $\mathbf{u}_{i,j}^\top\bm\zeta_j$s (predicted values from the model), $(y_i - \mathbf{u}_{i,j}^\top\bm\zeta_j)^2$s have large values.
The exponential function down-weighs such contaminated observations. The degree of down-weighing is adjusted by $\theta$: with $\theta$ getting smaller, contaminated observations have smaller influence. To accommodate censoring, $\omega_i$'s are imposed in a similar manner as the original Stute's approach. For low-dimensional linear regression model without censoring, the exponential squared loss has been examined in~\cite{wang2013robust}. Different from the existing studies, here we consider the more challenging high-dimensional genetic data especially interactions. In addition, the Kaplan-Meier weights are introduced to accommodate censoring. As to be shown in Appendix, such differences lead to significant differences in the statistical development.

\noindent{\bf Penalized estimation}
For gene $j(=1,\ldots, p)$, we propose the penalized objective function
\begin{eqnarray}
L_{\lambda,\theta}(\bm\zeta_j|\mathbf{y}, \mathbf{U}_j,\bm\omega) = Q_{ \theta}(\bm\zeta_j|\mathbf{y}, \mathbf{U}_j,\bm\omega)-\lambda||\bm\zeta_j||_1.
\label{eqn2}
\end{eqnarray}
$\lambda>0$ is the data-dependent tuning parameter, and $||\cdot||_1$ is the $\ell_1$ norm. Denote $\tilde{\bm\zeta_j}$ as the maximizer of $L_{\lambda,\theta}(\bm\zeta_j|\mathbf{y}, \mathbf{U}_j,\bm\omega)$. Interactions (and main effects) corresponding to the nonzero components of $\tilde{\bm\zeta_j}$ are identified as important.

The strategy of using penalization for identifying important interactions has been adopted in~\cite{shi2014gene} and other studies. Here a Lasso penalty is imposed, which is the most commonly adopted penalty. $\lambda$ controls the degree of sparsity and number of identified interactions. We impose the same $\lambda$ to all genes to ensure comparability.

Multiple penalties can take place of Lasso. In some recent studies such as~\cite{bien2013lasso} and~\cite{liu2013identification}, penalties have been developed to respect the ``main effects, interactions" hierarchy, which reinforces that the main effects corresponding to the identified interactions must be identified. In a large number of studies (\cite{zimmermann2011interaction},~\cite{Caspi2006gene}), it has been observed that genes can have important G$\times$E interactions but no main effects. In addition, if the hierarchy has to be reinforced, we can identify the important interactions first, and then add back the corresponding main effects. Computationally, Lasso is much simpler than the existing alternatives. Our limited experience suggests that with the complex robust loss function, complex penalties have a considerable probability of running into convergence problems. With the above considerations, the Lasso penalty is adopted here.


\noindent{\bf Consistency properties} A significant advantage of the proposed method is that it has the much deserved consistency properties under the ultrahigh dimensional setting. For many of the existing interaction analysis methods, there is a lack of theoretical development. The rigorously established consistency properties provide a solid statistical ground for the proposed method. Details are provided in Appendix.


\subsection*{Computation}

We propose a coordinate-wise updating procedure to compute the solution to~\eqref{eqn2}. For low-dimensional data under simpler settings, an iterative approach is suggested in~\cite{wang2013robust} to select the robust tuning parameter $\theta$. However, under the present high-dimensional settings and with the coordinate-wise updating procedure, such an approach is computationally infeasible. Alternatively, for each $(\lambda,\theta)$ pair, we compute the solution to each marginal model. This way, we can obtain a solution surface over the two-dimensional tuning parameter grid. Then a prediction error based method such as cross-validation can be used to select the robust and penalization tuning parameters. Computation is conducted for each gene separately and can be realized in a highly parallel manner to reduce computer time.

Consider gene $j(=1,\ldots, p)$. Let $r_i(\bm\zeta_j)=y_i-\mathbf{u}_{i,j}^\top\bm\zeta_j$. The first and second order derivatives of $Q(\bm\zeta_j)$ are
\begin{equation}
\begin{cases}
\dot{Q}_k(\bm\zeta_j)=\frac{\partial Q(\bm\zeta_j)}{\partial \zeta_{j,k}}=2\sum_{i=1}^n\omega_iu_{(i,j)_k}r_i(\bm\zeta_j)
\exp(-r_i^2(\bm\zeta_j)/\theta) /\theta, \\
\ddot{{Q}}_{kl}(\bm\zeta_j)=\frac{\partial^2 Q(\bm\zeta_j)}{\partial \zeta_{j,k}\partial \zeta_{j,l}}=2\sum_{i=1}^n\omega_iu_{(i,j)_k}u_{(i,j)_l}\exp(-r_i^2(\bm\zeta_j)/\theta) (2r_i^2(\bm\zeta_j)/\theta-1)/\theta.
\end{cases}
\label{eqn4}
\end{equation}
For a $\bm{\zeta}_j^m$ that is in a small neighborhood of $\bm{\zeta}_j$,  $Q(\bm\zeta_j)$ in~\eqref{eqn2} can be locally approximated by
\begin{equation}
Q(\bm\zeta_j) \approx Q(\bm\zeta^m_j) + \dot{Q}(\bm\zeta^m_j)^\top(\bm\zeta_j-\bm\zeta_j^m)+
\frac{1}{2}(\bm\zeta_j-\bm\zeta^m_j)^\top\ddot{Q}
(\bm\zeta^m_j)(\bm\zeta_j-\bm\zeta^m_j).
\label{eqn5}
\end{equation}
Replacing $Q(\bm\zeta_j)$ in~\eqref{eqn2} with this approximation and taking the first order derivative of $L(\bm\zeta_j)$ with respect to the $k$th element $\zeta_{j,k}$ give us
\begin{equation}
\zeta_{j,k}^{m+1}=\zeta_{j,k}^m -\ddot{Q}^{-1}_{kk}(\bm{\zeta}^m_j)\dot{Q}_k(\bm{\zeta}^m_j)+ \ddot{Q}_{kk}^{-1}(\bm{\zeta}^m_j)\lambda.
\end{equation}
Note that when $2r_i^2(\bm{\zeta}_j)/\theta>1$, $\ddot{Q}_{kk}(\bm{\zeta}_j)\geq0$. Then \eqref{eqn2} can be maximized at infinity, and the algorithm may fail to converge. To avoid this problem, notice that $$\ddot{Q}_{kk}(\bm{\zeta}_j)\geq -2\sum_{i=1}^n\omega_iu_{(i,j)_k}^2\exp(-r_i^2(\bm\zeta_j)/\theta)/\theta$$ for $\omega_i\geq 0$. The right hand side is always non-positive. We re-define $$\ddot{Q}_{kk}(\bm{\zeta}_j)\equiv -2\sum_{i=1}^n\omega_iu_{(i,j)_k}^2\exp(-r_i^2(\bm\zeta_j)/\theta)/\theta,$$ and use the minorization-maximization (MM) algorithm to ensure convergence. The algorithm that combines the coordinate descent and MM algorithm is summarized in Algorithm~\ref{alg1}.
\begin{algorithm}
\caption{The coordinate descent MM algorithm for marginal model $j$}
\label{alg1}
\SetKwInOut{Input}{input}\SetKwInOut{Output}{output}
\Input{$\mathbf{y}, \bm{\omega}, \mathbf{U}_j$, and tuning parameters $\lambda$ and $\theta$.}
\Output{the maximizer $\tilde{\bm{\zeta}}_j$ defined in~\eqref{eqn2}.}
\BlankLine
\textbf{initialization}: let $m = 0, \bm{\zeta}^0_j=\mathbf{0}$. Normalize $\mathbf{y}$ and $\mathbf{U}_j$ such that $\sum_{i=1}^n\omega_iy_i=0$, $\sum_{i=1}^n\omega_i u_{(i,j)_k}=0$ and $\sum_{i=1}^n\omega_iu_{(i,j)_k}^2=n$ for $k=1,\ldots,2q+2$.\;
\Repeat{$||\bm{\zeta}^{m_0}_j - \bm{\zeta}^{m}_j||_2\leq$ a predefined thresholding value}{
$m_0 = m$\;
\For{$k = 1, 2, \ldots, 2q+2$}{
\Repeat{$|{\zeta}^{m}_{j,k} - {\zeta}^{m-1}_{j,k}|\leq$ a predefined thresholding value. In our numerical study, we set the thresholding =$10^{-3}$}{
$\zeta_{j,k}^{m+1} = \zeta_{j,k}^m - \ddot{Q}^{-1}_{kk}(\bm{\zeta}^m_j)\dot{Q}_k(\bm{\zeta}^m_j)+ \ddot{Q}_{kk}^{-1}(\bm{\zeta}^m_j)\lambda$\;
$\zeta_{j,l}^{m+1} = \zeta_{j,l}^m$ for $l = 1, 2, \ldots, 2q+2$, $l\neq k$\;
$m = m + 1$\;
}
}
}
\Return{$\bm{\zeta}^m_j$} as the maximizer.
\end{algorithm}



As mentioned above, we search over a two-dimensional grid of $(\lambda,\theta)$ for the optimal. The range of $\lambda$ is determined as follows. First, its upper bound $\lambda_{\max}$ is selected such that $\tilde{\bm{\zeta}}_j=\mathbf{0}$ for $j=1,\ldots,p$. With the unrobust weighted least squared loss, $\lambda_{\max} = \max_{j=1}^p\{||\mathbf{U}_j^\top W \mathbf{y}||_{\infty}\}$ where $W$ is the diagonal matrix composed of $\omega_i$s. With the proposed robust method, the derivatives in~\eqref{eqn4} can be viewed as a weighted sum of $u_{(i,j)_k}r_i(\bm{\zeta}_j)$s. Because the weight for each subject changes with $\bm{\zeta}_j$, the previously defined $\lambda_{\max}$ may not guarantee that $\tilde{\bm{\zeta}}_j=\mathbf{0}$ for $j$. After some trials, we find that $\lambda_{\max} = 20\max_{j=1}^p\{||\mathbf{U}_j^\top W \mathbf{y}||_{\infty}\}$ is in general a safe upper bound for $\lambda$. The lower bound is chosen as $\lambda_{\min}=\lambda_{\max} / 1000.$ The range of $\theta$ depends on the degree of contamination in data, which is not known in practice. In our numerical studies, we find that the estimator is not very sensitive to the value of $\theta$. For cautions, we choose a relatively wide range for $\theta$. Specifically, after centralization, we choose $\theta \in (\min_{i=1}^n{y_i^2}/100, \max_{i=1}^ny_i^2 \times 100)$.

With the quadratic approximation and minorization, optimization in~\eqref{eqn2} is non-convex. Therefore, there is a risk of converging to local maximums. Thus when trying to solve the estimate at a new $(\lambda, \theta)$ point, instead of using $\bm{\zeta}^0_j=\mathbf{0}$ as the starting value, we use the estimate at a neighboring $(\lambda, \theta)$ point as the warm starting value. This modification speeds up the computation and also improves the convergence property.


\section*{Simulation}

In simulation, we set $n= 300$, $q=3$, and $p = 500, 1000$. There are a total of 18 nonzero effects: 3 main effects of E factors, 5 main effects of G factors, and 10 interactions. The positions of nonzero main G effects and interactions are uniformly placed. The nonzero regression coefficients are randomly generated from $uniform(0.5, 1.5)$. The E and G factors are generated from multivariate normal distributions with marginal means zero, marginal variances one, and the following variance matrix structures: Independent, AR(0.2), AR(0.8), Band(0.3), and Band(0.6). Under the AR($\rho$) correlation structure, for the $i$th and $j$th factors, $corr=\rho^{|i-j|}$. Under the Band($\rho$) correlation structure, for the $i$th and $j$th factors, $corr=\rho I(|i-j|\leq 2)$, where $I(\cdot)$ is the indicator function. Under each correlation scenario, consider seven different distributions for the random error $\epsilon$: N(0,1), 0.95N(0,1)+0.05Cauchy, 0.85N(0,1)+0.15Cauchy, 0.7N(0,1)+0.3Cauchy, 0.95N(0,1)+0.05t(3), 0.85N(0,1)+0.15t(3), and 0.7N(0,1)+0.3t(3). That is, we consider the scenarios with no contamination and three different levels of contamination. Two contamination distributions are considered with different thickness of tails. The log event times are generated from the AFT models. The censoring times are generated independently from exponential distributions. The parameters are adjusted so that the censoring rates are about 25\%.

Beyond the proposed method (referred to as ``Robust"), we also consider the following three alternatives: (a) the ``Unrobust" method, which adopts the weighted least squared loss function and applies the Lasso penalization for selecting important effects; (b) the ``Stute" method, which adopts the weighted least squared loss function, does not apply any penalization, and uses the significance level (p-value) as the criterion for quantifying the importance of effects; and (c) the ``Quantile" method, which adopts the quantile regression-based robust loss function and applies the Lasso penalization for selecting important effects. The Unrobust and Stute methods, as the proposed, are built on the weighted least squared estimation. Comparing with them can establish the advantage of robust loss and penalization-based identification, respectively. The Quantile method is also robust and has been recently developed. Comparing with it can establish the advantage of the weighted exponential squared loss. We acknowledge that multiple methods are potentially applicable to the simulated data. The three alternatives have frameworks closest to that of the proposed.

With the Robust, Unrobust, and Quantile methods, the numbers of selected interactions depend on the tuning parameter values. With the Stute method, the number depends on the p-value cutoff. To eliminate the effect of tuning parameter selection on identification, we examine a sequence of tuning parameter values, evaluate identification performance at each value, and use the ROC (receiver operating characteristics) curve to evaluate the interaction identification accuracy of different methods. A representative ROC plot is shown in Figure 1. In this plot, the proposed method has the dominatingly better accuracy.

Summary AUCs based on 100 replicates are shown in Tables 1 ($p=500$) and 2 ($p=1000$). When there is no contamination, performance of the proposed method is comparable to or slightly worse than that of the unrobust method. For example with $p=500$ and the independence correlation structure, the robust and unrobust methods have mean AUCs 0.861 and 0.889, respectively. And with the AR(0.2) correlation, they have mean AUCs 0.901 and 0.881, respectively. With contamination, the robust method outperforms the unrobust method. For example with $p=500$, the AR(0.2) correlation structure, and 0.7N(0,1)+0.3Cauchy error, the robust and unrobust methods have mean AUCs 0.886 and 0.751, respectively. Under all simulation scenarios, the Stute method, which adopts the robust loss function but significance level-based identification, has inferior identification accuracy. As has been suggested in the literature, with a moderate sample size, the unregularized estimates can be less reliable, leading to inaccurate identification. When comparing the proposed method with the quantile regression-based, we see that under the majority of the settings, the proposed method has superior performance. For example with $p=500$, AR(0.2) correlation, and 0.85N(0,1)+0.15Cauchy error, the proposed method and Quantile have mean AUCs 0.892 and 0.842, respectively. However, under a small number of settings, the Quantile method excels. For example with $p=500$, AR(0.8) correlation, and 0.7N(0,1)+0.3t(3) error, the two methods have mean AUCs 0.863 and 0.933, respectively. Such observations are also reasonable. No robust approach is expected to be able to dominate all others. The proposed method outperforms the quantile regression-based under most scenarios and provides a useful alternative.

\section*{Analysis of the lung squamous cell carcinoma data}

Lung squamous cell carcinoma is the second most common lung cancer and causes around 400,000 deaths each year worldwide~\citep{comprehensive2012TCGA}. Profiling studies have been extensively conducted searching for its prognostic factors. In this section, we analyze the TCGA (The Cancer Genome Atlas) data on the prognosis of lung squamous cell carcinoma. The TCGA data were recently collected and published by NCI and have a high quality. The dataset we analyze was downloaded in April of 2015.

The prognosis outcome of interest is overall survival. Multiple environmental and clinical variables are available for analysis. We first remove variables with a low quality of measurements (for example with a high missing rate). We then select the following variables for analysis: age, gender (female is coded as baseline), smoking pack years, and smoking status (non-smoker, reformed smoker for more than 15 years, reformed smoker for less than or equal to 15 years, current smoker; coded as 0, 1, 2 and 3). These variables have been suggested in the literature (\cite{miller2004Bronchioloalveolar} and~\cite{nakachi1991genetic}). A total of 422 samples have clinical and environmental measurements available.

For the genetic part, we analyze gene expression data. A total of 18,969 measurements are available for 502 samples. When matching the clinical/environmental data with genetic data, complete data are available for 404 samples. Among them, 129 died during following. The median followup was 30 months. 275 were censored, and the median followup was 18 months.

The tuning parameter can be chosen using data-driven methods for example cross validation. However, we note that the commonly used methods are mostly based on the notion of prediction. With a large number of marginal models, the goal is to identify markers top-ranked in a marginal sense,
not prediction. Thus, following published studies, we vary the tuning so that a prefixed number of interactions are selected. In Table~\ref{tab3}, we provide results on the 33 top-ranked interactions. Longer or shorter lists of identified interactions are available from the authors.
After the interactions are identified, we refit the marginal models without penalization on the main effects to satisfy the ``main effects, interactions" hierarchy. Beyond the proposed method, we also apply the three alternatives considered in simulation. Table~\ref{tab7} in Appendix suggests that different methods identify significantly different genes and interactions. Results under the proposed method are provided in Table~\ref{tab3}. Those under the alternative methods are provided in Appendix.

To assess the stability of our findings, we apply a cross validation-based approach. Specifically, one sample is removed from data, and the proposed method is applied. This step is repeated over all samples. For each identified interaction in Table 3, we compute its probability of being identified in the reduced datasets. It can be seen that three interactions have very small stability measures. All other interactions have stability measures close to 1. We have also examined those interactions not identified and found that their stability measures are all close to 0. This analysis suggests a certain degree of stability of the proposed method.

Literature search suggests that the identified interactions and corresponding genes may have important implications. Specifically, previous studies \citep{shi2014gene} have suggested that the major G$\times$E interactions occur between genes and smoking status. Among the identified gene-smoking interactions in the current study, the CEBPB (CCAAT/enhancer-binding proteins beta) protein is a transcription factor that works with other CCAAT/enhancer-binding protein family members in the regulation of cell cycle progress, differentiation and pro-inflammatory gene expression. CEBPB has been found to be upregulated by tobacco smoke in both human lung fibroblast~\citep{Miglino2012Cigarette} and mice emphysema~\citep{hirama2007increased}. Another gene that interacts with smoking is
EFNA1 (a.k.a. ephrin-A1). A recent {\em in vitro} study found that ephrin-A1 is overexpressed in tobacco smoke treated bronchial airway epithelial cells compared to control cells~\citep{nasreen2014tobacco}. In addition, the elevated expressions of ephrin-A1 are positively associated with tumor proliferative capacity in non-small cell lung carcinoma patients~\citep{giaginis2014ephrin}. The gene that shows both a strong main effect and interaction with duration of smoking is TERF1 (telomeric repeat binding factor 1). The TERF1 expression levels have been found to be decreased in lung cancer (\cite{lin2006expression},~\cite{hu2009analysis}), as well as other types of cancer in several studies (\cite{yamada2002down},~\cite{miyachi2002correlation}). It functions as an inhibitor of telomerase and is identified as a prognostic marker for overall survival in non-small cell lung cancer~\citep{hu2009analysis}. The molecular mechanism of why and how TERF1 decreases in the process of cancer is not clear. Our results suggest that smoking can be one of the factors.
Other than smoking duration, we also find that several genes interact with smoking intensity as well. The gene that interacts with both gender and smoking intensity is KATNB1. KATNB1 encodes protein katanin p80 subunit B1, which has been found to participate in cytokinesis by interacting with tumor suppressor gene LAPSER1. The disruption of cytokinesis process may potentially cause genetic instability and cancer~\citep{Sudo2007LAPSER1}.
In addition to smoking, we find eight genes interact with gender in Table~\ref{tab3}. Among these, STRADB and CA5BP1 draw our attention. STRADB is an important gene in lung cancer progression and metastasis through the activation of LKB1. LKB1 is essential for G1 cell cycle arrest, cell polarity and stress, cell detachment and adhesion~\citep{marcus2010LKB1}. The STRADB encoded protein also interacts with the X chromosome-linked inhibitor of apoptosis protein by enhancing its anti-apoptotic activity. In addition, gene CA5BP1 is located on X chromosome and found to be gender-associated.


\section*{Discussion}
G$\times$E interactions have important implications for the prognosis of a large number of complex diseases. In this study, we have developed a novel interaction analysis method. It is robust to contamination in the prognosis outcome, which is not rare in practice but cannot be accommodated using most of the existing methods. In addition, we adopt penalization for identifying important interactions. This strategy differs from the commonly adopted significance level-based identification and can have better marker identification accuracy. Significantly advancing from most of the existing interaction studies, the consistency property has been rigorously established under the ultrahigh dimensional setting, making the proposed method one of the very few with a strong theoretical basis. An effective computational algorithm has been developed. In simulation study, the proposed method is observed to have satisfactory performance. With contamination, robust methods are needed. Under the majority of the cases, the proposed method outperforms the quantile regression based method. Thus, it provides a useful alternative to the quantile regression-based and other existing methods. In the analysis of lung cancer data, the proposed method identifies meaningful genes and interactions, which differ from those using the alternatives, with satisfactory stability.

We choose the popular Lasso penalty because of its low computational cost and satisfactory performance. A limitation of Lasso is that its results may not respect the ``main effects, interactions" hierarchy. However, as shown in Appendix, the proposed method can consistently identify the important interactions. As demonstrated in data analysis, once the important interactions are identified, if the hierarchy is desired, the main effects can be added back. With a robust loss function, the proposed method is computationally more expensive. However, as the computing can be executed in a highly parallel manner, and with the effective coordinate descent algorithm, the computer time is very much affordable. A total of 70 scenarios are considered in simulation and show the satisfactory performance of proposed method. More comprehensive comparisons of different methods, especially robust methods, have not been conducted in the literature and are of interest to future studies. Interesting findings are made in the lung cancer prognosis data analysis. Unfortunately there is still a lack of objective way of determining which set of results is more sensible. Validation of the findings will be pursued in future studies.



