BACKGROUND
gene-set analysis is a microarray data analysis method that uses existing knowledge of biological pathways, or sets of individual genes that are linked via related biological functions.
gene-set analysis mainly aims to discover gene sets for which expression is associated with a phenotype of interest.
compared to single-gene analyses, gene-set analyses may lead to more interpretable results by yielding insights into biological mechanisms.
furthermore, considering gene sets rather than single genes can reduce problems associated with multiple testing, since there are typically far fewer gene pathways than individual genes.
initially, gene-set analyses focused on identifying biological pathways  that are differentially expressed between two classes of a phenotype such as tumor vs. normal cells.
mootha et al.
<cit>  proposed gene set enrichment analysis , based on the kolmogorov-smirnov statistic, which measures the maximum degree of differential gene expression in a gene set across a binary phenotype.
subramanian et al.
<cit>  improved gsea by weighting each gene according to its correlation with the phenotype, and calculating a running-sum statistic; in contrast, the original gsea uses equal weights regardless of the correlation between genes and the phenotype.
gsea calculates a p-value by permuting the original data set, which is computationally intensive for a large dataset.
kim and volsky  <cit>  proposed a parametric analysis of gene-set enrichment, which calculates a z-score for a given gene set from a parameter such as fold change value between two classes of a phenotype, and makes statistical inference from the asymptotic normal distribution of the z score.
dinu et al.
<cit>  described some critical problems with gsea, and proposed an alternative method by extending an individual gene analysis method, significance analysis of microarrays, to gene-set analysis .
in addition, they compared sam-gs to gsea using a mouse microarray dataset with simulated gene sets, and showed an advantage to sam-gs over gsea in the analysis of three real microarray datasets.
recently, gene-set analysis methods have been expanded to include other kinds of phenotypes, such as censored survival time and quantitative traits.
goeman et al.
<cit>  proposed the global test, a score statistic based on random-effects modeling of parameters corresponding to the coefficients of the individual genes in the pathway.
goeman et al.
<cit>  extended the global test to assess association of a set of genes with survival time based on a cox proportional hazards model.
in addition, binder and schumacher  <cit>  considered fitting high-dimensional survival models while allowing for mandatory clinical covariates using a boosting algorithm.
boulesteix and hothorn  <cit>  instead proposed the global boost test, which combines a cox model for modeling the clinical covariates with a boosting algorithm for modeling the additional predictive value of high-dimensional gene expression data.
furthermore, adewale et al.
<cit>  proposed a unified general analysis method for microarray data for identifying pathways whose expression is associated with a phenotype of any kind, adjusting for covariates that may also be associated with the phenotype of interest.
this unified pathway analysis method combines the regression-based test statistic for each individual gene in a pathway of interest into a real pathway-level test statistic.
the form of the test statistic is a sum of squares of the wald statistic for individual genes in the pathway of interest.
gene expression profiles have been used extensively in the prediction of tumor subtypes or patient survival .
initially, many studies focused on expression levels of single genes to predict tumor subtype or patient survival.
among many thousands of microarray measurements, each relating to the expression level of a single gene, a subset of significant genes can be identified by constructing a prediction model using the lasso , principal components analysis , supervised principal component analysis , support vector machines , and other methods.
however, single genes are often not of primary interest because the activities of entire pathways or genomic regions that are suspected to be more biologically relevant.
in addition, combining gene expression data with prior biological knowledge of groups of genes improves prediction accuracy and interpretability of survival models.
the prediction of patient survival may be improved by integrating gene expression data with prior biological knowledge such as gene sets and pathways, as well as by adjusting for covariates such as age, sex, and other clinical variables.
chen and wang  <cit>  proposed a general strategy for improving prediction accuracy and interpretability by constructing pathway-based prediction models for survival, which outperform the prediction models based on expression levels of single genes.
recently, liu et al.
<cit>  compared the statistical performance of three gene-set analysis methods - the global test, ancova global test, and sam-gs - for a binary phenotype based on simulated data and real microarray datasets.
they reported similar performances for all three methods after appropriate standardization, given the use of permutation-based inference.
they also showed the advantage of the global test and ancova global test, which are able to analyze survival phenotypes and adjust for covariates.
to our knowledge, however, the performances of gene-set analysis methods with a survival phenotype have rarely been studied via simulation.
in this paper, we compare the performances of five gene-set analysis tests: two different gsea tests  <cit> , global test   <cit> , wald-type test   <cit>  and global boost test   <cit>  for assessing differential expression associated with the survival phenotype based on a simulation dataset and a real dataset of ovarian cancer patients.
RESULTS
simulation experiment
we generated a simulation dataset by the following procedure to evaluate the performance of five gene-set analysis tests:
we randomly generated observations from a multivariate normal with a zero-mean vector and a variance-covariance matrix ∑, denoted by mvn.
we randomly generated a vector of regression coefficients, β, from either a uniform distribution or a normal distribution.
this represents the association between survival and gene expression.
using the observations generated in  and the vector of regression coefficients β generated in , we constructed a survival time from a cox model with a specified baseline hazard function.
censoring times were generated from an exponential distribution with a parameter λ. the parameter λ was determined by the censoring fraction.
in the simulation study, we considered various parameters: the total number of genes , the sample size , the fraction of censoring , the size of a gene set of interest , and the proportion of significant genes in the gene set .
to check the size of the five tests, we randomly generated gene expression variables from mvn, ∑ =  <dig> ip, where ip is an identity matrix of dimension p × p, and constructed survival times from a cox model with βj =  <dig>  j =  <dig> ..., p and the constant baseline hazard rate of  <dig> .
for the power calculation, we also randomly generated gene expression variables from mvn and constructed survival times from a cox model with the parameter β and the baseline hazard function has an exponential distribution with the hazard rate of  <dig> .
here ∑ and β were specified according to the following scenarios.
first, we considered four different correlation structures of gene expression variables as implemented in liu et al.
<cit>  and jung et al.
<cit> .
case  is that all gene expressions are independent, which assumes the correlation matrix as ∑ =  with σii =  <dig>  for i =  <dig> ..., p; σij =  <dig> for i ≠ j, with i, j =  <dig> ..., p. case  is that only significant genes are correlated within a gene set but non-significant genes are independent, that is, ∑ =  with σii =  <dig>  for i =  <dig> ..., p; σij =  <dig>  if two significant genes fall into the same gene set and σij =  <dig> otherwise.
case  is that there is an autoregressive correlation between significant genes, that is, ∑ =  with σii =  <dig>  for i =  <dig> ..., p; σij =  <dig>  ×  <dig> |i-j| if two significant genes fall into the same gene set and σij =  <dig> otherwise.
case  is that there is an unstructured correlation between significant genes, that is, ∑ =  with σii =  <dig>  for i =  <dig> ..., p; σij =  <dig>  × ρij, if two significant genes fall into the same gene set and σij =  <dig> otherwise, where ρij is a random variable generated from n. secondly, we considered two different ways of generating the regression coefficient, βj, for j =  <dig> ...,, to investigate how the association of survival with genes affects the power for detecting the significant gene sets.
here  represents the greatest integer less than or equal to a. case  is that the survival is positively associated with genes by generating only positive coefficients from a uniform distribution u. case  is that survival is randomly associated with genes by generating the coefficient from n. for the rest of the regression coefficients, we set βj =  <dig> for j = + <dig> ..., p.
to assess the size and power of the tests, we implemented the simulation procedure as follows:  we calculated the test statistic for each method,  permuted the samples  <dig> times, recalculated the test statistics, and used these permuted test statistics to estimate the p-value,  and then repeated procedures  and   <dig> times to estimate the size and  <dig> times to estimate the power, respectively.
the size is estimated as the observed proportion of replications with a p-value smaller than the nominal size α =  <dig> , and the power is estimated as the observed proportion of replications of in which the null hypothesis was correctly rejected at the nominal size α =  <dig> .
note: n = sample size; m = the size of gene set; cp = censoring proportion.
gsea <dig> = gsea with the equal weight; gsea <dig> = weighted gsea; gt = global test; wt = wald-type test; gbst = global boost test
note: m = the size of gene set; cp = censoring proportion; mp = the proportion of significant genes in a gene set.
gsea <dig> = gsea ith the equal weight; gsea <dig> = weighted gsea; gt = global test; wt = wald-type test; gbst = global boost test.
for case , all genes are independent; for case , there is within-correlation among significant genes; for case , there is an autoregressive correlation between significant genes; and for case , there is an unstructured correlation between significant genes
the results of tables  <dig> and  <dig> are depicted in figures  <dig> and  <dig> under the censoring fractions of cp =  <dig>  and  <dig> , respectively.
in each plot, the ten different lines represent the power of the five tests under the two different scenarios of cases  and .
the solid line represents the power of the five tests for case  and the dotted line for case .
the power is plotted with mp, the proportion of the significant genes in each gene set, set to be  <dig> , <dig> , <dig>  and  <dig> .
here mp is considered to be the effect size because the proportion of significant genes affects the association between the survival time and the gene set.
it is interesting that the power of the five tests differs depending on the combinations of the gene correlation structure and the association of genes with survival.
for case , the power of the five tests is higher for case  than for case .
on the other hand, in the plots for cases  and , the power of the five tests is much higher for case  than for case .
this result implies that there might be a synergistic effect in the power of detecting significant genes when the genes are correlated and the survival is positively associated with genes.
for case , the power of gt, wt and gbst is almost the same regardless of the association of genes with survival, whereas both gsea <dig> and gsea <dig> have higher power for case  than for case .
in general, as described in figure  <dig>  gt, wt and gbst consistently have higher power than gsea <dig> and gsea <dig>  comparing the two gsea tests, the power of gsea <dig> is equal to or slightly greater than the power of gsea <dig>  but these two tests have lower power than  <dig>  except for the combinations of cases  and  with case  when mp ≥  <dig> .
as shown in figure  <dig>  the same pattern of power is found when the censoring fraction increases from cp =  <dig>  to  <dig>  though power consistently decreases as censoring increases.
analysis of a real example of ovarian cancer data
we next evaluated the performance of the five tests using an ovarian cancer data set from dressman et al.
<cit> .
this dataset consists of  <dig> ovarian cancer samples that were obtained at the initial cytoreductive surgery from patients treated at duke university medical center and h. lee moffitt cancer center and research institute.
<dig> gene expression levels were used in this analysis, of which  <dig> pathways were identified by kegg.
among these  <dig> pathways, the smallest pathway consists of  <dig> genes while the largest one includes  <dig> genes, and the average number of genes across the  <dig> pathways is about  <dig>
crijns et al.
<cit>  identified survival-related profile, pathways and transcription factors by analyzing a dataset of  <dig> patients with advanced-stage serious ovarian cancer from the university medical center groningen in netherlands, collected in the time period 1990- <dig>  they used the dataset of  <dig> ovarian cancer samples in dressman et al.
<cit>  to validate the identified profile and pathways.
we compared the pathways identified by the five tests with those reported in dressman et al.
<cit>  and crijns et al.
<cit> .
in dressman et al.
<cit> , the significant pathways were found by a binary logistic regression model analysis and a stochastic regression model search, called shotgun stochastic search.
on the other hand, crijns et al.
<cit>  identified the significant pathways using the functional class scoring analysis, in which the p-value of a univariate cox proportional hazards is computed for all genes, and then p-values of each pathway are summarized by the mean negative logarithm of single gene p-values , as well as the kolmogorov-smirnov  statistic for testing if the p-values come from a uniform distribution.
the significance of each pathway is evaluated by computing the empirical distribution of these summary statistics in random samples of genes.
note: the underlined pathways are identified by crijns et al..
discussion
many gene-set analysis proposals have been compared in simulation studies and on real data sets.
however, most studies have dealt with a binary phenotype like the presence or absence of disease, or treatment vs. control.
in this paper, we focused on the survival phenotype and compared five different gene-set analysis tests, gsea <dig>  gsea <dig>  gt, wt and gbst, in a simulation study and on two ovarian cancer data sets.
from the simulation results, we found that gt, wt, and gbst are more powerful than gsea <dig> and gsea <dig>  furthermore, the power of the five tests is substantially affected by the correlation structure of genes and the association between survival and the genes.
the power of the five tests can approach  <dig>  when the genes are correlated, survival is positively associated with the gene expression values, and mp ≥  <dig> .
for more powerful test results, it might be desirable to check the correlation structure among genes within each pathway and the direction of association between genes and survival.
although gsea was originally based on the rank of the correlation coefficients between a binary phenotype and gene expression levels, we replaced the correlation coefficient by the regression coefficient from a cox model and took the rank of the absolute value of its standardized coefficient.
we then followed the rest of the gsea procedure in order to identify significant gene-sets or pathways.
this is a modification of gsea that takes into account the survival phenotype, which may in some cases be more informative than a binary phenotype.
in addition, the modified gsea allows covariance adjustment for age, sex and other clinical variables by taking the regression coefficient of genes from the cox model with those adjusting covariates.
however, the performance of the gsea tests, gsea <dig> and gsea <dig>  are not satisfactory except for a few cases.
this may be due to the fact that the gsea tests are nonparametric approaches using the rank-based statistic instead of using the value of the regression coefficient.
on the other hand, gt, wt and gbst have been proposed for regression-based models and can be extended to any phenotype such as binary, continuous, multi-class, or survival.
these three tests can be easily adjusted for covariates, in order to determine whether gene expression profiles have an association with survival beyond what is explained by the adjusting covariates.
the gt method assumes a random-effect model for the parameters corresponding to the coefficients of the individual genes in the pathway, in which the parameters are random variables and samples from n. here all parameters are assumed to have a common variance of τ <dig>  which can be extended to have a more complex covariance structure as mentioned in goeman et al.
<cit> .
the gt method uses a score test of τ <dig> =  <dig>  which is equivalent to the null hypothesis that there is no association between survival and a given set of genes.
therefore, this test does not depend on the number of genes in a given set and works under the assumption of a random-effect model with a common variance.
as discussed in goeman et al.
<cit> , the score test of gt can have the optimal power against alternatives with small values of parameter τ <dig>  in the simulation study, we generated the parameter from either a uniform or normal distribution with small variances.
this may have contributed to the high power that we observed for gt.
on the other hand, boulesteix and hothorn  <cit>  proposed the gbst for testing the additional predictive value of gene expression data while adjusting for clinical covariates.
they combined the standard regression models such as a logistic regression or a cox model with a boosting procedure and used a permutation-based testing scheme.
while goeman et al.
<cit>  reduced the multi-dimensional gene profiles into one-dimensional variance using a random-effect model, adewale et al.
<cit>  instead computed the sum of squares of the wald statistics for the genes in the pathway from the univariate regression model.
likewise, the wt does not depend on the multi-dimensionality of the gene sets since the test statistic can be summarized as the sum of individual association measures.
compared to the gsea tests, these three tests are based on parametric approaches since the values of the regression coefficients are taken into account in these tests via a score statistic or a sum of squares.
as pointed out by a referee, methodological issues should be considered when comparing gene-set analysis methods.
goeman and bühlmann  <cit>  addressed the definition of the null hypothesis, and in particular competitive versus self-contained tests, as well as the calculation of p-values, and the use of gene sampling versus subject sampling methods.
gt, wt and gbst are based on the classical statistical models which lead to test a self-contained null hypothesis.
furthermore, they calculate the p-value using a subject sampling model.
on the other hand, gsea is a hybrid method: a kolmogorov-smirnov test statistic is motivated by a gene-sampling model, whereas a subject-sampling model is used to calculate the p-value.
it was also pointed out that gsea sometimes shows low power since the model and null hypothesis used to motivate the test statistic are different from the model and null hypothesis used to calculate the p-value.
simulation results indicate that both gsea <dig> and gsea <dig> have lower power than gt, wt and gbst.
one of the important advantages of gene-set analysis over the single gene approach is that the multiple testing problems can be alleviated because the number of pathways is much smaller than that of genes.
in the ovarian cancer example, for example, there are  <dig> genes.
the number of pathways is  <dig>  which is more than a 100-fold reduction.
when a single pathway is of interest in pathway analysis, the issue of multiple testing is not a problem.
however, when multiple pathways are of interest, as in cases when pathway analysis is exploratory with many pathways of potential interest, then multiple comparison methods such as false-discovery rates  must be applied.
in the ovarian cancer example, we used fdr to control  <dig> pathway comparisons.
a few pathways were found to have q <  <dig> .
it is noted that all  <dig> pathways identified by crijns et al.
<cit>  contain more than  <dig> genes.
their sizes vary from  <dig> to  <dig>  and their median is  <dig>  however, gt detects  <dig> pathways whose sizes range from  <dig> to  <dig> with a median of  <dig>  wt detects  <dig> pathways whose sizes range from  <dig> to  <dig> with a median of  <dig>  and gbst detects  <dig> pathways whose sizes range from  <dig> to  <dig> with a median of  <dig>  in summary, gt, wt and gbst detected pathways with a variety of gene-set sizes, whereas crijns et al.
<cit>  detected only large gene sets.
to implement gt and gbst, we adapted the function 'gt' in the r package 'globaltest', and the function 'globalboosttest' in the r package 'globalboosttest'.
we implemented gsea <dig>  gsea <dig>  and wt from scratch in r, using the function 'coxph'.
for the readers' convenience, all codes are available as additional files.
the additional file  <dig> is r source codes for generating data, calling subroutines, and calculating permuted p-values and the additional file  <dig> is r source codes for calculating five test statistics, gsea <dig>  gsea <dig>  gt, wt, and gbst.
CONCLUSIONS
in conclusion, gt, wt and gbst have high power to identify gene sets whose expression is associated with survival.
survival is often strongly affected by well-known clinical variables such as stage, type of tumor, and tumor size, as well as demographic variables such as age, sex, and race.
therefore, it is very important to evaluate the effects of genes on survival while considering known covariates.
since these tests allow for adjustment for covariates, they assess whether expression of a given set of genes has an association with survival after controlling for confounders.
in general, gt, wt and gbst are more powerful than gsea <dig> and gsea <dig>  however, when genes are correlated within each pathway and survival is positively associated with genes, it does not matter which test is applied to detect significant gene sets, provided that the number of significant genes is moderately large within each pathway.

 gene-set analysis is a microarray data analysis method that uses existing knowledge of biological pathways, or sets of individual genes that are linked via related biological functions. gene-set analysis mainly aims to discover gene sets for which expression is associated with a phenotype of interest. compared to single-gene analyses, gene-set analyses may lead to more interpretable results by yielding insights into biological mechanisms. furthermore, considering gene sets rather than single genes can reduce problems associated with multiple testing, since there are typically far fewer gene pathways than individual genes. initially, gene-set analyses focused on identifying biological pathways that are differentially expressed between two classes of a phenotype such as tumor vs. normal cells. mootha et al. proposed gene set enrichment analysis, based on the kolmogorov-smirnov statistic, which measures the maximum degree of differential gene expression in a gene set across a binary phenotype. subramanian et al. improved gsea by weighting each gene according to its correlation with the phenotype, and calculating a running-sum statistic; in contrast, the original gsea uses equal weights regardless of the correlation between genes and the phenotype. gsea calculates a p-value by permuting the original data set, which is computationally intensive for a large dataset. kim and volsky proposed a parametric analysis of gene-set enrichment, which calculates a z-score for a given gene set from a parameter such as fold change value between two classes of a phenotype, and makes statistical inference from the asymptotic normal distribution of the z score. dinu et al. described some critical problems with gsea, and proposed an alternative method by extending an individual gene analysis method, significance analysis of microarrays, to gene-set analysis. in addition, they compared sam-gs to gsea using a mouse microarray dataset with simulated gene sets, and showed an advantage to sam-gs over gsea in the analysis of three real microarray datasets. recently, gene-set analysis methods have been expanded to include other kinds of phenotypes, such as censored survival time and quantitative traits. goeman et al. proposed the global test, a score statistic based on random-effects modeling of parameters corresponding to the coefficients of the individual genes in the pathway. goeman et al. extended the global test to assess association of a set of genes with survival time based on a cox proportional hazards model. in addition, binder and schumacher considered fitting high-dimensional survival models while allowing for mandatory clinical covariates using a boosting algorithm. boulesteix and hothorn instead proposed the global boost test, which combines a cox model for modeling the clinical covariates with a boosting algorithm for modeling the additional predictive value of high-dimensional gene expression data. furthermore, adewale et al. proposed a unified general analysis method for microarray data for identifying pathways whose expression is associated with a phenotype of any kind, adjusting for covariates that may also be associated with the phenotype of interest. this unified pathway analysis method combines the regression-based test statistic for each individual gene in a pathway of interest into a real pathway-level test statistic. the form of the test statistic is a sum of squares of the wald statistic for individual genes in the pathway of interest. gene expression profiles have been used extensively in the prediction of tumor subtypes or patient survival. initially, many studies focused on expression levels of single genes to predict tumor subtype or patient survival. among many thousands of microarray measurements, each relating to the expression level of a single gene, a subset of significant genes can be identified by constructing a prediction model using the lasso, principal components analysis, supervised principal component analysis, support vector machines, and other methods. however, single genes are often not of primary interest because the activities of entire pathways or genomic regions that are suspected to be more biologically relevant. in addition, combining gene expression data with prior biological knowledge of groups of genes improves prediction accuracy and interpretability of survival models. the prediction of patient survival may be improved by integrating gene expression data with prior biological knowledge such as gene sets and pathways, as well as by adjusting for covariates such as age, sex, and other clinical variables. chen and wang proposed a general strategy for improving prediction accuracy and interpretability by constructing pathway-based prediction models for survival, which outperform the prediction models based on expression levels of single genes. recently, liu et al. compared the statistical performance of three gene-set analysis methods - the global test, ancova global test, and sam-gs - for a binary phenotype based on simulated data and real microarray datasets. they reported similar performances for all three methods after appropriate standardization, given the use of permutation-based inference. they also showed the advantage of the global test and ancova global test, which are able to analyze survival phenotypes and adjust for covariates. to our knowledge, however, the performances of gene-set analysis methods with a survival phenotype have rarely been studied via simulation. in this paper, we compare the performances of five gene-set analysis tests: two different gsea tests, global test, wald-type test and global boost test for assessing differential expression associated with the survival phenotype based on a simulation dataset and a real dataset of ovarian cancer patients.  simulation experiment we generated a simulation dataset by the following procedure to evaluate the performance of five gene-set analysis tests: we randomly generated observations from a multivariate normal with a zero-mean vector and a variance-covariance matrix ∑, denoted by mvn. we randomly generated a vector of regression coefficients, β, from either a uniform distribution or a normal distribution. this represents the association between survival and gene expression. using the observations generated in and the vector of regression coefficients β generated in, we constructed a survival time from a cox model with a specified baseline hazard function. censoring times were generated from an exponential distribution with a parameter λ. the parameter λ was determined by the censoring fraction. in the simulation study, we considered various parameters: the total number of genes, the sample size, the fraction of censoring, the size of a gene set of interest, and the proportion of significant genes in the gene set. to check the size of the five tests, we randomly generated gene expression variables from mvn, ∑ = ip, where ip is an identity matrix of dimension p × p, and constructed survival times from a cox model with βj = j =..., p and the constant baseline hazard rate of. for the power calculation, we also randomly generated gene expression variables from mvn and constructed survival times from a cox model with the parameter β and the baseline hazard function has an exponential distribution with the hazard rate of. here ∑ and β were specified according to the following scenarios. first, we considered four different correlation structures of gene expression variables as implemented in liu et al. and jung et al.. case is that all gene expressions are independent, which assumes the correlation matrix as ∑ = with σii = for i =..., p; σij = for i ≠ j, with i, j =..., p. case is that only significant genes are correlated within a gene set but non-significant genes are independent, that is, ∑ = with σii = for i =..., p; σij = if two significant genes fall into the same gene set and σij = otherwise. case is that there is an autoregressive correlation between significant genes, that is, ∑ = with σii = for i =..., p; σij = × |i-j| if two significant genes fall into the same gene set and σij = otherwise. case is that there is an unstructured correlation between significant genes, that is, ∑ = with σii = for i =..., p; σij = × ρij, if two significant genes fall into the same gene set and σij = otherwise, where ρij is a random variable generated from n. secondly, we considered two different ways of generating the regression coefficient, βj, for j =...,, to investigate how the association of survival with genes affects the power for detecting the significant gene sets. here represents the greatest integer less than or equal to a. case is that the survival is positively associated with genes by generating only positive coefficients from a uniform distribution u. case is that survival is randomly associated with genes by generating the coefficient from n. for the rest of the regression coefficients, we set βj = for j = +..., p. to assess the size and power of the tests, we implemented the simulation procedure as follows: we calculated the test statistic for each method, permuted the samples times, recalculated the test statistics, and used these permuted test statistics to estimate the p-value, and then repeated procedures and times to estimate the size and times to estimate the power, respectively. the size is estimated as the observed proportion of replications with a p-value smaller than the nominal size α =, and the power is estimated as the observed proportion of replications of in which the null hypothesis was correctly rejected at the nominal size α =. note: n = sample size; m = the size of gene set; cp = censoring proportion. gsea = gsea with the equal weight; gsea = weighted gsea; gt = global test; wt = wald-type test; gbst = global boost test note: m = the size of gene set; cp = censoring proportion; mp = the proportion of significant genes in a gene set. gsea = gsea ith the equal weight; gsea = weighted gsea; gt = global test; wt = wald-type test; gbst = global boost test. for case, all genes are independent; for case, there is within-correlation among significant genes; for case, there is an autoregressive correlation between significant genes; and for case, there is an unstructured correlation between significant genes the results of tables and are depicted in figures and under the censoring fractions of cp = and, respectively. in each plot, the ten different lines represent the power of the five tests under the two different scenarios of cases and. the solid line represents the power of the five tests for case and the dotted line for case. the power is plotted with mp, the proportion of the significant genes in each gene set, set to be,, and. here mp is considered to be the effect size because the proportion of significant genes affects the association between the survival time and the gene set. it is interesting that the power of the five tests differs depending on the combinations of the gene correlation structure and the association of genes with survival. for case, the power of the five tests is higher for case than for case. on the other hand, in the plots for cases and, the power of the five tests is much higher for case than for case. this result implies that there might be a synergistic effect in the power of detecting significant genes when the genes are correlated and the survival is positively associated with genes. for case, the power of gt, wt and gbst is almost the same regardless of the association of genes with survival, whereas both gsea and gsea have higher power for case than for case. in general, as described in figure gt, wt and gbst consistently have higher power than gsea and gsea comparing the two gsea tests, the power of gsea is equal to or slightly greater than the power of gsea but these two tests have lower power than except for the combinations of cases and with case when mp ≥. as shown in figure the same pattern of power is found when the censoring fraction increases from cp = to though power consistently decreases as censoring increases. analysis of a real example of ovarian cancer data we next evaluated the performance of the five tests using an ovarian cancer data set from dressman et al.. this dataset consists of ovarian cancer samples that were obtained at the initial cytoreductive surgery from patients treated at duke university medical center and h. lee moffitt cancer center and research institute. gene expression levels were used in this analysis, of which pathways were identified by kegg. among these pathways, the smallest pathway consists of genes while the largest one includes genes, and the average number of genes across the pathways is about crijns et al. identified survival-related profile, pathways and transcription factors by analyzing a dataset of patients with advanced-stage serious ovarian cancer from the university medical center groningen in netherlands, collected in the time period 1990- they used the dataset of ovarian cancer samples in dressman et al. to validate the identified profile and pathways. we compared the pathways identified by the five tests with those reported in dressman et al. and crijns et al.. in dressman et al., the significant pathways were found by a binary logistic regression model analysis and a stochastic regression model search, called shotgun stochastic search. on the other hand, crijns et al. identified the significant pathways using the functional class scoring analysis, in which the p-value of a univariate cox proportional hazards is computed for all genes, and then p-values of each pathway are summarized by the mean negative logarithm of single gene p-values, as well as the kolmogorov-smirnov statistic for testing if the p-values come from a uniform distribution. the significance of each pathway is evaluated by computing the empirical distribution of these summary statistics in random samples of genes. note: the underlined pathways are identified by crijns et al.. discussion many gene-set analysis proposals have been compared in simulation studies and on real data sets. however, most studies have dealt with a binary phenotype like the presence or absence of disease, or treatment vs. control. in this paper, we focused on the survival phenotype and compared five different gene-set analysis tests, gsea gsea gt, wt and gbst, in a simulation study and on two ovarian cancer data sets. from the simulation results, we found that gt, wt, and gbst are more powerful than gsea and gsea furthermore, the power of the five tests is substantially affected by the correlation structure of genes and the association between survival and the genes. the power of the five tests can approach when the genes are correlated, survival is positively associated with the gene expression values, and mp ≥. for more powerful test results, it might be desirable to check the correlation structure among genes within each pathway and the direction of association between genes and survival. although gsea was originally based on the rank of the correlation coefficients between a binary phenotype and gene expression levels, we replaced the correlation coefficient by the regression coefficient from a cox model and took the rank of the absolute value of its standardized coefficient. we then followed the rest of the gsea procedure in order to identify significant gene-sets or pathways. this is a modification of gsea that takes into account the survival phenotype, which may in some cases be more informative than a binary phenotype. in addition, the modified gsea allows covariance adjustment for age, sex and other clinical variables by taking the regression coefficient of genes from the cox model with those adjusting covariates. however, the performance of the gsea tests, gsea and gsea are not satisfactory except for a few cases. this may be due to the fact that the gsea tests are nonparametric approaches using the rank-based statistic instead of using the value of the regression coefficient. on the other hand, gt, wt and gbst have been proposed for regression-based models and can be extended to any phenotype such as binary, continuous, multi-class, or survival. these three tests can be easily adjusted for covariates, in order to determine whether gene expression profiles have an association with survival beyond what is explained by the adjusting covariates. the gt method assumes a random-effect model for the parameters corresponding to the coefficients of the individual genes in the pathway, in which the parameters are random variables and samples from n. here all parameters are assumed to have a common variance of τ which can be extended to have a more complex covariance structure as mentioned in goeman et al.. the gt method uses a score test of τ = which is equivalent to the null hypothesis that there is no association between survival and a given set of genes. therefore, this test does not depend on the number of genes in a given set and works under the assumption of a random-effect model with a common variance. as discussed in goeman et al., the score test of gt can have the optimal power against alternatives with small values of parameter τ in the simulation study, we generated the parameter from either a uniform or normal distribution with small variances. this may have contributed to the high power that we observed for gt. on the other hand, boulesteix and hothorn proposed the gbst for testing the additional predictive value of gene expression data while adjusting for clinical covariates. they combined the standard regression models such as a logistic regression or a cox model with a boosting procedure and used a permutation-based testing scheme. while goeman et al. reduced the multi-dimensional gene profiles into one-dimensional variance using a random-effect model, adewale et al. instead computed the sum of squares of the wald statistics for the genes in the pathway from the univariate regression model. likewise, the wt does not depend on the multi-dimensionality of the gene sets since the test statistic can be summarized as the sum of individual association measures. compared to the gsea tests, these three tests are based on parametric approaches since the values of the regression coefficients are taken into account in these tests via a score statistic or a sum of squares. as pointed out by a referee, methodological issues should be considered when comparing gene-set analysis methods. goeman and bühlmann addressed the definition of the null hypothesis, and in particular competitive versus self-contained tests, as well as the calculation of p-values, and the use of gene sampling versus subject sampling methods. gt, wt and gbst are based on the classical statistical models which lead to test a self-contained null hypothesis. furthermore, they calculate the p-value using a subject sampling model. on the other hand, gsea is a hybrid method: a kolmogorov-smirnov test statistic is motivated by a gene-sampling model, whereas a subject-sampling model is used to calculate the p-value. it was also pointed out that gsea sometimes shows low power since the model and null hypothesis used to motivate the test statistic are different from the model and null hypothesis used to calculate the p-value. simulation results indicate that both gsea and gsea have lower power than gt, wt and gbst. one of the important advantages of gene-set analysis over the single gene approach is that the multiple testing problems can be alleviated because the number of pathways is much smaller than that of genes. in the ovarian cancer example, for example, there are genes. the number of pathways is which is more than a 100-fold reduction. when a single pathway is of interest in pathway analysis, the issue of multiple testing is not a problem. however, when multiple pathways are of interest, as in cases when pathway analysis is exploratory with many pathways of potential interest, then multiple comparison methods such as false-discovery rates must be applied. in the ovarian cancer example, we used fdr to control pathway comparisons. a few pathways were found to have q <. it is noted that all pathways identified by crijns et al. contain more than genes. their sizes vary from to and their median is however, gt detects pathways whose sizes range from to with a median of wt detects pathways whose sizes range from to with a median of and gbst detects pathways whose sizes range from to with a median of in summary, gt, wt and gbst detected pathways with a variety of gene-set sizes, whereas crijns et al. detected only large gene sets. to implement gt and gbst, we adapted the function'gt' in the r package'globaltest', and the function'globalboosttest' in the r package'globalboosttest'. we implemented gsea gsea and wt from scratch in r, using the function'coxph'. for the readers' convenience, all codes are available as additional files. the additional file is r source codes for generating data, calling subroutines, and calculating permuted p-values and the additional file is r source codes for calculating five test statistics, gsea gsea gt, wt, and gbst.  in conclusion, gt, wt and gbst have high power to identify gene sets whose expression is associated with survival. survival is often strongly affected by well-known clinical variables such as stage, type of tumor, and tumor size, as well as demographic variables such as age, sex, and race. therefore, it is very important to evaluate the effects of genes on survival while considering known covariates. since these tests allow for adjustment for covariates, they assess whether expression of a given set of genes has an association with survival after controlling for confounders. in general, gt, wt and gbst are more powerful than gsea and gsea however, when genes are correlated within each pathway and survival is positively associated with genes, it does not matter which test is applied to detect significant gene sets, provided that the number of significant genes is moderately large within each pathway. 
it seems that there is a synergistic effect in detecting significant gene sets when significant genes have within-class correlation and the association between survival and genes is positive or negative. we focused on the survival phenotype and compared the performances of gene set enrichment analysis, global test, wald-type test and global boost test methods in a simulation study and on two ovarian cancer data sets. for the ovarian cancer data set, using the fdr threshold of q <, the gt, wt and gbst detected and significant pathways, respectively, whereas neither gsea nor gsea detected any significant pathways. many gene-set analysis methods have been previously proposed and compared through simulation studies and analysis of real datasets for binary phenotypes. simulation studies and a real data example indicate that gt, wt and gbst tend to have high power, whereas gsea and gsea have lower power. simulation results indicated that gt, wt and gbst consistently have higher power than gsea and gsea across all scenarios. we considered two versions of gsea by allowing different weights: gsea uses equal weights, yielding results similar to the kolmogorov-smirnov test; while gsea2's weights are based on the correlation between genes and the phenotype. we compared gsea gsea gt, wt and gbst in a simulation study with various settings for the correlation structure of the genes and the association parameter between the survival outcome and the genes. however, the power of the five tests depends on the combination of correlation structure and association parameter.
 many gene-set analysis methods have been previously proposed and compared through simulation studies and analysis of real datasets for binary phenotypes. we focused on the survival phenotype and compared the performances of gene set enrichment analysis, global test, wald-type test and global boost test methods in a simulation study and on two ovarian cancer data sets. we considered two versions of gsea by allowing different weights: gsea uses equal weights, yielding results similar to the kolmogorov-smirnov test; while gsea2's weights are based on the correlation between genes and the phenotype.  we compared gsea gsea gt, wt and gbst in a simulation study with various settings for the correlation structure of the genes and the association parameter between the survival outcome and the genes. simulation results indicated that gt, wt and gbst consistently have higher power than gsea and gsea across all scenarios. however, the power of the five tests depends on the combination of correlation structure and association parameter. for the ovarian cancer data set, using the fdr threshold of q <, the gt, wt and gbst detected and significant pathways, respectively, whereas neither gsea nor gsea detected any significant pathways. in addition, among the pathways found significant by gt, wt, and gbst, three pathways - purine metabolism, leukocyte transendothelial migration and jak-stat signaling pathway - overlapped with those reported in previous ovarian cancer microarray studies.  simulation studies and a real data example indicate that gt, wt and gbst tend to have high power, whereas gsea and gsea have lower power. we also found that the power of the five tests is much higher when genes are correlated than when genes are independent, when survival is positively associated with genes. it seems that there is a synergistic effect in detecting significant gene sets when significant genes have within-class correlation and the association between survival and genes is positive or negative. 