in cancer diagnosis , cancer patients with the same diagnostic profile may have different clinical outcomes .
the difference probably lies in the limitation of the traditional classification of tumor types , based mainly on morphology .
a reliable and precise classification of tumors is essential for successful diagnosis.1 modern sequencing and microarray technology have enabled more detailed molecular characterization of cancer samples , leading to the discovery of many cancer subtypes .
in fact , cancer subtype identification has become an integral part of personalized medicine.1 traditionally , the problem of cancer subtype classification has been approached in many ways such as multinomial logit models,2 bayesian probit models,3,4 random forest,5 and support vector machine ( svm).69 other discriminatory methods , including linear discriminant analysis ( lda ) , k - nearest - neighbor ( knn ) classifier , and classification trees , were also investigated.10 among these , svm is a successful procedure applied to microarray - based cancer diagnosis problems.2,11 however , it suffers from the difficulty to interpret the resulted model as well as no direct outcome probability estimates produced.1214 multinomial logit model , on the other hand , is endowed with feature interpretability and probabilistic nature.15 classic multinomial logit model works well when the number of predictors is small . as the number of predictors increases , the generalization power of the model deteriorates because of model overfitting .
when the number of predictors exceeds the number of observations as is commonly seen in genomic studies , the method breaks down . to deal with the curse of high dimensionality as well as to increase the model interpretability ,
regularized procedures that incorporate a sparsity penalty have been proposed.1620 among these methods , group lasso is particularly appropriate for models with multiclass responses , in which all the coefficients linked to a common predictor form a group and are required to shrink to zero simultaneously in order to achieve the selection of predictors.16 although the sparse multinomial logit models are capable of selecting variables , they can not efficiently use underlying structure information such as a network of regulatory relationships between genes or gene - products .
such structure information could be obtained by a data - driven approach via clustering21 or other similarity - based methods,22,23 or be extracted from the external databases accumulated through years of biomedical research .
databases such as kegg , reactome , and mips have been developed to organize different types of biological network information .
cancer is a complex disease caused by dysregulation of pathways instead of a single gene.2426 thus , the incorporation of the network information can potentially increase the power of identifying cancer subtypes .
networks are often represented as graphs , with each vertex indicating a gene or a gene - product and each edge indexed by a relationship between two vertices .
a constraint , induced by the laplacian matrix of the graph,27 has been proposed to facilitate the selection of predictors in ordinary regression settings , enhancing both the global smoothness over the network and the interpretability of the association between selected genes and responses in the context of known biology . in this paper
, we propose a network - constrained sparse multinomial logit model for high - dimensional multinomial classification , aiming at improving prediction performance by utilizing the underlying network prior .
the remainder of this article is organized as follows . the model is first presented , followed by the algorithm to fit the model .
finally , the model is applied to a real data set of predicting the subtypes of glioblastoma multiforme ( gbm ) .
n with n observations and p predictors , yi denotes an observation of the categorical response variable y  { 1 ,  , k } and xi = ( xi1 , xi2 ,  , xip )  r indicates an observation of a p - dimensional vector of predictors . assuming that yi follows a multinomial distribution , a multinomial logit model is built with logit link , which is 
 ir = p(y = r|xi)=exp(r0+xir.t)s=1kexp(s0+xis.t)=exp(ir)s=1kexp(is)(1)where 
 r.=(r1,r2, 
we choose category k as the reference category by setting k0 = 0 and k .
k1 , correspond to the log odds ratio between category r and the reference category k. we regularize the multinomial logit model using a penalized likelihood approach , in which one maximizes the penalized log - likelihood 
 lp()=l()j(),(2)over a ( k  1 )  ( p + 1)-dimensional parameter vector  = ( 10 ,  , (k1)0 , 1 . ,  , k1 , . ) . in equation
( 2 ) , 
 l()=i=1nr=1kyirlogir=i=1n(r=1k1yirirlog(s=1kexp(is)))is the ordinary log - likelihood of a multinomial logit model , and j (  ) is a function penalizing the magnitude of the parameters and regularizing the structure of features .
 , the regulatory parameter , controls the strength of the regularization . assuming that all predictors are metric and standardized , that is , each predictor has one degree of freedom , and the differences in scale will not influence the penalty and , thus , the variable selection . in the multinomial logit model , we use a vector .j
= ( 1j , 2j ,  , k1j ) of parameters to capture the effect of variable xj , so that variable selection is achieved only when the k  1 parameters shrink to zero simultaneously .
since the ordinary lasso facilitates only parameter selection rather than predictor selection , a group lasso penalty is utilized to penalize the parameters at a group level , defined as 
 j1()=j=1pj.j=j=1pj(1j2+2j2++k1j2)1/2,(3)where j is a penalty weight , set as 1 by default . in group lasso
, all the parameters in a group .j would shrink to zero simultaneously . in an association study , the graphs or networks depicting relationships among predictors are informative , supplementary to numerical data . consider a network represented by a weighted graph g = ( v , e , w ) with the set of vertices v = 1 , 
p corresponding to p predictors , the set of edges e = { ( j , k ) : j and k are linked } , and the set of weights w = { wjk : ( j , k )  e}. wjk measures the level of concordance of predictors j and k , with 1 for identity and 0 for complete difference , if normalized to the scale of .
we then construct an adjacency matrix a by 
 ajk={wjk(j , k)e,0(j , k)e.and a degree matrix d = diag(d1 , d2, 
,dp ) , where dj = (j , k)e 
 jk is defined as the degree of vertex j. the laplacian matrix associated with graph g is l = d  a , which is always non - negative definite and can be factorized as l = ss . by simple algebra , 
 r.lrt can be written as 
 r.lr.t=(j,k)e(rjrk)2wjk .
thus , the network - constrained penalty,22,27,28 defined as 
 j2()=r=1k1r.lr.t,(4)induces a smooth solution of the vector r . with respect to the labeled weighted graph g. typically , the adjacency matrix can be constructed from external information using the abovestated method .
however , several data - driven methods are also applicable.22,23 adjacency coefficients can be defined based on similarity measures such as pearson correlation coefficient , which are transformed into adjacency by a monotonically increasing function .
the most widely used transformation functions include the signum function , the sigmoid function , and the power function .
detailed examples of adjacency measures are provided by huang and ma.22 to sum up , in our regularized model , the penalized log - likelihood function is given by 
 lp()=l()j()(5 ) 
 = i=1n(r=1k1yirirlogs=1keis)1j=0pj.j2r=1k1r.lr.t,(6)of which the second term , the sparse penalty , induces model sparsity and the third term , the network penalty , imposes smoothness over the network . when 2 = 0 , the model reduces to the group lasso multinomial logit model .
the incorporation of this extra tuning parameter expands the parameter search space and directs the search to more biological meaningful regions .
like ordinary lasso , group lasso also suffers from an issue of estimation bias , which is resulted from the fact that all predictors are penalized to the same degree . in order to reduce the bias
, we use adaptive group lasso , which penalizes predictors to different degrees by assigning a weight to each predictor . in our model ,
the weight is set to be the reciprocal of the l2 norm of the fitted coefficients in univariate analysis , where we fit the model with each individual predictor only . denoting 
 .j as the univariate estimate , the group lasso penalty term ( 3 ) becomes 
 j1()=j=1p1.j.j.
we use the proximal gradient - based fista ( fast iterative shrinkage - thresholding algorithm ) to fit the model.29,30 consider the optimization of the general penalized log- likelihood lp (  ) = l( )  1j1( ) , composed of a concave and continuously differentiable term l( ) , and a convex penalty term j1( ) . the penalized maximum likelihood ( ml )
estimator is defined by 
 ^=argmaxrdlp()=argmaxrd(l*()+1j1()),(7)where 
 l*()=i=1n(r=1k1yirirlogs=1keis)2r=1k1r.lr.tis a smooth function with respect to parameter . with a positive step size v , the quadratic approximation29 of lp( ) at a given point 0 is 
 qv(,0)=l*(0)l*(0)t(0)+12v02+1j1( ) .
l( ) , the first - order derivative of l( ) , is a ( k  1 )  ( p + 1)-dimensional vector , whose element corresponding to rj is 
 [ l*()]rj=l*rj=i=1n(yirir)xij+2lj.r.t .
the iterations of proximal gradient methods are defined by 
 ^(t+1)=argminrd{l*(^(t))l*(^(t))t(^(t))+12v^(t)2+1j1()}(8)which consists of a linear approximation of the negative modified log - likelihood at the current value 
 ^(t ) , a proximity term , and the penalty term .
first , we set 1 = 0 , and based on the standard formula for the iterates of gradient methods for smooth optimization , the unpenalized estimator 
 (t+1 ) has an explicit form 
 (t+1)=(t)+vl*((t ) ) .
then we move back to the optimization problem with an active penalty . via lagrange duality ,
equation ( 7 ) can be equivalently expressed by 
 ^=argminc(l*()),where 
 c={rd|j1()(1 ) } is the constraint region corresponding to j1( ) and (1 ) is a tuning parameter that is linked to 1 by a one - to - one mapping . given a search point u  r , the so - called proximal operator associated with j1( ) is defined as 
 p1(u)=argminrd(12u2+1j1(.(9)which is the projection of u onto region c. then the proximal gradient iterates defined in equation ( 8) can be equally expressed by the projection 
 ^(t+1)=p1v(^(t)+vl*(^(t ) ) ) .
next , consider the proximal operator ( 9 ) . owing to the block - separability of this specific penalty , the proximal operator can be written as 
 p1()=(p1(^.0),p1(.1), 
,p.(10 ) with ( u)+ = max ( u , 0 ) , the explicit solution to the proximal operator ( 10 ) can be derived from the karush 
set 
 (t+1)=^(t)+vl*(^(t ) ) , then the solution to the optimization problem ( 8) can be expressed as 
 ^(t+1)=p1v((t+1))=(p1v(.0(t+1)),p1v(.1(t+1)),  ,p1v(.p(t+1)))=(.0(t+1),(11j.1(t+1))+.1(t+1), 
,(11j.p(t+1))+.p(t+1 ) ) to summarize , the basic idea of proximal gradient methods is as follows : first , remove the l1 penalty of the objective function ( 6 ) and then optimize the smooth part by taking a step toward its ml estimator via first - order methods , which creates a search point .
second , project this search point onto the constraint region c in order to account for the non - smooth penalty term . to accelerate the convergence rate
, we extrapolate the current and the previous iterations with the help of deliberately chosen acceleration factors at,20 
 ^(t)=^(t)+at11at(^(t)^(t1 ) ) .
the extrapolate point 
 ^(t ) , instead of the current iterate 
 ^(t ) is used as a starting point to generate a search point , which is then projected on the penalty region . to select the tuning parameters 1 and 2
, we use cross - validation , where we divide the data set into training and test data sets .
the model is trained on the training data set , and prediction error is then assessed on the test data set .
we search on a grid of 1 , 2 values and choose the value of 1 , 2 that minimizes the cross - validated errors .
the prediction error is measured with the brier score , a measurement of the accuracy of probabilistic predictions defined as the euclidean distance between sample response and its estimated probabilities .
the purpose of the simulation is to show that the structure - constrained model dominates the alternative models that do not use such prior information in terms of parameter estimation and prediction . for each scenario presented ,
we first select the optimal tuning parameters through a five - fold cross - validation on the training set . with the selected tuning parameters ,
a final model is built on the whole training set and then tested on the test set . for each setting
, we run 50 simulations and calculate several criteria to evaluate the performance of the proposed model .
the numbers of total predictors are 20 for small and 200 for large models , and the numbers of relevant ones are 4 and 10 , respectively .
the predictors are continuous and follow a multivariate normal distribution with mean 0 and the p  p correlation matrix 
 =(12  p11  p221  p3p1p2p3  1),where  = 0.5 .
we divide the predictors into a few subsets ( subnetworks ) . in the small model ,
20 predictors are divided into 5 subnetworks evenly , and the 4 relevant predictors constitute the first subnetwork . in the large model ,
200 predictors are divided into 20 subnetworks evenly , and the 10 relevant predictors form the first subnetwork . ideally , we assume full connection within each subnetwork and no connection between them ; that is , the corresponding adjacency matrix is a block diagonal matrix , with the main diagonal blocks being all - ones square matrices and the off - diagonal blocks zero matrices .
we then construct the coefficient matrix , a 3  p matrix whose rows are indexed by all but the reference categories of the response and columns are indexed by predictors .
the columns corresponding to the irrelevant predictors are filled with zeros . for the relevant columns ,
the first setting requires equal coefficients in each category for relevant predictors . in the case of similar coefficients , entries in each category
share the same sign but have different values , indicating that all the relevant predictors impact the response in the same direction but with different magnitudes .
their absolute values are independently drawn from the set { 0.05 , 0.10 ,  , 0.50 } , and the sign for each category is random .
random coefficients are independently drawn from the set { 0.50 , 0.45 ,  , 0.05 , 0.05 ,  , 0.50}. in this case , the prior structure information is not useful , since we assume that the coefficients within each subnetwork should be at least similar .
this scenario serves as an example of model misspecification . to further study the effects of structure misspecification
, we also test our models in cases of incorrect networks and overlapping networks . in incorrect
network setting , the large and small adjacency coefficients are randomly drawn from ( 0.4 , 1 ) and from ( 0 , 0.6 ) instead of being constant 1 and 0 . in the overlapping scenario ,
based on the multinomial logistic model , the actual probabilities can be derived and the class label is then randomly drawn from a multinomial distribution for each observation .
in addition , the intercept is set to zero for the sake of a relatively balanced design .
to see the improved performance of using prior structure information in terms of parameter estimation and prediction accuracy , we compare the variants of the proposed model , network - constrained multinomial logit model with group lasso penalty ( ngl - mlm ) and the one with adaptive group lasso penalty ( ngl - mlma ) , to two traditional multinomial logit models with lasso ( l - mlm ) and group lasso ( gl - mlm ) , respectively , implemented in the package of glmnet in r.16 to measure the estimation accuracy , the mean - squared error ( mse ) between true parameter values and the estimated ones is used . in addition , the performance of prediction on test data is evaluated with brier score , the euclidean distance between sample response and the estimated probabilities , and prediction accuracy , the proportion of correctly predicted class labels .
we first simulate ideal network structure ; that is , all the relevant variables come from a fully connected subnetwork .
as expected , the structure information improves estimation significantly , especially for large models , which is particularly relevant for real applications .
the estimation of the adaptive method ( ngl - mlma ) outperforms others substantially . in case of random coefficients , where prior network does not provide any useful information ,
the proposed model is comparable to models without using the network information ( l - mlm , gl - mlm ) , and sometimes even better .
figure 2 shows that the prediction accuracy is also higher for the proposed model in almost all scenarios . when brier score is used ( fig .
3 ) , a similar trend follows : the network - constrained model always performs better when we simulate ideal and similar coefficients , and is comparable to traditional models without using structure information in case of random coefficients . to investigate the impact of structure misspecification
we simulate a medium - sized data set with 100 predictors , 10 being relevant .
the performance of our models is still satisfactory because of the flexible tuning parameter on the structure - constraint term ( fig .
, the prediction accuracy of ngl - mlm is comparable to that of gl - mlm in both situations , whereas , in terms of parameter estimation and brier score , the structure - constrained models ngl - mlm and ngl - mlma outperform the other two . in summary ,
our structure - constrained multinomial logit model has better performance in terms of parameter estimation and prediction when the prior network knowledge is at least partially correct , and the performance is comparable to traditional models when the network knowledge is incorrect .
this is because the gl - mlm is a special case of ngl - mlm with 2 = 0 .
cross - validation tends to select 2 = 0 when the prior assumption is not correct .
one important application of our method is cancer subtype prediction and relevant subnetwork identification on large - scale gene expression data .
we apply all four candidate methods , l - mlm , gl - mlm , ngl - mlm , and ngl - mlma , to a large - scale tcga ( the cancer genome atlas ) gbm subtype prediction problem , which contains the expression data of 11,861 genes across 202 samples and four categories , with 40 , 46 , 58 , and 58 samples in each category .
the network was built from a variety of sources , including reactome , kegg , as well as the inferred gene - interactions from protein interactions , gene co - expressions , protein domain interactions , and text - mined interactions .
the outcome is one of the four subtypes of gbm.31 the data set , the network information , as well as the subtype information were downloaded from http://bioen - compbio.bioen.illinois.edu / ncis/. since the number of genes in the gbm data set is much larger than the number of samples , which may lead to computation instability , we carry out gene screening before analysis .
starting with 11,861 genes , we screen genes based on the prior weights resulting from the ncis algorithm,31 by including the 599 most highly weighted genes . to construct the network smoother for model building , we tailor the original network subject to the remaining 599 genes .
then the laplacian matrix is constructed based on the tailored subnetwork . to compare the prediction performance of the four methods ,
202 samples are randomly divided into two subsets , a training set of 150 samples and a test set of 52 samples .
the random division is kept only when test samples have good representation of each category ( 1535% for each category ) . otherwise , we discard that division . feature selection and parameter estimation in model building
are completed strictly on the training set , and then the fitted models are tested on the test set . in practice , 50 valid divisions are obtained , and model building and testing are carried out on all pairs of data sets to assess variability , the results being summarized in table 1 .
the tuning parameter of the network - constraint controls the impact of the prior structure knowledge on model building .
the network information will have no effect if the tuning parameter is set to zero . among the 50 models built by ngl - mlm
, the network tuning parameter is chosen as zero in 28 models , whereas ngl - mlm is reduced to gl - mlm in this specific case .
in contrast , 48 ngl - mlma models choose non - zero tuning parameter for the network- constraint , which indicates that the structure knowledge is useful for prediction , explaining the higher prediction accuracy rate for ngl - mlma .
next , we apply ngl - mlma , the best model in both simulation and the application to gbm subtype analysis , on the whole sample set of gbm gene expression data and investigate the selected subnetwork .
it selects 35 predictors , among which 21 are non - singletons and form a subnetwork , shown in figure 5 .
for example , the most connected gene akt1 plays an important role in the pathogenicity of gbm .
akt1 is a downstream serine / threonine kinase in the rtk / pten / pi3k pathway , and large - scale genomic analysis of gbm has demonstrated that this pathway is mutated in many but not all gbms.32 therefore , the akt1 can be potentially used to define gbm subtypes .
cancer subtype prediction is of critical importance in the understanding , diagnosis , and treatment of cancer .
we introduced a classification model on the basis of multinomial logit regression to identify cancer subtypes from high - throughput gene expression data .
the group lasso penalizes the coefficients linked to a common predictor at a group level so that it facilitates variable selection by shrinking all elements within a group to zero simultaneously .
in addition , the network - constraint improves the smoothness of coefficients with respect to the prior structure information and results in more interpretable identification of genes and subnetworks .
the proposed model and its adaptive extension are compared to lasso and group lasso multinomial logit models with no network - constraint involved . from the results of simulation and the application to gbm gene expression data , we find that the proposed model is superior given correct prior network information and is comparable to traditional models given incorrect network information .
a key challenge to the future study is to correctly specify the networks . in the application to real data
, we may include too many misspecified edges on the network because of incomplete knowledge of pathways .
one possible solution is to use problem - specific network for a particular type of cancer , rather than using a general molecular interaction network .
the proposed method can be extended by using a non - convex sparsity penalty to reduce estimation bias .
scad ( smoothly clipped absolute deviation ) and mcp ( minimax concave penalty ) are two potential candidates.3335 the applications of the method also go beyond the cancer subtype prediction .
for example , it can also be used to predict the disease subtype based on the human microbiome data , where the phylogeny structure can be efficiently used.28,36,37
 classic multinomial logit model , commonly used in multiclass regression problem , is restricted to few predictors and does not take into account the relationship among variables . 
 it has limited use for genomic data , where the number of genomic features far exceeds the sample size . 
 genomic features such as gene expressions are usually related by an underlying biological network . 
 efficient use of the network information is important to improve classification performance as well as the biological interpretability . 
 we proposed a multinomial logit model that is capable of addressing both the high dimensionality of predictors and the underlying network information . 
 group lasso was used to induce model sparsity , and a network - constraint was imposed to induce the smoothness of the coefficients with respect to the underlying network structure . to deal with the non - smoothness of the objective function in optimization 
 , we developed a proximal gradient algorithm for efficient computation . 
 the proposed model was compared to models with no prior structure information in both simulations and a problem of cancer subtype prediction with real tcga ( the cancer genome atlas ) gene expression data . 
 the network - constrained mode outperformed the traditional ones in both cases . 
Introduction
Multinomial Logit Model and Penalized likelihood Approach
Proximal Gradient Method and Model Fitting
Simulation
Simulation Settings
Simulation Results
Application to the GBM Data Set
Conclusion and Discussion