{'fraction_non_alphanumeric': 0.03844180592201701, 'fraction_numerical': 0.015537965406039285, 'mean_word_length': 4.596595570139459, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Sparse group penalized integrative analysis of multiple cancer prognosis datasets

JIN LIU, J IAN HUANG, YANG XIE AND SHUANGGE MA* Division of Epidemiology and Biostatistics, UIC School of Public Health 1603 W. Taylor Street (MC 923), Chicago, IL 60612-4394, USA Department of Statistics and Actuarial Science, and Department of Biostatistics, University of Iowa, Iowa, USA Department of Clinical Sciences, UT Southwestern Medical Center, 5323 Harry Hines Blvd., Dallas, TX 75390, USA 4Department of Biostatistics, Yale University; and VA Cooperative Studies Program Coordinating Center, West Haven, CT, USA


Introduction
High-throughput genomic studies have been extensively conducted, searching for markers associated with the risk and prognosis of cancer. In this paper, we focus on cancer prognosis studies with gene expression measurements. Data generated in highthroughput studies have the 'large d, small n' characteristic, with the number of genes profiled d much larger than the sample size n. In addition, in whole-genome studies, only a subset of the profiled genes is associated with prognosis. Thus, analysis of cancer prognosis data with high-throughput genomic measurements demands simultaneous regularized estimation and selection.
It has been recognized that genomic markers identified from the analysis of single datasets can be unsatisfactory. Multiple factors may contribute to the unsatisfactory performance, including the highly noisy nature of cancer genomic data, technical variations of profiling techniques and, more importantly, small sample sizes of individual studies. Recent studies have shown that pooling and analysing multiple studies may effectively increase sample size and improve properties of the identified markers (Guerra & Goldsterin, 2009). For example, simulation and data analysis have shown that markers identified in multidataset analysis may have more true positives, fewer false positives and better prediction performance (Ma et al., 2009(Ma et al., , 2011. Multi-dataset methods include meta-analysis and integrative analysis methods. Integrative analysis pools and analyses raw data from multiple studies, and differs significantly from classic meta-analysis, which analyses multiple studies separately and then pools summary statistics (lists of identified genes, P-values, effect sizes, etc.).
In this paper, we conduct integrative analysis of multiple cancer prognosis studies, with note that the proposed method may also be applicable to other types of data (e.g. aetiology studies with categorical responses and treatment studies with continuous markers). In studies such as Ma et al. (2011 b), the homogeneity model has been assumed. Under this model, multiple datasets have exactly the same set of markers. This model has also been adopted with cancer diagnosis studies (Ma et al., 2011 a and references therein). The existing approaches have been designed to reinforce this homogeneity assumption (more details available in Section 2). In practice, different datasets may differ in terms of patient selection criteria, platforms and protocols used for profiling, normalization methods, and other aspects. Such differences may make the homogeneity assumption too restricted. In addition, data analyses in Ma et al. (2011 a, b) show that for some identified genes, the magnitudes of estimated regression coefficients may vary significantly across datasets. It is possible that the very small regression coefficients are actually zero. This can also be seen from data analysis presented in Section 4. Such an observation further suggests the necessity of relaxing the homogeneity assumption.
In what follows, we describe cancer survival using the accelerated failure time (AFT) model. Among the popular semiparametric survival models, the AFT model may have the lowest computational cost, which is especially desirable with high-dimensional data. In addition, its regression coefficients have more lucid interpretations. As an alternative to the homogeneity model, we consider the heterogeneity model. It includes the homogeneity model as a special case and is more flexible. For marker selection, we adopt a sparse group penalization method. The proposed penalization is intuitively reasonable and computationally feasible. This study complements the existing prognosis studies by conducting integrative survival analysis under the more flexible heterogeneity model and by adopting a penalization approach that has satisfactory performance under both the homogeneity and heterogeneity models. In a recent paper (Liu et al., 2013), we have studied the heterogeneity model and sparse penalization for binary data and logistic regression models. As the present data and model settings are significantly different from Liu et al. (2013), the performance of the sparse group penalization needs to be 're-examined ', the computational algorithm needs to be adjusted, and so the present study is warranted. In addition, this study may be the first to observe that for prognosis data under the homogeneity model, sparse penalization outperforms group penalization.
The rest of the paper is organized as follows. Data and model settings are described in Section 2. For integrity of this paper, we also briefly describe the homogeneity model and an existing penalization approach designed to reinforce the model assumption. The heterogeneity model and sparse group penalization approach are described in Section 3. The proposed estimate can be computed using an effective group coordinate descent (GCD) algorithm. Numerical studies, including simulation and data analysis, are conducted in Section 4. The paper concludes with the discussion in Section 5. Some additional analysis results are presented in the Appendix. In study m, denote T m as the logarithm (or another known monotone transformation) of the failure time. Denote X m as the length-d vector of gene expressions. For simplicity of notation, assume that the same set of genes is measured in all M studies. For the ith subject, the AFT model assumes that

Integrative analysis of cancer prognosis studies
where b m 0 is the intercept, b m is the length-d vector of regression coefficients, and " m i is the error term. When T m i is subject to right censoring, we observe is the logarithm of censoring time, and d m i =I{T m i fC m i } is the event indicator.
When the distribution of " m i is known, a parametric likelihood function can be constructed. Here we consider the more flexible case where this distribution is unknown (i.e., a semiparametric model). In the literature, multiple estimation approaches have been developed, including for example the Buckley-James and rank-based approaches. In this study, we adopt the weighted least-squares approach (Stute, 1996), which has the lowest computational cost. This property is especially desirable with high-dimensional data.
Let F m be the Kaplan-Meier estimator of the distribution function F m of T m . F m can be written as F m (y)=g , i=2, . . . , n m :

Penalized integrative prognosis analysis
Here Y m (1) f Á Á Á fY m (n m ) are the order statistics of Y m i s, and d m (1) , . . . , d m (n m ) are the associated event indicators. Similarly, let X m (1) , . . . , X m (n m ) be the associated gene expressions of the ordered Y m i s. Stute (1996) proposed the weighted least-squares estimator , respectively. With the weighted centred values, the intercept is zero. The weighted least-squares objective function can be written as Consider the overall objective function L(b)= g M m=1 L m (b m ). With this objective function, larger datasets have more contributions. When desirable, normalization by sample size can be applied.
That is, the M datasets have the same set of cancer-associated genes. This is a sensible model when multiple datasets have been generated under the same protocol.
For marker selection, Ma et al. (2011 b) propose using the group minimax concave penalty (GMCP) approach, where the estimate is defined aŝ (x=lc))+dx is the MCP penalty (Zhang, 2010). b j S j =kS 1=2 j b j k 2 , jj Á jj 2 is the L 2 norm, S j =n x1 X j kX j , and X j is the nrd j submatrix of X that corresponds to b j . In (4), d j is the size of the coefficient group corresponding to gene j. When the M datasets have matched gene sets, d j wM. We keep d j so that this formulation can easily accommodate partially matched gene sets. When gene j is not measured in dataset k, we take the convention b k j 0. l 1 is a tuning parameter, and c is a regularization parameter.
In this analysis, genes are the functional units. The overall penalty is the sum over d individual penalties, with one for each gene. For gene selection, MCP penalization is adopted. For a specific gene, its effects in the M studies are represented by a ' group ' of M regression coefficients. Under the homogeneity model, the M studies are expected to identify the same set of genes. Thus, within a group, no further selection is needed, and hence the L 2 norm is adopted. Note that here we adopt the jj Á jj S j norm, which rescales the regression coefficient vector by covariance matrix S j , so that computation is less ad hoc. This differs from the approach in Ma et al. (2011 b).

Heterogeneity model and penalized selection
(i) Heterogeneity model When multiple datasets are generated in independent studies, heterogeneity inevitably exists (Knudsen, 2006). The degree of heterogeneity depends on differences in study protocols, profiling techniques and many other factors. In cancer prognosis studies, the effort to unify the sets of identified markers across independent studies has not been very successful (Knudsen, 2006 ;Cheang et al., 2008). This can also be partly seen from Ma et al. (2011 b) and our data analysis in Section 4. Such observations raise the question whether the homogeneity model is too restricted and motivates the heterogeneity model. Under the heterogeneity model, one gene can be associated with prognosis in some studies but not others. This model includes the homogeneity model as a special case and can be more flexible.
There are also scenarios under which the homogeneity model is conceptually not sensible, but the heterogeneity model is. The first is where different studies are on different types of cancers (Ma et al., 2009). On the one hand, different cancers usually have different sets of markers. On the other hand, multiple pathways, such as apoptosis, cell cycle and signalling, are associated with prognosis of multiple cancers. The second scenario is the analysis of different subtypes of the same cancer. Different subtypes have different risks of occurrence and progression patterns, and it is not sensible to reinforce the same genomic basis. The third scenario is where subjects in different studies have different demographic measurements, clinical risk factors, environmental exposures and treatment regimens. For genes not intervened with those additional variables, their importance remains consistent across multiple studies. However, for other genes, they may be important in some studies but not others.
(ii) Penalized selection Here, notations have similar definitions as in Section (ii). l 2 is another tuning parameter. Under the heterogeneity model, marker selection needs to be achieved in two dimensions. The first is to determine whether a gene is associated with prognosis in any dataset at all. This step of selection is achieved using a GMCP penalty. For an important gene, the second dimension is to determine in which dataset(s) it is associated with prognosis. This step of selection is achieved using an MCP penalty. The sum of two penalties can achieve two-dimensional selection and fully recover the associations between multiple genes and prognosis in multiple studies.
In the single-dataset analysis with a continuous response and linear regression model, Friedman et al. (2010) proposed a sparse group Lasso (SGLasso) penalty. The penalty in (5) has been partly motivated by that study. Owing to the analogy, we refer to penalty (5) as 'sparse group MCP ' or SGMCP. Unlike in Friedman et al. (2010), multiple heterogeneous datasets are analysed in this study, and a group corresponds to one gene in multiple datasets, as opposed to multiple variables. In addition, the MCP penalty is adopted to replace the Lasso penalty, as in singledataset analysis, MCP has been shown to outperform Lasso (Zhang, 2010;Huang et al., 2012 a, b). In addition, censored survival data under the AFT model are analysed, which differ from the continuous response and linear model in Friedman et al. (2010). Note that MCP becomes Lasso when c='. Thus, SGMCP includes SGLasso as a special case.

(iii) Computation
For dataset m(=1, …, M), we standardize the gene expressions to have marginal means zero and satisfy X m j kX m j =n. Thus, n x1 X j kX j =I d j , where I d j is the d j rd j identity matrix. With slight abuse of notation, we denote jjb j jj for jjb j jj I d j . For computation, we consider a GCD algorithm. It optimizes the objective function with respect to one gene at a time, and iteratively cycles through all genes. The overall cycling is repeated multiple times until convergence. Consider the overall objective functioñ For j=1, …, d, given b k (klj) fixed at their current estimatesb b (s) k , we seek to minimizeL L(b, l 1 , l 2 , c) with respect to b j . Here, only the terms involving b j iñ L L(b, l 1 , l 2 , c) matter. This is equivalent to minimizing where r xj =Yxg klj X jb b (s) k . The first-order derivative of (7) is : where t= sgn(b 1 j ) With normalization, n x1 X k j X j =I d j . By setting expression (8) to be zero, we have : where z j =n x1 X j k r xj = (z 1 j , . . . , z M j )k, Penalized integrative prognosis analysis We first take g fixed at the current estimate of b j . The kth element in (9) can be rewritten as : The solution to equation (10) is , S 1 =sgn(u)(jujxl) + : For k=1, …M, we set u k = d gb k j gb k j and u=(u 1 , . . . , u M )k. Taking u back into its definition, we have Solving equation (11) where S 2 (u, t)= 1x t u k k + u: Solving equation (9) amounts to iteratively calculating equation (10) j . This is achieved by solving equation (9)-(12) iteratively until convergence. In our numerical study, convergence is achieved for all datasets within ten iterations; ; 3. Update sqs+1; 4. Repeat Steps 2-3 until convergence.
This algorithm starts with a null model. In each iteration, it cycles through all d genes. For each gene, it only involves simple calculations, and the update can be accomplished easily. We use the L 2 norm of the difference between two consecutive estimates smaller than 10 x5 as the convergence criterion. Convergence is achieved with all our simulated and real data.

(iv) Tuning parameter selection
With the MCP penalty, there is one tuning parameter l and one regularization parameter c. Generally speaking, smaller values of c are better at retaining the unbiasedness of the MCP penalty for large coefficients, but they also have the risk of creating objective functions with a non-convexity problem that are difficult to optimize and yield solutions that are discontinuous with respect to l. It is therefore advisable to choose a c value that is big enough to avoid this problem but not too big (Zhang, 2010 ;Breheny & Huang, 2011). Such results can be extended to the SGMCP.
In principle, it is possible to have different c values for the two penalties. In practice, adopting the same c value may reduce computational cost without affecting performance of the estimate much. In numerical study, we have experimented with a few values for c, particularly including 1 . 8, 3, 6 and 10 and ', as suggested in Zhang (2010) and Breheny & Huang (2011). l 1 and l 2 control shrinkage at the group and individual-coefficient levels, respectively. Let l 2max be the largest l 2 for which all regression coefficients are shrunk to zero. Let l 1max (l 2 ) be the largest l 1 under a fixed l 2 for which all regression coefficients are shrunk to zero. It can be shown that l 2 max = n x1 max j=1, ..., d, k=1, ..., M jX k j Y k j, where X k j is the component of X j that corresponds to the kth dataset. l 1 max (l 2 )= max j=1, ..., d S 1 (n x1 X j kY, l 2 ) = ffiffiffiffi d j p . This result puts constraints on the range of tuning parameters.
In this study, we search for c, l 1 and l 2 values jointly using V-fold cross-validation (V=5 in numerical study). This is computationally feasible as the GCD algorithm has low computational cost. It is expected that l 1 , l 2 cannot go down to very small values that correspond to regions not locally convex. The cross-validation criteria in such non-locally convex regions may not be monotonous. To avoid such regions, we select l 1 and l 2 where the cross-validation criterion first goes up.
For a simulated dataset, we analyse using SGMCP as well as the following alternatives : (i) MCP and Lasso (which is MCP with c='). Here, we analyse the three datasets separately and then combine analysis results across datasets. This is a meta-analysis approach. (ii) GMCP and GLasso (which is GMCP with c='). These two approaches reinforce that the same set of genes are identified across datasets, are sensible under the homogeneity model, but can be too restricted under the heterogeneity model ; and (iii) SGLasso (which is SGMCP with c='). With the above approaches, we consider three different c values for the MCP and GMCP penalties. For all approaches, tuning parameters are chosen using 5-fold cross-validation.
Summary statistics, including the number of true positives, number of false positives, and their standard deviations based on 100 replicates, are shown in Table 1 (homogeneity model, homogeneous errors) and 2 (heterogeneity model, homogeneous errors) and 5 (Appendix, homogeneity model, heterogeneous errors) and 6 (Appendix, heterogeneity model, heterogeneous errors), respectively. We observe that SGMCP has the best performance under all simulated scenarios. It can significantly outperform GMCP under the homogeneity model. SGMCP is able to identify most or all of the true positives. Under some simulated scenarios, there are a number of false positives. This is partly caused by correlations among genes. The false positive rate increases as genes become 'more correlated '. When c increases, then the number of false positives increases. This is as expected, as in single-dataset analysis, it has been suggested that Lasso-type penalization (c=') tends to overselect (Huang et al., 2012a). GMCP and MCP may identify a satisfactory number of true positives, however, at the price of a much larger number of false positives (Table 2).

(ii) Analysis of breast cancer studies
Worldwide, breast cancer is the commonest cancer death among women. In 2008, breast cancer caused 458 503 deaths worldwide (13 . 7 % of cancer deaths in women). Multiple high-throughput profiling studies have been conducted, searching for genomic markers associated with breast cancer prognosis. We collect data from three independent breast cancer prognosis studies, which were originally reported in Huang et al. (2003), Sotiriou et al. (2003) and van't Veer et al. (2002), respectively, and referred to as D1, D2 and D3 hereafter. In D1, Affymetrix chips were used to profile the expressions of 12 625 genes. There were a total of 71 subjects, among whom 35 died during follow-up. The median follow-up was 39 months. In D2, cDNA chips were used to profile the expressions of 7650 genes. There were a total of 98 subjects, among whom 45 died during follow-up. The median follow-up was 67 . 9 months. In D3, oligonucleotide arrays were used to profile 24 481 genes. There were a total of 78 patients, among whom 34 died during follow-up. The median follow-up was 64 . 2 months. We process each dataset separately as follows. With Affymetrix data, a floor and a ceiling are added, and then measurements are log 2 transformed. With both Affymetrix and cDNA data, we fill in missing expressions with means across samples. We then standardize each gene expression to have zero mean and unit variance.
In previous studies such as Ma et al. (2011 b), the homogeneity model is assumed. As the three datasets were generated in three independent studies, heterogeneity is expected to exist across datasets. This can be partly seen from the summary survival data, the profiling protocols, as well as results from analysis of each individual dataset using MCP (Table 7). We analyse the three datasets using the proposed approach as well as MCP, GMCP and SGLasso. The identified genes and corresponding estimates are shown in Tables 3 and 7-9 in the Appendix. Note that with all approaches, the small magnitudes of regression coefficients are caused by 'clustered ' log survival times. SGMCP identifies seven, eight and five genes for D1-D3, respectively. Estimates in Table 3 suggest that different datasets may have different prognosis-associated genes. This partly explains why it fails to unify the identified markers across different breast cancer prognosis studies (Cheang et al., 2008). As described in Section 1, multiple factors may contribute to this heterogeneity. Without having access to all the experimental details, we are not able to determine the exact cause of heterogeneity. Although there are overlaps, different approaches identify significantly different sets of genes. This observation is not surprising given the simulation results. We mine the published literature and find that genes identified by SGMCP may have important biological implications. More details are provided in the Appendix. Of note, as there is no objective way to determine which set of genes are 'biologically more meaningful ', we do not discuss the biological implications of genes identified by the alternative methods.
To provide a more comprehensive description of the three datasets and various approaches, we also conduct evaluation of prediction performance. Although, in principle, marker identification and prediction are two distinct objectives, evaluation of prediction performance can be informative for marker identification. In particular, if prediction is more accurate, then the identified markers are expected to be more meaningful. For prediction evaluation, we adopt a random sampling approach as in Ma et al. (2009). More specifically, we generate training sets and corresponding testing sets by random splitting data with sizes 3 : 1. Estimates are generated using the training sets only. We then make prediction for subjects in the testing sets. We dichotomize the predicted linear risk scores Xkb b at the median, create two risk groups, and compute the logrank statistics, which measure the difference in survival between the two groups. To avoid extreme splits, this procedure is repeated 100 times. The average logrank statistics are calculated as 7 . 30 (SGMCP), 3 . 43 (MCP), 2 . 77 (GMCP) and 3 . 33 (SGLasso). SGMCP is the only approach that can separate subjects into groups with significantly different survival risks (P-value=0 . 007). Such a result suggests that allowing for heterogeneity in markers across datasets can lead to better prediction performance.

(iii) Analysis of lung cancer studies
Lung cancer is the leading cause of death from cancer for both men and women in USA and in most other parts of the world. Non-small-cell lung cancer (NSCLC) is the most common cause of lung cancer death, accounting for up to 85 % of such deaths (Tsuboi et al., 2007). Gene-profiling studies have been extensively conducted on lung cancer, searching for markers associated with prognosis. Data were collected from three independent experiments (Shedden et al., 2008 ;Jeong et al., 2010 ;Xie et al., 2011). The UM (University of Michigan Cancer Center) study had a total of 175 patients, among whom 102 died during follow-up. The median follow-up was 53 months. The HLM (Moffitt Cancer Center) study had a total of 79 subjects, among whom 60 died during follow-up. The median follow-up was 39 months. The CAN/DF (Dana-Farber Cancer Institute) study had a total of 82 patients, among whom 35 died during follow-up. The median follow-up was 51 months. 22 283 genes were profiled in all three studies.
In studies such as Xie et al. (2011), the three datasets were combined and analysed. Such a strategy corresponds to the homogeneity model. Here, we analyse using SGMCP (Table 4), MCP (Table 10), GMCP (Table 11) and SGLasso (Table 12). The observed estimation patterns are similar to those for the breast cancer data. The SGMCP estimates again suggest that different datasets may have different prognosis-associated genes. In the Appendix, we show that the SGMCP identified genes may have important biological implications. We also evaluate prediction performance using the resampling approach described  in the last subsection. The average logrank statistics are calculated as 10 . 93 (SGMCP), 2 . 83 (MCP), 1 . 73 (GMCP) and 3 . 07 (SGLasso), respectively. The observation is again similar to that for the breast cancer data.

Discussion
In cancer genomic research, multi-dataset analysis provides an effective way to overcome certain drawbacks of single-dataset analysis. In most published studies, it has been reinforced that multiple datasets share the same set of prognosis-associated genes, that is, the homogeneity model. In this study, for multiple cancer prognosis datasets, we consider the heterogeneity model, which includes the homogeneity model as a special case and is less restricted. This model may provide a way to accommodate the failure to unify cancer prognosis markers across independent studies (Knudsen, 2006;Cheang, 2008). Under the heterogeneity model, we use sparse group penalization for marker identification. This penalization approach is intuitively reasonable and computationally feasible. Simulation study shows that it has satisfactory performance. It is interesting to note that SGMCP outperforms GMCP even under the homogeneity model. Thus, in practical data analysis, the heterogeneity model and SGMCP can be a 'safer ' choice. With the breast cancer and lung cancer prognosis datasets, existing analyses assume the same set of markers across datasets. Our analysis suggests that it may be more reasonable to allow for different sets of markers. Under the heterogeneity model, marker identification needs to be conducted in two dimensions. Beyond sparse group penalization, composite penalization and multi-step approaches may also be able to achieve such identification. In this paper, we focus on sparse group penalization. Comprehensive investigation and comparison of different approaches are beyond the scope of this paper. The proposed approach is based on the MCP penalty, which has been shown to have satisfactory performance in single-dataset analysis. We suspect that it is possible to develop similar approaches based on, for example, bridge and SCAD penalties. As in single-dataset analysis there is no evidence that such penalties are superior to MCP, such a development is not pursued. In our numerical study, we observe superior performance of the proposed approach. A limitation of this study is lack of theoretical support. Pursuit of theoretical properties is postponed to future studies.