BACKGROUND
in recent years, biomedical data have become complicated and high-dimensional  <cit> .
for example, a single human gene expression dataset contains tens of thousands of features, many of which are highly correlated  <cit> .
in addition, large mixed datasets are crucial for personalized treatment, in which the optimal treatment strategy is determined based on a dataset that combines a very large number of prognostic factors  <cit> .
from the viewpoint of statistical machine learning, supervised and unsupervised learning methods play central roles in such biomedical studies  <cit> .
in fact, shrinkage methods such as ridge and lasso are frequently used in the context of prediction  <cit> , and clustering methods are used in the context of interpretation  <cit> , potentially revealing novel findings.
supervised learning methods are often used to estimate risk scores when predicting dichotomous outcomes.
linear scores are among the most widely used forms because it is easy to learn the predictive score from a training dataset.
moreover, it is easy to understand the estimated score.
linear scores are often evaluated by linear discriminant or logistic regression analysis, and achieve not bad discriminative performance.
for example, the study of  <cit>  used a linear score, and their discoveries led to the development of mammaprint, a diagnostic kit for breast cancer metastasis.
unsupervised learning methods can also yield beneficial insights in high-dimensional data analysis.
for example,  <cit>  used biclustering to reveal more detailed subtypes in breast carcinomas with distinctive gene expression profiles from a group that was previously regarded as a homogenous unit.
their study revealed that in order to understand biodiversity, the heterogeneous structure of the targeted population must be considered, and that such heterogeneity can be clarified by the clustering method.
previously published reviews have described both clustering methods  <cit>  and biclustering algorithms  <cit> .
several studies have combined supervised and unsupervised learning methods.
for example,  <cit>  used clustering to discover different patterns of gene expression in different subgroups.
they then derived the respective scores for these groups and achieved good specificity without loss of sensitivity relative to existing diagnostic rules.
sample heterogeneity may result in marker heterogeneity.
as a result, different samples in different subgroups may have different intrinsic characteristics in their environmental and genetic factors as demonstrated by the motivated example in the “methods” section.
such heterogeneity may have unexpected effects on a therapy or treatment which is considered as best practice, and lead to an unfavorable risk in one part of the population.
bravo et al.
<cit>  focused on the marker heterogeneity by detecting the genes that showed different variation between healthy and disease samples.
they then defined an anti-profile score as the number of hyper-variable genes.
thus, more and more studies have considered heterogeneous structure and reflected this heterogeneity in their predictions.
however, the risk scores highlighted by published papers are linear, and heterogeneity is therefore not directly reflected in the score form.
in this study, we focused on heterogeneity and determined how to directly reflect this intrinsic characteristic in the score form.
we developed the quasi-linear score as a result, which combines linear scores as a kolmogorov-nagumo average  <cit> , enabling us to reflect the clustering result naturally, because it is based on separated feature vectors.
the rest of this paper is organized as follows.
in the “methods” section, we first present a motivated example of gene expression data and develop the quasi-linear score.
heterogeneity is observed via the clustering method, and we define the quasi-linear score to reflect gene clusters with a generalized average form.
we also formulate the quasi-linear logistic model and discuss the difference between the linear and quasi-linear scores.
we subsequently evaluated our method by numerical simulations and applications to real datasets.
we refer to the relationship between the quasi-linear score and traditional combined approaches in “discussion” section.
all technical details given as appendix are available in additional file  <dig>
methods
motivation and derivation
we studied the gene expression dataset from  <cit> .
this dataset is derived from  <dig> non-metastatic and  <dig> metastatic breast cancer patients.
in their study, the linear score was evaluated to discriminate metastatic events.
because estimation of the predictive linear score is often achieved by a diagonal fisher’s linear discriminant analysis   <cit> , we considered applying dlda to this dataset.
because the coefficients of the linear score estimated by dlda correspond to the t-statistic values, we checked the t-statistics directly for the purpose of visualization.
if the data have heterogeneous structure, it can be clarified by observing the difference between two divided, independent datasets.
therefore, we divided the full data into two independent sets, data <dig> and data <dig>  before calculating the t-statistics for each of them separately.
figure  <dig> shows the correspondence of the t-statistics.
some genes had no consistency in the signs of their t-values, indicating that some samples from the metastatic group had higher expression, whereas other samples had lower expression, relative to the non-metastatic group.
this phenomenon may be caused by heterogeneous factors  <cit> .
in fact, due to the existence of multiple subtypes of breast cancer, this disease is known to exhibit heterogeneity  <cit> .
for such heterogeneous data, clustering methods should work well, as shown in  <cit> .
we applied clustering according to the ward’s method  <cit> , as shown in fig.
<dig>  which highlights the results of clustering and the correlation matrix arranged by the clustered genes.
although biclustering result was not suggestive of the heterogeneity in appearance, it was observed via the correlation matrix.
some genes are strongly correlated with others in the same cluster.
thus, we observed the existence of heterogeneity using a t-statistics plot and trends in the expression patterns by clustering.
next, we developed an appropriate score form for discriminating such heterogeneous data based on clustering.
fig.
<dig> t-statistic values for two datasets from van’t veer et al.
.
the red points show the genes with sign mismatched t-values for these data
fig.
<dig> the hierarchical clustering and the correlation matrix of  <dig> genes for the dataset from van’t veer et al.
.
the figure shows the clustering result  and the correlation matrix .
there are  <dig> rows representing genes and  <dig> columns representing samples  and the gene expression data ranging from green  to red  are displayed.
elements of the correlation matrix  ranging from blue  to yellow  are displayed
we assume to know decomposition of p biomarkers into k groups by clustering.
based on these k sets of clustered markers, we define a quasi-linear score as
<dig> q=log∑k=1kexp,
where lk=αk+βk⊤x with the parameters α
k,β
k, and the marker vector x
for the k-th cluster of k= <dig> ,⋯,k.
when k= <dig>  the quasi-linear score q is reduced to the linear score,
<dig> l=α+β⊤x,
where α and β=β1⊤,⋯,βk⊤⊤ are parameters, and x is the full vector of x
,⋯,x
.
we note that the intercepts α
k’s in  are reduced to the single intercept α in .
additional file 1: appendix a gives another parameterization for q in which the intercepts α
k are uniquely decomposed to the overall intercept and weights of the k clusters.
when determining how to reflect cluster information in the score form, we had two main considerations: which scores should be integrated, and how integration should be performed.
all the linear scores l
k are integrated in the quasi-linear score q. we believe that this is reasonable because there are similar markers in each cluster, and we expect that heterogeneity would be caused by different mixed homogeneous features that are sufficiently described by the linear form.
although such an idea of combining the linear scores has already been proposed by  <cit>  as a composite link, it is different from the quasi-linear score in several ways.
one of the most significant differences between them is that the quasi-linear score is defined by disjointed sets of markers.
this results in a small number of parameters for the predictive score: the parsimonious expression.
the difference in these forms is mentioned in the discussion.
moreover, the quasi-linear score q summarizes l
k approximating the maximum function.
in fact,  is equal to a soft maximum function discussed by  <cit> , which can be approximated with
<dig> m=max1≤k≤klk.
therefore, the quasi-linear score q respects the maximum of k linear scores from all clusters.
see  <cit>  for a discussion of eq.
as maxout in neural networks.
the relationships among q, l, and m are clearly evaluated when a tuning parameter, τ, is introduced in the quasi-linear score q as
<dig> qτ=1τlog∑k=1kexp,
where τ is a positive parameter.
if l
k is fixed for all k, then the form of  is defined solely by the tuning parameter τ. when τ is equal to  <dig>  q
τ is equal to the quasi-linear score q by definition.
when τ goes to infinity, q
τ is simply the maximum score m. when τ goes to  <dig>  q
τ is equivalent to the linear score l. thus, these are unified by the hardness of approximation to the maximum function.
more details are provided in additional file 1: appendix b. the characteristics of the quasi-linear score q are understood by a more general expression of  : g=ϕ−1), where ϕ is an invertible function.
we define g as a generalized quasi-linear score because the form is the generalized mean called the kolmogorov-nagumo average.
if we take a simple average, ϕ=z, then the generalized quasi-linear score g corresponds to the linear score l. in this sense, the linear score l is a simple mean, and the quasi-linear score q is a generalized mean of linear scores l
k averaged by the exponential function.
although the simplest integration of clustered information is achieved by a simple average, resulting in the linear score form, it is intuitively unsatisfying because the predictive performance of these linear scores l
k differs among the clusters.
a cluster that strongly discriminates the outcome on its own should be respected in comparison with the other clusters.
if only the cluster with the highest linear score is reflected in the prediction, it is described by the maximum score m. however, this situation is still not ideal, because only one cluster is reflected in the prediction, and form  is difficult to handle mathematically because it is not differentiable when two or more linear scores are equal.
consequently, parameter estimation becomes impossible.
the quasi-linear score q is therefore naturally derived, and it is reasonable for discriminant analysis of the heterogeneous data because the quasi-linear score q plays an important role in cluster selection, as discussed later.
in the following sections, we let ϕ=exp, both because this form is approximated by the maximum function, and because it is optimal in the sense of bayes risk consistency when we consider the simple case in which the label conditional random variables follow a mixture of normal distribution and a normal distribution with equal variance, respectively.
additional file 1: appendix c provides more detail about the bayes risk consistency of the situation.
moreover, the exponential function gives us an understandable interpretation of the parameter estimation.
because we modify only the scoring form, the quasi-linear score q can be applied to all traditional settings in biostatistics, as the generalized linear model with the l <dig> and l <dig> shrinkage methods.
in particular, when we combine the quasi-linear score q with lasso shrinkage, the important clusters and variables in each cluster are determined simultaneously because of soft maximum property and l <dig> sparseness.
this property provides good performance for the discriminant problem when the data have a much larger number of correlated markers than the number of samples.
therefore, we derive the l <dig> and l <dig> shrinkage quasi-linear logistic model and display the performance of the quasi-linear score q when it is applied to gene expression data.
likelihood for logistic model and maximum likelihood estimation
consider the data {;i= <dig> ⋯,n}, where x
i is a covariate vector and y
i is a dichotomous outcome which takes  <dig> or  <dig> with the i-th individual.
assume that we know the decomposition of x
i as x
i,⋯,x
i with a fixed cluster size k, and that this is identical among individuals.
we denote the size of x
i as p
k, where ∑k=1kpk=p.
we note that the decomposition is given a priori by one of clustering methods for {x
i;i= <dig> ⋯,n}.
we derive the likelihood for a logistic model of the quasi-linear score because of versatility.
therefore, we assume that y
i is independently distributed according to the bernoulli distribution with a parameter π
i, and consider the logistic model.
below, we denote the quasi-linear score q based on x
i,⋯,x
i as q
i for simplicity.
the association between π
i and the quasi-linear score q
i is described by
<dig> logπi1−πi=qi.
in this setting, the unknown parameters are {α
k,β
k;k= <dig> ⋯,k} which specify q
i’s over individuals as in eq.
.
the log-likelihood function of parameter θ=α <dig> ⋯,αk,β1⊤,⋯,βk⊤⊤ is
<dig> l=∑i=1nyiqi−log).
the maximum likelihood estimator  of θ is therefore the solution of
<dig> ∂l∂θ=w⊤,
where w=⊤, y=⊤ and Π=⊤.
the solution is calculated by updating some initial value repeatedly by fisher’s scoring method as
<dig> θ=θ+w⊤vw−1w⊤y−Π,
where w=∂q1/∂θ,⋯,∂qn/∂θ⊤|θ=θ, v=diagπ11−π <dig> ⋯,πn1−πn and Π=π <dig> ⋯,πn⊤.
the r source code of the parameter estimation of the quasi-linear logistic model is available in additional file  <dig>  in the framework of a generalized linear model, z
=w
θ
+v
−1) is called the working response, and this algorithm is referred to as the iteratively reweighted least-square method  <cit>  because the eq.
is written as θ
=⊤
v
w
)−1
w
⊤
v
z
.
thus, the parameter estimation strategy is very similar to the linear-logistic model.
however, the estimation is not stable in a high dimensional setting.
in such situation, w
⊤
v
w
becomes a singular matrix.
it is thus difficult to compute the inverse matrix in eq.
for each step.
we can avoid the problem by regularization method, just as for the penalized linear logistic model  <cit> .
l <dig> and l <dig> regularization of the quasi-linear logistic model
the l <dig> penalized log-likelihood is described by
<dig> lridge=l−12λ0∑k=1kαk2−12∑k=1kλkβk⊤βk.
we note that we regularized α
k’s by λ
<dig> to avoid computational difficulty in calculating the inverse matrix, although the intercept parameters should not be regularized in the linear logistic model.
a mle with the ridge regularization of θ is calculated by fisher’s scoring method as
<dig> θ=w⊤vw+r−1w⊤v×w⊤θ+v−1y−Π.
here r=diag, where i
m denotes the identity matrix with size m. the derivation of the algorithm is described in additional file 1: appendix d in greater detail.
next we consider l <dig> regularization for the quasi-linear logistic model.
the l <dig> penalized log-likelihood is given by
<dig> llasso=l−∑k=1kλk|βk|.
this form is compatible with the group lasso  <cit> .
we note that the group lasso has a very similar concept in that regularizations are performed for each cluster.
however, the score forms are different between the two regularization methods.
the comparison of group lasso and quasi-linear score are performed in the “application” subsection of the “results”.
for the quasi-linear score q, it is computationally difficult to solve the problem of maximization with  by a method that involves the inverse matrix.
therefore, we applied the gradient ascent method of  <cit>  by using the directional derivative, which is a simple gradient ascent algorithm based on the components of a score function:
<dig> θ=θ+mintoptθ,tedgeθgθ,
where g=,⋯,g
p+k)⊤,
tedge=min1+k≤j≤p+k−θjgj:sign=−sign)≠ <dig>  and
topt=|g|g⊤∂2l∂θ∂θ⊤g.
here g
j=l
j for j= <dig> ⋯,k and
gj=lj−λksignifθj≠0lj−λksign)ifθj=0and|lj|>λk0otherwise
for j=k+ <dig> ⋯,p+k, where sign is a sign function, l
j is the j-th component of eq.
and k denotes the cluster number the j-th marker belongs to.
in each step, the t
opt provides the optimal solution of the gradient descent algorithm and t
edge controls the direction of the gradient so as to avoid changing the signs of the parameters.
the vector of the tuning parameters ⊤ is determined by a cross-validation method from candidate sets of parameters.
non-linearity of the quasi-linear score
the quasi-linear score q is non-linear by definition.
the non-linearity of the quasi-linear score q can be demonstrated by a simple illustration.
figure  <dig> shows the fitted curve of q when p=k= <dig>  in this figure, it looks as if two linear planes, specialized to each sub-space, are connected smoothly.
in this case, the linear surface is curved while still maintaining local linearity, thus forming a quasi-linear surface.
as an extreme case, let there be only one cluster with strong markers.
when all scores are integrated, the information from this cluster should not be affected by the others.
the quasi-linear score q makes up this nature because this approximates the maximum function.
if there is an ℓ,1≤ℓ≤k such that
<dig> lℓ≫lk
fig.
<dig> the boundary  and contour  of the quasi-linear score.
left and center; f=log+exp), the right; g=max.
the center panel is an expansion around the origin  of the left panel
for k≠ℓ, then ∑k=1kexp≈exp, so that q almost equals l
ℓ and the score is almost evaluated by the ℓ-th cluster.
in such a case, the quasi-linear score q achieves the cluster selection.
in the numerical sense, even if the inequality  is not very evident, selection is considered to be achieved because the exponential function inflates the input sufficiently.
for example, log{exp+ exp+ exp+ exp}= <dig> , which essentially means that only the first term is reflected in the construction of the quasi-linear score q. accordingly, q≈l
ℓ if x is in a set {x:l
ℓ=m}, say c
ℓ. we note that c
ℓ is expressed by the intersection of k− <dig> half planes, such that c
ℓ is a convex polyhedron.
thus the quasi-linear score q is locally linear over disjointed and exhaustive regions of the space of all biomarkers : ⋃ℓ=1kcℓ.
thus we observe that the quasi-linear score q is approximately equal to the linear score l that dominates over the other k− <dig> scores.
this property contrasts with the ordinary linear score, which is the sum of k linear scores.
in particular, the quasi-linear score q is advantageous in cases where there are predominant sets of separate biomarkers within the space of all biomarkers.
also, for both logistic models in the parameter estimation steps, we can see the difference between the linear and quasi-linear models reflected in the derivative term as:
<dig> ∂llθl∂θl= <dig> x⊤⊤,
<dig> ∂l∂θ=s <dig> ⋯,sk,s1x⊤,⋯,skx⊤⊤,
where l
l is a log-likelihood function of the linear logistic model with parameter θ
l=⊤ and sk=exp/∑k=1kexp.
a derivation of eq.
is given in additional file 1: appendix e. thus, the data space is decomposed by updated s
k and composed as one unit in each learning step.
this concept used in probabilistic models is referred to as the divide and conquer strategy, which is employed in many machine-learning studies as a mixture of expert models  <cit> .
RESULTS
simulation study
we examined the efficiency of the quasi-linear score q using logistic models , compared with the linear score l using logistic model .
we conducted simulations with five different settings.
for each dataset, the samples were divided between the disease group  and normal group .
first, to show the consistency of the quasi-linear logistic model without regularization, we used a simple setting that has an optimal solution of the quasi-linear form.
in this example, we simulated  <dig> random datasets.
each dataset was either small, containing  <dig> samples, or large, containing  <dig> samples.
next, we estimated the parameters using eq.
and checked the consistency.
second, we examined four high dimensional settings focusing on marker’s selection.
the divided populations were considered to have homogeneous or heterogeneous structure, which were described by normal or mixed normal distribution.
in these examples, we simulated  <dig> random datasets, each containing either  <dig> or  <dig> samples for training and test datasets, respectively.
for these settings, we use the l <dig> and l <dig> shrinkage method in order to avoid overfitting and hard computation.
below, we define rp=∈ℝp for the simple notations.
consistency
in this example, we assumed normality for the normal group and mixture normality for the disease group.
<dig> x|∼n02⊤,i <dig> x|∼∑g=12τgnμ1g,i <dig> ∑g=12τg= <dig>
we let μ
11=⊤ and μ
12=⊤.
in this setting, the bayes optimal form is log+exp).
figure  <dig> shows box plots of estimated parameters for  <dig> trials.
the optimal parameter derived from the true likelihood is =.
the means of the estimated parameters from  <dig> trials are = for the small datasets and = for the large datasets.
we observed that parameter estimation was more precise when the sample size was large, and that the estimated parameters were consistent.
fig.
<dig> box plot of the estimated parameters in the simple simulation.
left: small sample size setting ; right: large sample size setting .
the red lines show the optimal parameters derived from the true likelihood
high dimensional settings
: homo-homo
in this example we assumed normality for both groups.
<dig> x|∼nμy,ip.
we had three settings:  p= <dig> μ0=02⊤,μ1=12⊤,  p= <dig> μ0=0100⊤,μ1= <dig> ⊤,  p= <dig> μ0=0100⊤,μ1= <dig> ⊤.
for the quasi linear score q, we assumed the misspecification of heterogeneous structure, as k= <dig> and p
1=p
2= <dig> for  or p
1=p
2= <dig> for  and .
: homo-hetero
in this example, we assumed normality for the normal group and mixed normality for the disease group.
<dig> x|∼n,x|∼∑g=1gτgn,∑g=1gτg= <dig>
we had four settings.
in  and , we let g= <dig>  p= <dig>  τ
1=τ
2= <dig> , μ0=0100⊤.
in  and , we let g= <dig>  p= <dig>  τ
1=τ
2=τ
3=1/ <dig>  μ0=0100⊤.
the mean parameter for the disease group was set as  μ
11=⊤, μ
12=⊤,  μ
11=⊤, μ
12=⊤,  μ
11=⊤,  μ
11=⊤, μ
12=⊤, μ
13=⊤.
for the quasi-linear score q we assumed the correct specification of heterogeneous structure as k=g and p
1=p
2= <dig> or p
1= <dig> p
2=p
3= <dig>
: hetero-hetero
in this example, we assumed mixed normality for both groups.
<dig> x|∼∑g=1gτygn,∑g=1gτyg= <dig>
we used the following settings: g= <dig>  p= <dig>  τ
yg= <dig>  , μ01=0100⊤, μ
02=⊤, μ
11=⊤, μ
12=.
for the quasi-linear score q we assumed to specify there are heterogeneous structure as k= <dig> and p
1=p
2= <dig>
: correlated
in this example, we assumed normality for the normal group and mixed normality for the disease group.
<dig> x|∼n,x|∼∑g=1gτgn,∑g=1gτg= <dig>
the variance assumption was based on a real dataset, as shown in fig.
<dig>  we used the following settings:  g= <dig> p= <dig> τ1=τ2= <dig> ,μ0=070⊤,μ11=⊤, μ
12=⊤, Σ=Σ1Σ2Σ2⊤Σ <dig>  where Σ
1= <dig> i
35+ <dig> j
<dig> Σ
2=− <dig> j
<dig>  where j
m is a matrix of size m of which all components are  <dig>  for the quasi-linear score q we assumed to specify there are heterogeneous structure as k= <dig> and p
1=p
2= <dig>
table  <dig> estimated auc  of  <dig> repetitions
application
we applied our method for two datasets, namely breast cancer and prostate cancer data.
for both types of datasets, two independent datasets were used as training and testing to evaluate the predictive ability by test auc.
first, we compared the test auc among decision tree , random forest , support vector machine , naive bayes , group lasso , neural network , l <dig> or l <dig> penalized linear logistic  and l <dig> or l <dig> penalized quasi-linear logistic .
performance was evaluated by the test auc and the 95% cis of the test auc based on  <dig> bootstrapping sampling, as described in  <cit> .
the tuning parameters were determined with a grid search and resampling method as needed.
second, the stability for marker selection was compared among ll <dig>  ql <dig> and gl.
we used a similarity index proposed by  <cit>  defined by s=|a∩b|/|a∪b|, where a and b are subsets of marker index set, and |a| is a cardinality of the set a. s takes a value between  <dig> and  <dig> whose high value means high stability.
we evaluated the stability measure by 2r∑i=1r−1∑j=i+1rs, where m
<dig> ⋯,m
r are sets of the selected marker for r bootstrap sample sets from the training data set.
r was set to  <dig> below.
breast cancer data
the training dataset was taken from  <cit>  and the test datasets were from  <cit> .
yan et al.
<cit>  used these datasets and compared the aucs by the linear score l, which they evaluated by traditional methods as well as methods they proposed.
we focused on the  <dig> genes detected by  <cit> , as in  <cit> .
these datasets include  <dig> patients in one and  <dig> patients in the other.
for ql, grouping of  <dig> genes was based on the ward’s clustering method only by training dataset.
we had two options for dividing all the genes into clusters.
in the first option, the  <dig> genes were divided into two clusters, one with  <dig> and the other with  <dig> genes.
in the second option, the  <dig> genes were divided into three clusters of  <dig>   <dig>  and  <dig> genes.
for gl, we used two clusters option.
fig.
<dig> box plots of the test auc for all comparative methods by breast cancer data
when we use the linear score l, the absolute value of the coefficients of each marker reflects the order of importance of all markers for prediction.
therefore, the linear score is understandable in the sense that we can recognize strong markers.
this is no longer a consideration when we use a generalized non-linear score.
however, the quasi-linear score enables us to compare coefficients within the same cluster.
an example is shown in fig.
<dig>  which displays the ranking of the absolute values of the estimated coefficients by the ridge regularization method based on the existence of two clusters.
the gene labels are arranged in order of the rankings.
we observed that q and l gave quite different rankings.
this result shows that the quasi-linear score would produce different interpretations for the relationship between the markers.
fig.
<dig> ranking of the absolute values of the coefficients within the cluster with ridge regularization
fig.
<dig> learning and fitted plot for the training and test dataset when using the quasi-linear score of two clusters with lasso regularization.
the horizontal and vertical axis are the linear scores of the first and second cluster.
red points indicate the metastatic group and black points indicate the control group.
curve lines are contours of the quasi-linear score and blue line shows cut-off value based on youden-index
although the quasi-linear score q is approximately equivalent to the maximum function, the two are numerically different.
in fact, the test auc of the quasi-linear score with the lasso regularization method when we assumed two clusters was  <dig> , and the corresponding maximum score m is  <dig> , so that the smooth non-linearity of the quasi-linear form produced good predictive performance
the elastic net shrinkage method  <cit> , which combines the lasso and ridge shrinkage methods, is among the most frequently used.
when we combined the quasi-linear score and the elastic net regularization, the number of the tuning parameters was inflated.
although we used the elastic net experimentally for the application for some selected parameters, the predictive performance was not significantly different from the performance obtained with either the lasso or ridge.
detailed results are summarized in table  <dig>  moreover, to check the utility of the unsupervised clustering, we randomly divided the  <dig> genes into two subsets of  <dig> and  <dig> genes, and applied ql <dig> for the test dataset .
figure  <dig> shows that clustered subsets  performs better than randomly divided subsets.
thus, unsupervised clustering naturally benefits supervised learning via the quasi-linear form.
fig.
<dig> test aucs by the quasi-linear score for the dataset from buyse et al.
.
the score is learning by randomly divided genes subsets for the dataset from van’t veer et al .
the red line is the test auc by the quasi-linear score, which consists of subsets of genes clustered by unsupervised learning
a parameter ε denotes the proportion of ridge regularization to lasso regularization
prostate cancer data
the data set was taken from  <cit>  which contains expression data for  <dig> genes obtained from  <dig> prostate cancer tumors.
the tumors were from  <dig> subjects determined to be fusion status-positive and  <dig> subjects determined to be fusion status-negative.
we randomly divided the whole dataset into two independent datasets with the same number of tumor samples  while maintaining the ratio of positive to negative statuses.
first, we selected  <dig> relevant genes which had top  <dig> absolute value of t-statistic between the two statuses using only the training dataset.
such marker preselection has been performed in many studies  <cit> .
for ql, grouping of  <dig> genes was based on the ward’s clustering method only by training dataset.
we had two options for dividing all the genes into clusters.
in the first option, the  <dig> genes were divided into two clusters, one with  <dig> and the other with  <dig> genes.
in the second option, the  <dig> genes were divided into three clusters of  <dig>   <dig>  and  <dig> genes.
for gl, we used two clusters option.
we then compared the test auc among all comparative methods.
figure  <dig> displays the estimated auc for the test dataset.
as well as the application for breast cancer data, ql <dig> and ql <dig> performed better than any other comparative methods.
the numbers of selected markers in ll <dig>  ql <dig> , ql <dig>  and gl were  <dig>   <dig>   <dig> and  <dig>  respectively.
similarly, the stability measures were  <dig> ,  <dig> ,  <dig>  and  <dig> , respectively.
the stability of ql <dig> was higher than ll <dig>  we note that gl almost did not shrink any coefficients to zero as application for breast cancer data set.
fig.
<dig> box plots of the test auc for all comparative methods by prostate cancer data
discussion
we focused on heterogeneous structure and determined how to reflect such heterogeneity in the score function defined in .
for this purpose, the quasi-linear score was derived as the generalized mean called the kolmogorov-nagumo average.
the quasi-linear form is also called a soft maximum function or log-sum-exp function  <cit> .
in machine learning, the softmax function is often used as a differentiable approximation of the maximum.
in computer science, the log-sum-exp function is used to avoid computational problems such as overflow.
the non-linearity of the quasi-linear score is explained by the soft maximum function.
the quasi-linear score achieves cluster selection because of the soft maximum property as discussed in the subsection of “non-linearity of the quasi-linear score”.
this formulation does not require any prior information or assumption to separate markers into clusters, because this is achieved by the unsupervised learning step.
the quasi-linear score is based on the idea of combining predictors, which is related to several ideas in the literature.
for example, a mixture of expert models suggest the idea of decomposing input space  <cit> , in which the model divides the problem space probabilistically and the scores learned in all sub-spaces are combined.
the quasi-linear score utilizes the information given by the clustering method to reflect the heterogeneity of markers and combines the linear scores of all clusters.
hence, it relies on the disjointed decomposition of the markers.
the method of combining linear scores was also discussed in  <cit> , known as composite links, which assumes that the score is formed by a weighted sum of block-wise markers.
unlike the generalized linear model, the composite link model does not restrict the use of the single link function.
in a special case, the composite link logistic model corresponds to the quasi-linear logistic model.
however, these ideas differ in that the composite link considers the sum of the linked linear scores whereas the quasi-linear score considers a linkage of the summarization of linear scores in all clusters.
the key in our proposal is to model heterogeneity using information from the clustering method, thereby connecting supervised learning with unsupervised learning without any assumption via a change in the score form from the simple average to the kolmogorov-nagumo average.
for future work, we intend to extend some fixed settings presented in this report.
these include the choice of the clustering method, the size of the markers and clusters, the set of tuning parameters, the type of outcomes, and the format of the targeted data.
because the quasi-linear score can be defined by any decomposition ideas, the performance should be evaluated by clustering methods other than ward’s method, such as the k-means method  <cit> .
moreover, we need to investigate the sizes of markers and clusters, and the number of candidate sets of tuning parameters in addition to the parameter τ in , to obtain a more flexible form of the quasi-linear function.
although we applied and evaluated the proposed method after marker preselection in application, the performance should be evaluated in much higher dimensional setting.
an especially big concern is how to decide the cluster size for the quasi-linear score.
as described in the “application” subsection, the quasi-linear score by cluster size  <dig> gave the best predictive performance for breast cancer data, and adding more clusters yielded no improvement.
figure  <dig> supports this result because whole markers were divided into two primary clusters.
it is necessary to develop an objective index of definite cluster size selection for general applications.
the quasi-linear score would be also applicable in a case of the continuous outcomes and in a regression model, although we focused on binary outcomes and the logistic model in this study.
the performance of the quasi-linear score would be exhibited in the mixed large dataset, which would play an important role in biomedical studies in the near future, because such data must be heterogeneous.
furthermore, our method is not limited to biomedical data, and could also be beneficial for analyzing any data that have heterogeneous structure.
CONCLUSIONS
in this paper, we focused on heterogeneous structure.
such heterogeneity was captured well by a clustering method.
the quasi-linear score was naturally derived by bayes risk consistency between mixed and standard normal distributions.
moreover, the quasi-linear score approximates the maximum function and plays an important role in selecting the most effective cluster for prediction from given clusters.
the quasi-linear score has better predictive ability compared to linear score as shown in simulation studies and applications to real data.
additional files
additional file  <dig> technical derivations.
in this file, we perform some technical derivations and evaluations for quasi-linear score: parameterization; the relationship with linear and maximum score; the bayes risk consistency; l <dig> and l <dig> regularization methods; the derivatives.
additional file  <dig> r source code of the parameter estimation of the quasi-linear logistic model.
in this file, we introduce the r source code of the parameter estimation of the quasi-linear logistic model, which was used for simulation and application.
abbreviations
aucarea under curve
dldadiagonal fisher’s linear discriminant analysis
mlemaximum likelihood estimator
electronic supplementary material
the online version of this article  contains supplementary material, which is available to authorized users.

 in recent years, biomedical data have become complicated and high-dimensional. for example, a single human gene expression dataset contains tens of thousands of features, many of which are highly correlated. in addition, large mixed datasets are crucial for personalized treatment, in which the optimal treatment strategy is determined based on a dataset that combines a very large number of prognostic factors. from the viewpoint of statistical machine learning, supervised and unsupervised learning methods play central roles in such biomedical studies. in fact, shrinkage methods such as ridge and lasso are frequently used in the context of prediction, and clustering methods are used in the context of interpretation, potentially revealing novel findings. supervised learning methods are often used to estimate risk scores when predicting dichotomous outcomes. linear scores are among the most widely used forms because it is easy to learn the predictive score from a training dataset. moreover, it is easy to understand the estimated score. linear scores are often evaluated by linear discriminant or logistic regression analysis, and achieve not bad discriminative performance. for example, the study of used a linear score, and their discoveries led to the development of mammaprint, a diagnostic kit for breast cancer metastasis. unsupervised learning methods can also yield beneficial insights in high-dimensional data analysis. for example, used biclustering to reveal more detailed subtypes in breast carcinomas with distinctive gene expression profiles from a group that was previously regarded as a homogenous unit. their study revealed that in order to understand biodiversity, the heterogeneous structure of the targeted population must be considered, and that such heterogeneity can be clarified by the clustering method. previously published reviews have described both clustering methods and biclustering algorithms. several studies have combined supervised and unsupervised learning methods. for example, used clustering to discover different patterns of gene expression in different subgroups. they then derived the respective scores for these groups and achieved good specificity without loss of sensitivity relative to existing diagnostic rules. sample heterogeneity may result in marker heterogeneity. as a result, different samples in different subgroups may have different intrinsic characteristics in their environmental and genetic factors as demonstrated by the motivated example in the “methods” section. such heterogeneity may have unexpected effects on a therapy or treatment which is considered as best practice, and lead to an unfavorable risk in one part of the population. bravo et al. focused on the marker heterogeneity by detecting the genes that showed different variation between healthy and disease samples. they then defined an anti-profile score as the number of hyper-variable genes. thus, more and more studies have considered heterogeneous structure and reflected this heterogeneity in their predictions. however, the risk scores highlighted by published papers are linear, and heterogeneity is therefore not directly reflected in the score form. in this study, we focused on heterogeneity and determined how to directly reflect this intrinsic characteristic in the score form. we developed the quasi-linear score as a result, which combines linear scores as a kolmogorov-nagumo average, enabling us to reflect the clustering result naturally, because it is based on separated feature vectors. the rest of this paper is organized as follows. in the “methods” section, we first present a motivated example of gene expression data and develop the quasi-linear score. heterogeneity is observed via the clustering method, and we define the quasi-linear score to reflect gene clusters with a generalized average form. we also formulate the quasi-linear logistic model and discuss the difference between the linear and quasi-linear scores. we subsequently evaluated our method by numerical simulations and applications to real datasets. we refer to the relationship between the quasi-linear score and traditional combined approaches in “discussion” section. all technical details given as appendix are available in additional file methods motivation and derivation we studied the gene expression dataset from. this dataset is derived from non-metastatic and metastatic breast cancer patients. in their study, the linear score was evaluated to discriminate metastatic events. because estimation of the predictive linear score is often achieved by a diagonal fisher’s linear discriminant analysis, we considered applying dlda to this dataset. because the coefficients of the linear score estimated by dlda correspond to the t-statistic values, we checked the t-statistics directly for the purpose of visualization. if the data have heterogeneous structure, it can be clarified by observing the difference between two divided, independent datasets. therefore, we divided the full data into two independent sets, data and data before calculating the t-statistics for each of them separately. figure shows the correspondence of the t-statistics. some genes had no consistency in the signs of their t-values, indicating that some samples from the metastatic group had higher expression, whereas other samples had lower expression, relative to the non-metastatic group. this phenomenon may be caused by heterogeneous factors. in fact, due to the existence of multiple subtypes of breast cancer, this disease is known to exhibit heterogeneity. for such heterogeneous data, clustering methods should work well, as shown in. we applied clustering according to the ward’s method, as shown in fig. which highlights the results of clustering and the correlation matrix arranged by the clustered genes. although biclustering result was not suggestive of the heterogeneity in appearance, it was observed via the correlation matrix. some genes are strongly correlated with others in the same cluster. thus, we observed the existence of heterogeneity using a t-statistics plot and trends in the expression patterns by clustering. next, we developed an appropriate score form for discriminating such heterogeneous data based on clustering. fig. t-statistic values for two datasets from van’t veer et al.. the red points show the genes with sign mismatched t-values for these data fig. the hierarchical clustering and the correlation matrix of genes for the dataset from van’t veer et al.. the figure shows the clustering result and the correlation matrix. there are rows representing genes and columns representing samples and the gene expression data ranging from green to red are displayed. elements of the correlation matrix ranging from blue to yellow are displayed we assume to know decomposition of p biomarkers into k groups by clustering. based on these k sets of clustered markers, we define a quasi-linear score as q=log∑k=1kexp, where lk=αk+βk⊤x with the parameters α k,β k, and the marker vector x for the k-th cluster of k=,⋯,k. when k= the quasi-linear score q is reduced to the linear score, l=α+β⊤x, where α and β=β1⊤,⋯,βk⊤⊤ are parameters, and x is the full vector of x,⋯,x. we note that the intercepts α k’s in are reduced to the single intercept α in. additional file 1: appendix a gives another parameterization for q in which the intercepts α k are uniquely decomposed to the overall intercept and weights of the k clusters. when determining how to reflect cluster information in the score form, we had two main considerations: which scores should be integrated, and how integration should be performed. all the linear scores l k are integrated in the quasi-linear score q. we believe that this is reasonable because there are similar markers in each cluster, and we expect that heterogeneity would be caused by different mixed homogeneous features that are sufficiently described by the linear form. although such an idea of combining the linear scores has already been proposed by as a composite link, it is different from the quasi-linear score in several ways. one of the most significant differences between them is that the quasi-linear score is defined by disjointed sets of markers. this results in a small number of parameters for the predictive score: the parsimonious expression. the difference in these forms is mentioned in the discussion. moreover, the quasi-linear score q summarizes l k approximating the maximum function. in fact, is equal to a soft maximum function discussed by, which can be approximated with m=max1≤k≤klk. therefore, the quasi-linear score q respects the maximum of k linear scores from all clusters. see for a discussion of eq. as maxout in neural networks. the relationships among q, l, and m are clearly evaluated when a tuning parameter, τ, is introduced in the quasi-linear score q as qτ=1τlog∑k=1kexp, where τ is a positive parameter. if l k is fixed for all k, then the form of is defined solely by the tuning parameter τ. when τ is equal to q τ is equal to the quasi-linear score q by definition. when τ goes to infinity, q τ is simply the maximum score m. when τ goes to q τ is equivalent to the linear score l. thus, these are unified by the hardness of approximation to the maximum function. more details are provided in additional file 1: appendix b. the characteristics of the quasi-linear score q are understood by a more general expression of: g=ϕ−1), where ϕ is an invertible function. we define g as a generalized quasi-linear score because the form is the generalized mean called the kolmogorov-nagumo average. if we take a simple average, ϕ=z, then the generalized quasi-linear score g corresponds to the linear score l. in this sense, the linear score l is a simple mean, and the quasi-linear score q is a generalized mean of linear scores l k averaged by the exponential function. although the simplest integration of clustered information is achieved by a simple average, resulting in the linear score form, it is intuitively unsatisfying because the predictive performance of these linear scores l k differs among the clusters. a cluster that strongly discriminates the outcome on its own should be respected in comparison with the other clusters. if only the cluster with the highest linear score is reflected in the prediction, it is described by the maximum score m. however, this situation is still not ideal, because only one cluster is reflected in the prediction, and form is difficult to handle mathematically because it is not differentiable when two or more linear scores are equal. consequently, parameter estimation becomes impossible. the quasi-linear score q is therefore naturally derived, and it is reasonable for discriminant analysis of the heterogeneous data because the quasi-linear score q plays an important role in cluster selection, as discussed later. in the following sections, we let ϕ=exp, both because this form is approximated by the maximum function, and because it is optimal in the sense of bayes risk consistency when we consider the simple case in which the label conditional random variables follow a mixture of normal distribution and a normal distribution with equal variance, respectively. additional file 1: appendix c provides more detail about the bayes risk consistency of the situation. moreover, the exponential function gives us an understandable interpretation of the parameter estimation. because we modify only the scoring form, the quasi-linear score q can be applied to all traditional settings in biostatistics, as the generalized linear model with the l and l shrinkage methods. in particular, when we combine the quasi-linear score q with lasso shrinkage, the important clusters and variables in each cluster are determined simultaneously because of soft maximum property and l sparseness. this property provides good performance for the discriminant problem when the data have a much larger number of correlated markers than the number of samples. therefore, we derive the l and l shrinkage quasi-linear logistic model and display the performance of the quasi-linear score q when it is applied to gene expression data. likelihood for logistic model and maximum likelihood estimation consider the data {;i= ⋯,n}, where x i is a covariate vector and y i is a dichotomous outcome which takes or with the i-th individual. assume that we know the decomposition of x i as x i,⋯,x i with a fixed cluster size k, and that this is identical among individuals. we denote the size of x i as p k, where ∑k=1kpk=p. we note that the decomposition is given a priori by one of clustering methods for {x i;i= ⋯,n}. we derive the likelihood for a logistic model of the quasi-linear score because of versatility. therefore, we assume that y i is independently distributed according to the bernoulli distribution with a parameter π i, and consider the logistic model. below, we denote the quasi-linear score q based on x i,⋯,x i as q i for simplicity. the association between π i and the quasi-linear score q i is described by logπi1−πi=qi. in this setting, the unknown parameters are {α k,β k;k= ⋯,k} which specify q i’s over individuals as in eq.. the log-likelihood function of parameter θ=α ⋯,αk,β1⊤,⋯,βk⊤⊤ is l=∑i=1nyiqi−log). the maximum likelihood estimator of θ is therefore the solution of ∂l∂θ=w⊤, where w=⊤, y=⊤ and Π=⊤. the solution is calculated by updating some initial value repeatedly by fisher’s scoring method as θ=θ+w⊤vw−1w⊤y−Π, where w=∂q1/∂θ,⋯,∂qn/∂θ⊤|θ=θ, v=diagπ11−π ⋯,πn1−πn and Π=π ⋯,πn⊤. the r source code of the parameter estimation of the quasi-linear logistic model is available in additional file in the framework of a generalized linear model, z =w θ +v −1) is called the working response, and this algorithm is referred to as the iteratively reweighted least-square method because the eq. is written as θ =⊤ v w )−1 w ⊤ v z. thus, the parameter estimation strategy is very similar to the linear-logistic model. however, the estimation is not stable in a high dimensional setting. in such situation, w ⊤ v w becomes a singular matrix. it is thus difficult to compute the inverse matrix in eq. for each step. we can avoid the problem by regularization method, just as for the penalized linear logistic model. l and l regularization of the quasi-linear logistic model the l penalized log-likelihood is described by lridge=l−12λ0∑k=1kαk2−12∑k=1kλkβk⊤βk. we note that we regularized α k’s by λ to avoid computational difficulty in calculating the inverse matrix, although the intercept parameters should not be regularized in the linear logistic model. a mle with the ridge regularization of θ is calculated by fisher’s scoring method as θ=w⊤vw+r−1w⊤v×w⊤θ+v−1y−Π. here r=diag, where i m denotes the identity matrix with size m. the derivation of the algorithm is described in additional file 1: appendix d in greater detail. next we consider l regularization for the quasi-linear logistic model. the l penalized log-likelihood is given by llasso=l−∑k=1kλk|βk|. this form is compatible with the group lasso. we note that the group lasso has a very similar concept in that regularizations are performed for each cluster. however, the score forms are different between the two regularization methods. the comparison of group lasso and quasi-linear score are performed in the “application” subsection of the “results”. for the quasi-linear score q, it is computationally difficult to solve the problem of maximization with by a method that involves the inverse matrix. therefore, we applied the gradient ascent method of by using the directional derivative, which is a simple gradient ascent algorithm based on the components of a score function: θ=θ+mintoptθ,tedgeθgθ, where g=,⋯,g p+k)⊤, tedge=min1+k≤j≤p+k−θjgj:sign=−sign)≠ and topt=|g|g⊤∂2l∂θ∂θ⊤g. here g j=l j for j= ⋯,k and gj=lj−λksignifθj≠0lj−λksign)ifθj=0and|lj|>λk0otherwise for j=k+ ⋯,p+k, where sign is a sign function, l j is the j-th component of eq. and k denotes the cluster number the j-th marker belongs to. in each step, the t opt provides the optimal solution of the gradient descent algorithm and t edge controls the direction of the gradient so as to avoid changing the signs of the parameters. the vector of the tuning parameters ⊤ is determined by a cross-validation method from candidate sets of parameters. non-linearity of the quasi-linear score the quasi-linear score q is non-linear by definition. the non-linearity of the quasi-linear score q can be demonstrated by a simple illustration. figure shows the fitted curve of q when p=k= in this figure, it looks as if two linear planes, specialized to each sub-space, are connected smoothly. in this case, the linear surface is curved while still maintaining local linearity, thus forming a quasi-linear surface. as an extreme case, let there be only one cluster with strong markers. when all scores are integrated, the information from this cluster should not be affected by the others. the quasi-linear score q makes up this nature because this approximates the maximum function. if there is an ℓ,1≤ℓ≤k such that lℓ≫lk fig. the boundary and contour of the quasi-linear score. left and center; f=log+exp), the right; g=max. the center panel is an expansion around the origin of the left panel for k≠ℓ, then ∑k=1kexp≈exp, so that q almost equals l ℓ and the score is almost evaluated by the ℓ-th cluster. in such a case, the quasi-linear score q achieves the cluster selection. in the numerical sense, even if the inequality is not very evident, selection is considered to be achieved because the exponential function inflates the input sufficiently. for example, log{exp+ exp+ exp+ exp}=, which essentially means that only the first term is reflected in the construction of the quasi-linear score q. accordingly, q≈l ℓ if x is in a set {x:l ℓ=m}, say c ℓ. we note that c ℓ is expressed by the intersection of k− half planes, such that c ℓ is a convex polyhedron. thus the quasi-linear score q is locally linear over disjointed and exhaustive regions of the space of all biomarkers: ⋃ℓ=1kcℓ. thus we observe that the quasi-linear score q is approximately equal to the linear score l that dominates over the other k− scores. this property contrasts with the ordinary linear score, which is the sum of k linear scores. in particular, the quasi-linear score q is advantageous in cases where there are predominant sets of separate biomarkers within the space of all biomarkers. also, for both logistic models in the parameter estimation steps, we can see the difference between the linear and quasi-linear models reflected in the derivative term as: ∂llθl∂θl= x⊤⊤, ∂l∂θ=s ⋯,sk,s1x⊤,⋯,skx⊤⊤, where l l is a log-likelihood function of the linear logistic model with parameter θ l=⊤ and sk=exp/∑k=1kexp. a derivation of eq. is given in additional file 1: appendix e. thus, the data space is decomposed by updated s k and composed as one unit in each learning step. this concept used in probabilistic models is referred to as the divide and conquer strategy, which is employed in many machine-learning studies as a mixture of expert models.  simulation study we examined the efficiency of the quasi-linear score q using logistic models, compared with the linear score l using logistic model. we conducted simulations with five different settings. for each dataset, the samples were divided between the disease group and normal group. first, to show the consistency of the quasi-linear logistic model without regularization, we used a simple setting that has an optimal solution of the quasi-linear form. in this example, we simulated random datasets. each dataset was either small, containing samples, or large, containing samples. next, we estimated the parameters using eq. and checked the consistency. second, we examined four high dimensional settings focusing on marker’s selection. the divided populations were considered to have homogeneous or heterogeneous structure, which were described by normal or mixed normal distribution. in these examples, we simulated random datasets, each containing either or samples for training and test datasets, respectively. for these settings, we use the l and l shrinkage method in order to avoid overfitting and hard computation. below, we define rp=∈ℝp for the simple notations. consistency in this example, we assumed normality for the normal group and mixture normality for the disease group. x|∼n02⊤,i x|∼∑g=12τgnμ1g,i ∑g=12τg= we let μ 11=⊤ and μ 12=⊤. in this setting, the bayes optimal form is log+exp). figure shows box plots of estimated parameters for trials. the optimal parameter derived from the true likelihood is =. the means of the estimated parameters from trials are = for the small datasets and = for the large datasets. we observed that parameter estimation was more precise when the sample size was large, and that the estimated parameters were consistent. fig. box plot of the estimated parameters in the simple simulation. left: small sample size setting; right: large sample size setting. the red lines show the optimal parameters derived from the true likelihood high dimensional settings: homo-homo in this example we assumed normality for both groups. x|∼nμy,ip. we had three settings: p= μ0=02⊤,μ1=12⊤, p= μ0=0100⊤,μ1= ⊤, p= μ0=0100⊤,μ1= ⊤. for the quasi linear score q, we assumed the misspecification of heterogeneous structure, as k= and p 1=p 2= for or p 1=p 2= for and.: homo-hetero in this example, we assumed normality for the normal group and mixed normality for the disease group. x|∼n,x|∼∑g=1gτgn,∑g=1gτg= we had four settings. in and, we let g= p= τ 1=τ 2=, μ0=0100⊤. in and, we let g= p= τ 1=τ 2=τ 3=1/ μ0=0100⊤. the mean parameter for the disease group was set as μ 11=⊤, μ 12=⊤, μ 11=⊤, μ 12=⊤, μ 11=⊤, μ 11=⊤, μ 12=⊤, μ 13=⊤. for the quasi-linear score q we assumed the correct specification of heterogeneous structure as k=g and p 1=p 2= or p 1= p 2=p 3=: hetero-hetero in this example, we assumed mixed normality for both groups. x|∼∑g=1gτygn,∑g=1gτyg= we used the following settings: g= p= τ yg=, μ01=0100⊤, μ 02=⊤, μ 11=⊤, μ 12=. for the quasi-linear score q we assumed to specify there are heterogeneous structure as k= and p 1=p 2=: correlated in this example, we assumed normality for the normal group and mixed normality for the disease group. x|∼n,x|∼∑g=1gτgn,∑g=1gτg= the variance assumption was based on a real dataset, as shown in fig. we used the following settings: g= p= τ1=τ2=,μ0=070⊤,μ11=⊤, μ 12=⊤, Σ=Σ1Σ2Σ2⊤Σ where Σ 1= i 35+ j Σ 2=− j where j m is a matrix of size m of which all components are for the quasi-linear score q we assumed to specify there are heterogeneous structure as k= and p 1=p 2= table estimated auc of repetitions application we applied our method for two datasets, namely breast cancer and prostate cancer data. for both types of datasets, two independent datasets were used as training and testing to evaluate the predictive ability by test auc. first, we compared the test auc among decision tree, random forest, support vector machine, naive bayes, group lasso, neural network, l or l penalized linear logistic and l or l penalized quasi-linear logistic. performance was evaluated by the test auc and the 95% cis of the test auc based on bootstrapping sampling, as described in. the tuning parameters were determined with a grid search and resampling method as needed. second, the stability for marker selection was compared among ll ql and gl. we used a similarity index proposed by defined by s=|a∩b|/|a∪b|, where a and b are subsets of marker index set, and |a| is a cardinality of the set a. s takes a value between and whose high value means high stability. we evaluated the stability measure by 2r∑i=1r−1∑j=i+1rs, where m ⋯,m r are sets of the selected marker for r bootstrap sample sets from the training data set. r was set to below. breast cancer data the training dataset was taken from and the test datasets were from. yan et al. used these datasets and compared the aucs by the linear score l, which they evaluated by traditional methods as well as methods they proposed. we focused on the genes detected by, as in. these datasets include patients in one and patients in the other. for ql, grouping of genes was based on the ward’s clustering method only by training dataset. we had two options for dividing all the genes into clusters. in the first option, the genes were divided into two clusters, one with and the other with genes. in the second option, the genes were divided into three clusters of and genes. for gl, we used two clusters option. fig. box plots of the test auc for all comparative methods by breast cancer data when we use the linear score l, the absolute value of the coefficients of each marker reflects the order of importance of all markers for prediction. therefore, the linear score is understandable in the sense that we can recognize strong markers. this is no longer a consideration when we use a generalized non-linear score. however, the quasi-linear score enables us to compare coefficients within the same cluster. an example is shown in fig. which displays the ranking of the absolute values of the estimated coefficients by the ridge regularization method based on the existence of two clusters. the gene labels are arranged in order of the rankings. we observed that q and l gave quite different rankings. this result shows that the quasi-linear score would produce different interpretations for the relationship between the markers. fig. ranking of the absolute values of the coefficients within the cluster with ridge regularization fig. learning and fitted plot for the training and test dataset when using the quasi-linear score of two clusters with lasso regularization. the horizontal and vertical axis are the linear scores of the first and second cluster. red points indicate the metastatic group and black points indicate the control group. curve lines are contours of the quasi-linear score and blue line shows cut-off value based on youden-index although the quasi-linear score q is approximately equivalent to the maximum function, the two are numerically different. in fact, the test auc of the quasi-linear score with the lasso regularization method when we assumed two clusters was, and the corresponding maximum score m is, so that the smooth non-linearity of the quasi-linear form produced good predictive performance the elastic net shrinkage method, which combines the lasso and ridge shrinkage methods, is among the most frequently used. when we combined the quasi-linear score and the elastic net regularization, the number of the tuning parameters was inflated. although we used the elastic net experimentally for the application for some selected parameters, the predictive performance was not significantly different from the performance obtained with either the lasso or ridge. detailed results are summarized in table moreover, to check the utility of the unsupervised clustering, we randomly divided the genes into two subsets of and genes, and applied ql for the test dataset. figure shows that clustered subsets performs better than randomly divided subsets. thus, unsupervised clustering naturally benefits supervised learning via the quasi-linear form. fig. test aucs by the quasi-linear score for the dataset from buyse et al.. the score is learning by randomly divided genes subsets for the dataset from van’t veer et al. the red line is the test auc by the quasi-linear score, which consists of subsets of genes clustered by unsupervised learning a parameter ε denotes the proportion of ridge regularization to lasso regularization prostate cancer data the data set was taken from which contains expression data for genes obtained from prostate cancer tumors. the tumors were from subjects determined to be fusion status-positive and subjects determined to be fusion status-negative. we randomly divided the whole dataset into two independent datasets with the same number of tumor samples while maintaining the ratio of positive to negative statuses. first, we selected relevant genes which had top absolute value of t-statistic between the two statuses using only the training dataset. such marker preselection has been performed in many studies. for ql, grouping of genes was based on the ward’s clustering method only by training dataset. we had two options for dividing all the genes into clusters. in the first option, the genes were divided into two clusters, one with and the other with genes. in the second option, the genes were divided into three clusters of and genes. for gl, we used two clusters option. we then compared the test auc among all comparative methods. figure displays the estimated auc for the test dataset. as well as the application for breast cancer data, ql and ql performed better than any other comparative methods. the numbers of selected markers in ll ql, ql and gl were and respectively. similarly, the stability measures were,, and, respectively. the stability of ql was higher than ll we note that gl almost did not shrink any coefficients to zero as application for breast cancer data set. fig. box plots of the test auc for all comparative methods by prostate cancer data discussion we focused on heterogeneous structure and determined how to reflect such heterogeneity in the score function defined in. for this purpose, the quasi-linear score was derived as the generalized mean called the kolmogorov-nagumo average. the quasi-linear form is also called a soft maximum function or log-sum-exp function. in machine learning, the softmax function is often used as a differentiable approximation of the maximum. in computer science, the log-sum-exp function is used to avoid computational problems such as overflow. the non-linearity of the quasi-linear score is explained by the soft maximum function. the quasi-linear score achieves cluster selection because of the soft maximum property as discussed in the subsection of “non-linearity of the quasi-linear score”. this formulation does not require any prior information or assumption to separate markers into clusters, because this is achieved by the unsupervised learning step. the quasi-linear score is based on the idea of combining predictors, which is related to several ideas in the literature. for example, a mixture of expert models suggest the idea of decomposing input space, in which the model divides the problem space probabilistically and the scores learned in all sub-spaces are combined. the quasi-linear score utilizes the information given by the clustering method to reflect the heterogeneity of markers and combines the linear scores of all clusters. hence, it relies on the disjointed decomposition of the markers. the method of combining linear scores was also discussed in, known as composite links, which assumes that the score is formed by a weighted sum of block-wise markers. unlike the generalized linear model, the composite link model does not restrict the use of the single link function. in a special case, the composite link logistic model corresponds to the quasi-linear logistic model. however, these ideas differ in that the composite link considers the sum of the linked linear scores whereas the quasi-linear score considers a linkage of the summarization of linear scores in all clusters. the key in our proposal is to model heterogeneity using information from the clustering method, thereby connecting supervised learning with unsupervised learning without any assumption via a change in the score form from the simple average to the kolmogorov-nagumo average. for future work, we intend to extend some fixed settings presented in this report. these include the choice of the clustering method, the size of the markers and clusters, the set of tuning parameters, the type of outcomes, and the format of the targeted data. because the quasi-linear score can be defined by any decomposition ideas, the performance should be evaluated by clustering methods other than ward’s method, such as the k-means method. moreover, we need to investigate the sizes of markers and clusters, and the number of candidate sets of tuning parameters in addition to the parameter τ in, to obtain a more flexible form of the quasi-linear function. although we applied and evaluated the proposed method after marker preselection in application, the performance should be evaluated in much higher dimensional setting. an especially big concern is how to decide the cluster size for the quasi-linear score. as described in the “application” subsection, the quasi-linear score by cluster size gave the best predictive performance for breast cancer data, and adding more clusters yielded no improvement. figure supports this result because whole markers were divided into two primary clusters. it is necessary to develop an objective index of definite cluster size selection for general applications. the quasi-linear score would be also applicable in a case of the continuous outcomes and in a regression model, although we focused on binary outcomes and the logistic model in this study. the performance of the quasi-linear score would be exhibited in the mixed large dataset, which would play an important role in biomedical studies in the near future, because such data must be heterogeneous. furthermore, our method is not limited to biomedical data, and could also be beneficial for analyzing any data that have heterogeneous structure.  in this paper, we focused on heterogeneous structure. such heterogeneity was captured well by a clustering method. the quasi-linear score was naturally derived by bayes risk consistency between mixed and standard normal distributions. moreover, the quasi-linear score approximates the maximum function and plays an important role in selecting the most effective cluster for prediction from given clusters. the quasi-linear score has better predictive ability compared to linear score as shown in simulation studies and applications to real data. additional files additional file technical derivations. in this file, we perform some technical derivations and evaluations for quasi-linear score: parameterization; the relationship with linear and maximum score; the bayes risk consistency; l and l regularization methods; the derivatives. additional file r source code of the parameter estimation of the quasi-linear logistic model. in this file, we introduce the r source code of the parameter estimation of the quasi-linear logistic model, which was used for simulation and application. abbreviations aucarea under curve dldadiagonal fisher’s linear discriminant analysis mlemaximum likelihood estimator electronic supplementary material the online version of this article contains supplementary material, which is available to authorized users. 
we observed that two shrinkage methods worked well: ridge shrinkage for estimating the quasi-linear score, and lasso shrinkage for selecting markers within each cluster. we extended the linear score to the quasi-linear score by a generalized average form, the kolmogorov-nagumo average. simulation studies and applications to real data show that the proposed method has good predictive performance compared with existing methods. linear scores are widely used to predict dichotomous outcomes in biomedical studies because of their learnability and understandability. accordingly, we developed a quasi-linear score, which effectively combines the linear scores of clustered markers. quasi-linear scores combine such heterogeneity and have a better predictive ability compared with linear scores. the online version of this article contains supplementary material, which is available to authorized users. discriminant analysisheterogeneitykolmogorov-nagumo averagepredictionjapan society for the promotion of science 25280008eguchi shinto issue-copyright-statement© the author 2017 such approaches, however, cannot be used to elucidate biodiversity when there is heterogeneous structure in target population. heterogeneous structure is captured by a clustering method.
 linear scores are widely used to predict dichotomous outcomes in biomedical studies because of their learnability and understandability. such approaches, however, cannot be used to elucidate biodiversity when there is heterogeneous structure in target population.  our study was focused on describing intrinsic heterogeneity in predictions. because heterogeneity can be captured by a clustering method, integrating different information from different clusters should yield better predictions. accordingly, we developed a quasi-linear score, which effectively combines the linear scores of clustered markers. we extended the linear score to the quasi-linear score by a generalized average form, the kolmogorov-nagumo average. we observed that two shrinkage methods worked well: ridge shrinkage for estimating the quasi-linear score, and lasso shrinkage for selecting markers within each cluster. simulation studies and applications to real data show that the proposed method has good predictive performance compared with existing methods.  heterogeneous structure is captured by a clustering method. quasi-linear scores combine such heterogeneity and have a better predictive ability compared with linear scores. electronic supplementary material the online version of this article contains supplementary material, which is available to authorized users. keywords discriminant analysisheterogeneitykolmogorov-nagumo averagepredictionjapan society for the promotion of science 25280008eguchi shinto issue-copyright-statement© the author 2017 