many of the advances in genome - wide association studies ( gwass ) are based on the analysis of case 
the goal of such studies is usually to correlate a phenotype , often coded as a binary variable , with genetic markers such as single nucleotide polymorphisms ( snps ) .
the outcomes of gwass are increasingly integrated with next - generation sequencing results for snp validation and disease marker extraction .
therefore , gwass are powerful tools for the identification of interesting regions at the genome - wide level for further investigation by high - definition procedures.1 regrettably , achieving the promising goals mentioned above is hindered by the technical difficulties accompanying gwass .
recent genotyping technologies allow the extraction of millions of genetic measurements , but the number of subjects is smaller by several orders of magnitude . in the literature , this is known as a  p > n problem . the key task is to reduce the dimensionality of the problem by identifying the optimal subset of predictors that are informative with respect to the outcome variable .
widely used techniques such as univariate or stepwise regression , are recognized as insufficient to fully address the above issues .
advanced multivariate methods are therefore required in order to determine the biological truth.24 recently , a number of model search algorithms have been proposed , with particular emphasis on sparse linear models based on regularization methods5 
 6 and naive - bayes techniques.7 furthermore , bottolo and richardson designed a new method for bayesian model selection for linear regression with continuous outcomes.8 unlike regularization analysis , this method performs a model search by sampling the predictors on the basis of a general and efficient markov chain monte carlo ( mcmc ) technique that exploits the conjugacy structure of data and parameters . in the case of binary outcomes , however , easy analytical tractability is only partially possible , mainly because of the non - linear link between the observed variables and the underlying predictors .
although a detailed bayesian analysis for binary and categorical responses is available,9 model search techniques for binary data have only been explored in a limited number of cases.1012 in the present paper , we propose an innovative stochastic model search technique for binary outcomes .
the relationship between the observed responses and the available predictors is described by a latent variable model with a probit link,9 rather than resorting to gaussian approximations.12 prior distributions are assigned both to the regression coefficients and the model size , therefore allowing the user to specify a prior belief on the model complexity . a computationally efficient metropolis - hastings sampling algorithm8 is adopted .
in the present paper , we propose an innovative stochastic model search technique for binary outcomes .
the relationship between the observed responses and the available predictors is described by a latent variable model with a probit link,9 rather than resorting to gaussian approximations.12 prior distributions are assigned both to the regression coefficients and the model size , therefore allowing the user to specify a prior belief on the model complexity . a computationally efficient metropolis - hastings sampling algorithm8 is adopted .
denote with y the vector of the observed individual binary responses yi , where i=1  n and n is the number of subjects .
the yi are assumed to be the results of independent bernoulli trials with success probability i .
let xi=(xi1 ,  , xip ) be the covariate vector associated with response yi , where p is the total number of predictors .
additionally , let =(1 ,  , p ) be the vector of regression coefficients for the p predictors .
the relationship between observed binary responses and covariates is usually modeled with a link function g , which expresses the dependency of the success probabilities i on the linear regressor xi. taking g as the standard normal distribution  yields the probit model:(1)i=(xit)=xitn(0,1)duwhere  , in a bayesian perspective , can be assigned a multivariate normal prior .
unfortunately , the resulting posterior density of  is not analytically tractable,9 although approximations are available in the literature . in order to circumvent the above issue and
obtain the posterior of  , one can resort to a data augmentation approach as described by albert and chib.9 the idea is to introduce n independent latent variables z=(z1 ,  , zn ) such that each zi  n(xi , 1 ) , letting yi=0 if zi0 and yi=1 if zi>0 .
equation 2 summarizes the latent variable model:(2)yi={0zi01zi>0z = x+where x is the np design matrix with i - th row equal to xi and  is a vector of independent normal errors with zero mean and unit variance . of course , if the zi were known , and if  were assigned a multivariate normal prior , the posterior of  would be obtained through usual bayesian linear regression results .
although the zi are actually not known , their distribution conditional on the data yi is a truncated normal with mean xi and unit variance , the side of truncation depending on the value of yi .
one could therefore resort to gibbs sampling to iteratively sample z and  , thus obtaining the required posteriors .
the goal of variable selection is to model the dependency of y ( or equivalently z ) on a subset of the predictors x1 ,  , xp .
we follow the approach of bottolo and richardson8 and introduce a latent binary vector =(1 ,  , p ) such that j=0 if j=0 and j=1 if j0 .
the model space is therefore given by the 2 possible combinations of the indicator variables j .
given  , the gaussian linear model in equation 2 can therefore be modified as:(3)z=1+x+that also accounts for the intercept . the symbol 1 denotes a column vector of ones .
its length is denoted by p. similarly , the np matrix x is obtained from x by taking the columns for which j=1 .
the intercept  is assigned a flat prior , f()1 , whereas the prior distribution for  is taken as a multivariate normal : moreover , the error variance  is equal to 1 by definition.9 the matrix  can be expressed as i , where  controls the degree of coefficient shrinkage and i is the identity matrix .
other specifications of the prior covariance matrix in equation 4 are covered by bottolo and richardson.8 the prior for the indicator vector  is defined as in kohn et al13:(5)f()=b(p+a , pp+b)b(a , b)which gives rise to a beta - binomial distribution on the model size p. parameters a and b can be elicited based on prior knowledge on the model size , for example , expected value and variance of p.8 
 13 in order to obtain  , observe that , if z were known , one could write its associated marginal likelihood by integrating out  and :(6)f(z|)=f(z|,,)f(|)f()ddwhere the dependence on  has been dropped since it is fixed to 1 . the term f(z |  ,  ,  ) is the multivariate normal density arising from the distribution of  . note that , up to a normalizing factor , equations 5 and 6 together yield the conditional density f( | z ) . also observe that such density does not depend on y since z is given . additionally , observe that , if z were known , the posterior of  would follow from the usual bayes formula for linear gaussian regression .
the model specification is completed by defining the full conditional for z. as stated in the previous paragraph , conditional on yi and  , the densities of zi are independent truncated normals:(7)f(zi|yi,)={tn(xi,t,1)yi=0tn+(xi,t,1)yi=1where tn and tn denote right- and left - truncated normal densities respectively , and xi, is the i - th row of x. equations 47 allow devising a sampling scheme to obtain the posterior density of the latent binary vector  , which is the main goal of the model search procedure . our algorithm is based on an efficient fast scan metropolis - hastings ( fsmh ) sampler introduced by bottolo and richardson8 and generalized here to account for binary responses .
the key idea of our algorithm is that , if z is assumed known , the variable selection problem for binary outcomes reduces to one for continuous outcomes .
in such a case , efficient algorithms are available to sample from the conditional density of .8 
 12 as shown by bottolo and richardson,8 sampling  can be achieved through the fsmh sampler using equations 5 and 6 .
once a sample of  is available , a new value of z is then sampled from its full conditional distribution and the process is iterated . the sampling scheme , depicted in figure 1 ,
squares denote constants ( eg , the design matrix ) , and circles denote random variables .
choose a random initialization for z ( positive or negative according to its corresponding yi , see equation 7 ) ; sample  through the fsmh sampler , using the prior model size in equation 5 and the marginal likelihood of z in equation 6 ; given z and  , sample  through the standard bayesian linear regression .
note that , given z , such density does not depend on y ; given  and y , sample z from equation 7 ; restart from step 2 until a pre - determined number of iterations is reached .
the proposed approach has been applied to a first gwas on the genetic bases of longevity,14 and a second gwas aimed at the genetic dissection of the type 2 diabetes ( t2d ) trait.15 the latter dataset has been provided by the wellcome trust case control consortium ( wtccc ) after data access approval .
individual and marker - level data quality control , inbred patients , and genetic outliers removal have been described by malovini et al14 and the wtccc.15 a preliminary feature selection has been performed based on the results from univariate pearson 's  tests with 2 degrees of freedom comparing genotype distributions between cases and controls . only snps characterized by complete genotype data , passing the significance threshold ( p<0.40 for the longevity dataset and p<0.01 for the t2d dataset ) and featuring at least five counts per cell have been retained .
missing values in the t2d dataset were imputed to the modal value for each snp .
the above statistical analyses have been performed with plink.16 boss was implemented in matlab 7.5.0 r2007b17 together with the randraw routine18 to sample from the truncated normal density .
denote with y the vector of the observed individual binary responses yi , where i=1  n and n is the number of subjects .
the yi are assumed to be the results of independent bernoulli trials with success probability i .
let xi=(xi1 ,  , xip ) be the covariate vector associated with response yi , where p is the total number of predictors .
additionally , let =(1 ,  , p ) be the vector of regression coefficients for the p predictors .
the relationship between observed binary responses and covariates is usually modeled with a link function g , which expresses the dependency of the success probabilities i on the linear regressor xi. taking g as the standard normal distribution  yields the probit model:(1)i=(xit)=xitn(0,1)duwhere  , in a bayesian perspective , can be assigned a multivariate normal prior .
unfortunately , the resulting posterior density of  is not analytically tractable,9 although approximations are available in the literature . in order to circumvent the above issue and
obtain the posterior of  , one can resort to a data augmentation approach as described by albert and chib.9 the idea is to introduce n independent latent variables z=(z1 ,  , zn ) such that each zi  n(xi , 1 ) , letting yi=0 if zi0 and yi=1 if zi>0 .
equation 2 summarizes the latent variable model:(2)yi={0zi01zi>0z = x+where x is the np design matrix with i - th row equal to xi and  is a vector of independent normal errors with zero mean and unit variance . of course , if the zi were known , and if  were assigned a multivariate normal prior , the posterior of  would be obtained through usual bayesian linear regression results .
although the zi are actually not known , their distribution conditional on the data yi is a truncated normal with mean xi and unit variance , the side of truncation depending on the value of yi .
one could therefore resort to gibbs sampling to iteratively sample z and  , thus obtaining the required posteriors .
the goal of variable selection is to model the dependency of y ( or equivalently z ) on a subset of the predictors x1 ,  , xp . however , there is uncertainty about which subset to use . in this work ,
we follow the approach of bottolo and richardson8 and introduce a latent binary vector =(1 ,  , p ) such that j=0 if j=0 and j=1 if j0 .
the model space is therefore given by the 2 possible combinations of the indicator variables j .
given  , the gaussian linear model in equation 2 can therefore be modified as:(3)z=1+x+that also accounts for the intercept . the symbol 1 denotes a column vector of ones .
its length is denoted by p. similarly , the np matrix x is obtained from x by taking the columns for which j=1 .
the intercept  is assigned a flat prior , f()1 , whereas the prior distribution for  is taken as a multivariate normal : moreover , the error variance  is equal to 1 by definition.9 the matrix  can be expressed as i , where  controls the degree of coefficient shrinkage and i is the identity matrix .
other specifications of the prior covariance matrix in equation 4 are covered by bottolo and richardson.8 the prior for the indicator vector  is defined as in kohn et al13:(5)f()=b(p+a , pp+b)b(a , b)which gives rise to a beta - binomial distribution on the model size p. parameters a and b can be elicited based on prior knowledge on the model size , for example , expected value and variance of p.8 
 13 in order to obtain  , observe that , if z were known , one could write its associated marginal likelihood by integrating out  and :(6)f(z|)=f(z|,,)f(|)f()ddwhere the dependence on  has been dropped since it is fixed to 1 . the term f(z |  ,  ,  ) is the multivariate normal density arising from the distribution of  .
note that , up to a normalizing factor , equations 5 and 6 together yield the conditional density f( | z ) . also observe that such density does not depend on y since z is given .
additionally , observe that , if z were known , the posterior of  would follow from the usual bayes formula for linear gaussian regression .
the model specification is completed by defining the full conditional for z. as stated in the previous paragraph , conditional on yi and  , the densities of zi are independent truncated normals:(7)f(zi|yi,)={tn(xi,t,1)yi=0tn+(xi,t,1)yi=1where tn and tn denote right- and left - truncated normal densities respectively , and xi, is the i - th row of x.
equations 47 allow devising a sampling scheme to obtain the posterior density of the latent binary vector  , which is the main goal of the model search procedure .
our algorithm is based on an efficient fast scan metropolis - hastings ( fsmh ) sampler introduced by bottolo and richardson8 and generalized here to account for binary responses .
the key idea of our algorithm is that , if z is assumed known , the variable selection problem for binary outcomes reduces to one for continuous outcomes .
in such a case , efficient algorithms are available to sample from the conditional density of .8 
 12 as shown by bottolo and richardson,8 sampling  can be achieved through the fsmh sampler using equations 5 and 6 .
once a sample of  is available , a new value of z is then sampled from its full conditional distribution and the process is iterated .
the sampling scheme , depicted in figure 1 , is implemented as follows : relationships among the stochastic nodes sampled by boss .
squares denote constants ( eg , the design matrix ) , and circles denote random variables .
choose a random initialization for z ( positive or negative according to its corresponding yi , see equation 7 ) ; sample  through the fsmh sampler , using the prior model size in equation 5 and the marginal likelihood of z in equation 6 ; given z and  , sample  through the standard bayesian linear regression .
note that , given z , such density does not depend on y ; given  and y , sample z from equation 7 ; restart from step 2 until a pre - determined number of iterations is reached .
the proposed approach has been applied to a first gwas on the genetic bases of longevity,14 and a second gwas aimed at the genetic dissection of the type 2 diabetes ( t2d ) trait.15 the latter dataset has been provided by the wellcome trust case control consortium ( wtccc ) after data access approval .
individual and marker - level data quality control , inbred patients , and genetic outliers removal have been described by malovini et al14 and the wtccc.15 a preliminary feature selection has been performed based on the results from univariate pearson 's  tests with 2 degrees of freedom comparing genotype distributions between cases and controls . only snps characterized by complete genotype data , passing the significance threshold ( p<0.40 for the longevity dataset and p<0.01 for the t2d dataset ) and featuring at least five counts per cell have been retained .
missing values in the t2d dataset were imputed to the modal value for each snp .
the above statistical analyses have been performed with plink.16 boss was implemented in matlab 7.5.0 r2007b17 together with the randraw routine18 to sample from the truncated normal density .
a simulation study was performed to assess the performances of boss , stepwise regression ( stepwisefit function in matlab17 ) , and two state - of - the - art methods , logistic lasso and elastic net11 ( package glmnet in r v.2.13.119 ) .
our benchmark featured simulated scenarios of varying complexity , according to the following characteristics of the design matrix x : uncorrelated or correlated columns ( pairwise correlation of 0.5 ) ; percentage of explained variance ( ie , proportion of variance of z explained by the predictors ) .
therefore , we simulated 23=6 scenarios , labeled as follows : a90 : uncorrelated x , 90% explained variance ; a70 : uncorrelated x , 70% explained variance ; a50 : uncorrelated x , 50% explained variance ; b90 : correlated x , 90% explained variance ; b70 : correlated x , 70% explained variance ; b50 : correlated x , 50% explained variance .
synthetic datasets were generated by simulating a design matrix x according to the characteristics described above , and a vector of independent normally distributed errors  with zero mean and unit variance .
only the first 10 parameters of the vector  were set to a non - zero value , so as to yield 10 true predictors .
the module of such parameters was set so as to obtain the given percentage of explained variance .
the values of zi so generated were then binarized according to their sign to obtain yi ( equation 2 ) .
the prior model size p was specified by e[p]=20 and var[p]=100 , which elicited the values a=4.52 and b=63.26 .
parameter  was fixed to 1 , as done by bottolo and richardson8 and hans et al,12 except in scenarios a50 and b50 where a value of 0.1 proved more suitable . in principle
the prior mean and variance were chosen so as to cover a reasonable range of plausible model sizes .
logistic lasso and elastic net were run with the  one - standard - error rule to choose the penalty parameter.11 the elastic net method was executed with equal lasso and ridge - regression penalties .
the goal of the simulated benchmark was to assess the capabilities of boss , logistic lasso , elastic net , and stepwise regression to capture the true predictors used to generate the data , while discarding useless predictors .
for this purpose , we evaluated precision , recall , and f - measure relative to the number of predictors correctly / incorrectly identified as true / false . in the case of boss , a predictor was deemed 
results of the comparison are reported in figure 3 for uncorrelated predictors ( a90 , a70 , and a50 ) and in figure 4 for correlated predictors ( b90 , b70 , and b50 ) .
the two figures show the distribution of f - measure values obtained with boss , logistic lasso , elastic net , and stepwise regression in each scenario .
f - measure comparison in scenarios a90 , a70 , and a50 ( uncorrelated predictors ) .
f - measure comparison in scenarios b90 , b70 , and b50 ( correlated predictors ) .
en , elastic net ; ll , logistic lasso ; sw , stepwise regression . comparison of precision , recall , and f - measure obtained with stepwise regression ( sw ) , elastic net ( en ) , logistic lasso ( ll ) and boss each value is the average across the 50 datasets simulated in each scenario . in terms of f - measure ,
boss performed better than the three reference methods ( p<0.01 , one - sided paired t test ) . in most scenarios considered here
although boss achieved generally lower recall rates , it attained higher values of precision and a better overall tradeoff .
larger priors on the model size and a larger number of iterations did not substantially affect the results obtained with boss .
interestingly , the performances of the four methods did not degrade if predictors were correlated .
overall , logistic lasso proved superior to the traditional stepwise regression and comparable to the elastic net .
furthermore , we assessed the ability of the four methods to correctly predict newly observed binary responses that were not used for parameter estimation .
prediction performances were evaluated in terms of correctly / incorrectly predicted positive / negative responses . a new outcome was predicted as positive ( y=1 ) if the corresponding probability , reconstructed with the selected features , was greater than a given threshold .
we then evaluated specificity , sensitivity , and prediction accuracy for thresholds from 0.05 to 0.95 , with steps of 0.05 . in order to obtain summary statistics , we computed the area under the receiver operating characteristic ( roc auc ) using mean specificity and sensitivity for a given threshold . additionally , we calculated the average prediction accuracy across all thresholds .
the auc obtained with boss is only slightly greater than with the other three methods . in terms of accuracy
comparison of prediction performances on new data using the selected features except where indicated ( * ) , boss improved on the three reference methods ( p<0.05 , one - sided paired t test ) .
en , elastic net ; ll , logistic lasso ; roc auc , area under the receiver operating characteristic ; sw , stepwise regression .
this section reports the results obtained by boss and the three reference methods on the two experimental gwass .
a total of 410 long aged individuals ( cases ) , 553 average living controls and 290  364 autosomal snps passing data quality control underwent a feature selection pre - processing phase . a subset of 5707 snps , selected according to the inclusion criteria reported in the materials and methods section , were then considered for subsequent analyses .
boss was run for 20  000 iterations ( first 5000 of burn - in ; running time of about 4.5  days ) .
the prior model size p was specified by e[p]=20 and var[p]=100 , so as to encompass a range of plausible model sizes without being overly restrictive .
in particular , the intergenic snp rs2147556 achieved a 100% probability of inclusion , while rs10491334 , mapping to camk4 and previously identified as strongly associated,14 achieved a probability of about 76% .
moreover , boss was able to identify rs522796 ( kl gene ) as mildly associated with the outcome ( probability of 14% ) .
this snp was reported to be related to non - diabetic end - stage renal disease20 and preterm birth.21 snp rs800292 , mapping to cfh and known to be associated with age - related macular degeneration,22 was also identified , although with a relatively low probability of inclusion ( 10% ) .
we limited the plot to the top 10 snps in order to summarize the most relevant findings only .
marginal posterior probabilities of the 10 top single nucleotide polymorphisms ( snps ) identified by boss in the longevity dataset .
darker bars denote snps previously reported in the literature . in order to further evaluate the results obtained by the proposed model search technique
, we additionally fitted a multivariate logistic regression model on the top 10 predictors of figure 5 .
summary of results from logistic regression on the top 10 predictors identified by boss in the longevity dataset chr , chromosome ; position ( bp ) , physical position on the chromosome expressed in terms of base pairs ; snp , single nucleotide polymorphism .
it is interesting to observe that the minor allele of rs800292 has a predisposing additive impact on longevity ( or=1.53 ) : this is in accordance with the results obtained by yang et al , who showed a significant association of rs800292 with a reduced risk for exudative age - related macular degeneration.22 as a further step , we ran logistic lasso , elastic net , and stepwise regression on the longevity data to assess possible differences with the results obtained by boss .
the three algorithms were run with the same settings described in the simulated benchmark section of the results section . in our analysis ,
the logistic lasso included most of the features reported in figure 5 , but failed to detect rs800292 .
the elastic net yielded a more parsimonious model , which included four snps ( rs2147556 , rs10491334 , rs621011 , and rs522796 ) .
stepwise regression selected a large number of features ( 891 ) : all snps in figure 5 were included except rs522796 .
a total of 1924 t2d patients , 1458 uk blood service ( ukbs ) controls and 456  856 autosomal snps ( mapping to chromosomes 122 ) passing data quality control15 underwent a features selection phase ( see materials and methods section ) .
we ran boss for 700 iterations ( first 200 of burn - in ; running time of about 1  day ) . increasing the number of iterations did not yield appreciably different results .
the prior model size p was specified by e[p]=4 and var[p]=4 , in order to allow boss identify only the most relevant genetic features .
our analysis detected several snps mapping to gene regions which have been previously associated with metabolic syndromes or other disease classes .
in particular , the top ranked snp rs7193144 is an intronic marker mapping to fto , a gene that has been already associated with t2d and obesity conditions in several european cohorts.15 
 2325 moreover , rs4132670 and rs10885409 map to tcf7l2 , previously reported to be associated with t2d.15 
 26 
 27 last , rs11693602 is localized within the rbms1 gene , which has been associated with t2d by qi et al.28 posterior probabilities of snps obtained with boss are shown in figure 6 .
we limited the plot to the top 20 snps in order to summarize the most relevant findings only .
marginal posterior probabilities of the 20 top single nucleotide polymorphisms ( snps ) identified by boss in the type 2 diabetes ( t2d ) dataset .
darker bars denote snps previously reported in the literature . additionally , we performed further model evaluation by fitting a multivariate logistic regression on the genetic markers reported in figure 6 .
summary of results from logistic regression on the top 20 predictors identified by boss in the t2d dataset chr , chromosome ; position ( bp ) , physical position on the chromosome expressed in terms of base pairs ; snp , single nucleotide polymorphism .
we also performed logistic lasso , elastic net , and stepwise regression on the t2d data .
however , stepwise regression and logistic lasso detected none of the four snps introduced above . the elastic net missed rs10885409 , but captured rs7193144 , rs11693602 , and rs4132670 .
a simulation study was performed to assess the performances of boss , stepwise regression ( stepwisefit function in matlab17 ) , and two state - of - the - art methods , logistic lasso and elastic net11 ( package glmnet in r v.2.13.119 ) .
our benchmark featured simulated scenarios of varying complexity , according to the following characteristics of the design matrix x : uncorrelated or correlated columns ( pairwise correlation of 0.5 ) ; percentage of explained variance ( ie , proportion of variance of z explained by the predictors ) .
therefore , we simulated 23=6 scenarios , labeled as follows : a90 : uncorrelated x , 90% explained variance ; a70 : uncorrelated x , 70% explained variance ; a50 : uncorrelated x , 50% explained variance ; b90 : correlated x , 90% explained variance ; b70 : correlated x , 70% explained variance ; b50 : correlated x , 50% explained variance .
synthetic datasets were generated by simulating a design matrix x according to the characteristics described above , and a vector of independent normally distributed errors  with zero mean and unit variance .
only the first 10 parameters of the vector  were set to a non - zero value , so as to yield 10 true predictors .
the module of such parameters was set so as to obtain the given percentage of explained variance .
the values of zi so generated were then binarized according to their sign to obtain yi ( equation 2 ) .
the prior model size p was specified by e[p]=20 and var[p]=100 , which elicited the values a=4.52 and b=63.26 .
parameter  was fixed to 1 , as done by bottolo and richardson8 and hans et al,12 except in scenarios a50 and b50 where a value of 0.1 proved more suitable . in principle
the prior mean and variance were chosen so as to cover a reasonable range of plausible model sizes .
logistic lasso and elastic net were run with the  one - standard - error rule to choose the penalty parameter.11 the elastic net method was executed with equal lasso and ridge - regression penalties .
the goal of the simulated benchmark was to assess the capabilities of boss , logistic lasso , elastic net , and stepwise regression to capture the true predictors used to generate the data , while discarding useless predictors .
for this purpose , we evaluated precision , recall , and f - measure relative to the number of predictors correctly / incorrectly identified as true / false . in the case of boss , a predictor was deemed 
results of the comparison are reported in figure 3 for uncorrelated predictors ( a90 , a70 , and a50 ) and in figure 4 for correlated predictors ( b90 , b70 , and b50 ) .
the two figures show the distribution of f - measure values obtained with boss , logistic lasso , elastic net , and stepwise regression in each scenario .
f - measure comparison in scenarios a90 , a70 , and a50 ( uncorrelated predictors ) .
f - measure comparison in scenarios b90 , b70 , and b50 ( correlated predictors ) .
en , elastic net ; ll , logistic lasso ; sw , stepwise regression . comparison of precision , recall , and f - measure obtained with stepwise regression ( sw ) , elastic net ( en ) , logistic lasso ( ll ) and boss each value is the average across the 50 datasets simulated in each scenario . in terms of f - measure ,
boss performed better than the three reference methods ( p<0.01 , one - sided paired t test ) . in most scenarios considered here
although boss achieved generally lower recall rates , it attained higher values of precision and a better overall tradeoff .
larger priors on the model size and a larger number of iterations did not substantially affect the results obtained with boss .
interestingly , the performances of the four methods did not degrade if predictors were correlated .
overall , logistic lasso proved superior to the traditional stepwise regression and comparable to the elastic net .
furthermore , we assessed the ability of the four methods to correctly predict newly observed binary responses that were not used for parameter estimation .
prediction performances were evaluated in terms of correctly / incorrectly predicted positive / negative responses . a new outcome was predicted as positive ( y=1 ) if the corresponding probability , reconstructed with the selected features , was greater than a given threshold .
we then evaluated specificity , sensitivity , and prediction accuracy for thresholds from 0.05 to 0.95 , with steps of 0.05 . in order to obtain summary statistics , we computed the area under the receiver operating characteristic ( roc auc ) using mean specificity and sensitivity for a given threshold . additionally , we calculated the average prediction accuracy across all thresholds .
the auc obtained with boss is only slightly greater than with the other three methods . in terms of accuracy
comparison of prediction performances on new data using the selected features except where indicated ( * ) , boss improved on the three reference methods ( p<0.05 , one - sided paired t test ) .
en , elastic net ; ll , logistic lasso ; roc auc , area under the receiver operating characteristic ; sw , stepwise regression .
this section reports the results obtained by boss and the three reference methods on the two experimental gwass .
a total of 410 long aged individuals ( cases ) , 553 average living controls and 290  364 autosomal snps passing data quality control underwent a feature selection pre - processing phase . a subset of 5707 snps , selected according to the inclusion criteria reported in the materials and methods section , were then considered for subsequent analyses .
boss was run for 20  000 iterations ( first 5000 of burn - in ; running time of about 4.5  days ) .
the prior model size p was specified by e[p]=20 and var[p]=100 , so as to encompass a range of plausible model sizes without being overly restrictive .
in particular , the intergenic snp rs2147556 achieved a 100% probability of inclusion , while rs10491334 , mapping to camk4 and previously identified as strongly associated,14 achieved a probability of about 76% .
moreover , boss was able to identify rs522796 ( kl gene ) as mildly associated with the outcome ( probability of 14% ) .
this snp was reported to be related to non - diabetic end - stage renal disease20 and preterm birth.21 snp rs800292 , mapping to cfh and known to be associated with age - related macular degeneration,22 was also identified , although with a relatively low probability of inclusion ( 10% ) .
we limited the plot to the top 10 snps in order to summarize the most relevant findings only .
marginal posterior probabilities of the 10 top single nucleotide polymorphisms ( snps ) identified by boss in the longevity dataset .
darker bars denote snps previously reported in the literature . in order to further evaluate the results obtained by the proposed model search technique
, we additionally fitted a multivariate logistic regression model on the top 10 predictors of figure 5 .
summary of results from logistic regression on the top 10 predictors identified by boss in the longevity dataset chr , chromosome ; position ( bp ) , physical position on the chromosome expressed in terms of base pairs ; snp , single nucleotide polymorphism .
it is interesting to observe that the minor allele of rs800292 has a predisposing additive impact on longevity ( or=1.53 ) : this is in accordance with the results obtained by yang et al , who showed a significant association of rs800292 with a reduced risk for exudative age - related macular degeneration.22 as a further step , we ran logistic lasso , elastic net , and stepwise regression on the longevity data to assess possible differences with the results obtained by boss .
the three algorithms were run with the same settings described in the simulated benchmark section of the results section . in our analysis ,
the logistic lasso included most of the features reported in figure 5 , but failed to detect rs800292 .
the elastic net yielded a more parsimonious model , which included four snps ( rs2147556 , rs10491334 , rs621011 , and rs522796 ) .
stepwise regression selected a large number of features ( 891 ) : all snps in figure 5 were included except rs522796 .
a total of 1924 t2d patients , 1458 uk blood service ( ukbs ) controls and 456  856 autosomal snps ( mapping to chromosomes 122 ) passing data quality control15 underwent a features selection phase ( see materials and methods section ) .
we ran boss for 700 iterations ( first 200 of burn - in ; running time of about 1  day ) . increasing the number of iterations did not yield appreciably different results .
the prior model size p was specified by e[p]=4 and var[p]=4 , in order to allow boss identify only the most relevant genetic features .
our analysis detected several snps mapping to gene regions which have been previously associated with metabolic syndromes or other disease classes .
in particular , the top ranked snp rs7193144 is an intronic marker mapping to fto , a gene that has been already associated with t2d and obesity conditions in several european cohorts.15 
 2325 moreover , rs4132670 and rs10885409 map to tcf7l2 , previously reported to be associated with t2d.15 
 26 
 27 last , rs11693602 is localized within the rbms1 gene , which has been associated with t2d by qi et al.28 posterior probabilities of snps obtained with boss are shown in figure 6 .
we limited the plot to the top 20 snps in order to summarize the most relevant findings only .
marginal posterior probabilities of the 20 top single nucleotide polymorphisms ( snps ) identified by boss in the type 2 diabetes ( t2d ) dataset .
darker bars denote snps previously reported in the literature . additionally , we performed further model evaluation by fitting a multivariate logistic regression on the genetic markers reported in figure 6 .
summary of results from logistic regression on the top 20 predictors identified by boss in the t2d dataset chr , chromosome ; position ( bp ) , physical position on the chromosome expressed in terms of base pairs ; snp , single nucleotide polymorphism .
we also performed logistic lasso , elastic net , and stepwise regression on the t2d data .
however , stepwise regression and logistic lasso detected none of the four snps introduced above . the elastic net missed rs10885409 , but captured rs7193144 , rs11693602 , and rs4132670 .
a total of 410 long aged individuals ( cases ) , 553 average living controls and 290  364 autosomal snps passing data quality control underwent a feature selection pre - processing phase . a subset of 5707 snps , selected according to the inclusion criteria reported in the materials and methods section , were then considered for subsequent analyses .
boss was run for 20  000 iterations ( first 5000 of burn - in ; running time of about 4.5  days ) .
the prior model size p was specified by e[p]=20 and var[p]=100 , so as to encompass a range of plausible model sizes without being overly restrictive .
in particular , the intergenic snp rs2147556 achieved a 100% probability of inclusion , while rs10491334 , mapping to camk4 and previously identified as strongly associated,14 achieved a probability of about 76% .
moreover , boss was able to identify rs522796 ( kl gene ) as mildly associated with the outcome ( probability of 14% ) .
this snp was reported to be related to non - diabetic end - stage renal disease20 and preterm birth.21 snp rs800292 , mapping to cfh and known to be associated with age - related macular degeneration,22 was also identified , although with a relatively low probability of inclusion ( 10% ) .
we limited the plot to the top 10 snps in order to summarize the most relevant findings only .
marginal posterior probabilities of the 10 top single nucleotide polymorphisms ( snps ) identified by boss in the longevity dataset .
darker bars denote snps previously reported in the literature . in order to further evaluate the results obtained by the proposed model search technique
, we additionally fitted a multivariate logistic regression model on the top 10 predictors of figure 5 .
summary of results from logistic regression on the top 10 predictors identified by boss in the longevity dataset chr , chromosome ; position ( bp ) , physical position on the chromosome expressed in terms of base pairs ; snp , single nucleotide polymorphism .
it is interesting to observe that the minor allele of rs800292 has a predisposing additive impact on longevity ( or=1.53 ) : this is in accordance with the results obtained by yang et al , who showed a significant association of rs800292 with a reduced risk for exudative age - related macular degeneration.22 as a further step , we ran logistic lasso , elastic net , and stepwise regression on the longevity data to assess possible differences with the results obtained by boss .
the three algorithms were run with the same settings described in the simulated benchmark section of the results section . in our analysis ,
the logistic lasso included most of the features reported in figure 5 , but failed to detect rs800292 .
the elastic net yielded a more parsimonious model , which included four snps ( rs2147556 , rs10491334 , rs621011 , and rs522796 ) .
stepwise regression selected a large number of features ( 891 ) : all snps in figure 5 were included except rs522796 .
a total of 1924 t2d patients , 1458 uk blood service ( ukbs ) controls and 456  856 autosomal snps ( mapping to chromosomes 122 ) passing data quality control15 underwent a features selection phase ( see materials and methods section ) .
we ran boss for 700 iterations ( first 200 of burn - in ; running time of about 1  day ) . increasing the number of iterations did not yield appreciably different results .
the prior model size p was specified by e[p]=4 and var[p]=4 , in order to allow boss identify only the most relevant genetic features .
our analysis detected several snps mapping to gene regions which have been previously associated with metabolic syndromes or other disease classes .
in particular , the top ranked snp rs7193144 is an intronic marker mapping to fto , a gene that has been already associated with t2d and obesity conditions in several european cohorts.15 
 2325 moreover , rs4132670 and rs10885409 map to tcf7l2 , previously reported to be associated with t2d.15 
 26 
 27 last , rs11693602 is localized within the rbms1 gene , which has been associated with t2d by qi et al.28 posterior probabilities of snps obtained with boss are shown in figure 6 .
we limited the plot to the top 20 snps in order to summarize the most relevant findings only .
marginal posterior probabilities of the 20 top single nucleotide polymorphisms ( snps ) identified by boss in the type 2 diabetes ( t2d ) dataset .
darker bars denote snps previously reported in the literature . additionally , we performed further model evaluation by fitting a multivariate logistic regression on the genetic markers reported in figure 6 .
summary of results from logistic regression on the top 20 predictors identified by boss in the t2d dataset chr , chromosome ; position ( bp ) , physical position on the chromosome expressed in terms of base pairs ; snp , single nucleotide polymorphism .
we also performed logistic lasso , elastic net , and stepwise regression on the t2d data .
however , stepwise regression and logistic lasso detected none of the four snps introduced above . the elastic net missed rs10885409 , but captured rs7193144 , rs11693602 , and rs4132670 .
the results obtained offer an interesting insight into our new model search technique and its performances relative to three established methods for feature selection . in the simulated benchmark ,
stepwise regression captured many features , including important ones , although at the cost of reduced precision .
logistic lasso performed better than stepwise regression and comparably to the elastic net , although both techniques attained relatively low precisions ( 0.260.59 ; see table 1 ) . by contrast , as far as this simulation study is concerned , boss was able to discover true signals even in difficult scenarios with as low as 50% of explained variance .
note that , even if the prior model size allowed for possibly large models ( up to about 60 predictors , see figure 2 ) , boss was able to identify a small subset of key features that correlated well with the data , while discarding unimportant predictors .
our analysis suggests that the greater ability of boss to detect the true predictors entails increased prediction accuracy on newly observed data .
it is worth remarking that , in principle , boss is also capable of handling interactions between features .
although this was not extensively explored , a preliminary simulated analysis ( similar to the a90 scenario ) showed that boss successfully captured both the main and interaction effects between pairs of features . in the longevity study ,
boss effectively accomplished variable selection with many possible predictors ( about 5700 snps ) and a smaller number of observations ( n=963 ) .
boss successfully identified relevant genetic markers and yielded a parsimonious model , which can be used for predictive purposes .
in particular , rs10491334 has been previously associated with longevity in the same population,14 while rs800292 has been previously associated with age - related macular degeneration in a cohort of han chinese individuals.22 moreover , rs522796 has been reported both in relation to non - diabetic end - stage renal disease20 and preterm birth.21 the results from logistic lasso , elastic net , and stepwise regression obtained on the longevity data were comparable , although neither the lasso nor the elastic net were able to detect rs800292 .
stepwise regression yielded an overly parametrized model , with almost as many features ( 891 ) as subjects ( 963 ) , but missed the association of rs522796 . in the t2d study ,
boss identified snps mapping to genes known to be associated with diabetes or other metabolic syndromes.15 
 2328 further , boss identified functionally relevant variants that would have been discarded by standard univariate association tests : for example , the top snp detected by boss , rs7193144 , achieved a univariate p value of 2.6110 (  distribution with 2 degrees of freedom ) , which is greater than usual significance thresholds , for example , 110 or less .
a direct comparison of our findings with those reported by the wtccc15 is only partially feasible , since our research group has no access to the wtccc 1958 british birth cohort , which was included in the original analyses as the control group .
locus - based analyses will allow identification of structural correlations and/or potential interactions between the selected snps and other functionally relevant unobserved variants .
interestingly , stepwise regression and logistic lasso were able to detect a large set of features , but failed to identify the four snps reported in the results section .
this may be explained by convergence to a local optimum for the stepwise regression , and excessive coefficient shrinkage for the lasso .
the elastic net did not detect rs10885409 but found the other three markers described above .
it is a known in the literature that advanced sampling - based techniques for model selection may be computationally demanding ( see , for example , o'hara and sillanp29 ) , especially if compared to lasso - based methods.11 in the two experimental studies , we preliminarily filtered the genetic features to approximately the same number ( about 4000 and 6000 ) , so as to demonstrate the effectiveness of our approach while keeping the computational complexity reasonable .
of note , variational bayes approaches ( see , for example , girolami et al30 ) could be adapted for model selection to further speed up the model search .
the algorithm presented in this paper represents a step forward in the statistical methodology for variable selection .
binary outcomes have traditionally received less attention than continuous outcomes , mainly because of their lower analytical tractability . moreover ,
n problem requires advanced techniques that allow discovery of important features without yielding overparametrized models .
we tackled such challenges using a bayesian paradigm and developed a general mcmc framework for variable selection in the presence of binary responses .
the exact formulation of the likelihood paves the way to generalizations such as multinomial responses that can not be handled simply by gaussian approximation .
the results obtained from our simulated benchmark and the two experimental case studies suggest that boss is an effective and reliable tool for model selection when the number of features is very large .
 objectivethe spread of case  control genome - wide association studies ( gwass ) has stimulated the development of new variable selection methods and predictive models . 
 we introduce a novel bayesian model search algorithm , binary outcome stochastic search ( boss ) , which addresses the model selection problem when the number of predictors far exceeds the number of binary responses.materials and methodsour method is based on a latent variable model that links the observed outcomes to the underlying genetic variables . 
 a markov chain monte carlo approach is used for model search and to evaluate the posterior probability of each predictor.resultsboss is compared with three established methods ( stepwise regression , logistic lasso , and elastic net ) in a simulated benchmark . 
 two real case studies are also investigated : a gwas on the genetic bases of longevity , and the type 2 diabetes study from the wellcome trust case control consortium . 
 simulations show that boss achieves higher precisions than the reference methods while preserving good recall rates . in both experimental studies , 
 boss successfully detects genetic polymorphisms previously reported to be associated with the analyzed phenotypes.discussionboss outperforms the other methods in terms of f - measure on simulated data . in the two real studies , 
 boss successfully detects biologically relevant features , some of which are missed by univariate analysis and the three reference techniques.conclusionthe proposed algorithm is an advance in the methodology for model selection with a large number of features . 
 our simulated and experimental results showed that boss proves effective in detecting relevant markers while providing a parsimonious model . 
Introduction and background
Objective
Materials and methods
Latent variable model for binary data
Variable selection
Sampling algorithm: Binary Outcome Stochastic Search (BOSS)
Experimental data and software implementation
Results
Simulated benchmark
Experimental case studies
GWAS on exceptional longevity
GWAS on T2D
Discussion
Conclusion