BACKGROUND
microarray experiments that monitor gene expression profiles associated with different disease phenotypes have become commonplace in biomedical research.
each dna sequence represented in microarrays can be considered a potential biomarker.
classification and prediction using such genomic measurements may contribute to a better understanding of the genetic pathways involved in diseases, and hence may lead to better diagnosis of disease and better prescription of treatment.
analysis of genomic data is challenging due to high dimensionality of data and low sample size.
although the number of genes assayed is large, but there may be only a small number of genes that are associated with variations of phenotypes.
by employing standard methods directly, we usually obtain estimates that are not "regular", i.e., estimates are not unique or ill-behaved.
regularization, through which we achieve unique and well-behaved estimates, is usually needed.
regularization can be achieved via model reduction or variable selection methods.
several dimension reduction techniques have been employed for classification using genomic data.
examples include the partial least squares  <cit> , the principal component regression  <cit> , and the singular value decomposition under the bayesian framework  <cit> , among others.
by using low dimensional projections of covariates as features in model estimation, one may obtain estimators with better prediction performance due to the bias-variance tradeoff.
one drawback of such dimension reduction techniques is that all genes are used in estimation and prediction.
biological interpretation of such classifiers are usually not straightforward.
moreover, if certain genes are not associated with the clinical outcome, it is important to exclude them from the predictive model.
an alternative approach to the dimension reduction techniques is to use methods that are capable of simultaneous biomarker selection and model fitting, which can be realized by penalization.
such methods include the least absolute shrinkage and selection operator-lasso  <cit> , the least angle regression-lars  <cit> , and the threshold gradient directed regularization method-tgdr  <cit> .
these methods can produce parsimonious models with a small number of biomarkers and hence more lucent biological interpretations.
in this article, we propose an approach for simultaneous estimation and biomarker selection using a scaled tgdr method.
it is important to assess both false-positive and false-negative errors, since the clinical and financial consequences of the two types of errors can differ significantly.
a common practice is to use the receiver operating characteristic  curve  <cit> , where the classification performance can be measured by the area under the roc curve .
advantages of the roc method include:  it does not assume a parametric form of the class probability.
this is different from the logistic regression method.
although one can construct an roc curve from logistic regression, it assumes a parametric form of the class probability.
this parametric assumption may be not be satisfied;  it is adaptable to outcome-dependent samplings, for example the case control design; and  it is capable of penalizing false positives and false negatives differently.
therefore, the roc method may be preferable in biomarker selection and classification using genomic measurements.
pepe  <cit>  proposed the empirical auc as the objective function for combining multiple biomarkers in a low dimensional setting.
ma and huang  <cit>  proposed a smooth sigmoid approximation of the empirical auc for high-dimensional data.
an alternative to the empirical auc is the binormal auc  <cit> .
the binormal auc technique was developed parallel to, but separated from the empirical auc method.
for small sample sizes, the empirical auc may change dramatically due to small perturbations and differ significantly from the expected auc, whereas the binormal auc is more stable.
studies with low-dimensional biomarkers show that the binormal auc may provide valuable information beyond the empirical auc  <cit> .
for data with high dimensional covariates or large sample sizes, the binormal auc is computationally much more affordable than the empirical auc.
since both the empirical auc and the binormal auc are extensively used in low dimensional settings, it is of great interest to extend the study of  <cit>  and explore use of the binormal auc for disease classification with microarray data.
in this article, we proposed an approach for biomarker selection and classification with microarray data by optimizing the binormal auc.
the scaled tgdr method, which is a modified version of the tgdr, is adopted for estimation and biomarker selection.
tuning parameters are selected using v-fold cross validation, and monte carlo based methods are proposed for evaluation purposes.
we assess the proposed approach by extensive simulation studies and demonstrate it on two cancer datasets.
comparing to the method in  <cit> , which uses a smoothed version of the empirical auc as the objective function, the contributions of this paper are as follows.
first, the binormal auc, which is at least as important as the empirical auc, is used as the objective function.
second, we propose using the scaled tgdr, which can significantly reduce the computational cost.
moreover, the occurrence index, which is a way to rank the selected biomarkers and measure their relative stability in the presence of sampling variation, is proposed in this paper.
RESULTS
simulation studies
simulation studies are used to examine the finite sample performance of the proposed regularized roc method.
in simulation study i, we look at the case where the genes are assumed to be independently normally distributed.
this setting satisfies all the assumptions for the binormal roc method  <cit>  and the underlying model is the logistic regression model.
in simulation study ii, we consider independent genes with marginal distributions other than the normal.
the corresponding rocs are still binormal since the linear predictors for both the diseased and non-diseased groups are normally distributed after a common monotone transformation.
we investigate the performance of the proposed method when the transformation is intentionally neglected and the marginal normality assumption is violated.
in simulation study iii, we study the effect of correlations among genes.
in simulation study iv, we investigate the case where marginal normality cannot be obtained with simple transformations.
the simulation settings are related to but much more comprehensive than those in  <cit> .
using simulation settings similar to those in  <cit> , we first generate d =  <dig> dimensional vectors from the diseased and non-diseased populations.
we consider the following sample size combinations  = , ,  and , where nd and nh denote the number of samples in the groups with y =  <dig> and y =  <dig> , respectively.
different marginal and joint distributions of the genes are considered, representing a variety of data structures encountered in real medical studies.
we assume a model in which a fraction π of the potential biomarkers are differentially expressed and hence associated with the binary outcome of interest.
simulation study i
all the  <dig> biomarkers are assumed to be independent and normally distributed with variance  <dig>  we consider π =  <dig> , π =  <dig>  and π =  <dig> , which represent the cases where a small number, a moderate number and a large number of biomarkers are associated with the outcome.
moreover, we investigate the following three scenarios:
<dig>   there is a shift of  <dig>  in the mean expression levels for the differentially expressed biomarkers;
<dig>   the fold change is  <dig>  unit difference in mean;
<dig>   a shift of  <dig> units in the mean is considered.
the underlying model is the logistic regression model with regression coefficients equal to two times the mean shifts for genes differentially expressed and zero for the rest.
for each simulated dataset, we first randomly sample a training set of size 2/ <dig> ×  without consideration of the stratification of the binary outcome.
the testing set is composed of the rest observations.
we construct the proposed roc estimate based on the training set only.
tuning parameters are selected using 3-fold cross validation.
the predictive classification error and binormal auc are computed for the testing set based on the estimate obtained from the training set.
summary statistics based on  <dig> simulated datasets are given in table  <dig>
it can be seen that for moderate and large effect cases, the binormal aucs are close to  <dig> with very small classification errors.
the classification performance for small mean shift cases is less satisfactory, although still acceptable when the sample size is not too small, due to the relatively small signal to noise ratios.
as sample size increases, the classification performance improves.
the classification performance also improves with more differentially expressed genes.
the empirical aucs are very close to the binormal aucs.
results for the empirical aucs are omitted here.
simulation study ii
we consider data settings similar to those in simulation i, with π =  <dig>  and  <dig>  unit shift in mean for differentially expressed genes.
we assume all biomarkers are independent and marginally distributed as uniform[-3
mathtype@mtef@5@5@+=feaafiart1ev1aaatcvaufkttlearuwrp9mdh5mbpbiqv92aaexatlxbi9gbaebbnrfifhhdyfgasaach8aky=wiffydh8gipec8eeeu0xxdbba9frfj0=oqffea0dxdd9vqai=hguq8kuc9pgc9s8qqaq=dirpe0xb9q8qilsfr0=vr0=vr0dc8meaabaqaciaacagaaeqabaqabegadaaakeaadagcaaqaaiabiodazawcbeaaaaa@2dbb@, 3
mathtype@mtef@5@5@+=feaafiart1ev1aaatcvaufkttlearuwrp9mdh5mbpbiqv92aaexatlxbi9gbaebbnrfifhhdyfgasaach8aky=wiffydh8gipec8eeeu0xxdbba9frfj0=oqffea0dxdd9vqai=hguq8kuc9pgc9s8qqaq=dirpe0xb9q8qilsfr0=vr0=vr0dc8meaabaqaciaacagaaeqabaqabegadaaakeaadagcaaqaaiabiodazawcbeaaaaa@2dbb@] gamma, and 22
mathtype@mtef@5@5@+=feaafiart1ev1aaatcvaufkttlearuwrp9mdh5mbpbiqv92aaexatlxbi9gbaebbnrfifhhdyfgasaach8aky=wiffydh8gipec8eeeu0xxdbba9frfj0=oqffea0dxdd9vqai=hguq8kuc9pgc9s8qqaq=dirpe0xb9q8qilsfr0=vr0=vr0dc8meaabaqaciaacagaaeqabaqabegadaaakeaadagcaaqaaiabikdayawcbeaaaaa@2db9@beta for the non-diseased group.
the diseased group has the same marginal distributions with only mean shifts.
all three marginal distributions have variance  <dig> and represent three different distribution scenarios: uniformly distributed gene expressions; gene expressions with a skewed unimodal distribution ; and gene expressions clustering around the upper and lower limits .
proper transformations can lead to marginal normal distributions.
we intentionally neglect the transformations and employ the proposed approach directly.
the purpose is to investigate the performance of the proposed approach when marginal distribution assumption is violated.
summary statistics based on  <dig> simulated datasets are shown in table  <dig>  for comparison, we also reproduce the results from the normal distribution case in simulation i. it can be seen that the proposed approach has satisfactory classification performance under all scenarios.
the effect of different marginal distributions is not significant.
the empirical aucs are also reasonably close to the binormal aucs.
simulation studies with different mean shifts and different choices of π lead to conclusions similar to those for simulation study i. those results are omitted here.
simulation study iii
in simulation studies i and ii, possible correlations among genes have been neglected.
in simulation study iii, we investigate the case where the genes are correlated.
we still assume π =  <dig>  and  <dig>  unit shift in mean for differentially expressed genes.
for the non-diseased subjects, let zi <dig>  zi <dig>  ..., i =  <dig>  ..., nh, be a sequence of i.i.d.
standard normal random variables.
denote the d × nh dimensional gene expression matrix for healthy subjects as x. for the jth gene of the ith subject, we consider the following models:
<dig>  model 1: xij=h×m+1×m+kzil/k)
mathtype@mtef@5@5@+=feaafiart1ev1aaatcvaufkttlearuwrp9mdh5mbpbiqv92aaexatlxbi9gbaebbnrfifhhdyfgasaach8aky=wiffydh8gipec8eeeu0xxdbba9frfj0=oqffea0dxdd9vqai=hguq8kuc9pgc9s8qqaq=dirpe0xb9q8qilsfr0=vr0=vr0dc8meaabaqaciaacagaaeqabaqabegadaaakeaacqwgybawdawgaawcbagaemyaakmaemoaaogabeaakiabg2da9iabdieainaabmaabawaaabmaeaacqwgabgwdawgaawcbagaemyaakmaemibawgabeaakiabc+cavmaakaaabagaem4aasgaleqaaaqaaiabdygasjabg2da9maabmaabagaemoaaomaeyoei0iaegymaedacagloagaayzkaagaey41aqraemyba0maey4kasiaegymaedabawaaewaaeaacqwgqbgacqghsislcqaixaqmaiaawicacaglpaaacqghxdatcqwgtbqbcqghrawkcqwgrbwaa0gaeyyeiuoaaogaayjkaiaawmcaaaaa@53ff@ with k =  <dig>  m = 7; h = 3
mathtype@mtef@5@5@+=feaafiart1ev1aaatcvaufkttlearuwrp9mdh5mbpbiqv92aaexatlxbi9gbaebbnrfifhhdyfgasaach8aky=wiffydh8gipec8eeeu0xxdbba9frfj0=oqffea0dxdd9vqai=hguq8kuc9pgc9s8qqaq=dirpe0xb9q8qilsfr0=vr0=vr0dc8meaabaqaciaacagaaeqabaqabegadaaakeaadagcaaqaaiabiodazawcbeaaaaa@2dbb@, where Φ is the cumulative distribution for the standard normal.
<dig>  model 2: same as model  <dig>  but with k =  <dig>  m =  <dig>
<dig>  model 3: same as model  <dig>  with h as the identify function.
<dig>  model 4: same as model  <dig>  with h as the identify function.
this simulation setting has been considered in  <cit> .
for models  <dig> and  <dig>  the genes are marginally uniformly distributed.
for models  <dig> and  <dig>  the genes are marginally normally distributed.
for genes in models  <dig> and  <dig>  the correlations between all genes  are weak; whereas models  <dig> and  <dig> are examples of strongly correlated genes.
note that we not only assume differential genes are correlated, we also assume differentially expressed and non-differentially expressed genes are correlated.
summary statistics based on  <dig> replicates are shown in table  <dig>  for comparison, we also reproduce the results with independently, marginally uniformly  distributed genes.
we can see from table  <dig> that the effects of correlation are negligible.
we looked at other correlation structures and obtained similar results.
moreover the effect of different marginal distributions are not significant.
simulation study iv
the binormal auc is theoretically constructed based on the joint normal distribution assumption on the linear combination of the covariates in the diseased and non-diseased population after a common monotone transformation.
we now consider the case when marginal normality cannot be achieved with simple transformation.
specifically, we investigate the case when the marginal gene distributions for the diseased class and the non-diseased class are different.
we assume the number of genes is  <dig> and there are  <dig> subjects in each group.
we assume that π =  <dig>  percent genes are differentially expressed.
moreover, we assume the genes in the non-diseased group are independent and have standard normal distribution.
for subjects in the diseased class, we assume the genes are independently distributed with:
<dig>  scenario 1: uniform [ <dig>  - 3
mathtype@mtef@5@5@+=feaafiart1ev1aaatcvaufkttlearuwrp9mdh5mbpbiqv92aaexatlxbi9gbaebbnrfifhhdyfgasaach8aky=wiffydh8gipec8eeeu0xxdbba9frfj0=oqffea0dxdd9vqai=hguq8kuc9pgc9s8qqaq=dirpe0xb9q8qilsfr0=vr0=vr0dc8meaabaqaciaacagaaeqabaqabegadaaakeaadagcaaqaaiabiodazawcbeaaaaa@2dbb@,  <dig>  + 3
mathtype@mtef@5@5@+=feaafiart1ev1aaatcvaufkttlearuwrp9mdh5mbpbiqv92aaexatlxbi9gbaebbnrfifhhdyfgasaach8aky=wiffydh8gipec8eeeu0xxdbba9frfj0=oqffea0dxdd9vqai=hguq8kuc9pgc9s8qqaq=dirpe0xb9q8qilsfr0=vr0=vr0dc8meaabaqaciaacagaaeqabaqabegadaaakeaadagcaaqaaiabiodazawcbeaaaaa@2dbb@] for differential genes and uniform [-3
mathtype@mtef@5@5@+=feaafiart1ev1aaatcvaufkttlearuwrp9mdh5mbpbiqv92aaexatlxbi9gbaebbnrfifhhdyfgasaach8aky=wiffydh8gipec8eeeu0xxdbba9frfj0=oqffea0dxdd9vqai=hguq8kuc9pgc9s8qqaq=dirpe0xb9q8qilsfr0=vr0=vr0dc8meaabaqaciaacagaaeqabaqabegadaaakeaadagcaaqaaiabiodazawcbeaaaaa@2dbb@, 3
mathtype@mtef@5@5@+=feaafiart1ev1aaatcvaufkttlearuwrp9mdh5mbpbiqv92aaexatlxbi9gbaebbnrfifhhdyfgasaach8aky=wiffydh8gipec8eeeu0xxdbba9frfj0=oqffea0dxdd9vqai=hguq8kuc9pgc9s8qqaq=dirpe0xb9q8qilsfr0=vr0=vr0dc8meaabaqaciaacagaaeqabaqabegadaaakeaadagcaaqaaiabiodazawcbeaaaaa@2dbb@]otherwise;
<dig>  scenario 2: 22
mathtype@mtef@5@5@+=feaafiart1ev1aaatcvaufkttlearuwrp9mdh5mbpbiqv92aaexatlxbi9gbaebbnrfifhhdyfgasaach8aky=wiffydh8gipec8eeeu0xxdbba9frfj0=oqffea0dxdd9vqai=hguq8kuc9pgc9s8qqaq=dirpe0xb9q8qilsfr0=vr0=vr0dc8meaabaqaciaacagaaeqabaqabegadaaakeaadagcaaqaaiabikdayawcbeaaaaa@2db9@beta - 2
mathtype@mtef@5@5@+=feaafiart1ev1aaatcvaufkttlearuwrp9mdh5mbpbiqv92aaexatlxbi9gbaebbnrfifhhdyfgasaach8aky=wiffydh8gipec8eeeu0xxdbba9frfj0=oqffea0dxdd9vqai=hguq8kuc9pgc9s8qqaq=dirpe0xb9q8qilsfr0=vr0=vr0dc8meaabaqaciaacagaaeqabaqabegadaaakeaadagcaaqaaiabikdayawcbeaaaaa@2db9@ +  <dig>  for differential genes and 22
mathtype@mtef@5@5@+=feaafiart1ev1aaatcvaufkttlearuwrp9mdh5mbpbiqv92aaexatlxbi9gbaebbnrfifhhdyfgasaach8aky=wiffydh8gipec8eeeu0xxdbba9frfj0=oqffea0dxdd9vqai=hguq8kuc9pgc9s8qqaq=dirpe0xb9q8qilsfr0=vr0=vr0dc8meaabaqaciaacagaaeqabaqabegadaaakeaadagcaaqaaiabikdayawcbeaaaaa@2db9@beta - 2
mathtype@mtef@5@5@+=feaafiart1ev1aaatcvaufkttlearuwrp9mdh5mbpbiqv92aaexatlxbi9gbaebbnrfifhhdyfgasaach8aky=wiffydh8gipec8eeeu0xxdbba9frfj0=oqffea0dxdd9vqai=hguq8kuc9pgc9s8qqaq=dirpe0xb9q8qilsfr0=vr0=vr0dc8meaabaqaciaacagaaeqabaqabegadaaakeaadagcaaqaaiabikdayawcbeaaaaa@2db9@ otherwise;
<dig>  scenario 3: gamma +  <dig> for differential genes and gamma -  <dig>  otherwise.
under all three scenarios, the differentially expressed genes for the diseased subjects have mean  <dig>  and variance  <dig>  and non-differential genes have mean  <dig> and variance  <dig>  based on  <dig> replicates, the mean aucs are  <dig>  ,  <dig>   and  <dig>   under scenarios 1– <dig>  respectively, where the values in the parentheses are the corresponding standard errors.
the mean classification errors are  <dig>  ,  <dig>   and  <dig>  , respectively.
colon and estrogen data
data description
colon data
in this dataset, expression levels of  <dig> tumor and  <dig> normal colon tissues for  <dig> human genes are measured using the affymetrix gene chip.
a selection of  <dig> genes with the highest minimal intensity across the samples has been made by  <cit> , and these data are publicly available at  <cit> .
the colon data have been analyzed in several previous studies using other statistical approaches, see for example  <cit> .
estrogen data
this dataset was first presented by  <cit> .
it includes expression values of  <dig> genes of  <dig> breast tumor samples.
the expression data were obtained using the affymetrix gene chip technology and are available at  <cit> .
the response describes the status of the estrogen receptor .
among the  <dig> samples,  <dig> are positive  and  <dig> are negative .
we threshold the raw data with a floor of  <dig> and a ceiling of  <dig>  genes with max/min <  <dig> and/or max - min <  <dig> are also excluded  <cit> .
<dig>  genes pass the first step screening.
a base  <dig> logarithmic transformation is then applied.
the estrogen data have also been studied by  <cit> .
estimation
the roc estimate is identifiable only up to a scale constant  <cit> .
so prior to the analysis, we need to identify the "anchor biomarker," i.e, the gene whose estimated coefficient will be set as a constant for identifiability purposes.
we select the anchor biomarker as follows.
compute the sample standard errors of the d biomarkers se, ..., se and denote their median as med.se.
compute the adjusted standard errors as  <dig>  + med.se), ...,  <dig>  + med.se).
then the biomarkers are ranked based on the t-statistics computed with the adjusted standard errors.
this adjusted t-statistic is similar to a simple shrinkage method discussed in  <cit> .
it has been observed in microarray studies that the standard t-statistic can be very large even if there is very little difference between the mean expressions , due to the small variance estimate caused by small sample size  <cit> .
so compared with the standard t-statistic, the adjusted t-statistic takes into account variability as well as fold change.
the biomarker with the largest absolute value of the adjusted t-statistic is chosen as the anchor biomarker.
for the anchor biomarker, if the sample mean of the diseased class is larger, β =  <dig> , otherwise β = - <dig> .
further investigation of the genes that have passed the first step screening shows that quite a few of those genes have very little variations across sample.
including such genes in the joint modeling may create unreliable estimates.
moreover, it is unlikely all the  <dig> genes for the colon data  will be associated with the outcome.
so we only use the  <dig> genes with the largest absolute values of the adjusted t-statistics for classification.
the genes are then standardized to have zero means and unit variances.
the same anchor marker detection and pre-processing have been considered in  <cit>  and are shown to be well-behaved.
note that there is no computational limitation on how many genes can be used in the roc classification.
empirical studies in  <cit>  show that the performance of certain statistical approaches decreases dramatically as the number of genes increases.
so it is believed that identification of the top  <dig> genes will help to generate more reliable estimate.
we use three-fold cross validation for tuning parameter selection.
we show in table  <dig> the number of genes with nonzero coefficients and the cv scores for each fixed threshold value τ. it can be seen that generally the number of nonzero coefficients decreases as τ increases.
however, the change of cv score is very small.
our extensive empirical studies show that in general a larger τ will lead to a more parsimonious model.
parsimonious models are preferred when the cv scores are comparable.
so we choose τ =  <dig>  for both datasets.
in the tgdr, variable selection is achieved by thresholding some estimated coefficients to be exactly zero based on their derivatives, so the corresponding covariates are not selected in the final model.
the genes in tables  <dig> and  <dig> have nonzero coefficients and other genes have estimated coefficients zero and hence not selected.
since the gene expressions have been normalized to have unit variances, the estimated coefficients are directly comparable.
larger absolute values of coefficients indicate stronger influences.
more detailed information about these genes can be found on the ncbi website.
for example, for the colon data, hsa.
<dig> is from the interleukin  <dig> gene .
this gene encodes a member of the cytokine family.
expression of this protein is increased after the activation of t-cells by mitogens or the activation of nk cells by il- <dig>  this protein induces the production of tnfalpha from macrophage cells.
hsa.
<dig>  encodes a potent cyclin-dependent kinase inhibitor.
the encoded protein binds to and inhibits the activity of cyclin-cdk <dig> or -cdk <dig> complexes, and thus functions as a regulator of cell cycle progression at gl.
the expression of this gene is tightly controlled by the tumor suppressor protein p <dig>  through which this protein mediates the p53-dependent cell cycle g l phase arrest in response to a variety of stress stimuli.
this protein can interact with proliferating cell nuclear antigen , a dna polymerase accessory factor, and plays a regulatory role in s phase dna replication and dna damage repair.
hsa.
<dig> encodes a protein that is a member of the cxc chemokine family.
this chemokine is one of the major mediators of the inflammatory response and is secreted by several cell types.
it functions as a chemoattractant, and is also a potent angiogenic factor.
this gene is believed to play a role in the pathogenesis of bronchiolitis, a common respiratory tract disease caused by viral infection.
this gene and the other ten members of the cxc chemokine gene family form a chemokine gene cluster in a region mapped to chromosome 4q.
for the estrogen data, d <dig>  encodes a protein belonging to a small family of adapter proteins that are known to interact with a number of receptor tyrosine kinases and signaling molecules.
this gene encodes a growth factor receptor-binding protein that interacts with epidermal growth factor receptor  and ephrin receptors.
the protein plays a role in the integrin signaling pathway and cell migration by binding with focal adhesion kinase .
k <dig>  encodes a member of the bombesin-like family of gastrin-releasing peptides.
its preproprotein, following cleavage of a signal peptide, is further processed to produce either the  <dig> aa gastrin-releasing peptide or the  <dig> aa neuromedin c. these smaller peptides regulate numerous functions of the gastrointestinal and central nervous systems, including release of gastrointestinal hormones, smooth muscle cell contraction, and epithelial cell proliferation.
these peptides are also likely to play a role in human cancers of the lung, colon, stomach, pancreas, breast, and prostate.
u <dig>  plays a role in a number of diseases.
for example, it impacts the progression of hormone receptor-negative breast cancer.
it is also expressed in cutaneous cd30+ lymphoproliferative disorders and is a candidate marker of disease progression.
x <dig>  is the estrogen receptor  that is a ligand-activated transcription factor composed of several domains important for hormone binding, dna binding, and activation of transcription.
evaluation
we first randomly partition the data into training and testing sets.
we then identify the top  <dig> genes with the largest adjusted t-statistics based on the training sets only.
so the  <dig> genes used in each individual partition and computation may be different, although there is considerable overlap.
the occurrence index, which measures the relative stability of individual selected genes in the joint model, is shown in figures  <dig> and  <dig> for the colon and estrogen data, respectively.
we observe that many genes selected with the proposed approach have relatively high occurrence index compared to genes not selected, which suggests that the roc estimates are relatively stable for these two specific datasets.
on the other hand, we also see that the occurrence index, even for the genes selected, can be low in absolute terms.
most selected genes have occurrence index less than  <dig> .
similar phenomenon has been observed in  <cit> .
the genes identified in the final model do not necessarily have the largest occurrence index.
for example, in the colon data the gene with the highest occurrence index was not selected in the model.
this is not surprising since  the objective function is the binormal auc and the genes were selected to optimize this objective function;  the genes were selected based on all the available data, whereas the occurrence index was calculated using partitioned data.
the occurrence index is one way to evaluate the relative importance of the genes after they are selected and it provides a measure of stability of the genes in the presence of sampling variation, but it does not measure the prediction power of the genes.
furthermore, the genes significantly associated with the outcome in different subsets are not necessarily important when the whole dataset is used.
thus we do not select genes based on the occurrence index.
the classification performance of the proposed method for the two datasets is evaluated using the monte carlo approach described in the methods section.
both the number of random partitions for the observed data and the number of permutations are  <dig>  for the colon data, the roc estimate has mean aucs  <dig>   and  <dig>   for the opd and ppd, respectively, where the values in the parentheses are the standard errors.
corresponding mean classification errors are  <dig>   for observed data and  <dig>   for randomly permuted data, respectively.
in the upper-right panel of figures  <dig> and  <dig>  we also show the kernel density estimates of the opd and ppd of the auc.
we can see that the two distributions are well separated.
a simple wilcoxon test of the difference of these two distributions gives p-value <  <dig> .
for the same dataset and using the top ranked  <dig> genes based on the marginal wilcoxon rank test, table  <dig> in  <cit>  shows that the mean classification errors are  <dig>  ,  <dig>  ,  <dig>   and  <dig>  .
in  <cit> ,  <dig> genes are used and the svm based methods have mean aucs  <dig>  and mean classification errors  <dig> , whereas the principal component analysis based approaches have even smaller aucs and larger classification errors.
using the same set of  <dig> genes, the logistic-lasso method in  <cit>  has mean aucs  <dig>  and mean classification errors  <dig> ; the smrc method in  <cit>  has mean aucs  <dig>  and mean classification errors  <dig> .
so for the colon data, the proposed binormal roc method performs better than other approaches in terms of auc and classification error, although the improvement is not dramatic since most other approaches behave reasonably well.
we also note that since different sets of genes are used in  <cit> , the above results can only provide a rough comparison.
for the estrogen data, the proposed roc estimate has mean aucs  <dig>   and  <dig>   for the opd and ppd, respectively, with corresponding mean classification errors  <dig>   and  <dig>  .
for the same dataset, dettling and buhlmann  use  <dig> genes selected based on the marginal wilcoxon tests and yield classification errors  <dig> ,  <dig>  and  <dig>  with different logitboost methods, classification error  <dig>  for adaboost,  <dig>  for 1-nearest-neighbor approach and  <dig>  by using the classification tree.
using the same set of genes, the logistic-lasso  <cit>  has mean aucs  <dig>  and mean classification errors  <dig> ; whereas the smrc  <cit>  has mean aucs  <dig>  and mean classification errors  <dig> .
note that in our study, joint models are considered, i.e, we consider the joint effects instead of marginal effects of genes.
identification of differentially expressed genes associated with clinical outcomes based on marginal significance has been well investigated and accepted.
so as a simple evaluation, we show the marginal binormal auc and t-statistic for each individual gene in figures  <dig> and  <dig>  similar evaluation has been considered in  <cit>  as a gene ranking method.
we can see from figures  <dig> and  <dig> that several genes we identify do have large marginal binormal aucs and t-statistics; both up- and down- regulated genes are identified in the final models.
however, we also identify some genes with small binormal aucs or small t-statistics.
it is believed that these genes are important in the joint models, but not in the marginal models.
variable selection using tgdr looks at all the variables and considers the relative importance of each variable in terms of its contribution to change in the objective function, but is not based on marginal significance.
so the results do not necessarily correspond to those based on marginal significance.
when binormal roc method is used, it is of interest to examine whether the roc curve is concave.
pepe  <cit>  pointed out the binormal roc is not concave when the variances of the two groups are not equal.
however, she also pointed that "swets  and hanley  conclude that the binomal roc curve provides a good approximation to a wide range of roc curves that occur in practice  <cit> .
this may be because the undesirable behavior of the approximating binormal roc usually occurs over a very small part of the roc curve."
as a simple check, we show in figure  <dig> the plots of binormal roc for a randomly selected training set and a testing set.
we can clearly see the concavity of the binormal roc.
examination of other training and testing sets shows similar results.
we note that although the genes identified using the two datasets exhibit strong predictive power for disease statuses as demonstrated above, the analysis here does not provide information on whether they are just genomic markers correlated with the disease statuses or are actually in the pathways leading to the diseases.
indeed, as is typical with microarray results, further investigation based on independent assays and/or independent samples is desirable.
CONCLUSIONS
the roc method has been successfully used for disease classification using low- dimensional biomarkers.
the empirical auc method has been successfully implemented in  <cit> .
in this article, we propose using the binormal auc, which is a useful alternative to the empirical auc, as the objective function for disease classification using high dimensional genomic data.
the scaled tgdr is employed for regularization and biomarker selection.
extensive simulation studies and applications to two cancer data show that the proposed approach can generate satisfactory classification results with a small number of identified biomarkers with biologically sound implications.
although theoretically the binormal auc is constructed under the normal distribution assumption, the proposed approach still behaves satisfactorily when this assumption is seriously violated.
beyond the merits carried by the roc method, the proposed approach has the advantage of little computational cost.
using the scaled gradient significantly reduces the computational complexity.
it takes only a few seconds to analyze dataset with  <dig> genes and  <dig> subjects using r  <cit>  code, whereas the empirical auc method takes several minutes or longer.
in our data analyses, we selected  <dig> genes with the largest adjusted t-statistics.
the proposed approach has no limitation on the number of genes that can be used in the final model.
we note that use of  <dig> genes is somewhat subjective.
same analyses with  <dig> genes and all the genes have also been carried out , and there are considerable overlaps among genes identified.
we prefer using  <dig> genes because this set of genes should include most  important genes without including too much noise.
we used the same datasets as in  <cit> .
although both methods are roc based and the same sets of genes are used, the identified biomarkers are different.
considering that both the empirical auc and the binormal auc methods are extensively used, such difference warrants further investigation and confirmation by independent studies.
in our study, the adjusted t-statistic is used instead of the standard t-statistic.
the reason we used the adjusted t-statistic is that the unadjusted t-statistics can be large for genes with small mean differences in expression levels, if the standard error is very small.
because of a large number of genes on each array, there will usually be genes with very small standard errors, and some of these genes will have small means as well.
this phenomenon has been discussed by  <cit>  in motivating their adjusted t-statistic for significance analysis of microarrays .
similar variance adjust has also been considered in  <cit> .
however, a recent study  <cit>  points out that such variance adjust may not be necessary.
since current study is not conclusive, we use the adjusted t-statistic and leave the choice of t-statistic versus adjusted t-statistic to a future study.
in this paper, we proposed monte carlo based methods for evaluation of prediction and stability of individual genes.
the proposed methods are closely related to cross validation and permutation methods.
the problem of evaluation is extremely important for high dimensional classification, but it has not been extensively studied.
although the proposed methods have sound, intuitive implications, further studies of their theoretical properties would be of interest.
we focus on the linear combination of genes that maximizes the binormal-aug. in  <cit> , linear combinations of genes favoring certain specificity or other criterion are investigated.
if maximization of specificity is of interest, we can then combine liu's methods with certain variable selection techniques, for example lasso or the proposed tgdr.
we leave this to a future study.
we only discussed two-sample classification in this article.
when extension of the binormal auc to multi-class study exists, the proposed estimation and regularization method can be used.
we postpone multi-class study to a future report.

 microarray experiments that monitor gene expression profiles associated with different disease phenotypes have become commonplace in biomedical research. each dna sequence represented in microarrays can be considered a potential biomarker. classification and prediction using such genomic measurements may contribute to a better understanding of the genetic pathways involved in diseases, and hence may lead to better diagnosis of disease and better prescription of treatment. analysis of genomic data is challenging due to high dimensionality of data and low sample size. although the number of genes assayed is large, but there may be only a small number of genes that are associated with variations of phenotypes. by employing standard methods directly, we usually obtain estimates that are not "regular", i.e., estimates are not unique or ill-behaved. regularization, through which we achieve unique and well-behaved estimates, is usually needed. regularization can be achieved via model reduction or variable selection methods. several dimension reduction techniques have been employed for classification using genomic data. examples include the partial least squares, the principal component regression, and the singular value decomposition under the bayesian framework, among others. by using low dimensional projections of covariates as features in model estimation, one may obtain estimators with better prediction performance due to the bias-variance tradeoff. one drawback of such dimension reduction techniques is that all genes are used in estimation and prediction. biological interpretation of such classifiers are usually not straightforward. moreover, if certain genes are not associated with the clinical outcome, it is important to exclude them from the predictive model. an alternative approach to the dimension reduction techniques is to use methods that are capable of simultaneous biomarker selection and model fitting, which can be realized by penalization. such methods include the least absolute shrinkage and selection operator-lasso, the least angle regression-lars, and the threshold gradient directed regularization method-tgdr. these methods can produce parsimonious models with a small number of biomarkers and hence more lucent biological interpretations. in this article, we propose an approach for simultaneous estimation and biomarker selection using a scaled tgdr method. it is important to assess both false-positive and false-negative errors, since the clinical and financial consequences of the two types of errors can differ significantly. a common practice is to use the receiver operating characteristic curve, where the classification performance can be measured by the area under the roc curve. advantages of the roc method include: it does not assume a parametric form of the class probability. this is different from the logistic regression method. although one can construct an roc curve from logistic regression, it assumes a parametric form of the class probability. this parametric assumption may be not be satisfied; it is adaptable to outcome-dependent samplings, for example the case control design; and it is capable of penalizing false positives and false negatives differently. therefore, the roc method may be preferable in biomarker selection and classification using genomic measurements. pepe proposed the empirical auc as the objective function for combining multiple biomarkers in a low dimensional setting. ma and huang proposed a smooth sigmoid approximation of the empirical auc for high-dimensional data. an alternative to the empirical auc is the binormal auc. the binormal auc technique was developed parallel to, but separated from the empirical auc method. for small sample sizes, the empirical auc may change dramatically due to small perturbations and differ significantly from the expected auc, whereas the binormal auc is more stable. studies with low-dimensional biomarkers show that the binormal auc may provide valuable information beyond the empirical auc. for data with high dimensional covariates or large sample sizes, the binormal auc is computationally much more affordable than the empirical auc. since both the empirical auc and the binormal auc are extensively used in low dimensional settings, it is of great interest to extend the study of and explore use of the binormal auc for disease classification with microarray data. in this article, we proposed an approach for biomarker selection and classification with microarray data by optimizing the binormal auc. the scaled tgdr method, which is a modified version of the tgdr, is adopted for estimation and biomarker selection. tuning parameters are selected using v-fold cross validation, and monte carlo based methods are proposed for evaluation purposes. we assess the proposed approach by extensive simulation studies and demonstrate it on two cancer datasets. comparing to the method in, which uses a smoothed version of the empirical auc as the objective function, the contributions of this paper are as follows. first, the binormal auc, which is at least as important as the empirical auc, is used as the objective function. second, we propose using the scaled tgdr, which can significantly reduce the computational cost. moreover, the occurrence index, which is a way to rank the selected biomarkers and measure their relative stability in the presence of sampling variation, is proposed in this paper.  simulation studies simulation studies are used to examine the finite sample performance of the proposed regularized roc method. in simulation study i, we look at the case where the genes are assumed to be independently normally distributed. this setting satisfies all the assumptions for the binormal roc method and the underlying model is the logistic regression model. in simulation study ii, we consider independent genes with marginal distributions other than the normal. the corresponding rocs are still binormal since the linear predictors for both the diseased and non-diseased groups are normally distributed after a common monotone transformation. we investigate the performance of the proposed method when the transformation is intentionally neglected and the marginal normality assumption is violated. in simulation study iii, we study the effect of correlations among genes. in simulation study iv, we investigate the case where marginal normality cannot be obtained with simple transformations. the simulation settings are related to but much more comprehensive than those in. using simulation settings similar to those in, we first generate d = dimensional vectors from the diseased and non-diseased populations. we consider the following sample size combinations =,, and, where nd and nh denote the number of samples in the groups with y = and y =, respectively. different marginal and joint distributions of the genes are considered, representing a variety of data structures encountered in real medical studies. we assume a model in which a fraction π of the potential biomarkers are differentially expressed and hence associated with the binary outcome of interest. simulation study i all the biomarkers are assumed to be independent and normally distributed with variance we consider π =, π = and π =, which represent the cases where a small number, a moderate number and a large number of biomarkers are associated with the outcome. moreover, we investigate the following three scenarios: there is a shift of in the mean expression levels for the differentially expressed biomarkers; the fold change is unit difference in mean; a shift of units in the mean is considered. the underlying model is the logistic regression model with regression coefficients equal to two times the mean shifts for genes differentially expressed and zero for the rest. for each simulated dataset, we first randomly sample a training set of size 2/ × without consideration of the stratification of the binary outcome. the testing set is composed of the rest observations. we construct the proposed roc estimate based on the training set only. tuning parameters are selected using 3-fold cross validation. the predictive classification error and binormal auc are computed for the testing set based on the estimate obtained from the training set. summary statistics based on simulated datasets are given in table it can be seen that for moderate and large effect cases, the binormal aucs are close to with very small classification errors. the classification performance for small mean shift cases is less satisfactory, although still acceptable when the sample size is not too small, due to the relatively small signal to noise ratios. as sample size increases, the classification performance improves. the classification performance also improves with more differentially expressed genes. the empirical aucs are very close to the binormal aucs. results for the empirical aucs are omitted here. simulation study ii we consider data settings similar to those in simulation i, with π = and unit shift in mean for differentially expressed genes. we assume all biomarkers are independent and marginally distributed as uniform[-3 mathtype@mtef@5@5@+=feaafiart1ev1aaatcvaufkttlearuwrp9mdh5mbpbiqv92aaexatlxbi9gbaebbnrfifhhdyfgasaach8aky=wiffydh8gipec8eeeu0xxdbba9frfj0=oqffea0dxdd9vqai=hguq8kuc9pgc9s8qqaq=dirpe0xb9q8qilsfr0=vr0=vr0dc8meaabaqaciaacagaaeqabaqabegadaaakeaadagcaaqaaiabiodazawcbeaaaaa@2dbb@, 3 mathtype@mtef@5@5@+=feaafiart1ev1aaatcvaufkttlearuwrp9mdh5mbpbiqv92aaexatlxbi9gbaebbnrfifhhdyfgasaach8aky=wiffydh8gipec8eeeu0xxdbba9frfj0=oqffea0dxdd9vqai=hguq8kuc9pgc9s8qqaq=dirpe0xb9q8qilsfr0=vr0=vr0dc8meaabaqaciaacagaaeqabaqabegadaaakeaadagcaaqaaiabiodazawcbeaaaaa@2dbb@] gamma, and 22 mathtype@mtef@5@5@+=feaafiart1ev1aaatcvaufkttlearuwrp9mdh5mbpbiqv92aaexatlxbi9gbaebbnrfifhhdyfgasaach8aky=wiffydh8gipec8eeeu0xxdbba9frfj0=oqffea0dxdd9vqai=hguq8kuc9pgc9s8qqaq=dirpe0xb9q8qilsfr0=vr0=vr0dc8meaabaqaciaacagaaeqabaqabegadaaakeaadagcaaqaaiabikdayawcbeaaaaa@2db9@beta for the non-diseased group. the diseased group has the same marginal distributions with only mean shifts. all three marginal distributions have variance and represent three different distribution scenarios: uniformly distributed gene expressions; gene expressions with a skewed unimodal distribution; and gene expressions clustering around the upper and lower limits. proper transformations can lead to marginal normal distributions. we intentionally neglect the transformations and employ the proposed approach directly. the purpose is to investigate the performance of the proposed approach when marginal distribution assumption is violated. summary statistics based on simulated datasets are shown in table for comparison, we also reproduce the results from the normal distribution case in simulation i. it can be seen that the proposed approach has satisfactory classification performance under all scenarios. the effect of different marginal distributions is not significant. the empirical aucs are also reasonably close to the binormal aucs. simulation studies with different mean shifts and different choices of π lead to conclusions similar to those for simulation study i. those results are omitted here. simulation study iii in simulation studies i and ii, possible correlations among genes have been neglected. in simulation study iii, we investigate the case where the genes are correlated. we still assume π = and unit shift in mean for differentially expressed genes. for the non-diseased subjects, let zi zi..., i =..., nh, be a sequence of i.i.d. standard normal random variables. denote the d × nh dimensional gene expression matrix for healthy subjects as x. for the jth gene of the ith subject, we consider the following models: model 1: xij=h×m+1×m+kzil/k) mathtype@mtef@5@5@+=feaafiart1ev1aaatcvaufkttlearuwrp9mdh5mbpbiqv92aaexatlxbi9gbaebbnrfifhhdyfgasaach8aky=wiffydh8gipec8eeeu0xxdbba9frfj0=oqffea0dxdd9vqai=hguq8kuc9pgc9s8qqaq=dirpe0xb9q8qilsfr0=vr0=vr0dc8meaabaqaciaacagaaeqabaqabegadaaakeaacqwgybawdawgaawcbagaemyaakmaemoaaogabeaakiabg2da9iabdieainaabmaabawaaabmaeaacqwgabgwdawgaawcbagaemyaakmaemibawgabeaakiabc+cavmaakaaabagaem4aasgaleqaaaqaaiabdygasjabg2da9maabmaabagaemoaaomaeyoei0iaegymaedacagloagaayzkaagaey41aqraemyba0maey4kasiaegymaedabawaaewaaeaacqwgqbgacqghsislcqaixaqmaiaawicacaglpaaacqghxdatcqwgtbqbcqghrawkcqwgrbwaa0gaeyyeiuoaaogaayjkaiaawmcaaaaa@53ff@ with k = m = 7; h = 3 mathtype@mtef@5@5@+=feaafiart1ev1aaatcvaufkttlearuwrp9mdh5mbpbiqv92aaexatlxbi9gbaebbnrfifhhdyfgasaach8aky=wiffydh8gipec8eeeu0xxdbba9frfj0=oqffea0dxdd9vqai=hguq8kuc9pgc9s8qqaq=dirpe0xb9q8qilsfr0=vr0=vr0dc8meaabaqaciaacagaaeqabaqabegadaaakeaadagcaaqaaiabiodazawcbeaaaaa@2dbb@, where Φ is the cumulative distribution for the standard normal. model 2: same as model but with k = m = model 3: same as model with h as the identify function. model 4: same as model with h as the identify function. this simulation setting has been considered in. for models and the genes are marginally uniformly distributed. for models and the genes are marginally normally distributed. for genes in models and the correlations between all genes are weak; whereas models and are examples of strongly correlated genes. note that we not only assume differential genes are correlated, we also assume differentially expressed and non-differentially expressed genes are correlated. summary statistics based on replicates are shown in table for comparison, we also reproduce the results with independently, marginally uniformly distributed genes. we can see from table that the effects of correlation are negligible. we looked at other correlation structures and obtained similar results. moreover the effect of different marginal distributions are not significant. simulation study iv the binormal auc is theoretically constructed based on the joint normal distribution assumption on the linear combination of the covariates in the diseased and non-diseased population after a common monotone transformation. we now consider the case when marginal normality cannot be achieved with simple transformation. specifically, we investigate the case when the marginal gene distributions for the diseased class and the non-diseased class are different. we assume the number of genes is and there are subjects in each group. we assume that π = percent genes are differentially expressed. moreover, we assume the genes in the non-diseased group are independent and have standard normal distribution. for subjects in the diseased class, we assume the genes are independently distributed with: scenario 1: uniform [ - 3 mathtype@mtef@5@5@+=feaafiart1ev1aaatcvaufkttlearuwrp9mdh5mbpbiqv92aaexatlxbi9gbaebbnrfifhhdyfgasaach8aky=wiffydh8gipec8eeeu0xxdbba9frfj0=oqffea0dxdd9vqai=hguq8kuc9pgc9s8qqaq=dirpe0xb9q8qilsfr0=vr0=vr0dc8meaabaqaciaacagaaeqabaqabegadaaakeaadagcaaqaaiabiodazawcbeaaaaa@2dbb@, + 3 mathtype@mtef@5@5@+=feaafiart1ev1aaatcvaufkttlearuwrp9mdh5mbpbiqv92aaexatlxbi9gbaebbnrfifhhdyfgasaach8aky=wiffydh8gipec8eeeu0xxdbba9frfj0=oqffea0dxdd9vqai=hguq8kuc9pgc9s8qqaq=dirpe0xb9q8qilsfr0=vr0=vr0dc8meaabaqaciaacagaaeqabaqabegadaaakeaadagcaaqaaiabiodazawcbeaaaaa@2dbb@] for differential genes and uniform [-3 mathtype@mtef@5@5@+=feaafiart1ev1aaatcvaufkttlearuwrp9mdh5mbpbiqv92aaexatlxbi9gbaebbnrfifhhdyfgasaach8aky=wiffydh8gipec8eeeu0xxdbba9frfj0=oqffea0dxdd9vqai=hguq8kuc9pgc9s8qqaq=dirpe0xb9q8qilsfr0=vr0=vr0dc8meaabaqaciaacagaaeqabaqabegadaaakeaadagcaaqaaiabiodazawcbeaaaaa@2dbb@, 3 mathtype@mtef@5@5@+=feaafiart1ev1aaatcvaufkttlearuwrp9mdh5mbpbiqv92aaexatlxbi9gbaebbnrfifhhdyfgasaach8aky=wiffydh8gipec8eeeu0xxdbba9frfj0=oqffea0dxdd9vqai=hguq8kuc9pgc9s8qqaq=dirpe0xb9q8qilsfr0=vr0=vr0dc8meaabaqaciaacagaaeqabaqabegadaaakeaadagcaaqaaiabiodazawcbeaaaaa@2dbb@]otherwise; scenario 2: 22 mathtype@mtef@5@5@+=feaafiart1ev1aaatcvaufkttlearuwrp9mdh5mbpbiqv92aaexatlxbi9gbaebbnrfifhhdyfgasaach8aky=wiffydh8gipec8eeeu0xxdbba9frfj0=oqffea0dxdd9vqai=hguq8kuc9pgc9s8qqaq=dirpe0xb9q8qilsfr0=vr0=vr0dc8meaabaqaciaacagaaeqabaqabegadaaakeaadagcaaqaaiabikdayawcbeaaaaa@2db9@beta - 2 mathtype@mtef@5@5@+=feaafiart1ev1aaatcvaufkttlearuwrp9mdh5mbpbiqv92aaexatlxbi9gbaebbnrfifhhdyfgasaach8aky=wiffydh8gipec8eeeu0xxdbba9frfj0=oqffea0dxdd9vqai=hguq8kuc9pgc9s8qqaq=dirpe0xb9q8qilsfr0=vr0=vr0dc8meaabaqaciaacagaaeqabaqabegadaaakeaadagcaaqaaiabikdayawcbeaaaaa@2db9@ + for differential genes and 22 mathtype@mtef@5@5@+=feaafiart1ev1aaatcvaufkttlearuwrp9mdh5mbpbiqv92aaexatlxbi9gbaebbnrfifhhdyfgasaach8aky=wiffydh8gipec8eeeu0xxdbba9frfj0=oqffea0dxdd9vqai=hguq8kuc9pgc9s8qqaq=dirpe0xb9q8qilsfr0=vr0=vr0dc8meaabaqaciaacagaaeqabaqabegadaaakeaadagcaaqaaiabikdayawcbeaaaaa@2db9@beta - 2 mathtype@mtef@5@5@+=feaafiart1ev1aaatcvaufkttlearuwrp9mdh5mbpbiqv92aaexatlxbi9gbaebbnrfifhhdyfgasaach8aky=wiffydh8gipec8eeeu0xxdbba9frfj0=oqffea0dxdd9vqai=hguq8kuc9pgc9s8qqaq=dirpe0xb9q8qilsfr0=vr0=vr0dc8meaabaqaciaacagaaeqabaqabegadaaakeaadagcaaqaaiabikdayawcbeaaaaa@2db9@ otherwise; scenario 3: gamma + for differential genes and gamma - otherwise. under all three scenarios, the differentially expressed genes for the diseased subjects have mean and variance and non-differential genes have mean and variance based on replicates, the mean aucs are, and under scenarios 1– respectively, where the values in the parentheses are the corresponding standard errors. the mean classification errors are, and, respectively. colon and estrogen data data description colon data in this dataset, expression levels of tumor and normal colon tissues for human genes are measured using the affymetrix gene chip. a selection of genes with the highest minimal intensity across the samples has been made by, and these data are publicly available at. the colon data have been analyzed in several previous studies using other statistical approaches, see for example. estrogen data this dataset was first presented by. it includes expression values of genes of breast tumor samples. the expression data were obtained using the affymetrix gene chip technology and are available at. the response describes the status of the estrogen receptor. among the samples, are positive and are negative. we threshold the raw data with a floor of and a ceiling of genes with max/min < and/or max - min < are also excluded. genes pass the first step screening. a base logarithmic transformation is then applied. the estrogen data have also been studied by. estimation the roc estimate is identifiable only up to a scale constant. so prior to the analysis, we need to identify the "anchor biomarker," i.e, the gene whose estimated coefficient will be set as a constant for identifiability purposes. we select the anchor biomarker as follows. compute the sample standard errors of the d biomarkers se,..., se and denote their median as med.se. compute the adjusted standard errors as + med.se),..., + med.se). then the biomarkers are ranked based on the t-statistics computed with the adjusted standard errors. this adjusted t-statistic is similar to a simple shrinkage method discussed in. it has been observed in microarray studies that the standard t-statistic can be very large even if there is very little difference between the mean expressions, due to the small variance estimate caused by small sample size. so compared with the standard t-statistic, the adjusted t-statistic takes into account variability as well as fold change. the biomarker with the largest absolute value of the adjusted t-statistic is chosen as the anchor biomarker. for the anchor biomarker, if the sample mean of the diseased class is larger, β =, otherwise β = -. further investigation of the genes that have passed the first step screening shows that quite a few of those genes have very little variations across sample. including such genes in the joint modeling may create unreliable estimates. moreover, it is unlikely all the genes for the colon data will be associated with the outcome. so we only use the genes with the largest absolute values of the adjusted t-statistics for classification. the genes are then standardized to have zero means and unit variances. the same anchor marker detection and pre-processing have been considered in and are shown to be well-behaved. note that there is no computational limitation on how many genes can be used in the roc classification. empirical studies in show that the performance of certain statistical approaches decreases dramatically as the number of genes increases. so it is believed that identification of the top genes will help to generate more reliable estimate. we use three-fold cross validation for tuning parameter selection. we show in table the number of genes with nonzero coefficients and the cv scores for each fixed threshold value τ. it can be seen that generally the number of nonzero coefficients decreases as τ increases. however, the change of cv score is very small. our extensive empirical studies show that in general a larger τ will lead to a more parsimonious model. parsimonious models are preferred when the cv scores are comparable. so we choose τ = for both datasets. in the tgdr, variable selection is achieved by thresholding some estimated coefficients to be exactly zero based on their derivatives, so the corresponding covariates are not selected in the final model. the genes in tables and have nonzero coefficients and other genes have estimated coefficients zero and hence not selected. since the gene expressions have been normalized to have unit variances, the estimated coefficients are directly comparable. larger absolute values of coefficients indicate stronger influences. more detailed information about these genes can be found on the ncbi website. for example, for the colon data, hsa. is from the interleukin gene. this gene encodes a member of the cytokine family. expression of this protein is increased after the activation of t-cells by mitogens or the activation of nk cells by il- this protein induces the production of tnfalpha from macrophage cells. hsa. encodes a potent cyclin-dependent kinase inhibitor. the encoded protein binds to and inhibits the activity of cyclin-cdk or -cdk complexes, and thus functions as a regulator of cell cycle progression at gl. the expression of this gene is tightly controlled by the tumor suppressor protein p through which this protein mediates the p53-dependent cell cycle g l phase arrest in response to a variety of stress stimuli. this protein can interact with proliferating cell nuclear antigen, a dna polymerase accessory factor, and plays a regulatory role in s phase dna replication and dna damage repair. hsa. encodes a protein that is a member of the cxc chemokine family. this chemokine is one of the major mediators of the inflammatory response and is secreted by several cell types. it functions as a chemoattractant, and is also a potent angiogenic factor. this gene is believed to play a role in the pathogenesis of bronchiolitis, a common respiratory tract disease caused by viral infection. this gene and the other ten members of the cxc chemokine gene family form a chemokine gene cluster in a region mapped to chromosome 4q. for the estrogen data, d encodes a protein belonging to a small family of adapter proteins that are known to interact with a number of receptor tyrosine kinases and signaling molecules. this gene encodes a growth factor receptor-binding protein that interacts with epidermal growth factor receptor and ephrin receptors. the protein plays a role in the integrin signaling pathway and cell migration by binding with focal adhesion kinase. k encodes a member of the bombesin-like family of gastrin-releasing peptides. its preproprotein, following cleavage of a signal peptide, is further processed to produce either the aa gastrin-releasing peptide or the aa neuromedin c. these smaller peptides regulate numerous functions of the gastrointestinal and central nervous systems, including release of gastrointestinal hormones, smooth muscle cell contraction, and epithelial cell proliferation. these peptides are also likely to play a role in human cancers of the lung, colon, stomach, pancreas, breast, and prostate. u plays a role in a number of diseases. for example, it impacts the progression of hormone receptor-negative breast cancer. it is also expressed in cutaneous cd30+ lymphoproliferative disorders and is a candidate marker of disease progression. x is the estrogen receptor that is a ligand-activated transcription factor composed of several domains important for hormone binding, dna binding, and activation of transcription. evaluation we first randomly partition the data into training and testing sets. we then identify the top genes with the largest adjusted t-statistics based on the training sets only. so the genes used in each individual partition and computation may be different, although there is considerable overlap. the occurrence index, which measures the relative stability of individual selected genes in the joint model, is shown in figures and for the colon and estrogen data, respectively. we observe that many genes selected with the proposed approach have relatively high occurrence index compared to genes not selected, which suggests that the roc estimates are relatively stable for these two specific datasets. on the other hand, we also see that the occurrence index, even for the genes selected, can be low in absolute terms. most selected genes have occurrence index less than. similar phenomenon has been observed in. the genes identified in the final model do not necessarily have the largest occurrence index. for example, in the colon data the gene with the highest occurrence index was not selected in the model. this is not surprising since the objective function is the binormal auc and the genes were selected to optimize this objective function; the genes were selected based on all the available data, whereas the occurrence index was calculated using partitioned data. the occurrence index is one way to evaluate the relative importance of the genes after they are selected and it provides a measure of stability of the genes in the presence of sampling variation, but it does not measure the prediction power of the genes. furthermore, the genes significantly associated with the outcome in different subsets are not necessarily important when the whole dataset is used. thus we do not select genes based on the occurrence index. the classification performance of the proposed method for the two datasets is evaluated using the monte carlo approach described in the methods section. both the number of random partitions for the observed data and the number of permutations are for the colon data, the roc estimate has mean aucs and for the opd and ppd, respectively, where the values in the parentheses are the standard errors. corresponding mean classification errors are for observed data and for randomly permuted data, respectively. in the upper-right panel of figures and we also show the kernel density estimates of the opd and ppd of the auc. we can see that the two distributions are well separated. a simple wilcoxon test of the difference of these two distributions gives p-value <. for the same dataset and using the top ranked genes based on the marginal wilcoxon rank test, table in shows that the mean classification errors are,, and. in, genes are used and the svm based methods have mean aucs and mean classification errors, whereas the principal component analysis based approaches have even smaller aucs and larger classification errors. using the same set of genes, the logistic-lasso method in has mean aucs and mean classification errors; the smrc method in has mean aucs and mean classification errors. so for the colon data, the proposed binormal roc method performs better than other approaches in terms of auc and classification error, although the improvement is not dramatic since most other approaches behave reasonably well. we also note that since different sets of genes are used in, the above results can only provide a rough comparison. for the estrogen data, the proposed roc estimate has mean aucs and for the opd and ppd, respectively, with corresponding mean classification errors and. for the same dataset, dettling and buhlmann use genes selected based on the marginal wilcoxon tests and yield classification errors, and with different logitboost methods, classification error for adaboost, for 1-nearest-neighbor approach and by using the classification tree. using the same set of genes, the logistic-lasso has mean aucs and mean classification errors; whereas the smrc has mean aucs and mean classification errors. note that in our study, joint models are considered, i.e, we consider the joint effects instead of marginal effects of genes. identification of differentially expressed genes associated with clinical outcomes based on marginal significance has been well investigated and accepted. so as a simple evaluation, we show the marginal binormal auc and t-statistic for each individual gene in figures and similar evaluation has been considered in as a gene ranking method. we can see from figures and that several genes we identify do have large marginal binormal aucs and t-statistics; both up- and down- regulated genes are identified in the final models. however, we also identify some genes with small binormal aucs or small t-statistics. it is believed that these genes are important in the joint models, but not in the marginal models. variable selection using tgdr looks at all the variables and considers the relative importance of each variable in terms of its contribution to change in the objective function, but is not based on marginal significance. so the results do not necessarily correspond to those based on marginal significance. when binormal roc method is used, it is of interest to examine whether the roc curve is concave. pepe pointed out the binormal roc is not concave when the variances of the two groups are not equal. however, she also pointed that "swets and hanley conclude that the binomal roc curve provides a good approximation to a wide range of roc curves that occur in practice. this may be because the undesirable behavior of the approximating binormal roc usually occurs over a very small part of the roc curve." as a simple check, we show in figure the plots of binormal roc for a randomly selected training set and a testing set. we can clearly see the concavity of the binormal roc. examination of other training and testing sets shows similar results. we note that although the genes identified using the two datasets exhibit strong predictive power for disease statuses as demonstrated above, the analysis here does not provide information on whether they are just genomic markers correlated with the disease statuses or are actually in the pathways leading to the diseases. indeed, as is typical with microarray results, further investigation based on independent assays and/or independent samples is desirable.  the roc method has been successfully used for disease classification using low- dimensional biomarkers. the empirical auc method has been successfully implemented in. in this article, we propose using the binormal auc, which is a useful alternative to the empirical auc, as the objective function for disease classification using high dimensional genomic data. the scaled tgdr is employed for regularization and biomarker selection. extensive simulation studies and applications to two cancer data show that the proposed approach can generate satisfactory classification results with a small number of identified biomarkers with biologically sound implications. although theoretically the binormal auc is constructed under the normal distribution assumption, the proposed approach still behaves satisfactorily when this assumption is seriously violated. beyond the merits carried by the roc method, the proposed approach has the advantage of little computational cost. using the scaled gradient significantly reduces the computational complexity. it takes only a few seconds to analyze dataset with genes and subjects using r code, whereas the empirical auc method takes several minutes or longer. in our data analyses, we selected genes with the largest adjusted t-statistics. the proposed approach has no limitation on the number of genes that can be used in the final model. we note that use of genes is somewhat subjective. same analyses with genes and all the genes have also been carried out, and there are considerable overlaps among genes identified. we prefer using genes because this set of genes should include most important genes without including too much noise. we used the same datasets as in. although both methods are roc based and the same sets of genes are used, the identified biomarkers are different. considering that both the empirical auc and the binormal auc methods are extensively used, such difference warrants further investigation and confirmation by independent studies. in our study, the adjusted t-statistic is used instead of the standard t-statistic. the reason we used the adjusted t-statistic is that the unadjusted t-statistics can be large for genes with small mean differences in expression levels, if the standard error is very small. because of a large number of genes on each array, there will usually be genes with very small standard errors, and some of these genes will have small means as well. this phenomenon has been discussed by in motivating their adjusted t-statistic for significance analysis of microarrays. similar variance adjust has also been considered in. however, a recent study points out that such variance adjust may not be necessary. since current study is not conclusive, we use the adjusted t-statistic and leave the choice of t-statistic versus adjusted t-statistic to a future study. in this paper, we proposed monte carlo based methods for evaluation of prediction and stability of individual genes. the proposed methods are closely related to cross validation and permutation methods. the problem of evaluation is extremely important for high dimensional classification, but it has not been extensively studied. although the proposed methods have sound, intuitive implications, further studies of their theoretical properties would be of interest. we focus on the linear combination of genes that maximizes the binormal-aug. in, linear combinations of genes favoring certain specificity or other criterion are investigated. if maximization of specificity is of interest, we can then combine liu's methods with certain variable selection techniques, for example lasso or the proposed tgdr. we leave this to a future study. we only discussed two-sample classification in this article. when extension of the binormal auc to multi-class study exists, the proposed estimation and regularization method can be used. we postpone multi-class study to a future report. 
we propose using the binormal auc as the objective function for two-sample classification, and the scaled threshold gradient directed regularization method for regularized estimation and biomarker selection. the roc technique has been widely used in disease classification with low dimensional biomarkers. we develop monte carlo based methods for evaluating the stability of individual biomarkers and overall prediction performance. application of the method to two cancer studies shows that the identified genes are reasonably stable with satisfactory prediction performance and biologically sound implications. extensive simulation studies show that the proposed approach can generate parsimonious models with excellent classification and prediction performance, under most simulated scenarios including model mis-specification. tuning parameter selection is based on v-fold cross validation. it is also of interest to develop methods for evaluation of classification performance and ranking of identified biomarkers. the overall classification performance is satisfactory, with small classification errors and large aucs. compared with the empirical roc approach, the binormal roc is computationally more affordable and robust in small sample size cases.
 an important application of microarrays is to discover genomic biomarkers, among tens of thousands of genes assayed, for disease diagnosis and prognosis. thus it is of interest to develop efficient statistical methods that can simultaneously identify important biomarkers from such high-throughput genomic data and construct appropriate classification rules. it is also of interest to develop methods for evaluation of classification performance and ranking of identified biomarkers.  the roc technique has been widely used in disease classification with low dimensional biomarkers. compared with the empirical roc approach, the binormal roc is computationally more affordable and robust in small sample size cases. we propose using the binormal auc as the objective function for two-sample classification, and the scaled threshold gradient directed regularization method for regularized estimation and biomarker selection. tuning parameter selection is based on v-fold cross validation. we develop monte carlo based methods for evaluating the stability of individual biomarkers and overall prediction performance. extensive simulation studies show that the proposed approach can generate parsimonious models with excellent classification and prediction performance, under most simulated scenarios including model mis-specification. application of the method to two cancer studies shows that the identified genes are reasonably stable with satisfactory prediction performance and biologically sound implications. the overall classification performance is satisfactory, with small classification errors and large aucs.  in comparison to existing methods, the proposed approach is computationally more affordable without losing the optimality possessed by the standard roc method. 