A Scalable Empirical Bayes Approach to Variable Selection in Generalized Linear Models
Haim Bar James Booth Martin T. Wells


#1 1




0
  
  
  A Scalable Empirical Bayes Approach to Variable
Selection in  Generalized Linear Models
    Haim Y. Bar
    The authors gratefully acknowledge the following funding support: Prof. Bar's research was supported by NSF-DMS
1612625. Professor Booth's research was partially supported
by an NSF grant, NSF-DMS 1208488, and by NSF-DMS 1611893. Professor Wells' research was partially supported by
NSF-DMS 1208488, NSF-DMS 1611893, and NIH grant U19 AI111143.

    Department of Statistics University of Connecticut,
 Storrs CT, USA.

    and 

    James G. Booth and Martin T. Wells

    Department of Statistics and Data Science, Cornell University, 
 Ithaca NY, 14853, USA.
    December 30, 2023
=====================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================

 

1
  

    A Scalable Empirical Bayes Approach to Variable
Selection in  Generalized Linear Models


   




A new empirical Bayes approach to variable selection in the context of
generalized linear models is developed. The proposed algorithm scales to
situations in which the number of putative explanatory variables is very large, possibly
much larger than the number of responses. The coefficients
in the linear predictor are modeled as a three-component mixture
allowing the explanatory variables to have a random positive effect on the response, a
random negative effect, or no effect. A key assumption is that only a
small (but unknown) fraction of the candidate variables have a
non-zero effect. This assumption, in addition to treating the
coefficients as random effects facilitates an approach that is computationally
efficient. In particular, the number of parameters that have to be estimated is small, and
remains constant regardless of the number of explanatory variables.
The model parameters are estimated using a Generalized Alternating Maximization
algorithm which
is scalable, and leads to significantly faster convergence compared
with simulation-based fully Bayesian methods.



Keywords:  Feature selection; Generalized linear mixed model;
 High dimensional data; EM algorithm; Mixture model; Sparsity



1.45 



§ INTRODUCTION

This paper concerns variable selection in generalized linear regression models
when there are a large number of candidate explanatory variables
(putative predictors), most of which
have little or no effect on the dependent variable. An empirical
Bayes, model-based approach is proposed that is implemented via a fast
and scalable Generalized Alternating Maximization algorithm.

The new age of high speed computing and technological advances in
genetics and molecular biology, for example, have dramatically changed
modeling and computation needs. It is now common for researchers
to want to estimate the effects of hundreds or
even thousands of predictors (K) on a given response, often with a far
smaller sample size (N).  In such cases, traditional fitting methods such
as least squares break down. In addition, even with a relatively modest
number of predictors, the model space can be large enough to render exhaustive
search-based algorithms impractical.

Automated methods for variable selection in normal linear regression models have
long been studied in the literature; see, for example
<cit.>. Nowadays virtually every
statistical package contains an implementation of standard stepwise
methods that typically add or remove one variable from the model in each
iteration, based on sequential F-tests and a threshold, or a well-known
selection criterion such as AIC , BIC, or Mallow's-C_p.
A modern alternative is to use false discovery rate (FDR) for
stepwise model selection <cit.>.

Much of the recent literature has focused on variations of penalized likelihood approaches
in which coefficient estimation and variable selection are done simultaneously.
The most well-known method of this type is the LASSO <cit.>
which minimizes the residuals sum of squares subject to an ℓ_1 constraint.
This constraint allows the number of non-zero parameter
estimates to be controlled and adapt to sparsity.
Other related methods that are based on a minimizing
a loss function, subject to a constraint on the complexity of the
model, include SCAD <cit.>, the adaptive LASSO <cit.>,
LARS <cit.>, and more recent proposals by
<cit.>, <cit.> and <cit.>.

Bayesian approaches are another important direction in model selection
research. Significant contributions include <cit.>,
<cit.>, and the spike-and-slab method in
<cit.>.  The model proposed here is similar to
<cit.>, and <cit.> whose work is motivated by QTL
and genome-wide association studies (GWAS).  Our model allows for a
fully-Bayesian implementation, but an empirical Bayes analysis via the
Generalized Alternating Maximization algorithm <cit.> is
proposed instead because the running time of an MCMC sampler is too
long for many data sets in modern applications.
For example, in our simulations (Section <ref>) we find that the faster of
two MCMC-based variable selection methods required  40 minutes to
complete 1000 MCMC iterations. In contrast, non-MCMC methods required
only a few seconds to perform variable selection, and generally, gave better results.
In this sense the
algorithm is a close competitor to the exact EM algorithm of
<cit.> based on a Bayesian spike-and-slab model. A key
advantage of the approach proposed in this paper is that it extends in a straightforward
manner to the generalized linear model framework.

Recently developed continuous prior distributions have proven more effective in sparse regression than the a Bayesian lasso <cit.> that uses a Laplace prior, which fails to simultaneously induce sparsity while efficiently recovering non-null parameters <cit.>. Among these are “global-local priors” including the Horseshoe-type priors <cit.> and the Gamma Gamma prior <cit.>. These approaches offer computational advantages relative to two-component spike-and-slab mixture priors since one does not need to explore a complex discrete model space of size 2^K. A limitation of the continuous shrinkage prior approach is that it fails to provide a sparse solution. To address the problem, ad hoc post-processing methods for producing sparse estimates from posterior samples that decoupled shrinkage and selection   have been developed <cit.>.

The remainder of the article is organized as follows.  The model and
notation are introduced in Section <ref>.  A Generalized Alternating Maximization
fitting algorithm and selection procedure is described
in Section <ref>. Section <ref> describes results
of simulation studies in which the proposed procedure is compared with
several variable selection packages, including
 <cit.>
and  <cit.>, both of which implement three types of penalties:
LASSO <cit.>, SCAD <cit.>, and MCP
<cit.>. Other packages we included in our simulations are:  <cit.>,
 <cit.>,   <cit.>,
 <cit.>,  <cit.>, and
 <cit.>. We also compared the performance of our method
with a one-predictor-at-a-time approach, controlling the false discovery rate <cit.>.
In Section <ref> we discuss some applications. First,
we demonstrate our method when the response is normal (the logarithm of vitamin B12
production rate). In this case the sample size is N=71 and the number of predictors is 4,088.
We compare the model obtained from our method with ones obtained from TREX <cit.>,
and the results obtained by <cit.>. The second example demonstrates an
application to a binary response (categorized Body Mass Index), which we compare to an analysis
using the continuous outcome.
In this case, N=96 and there are
45 (compositional) predictors. Finally, we explain how to perform variable selection
in survival analysis by treating the number of deaths in a sequence of
non-overlapping time intervals as a Poisson counts.
The paper concludes with with some discussion in Section <ref>.
Important implementation considerations and additional examples are
discussed in the supplementary material.



§ A STATISTICAL MODEL FOR AUTOMATIC VARIABLE
SELECTION

Consider  responses y_i, i=1,…,N, and assume that the
mean of the ith response, λ_i, is linked to a linear predictor, η_i, as follows:

    g(λ_i)≡η_i=∑_j=1^Jx_ijβ _j+
      ∑_k=1^Kz_ikγ _ku_k .

The model (<ref>) allows for J≥ 0 predictors, x_ij, that are
always included in the model and a set of K>0 `putative' predictors, z_ik,
from which it is expected only a small subset are to be
included. Here, β_j is the coefficient associated with the jth
'locked-in' predictor, u_k is a random coefficient associated with
the kth putative predictor. We assume that each γ _ku_k
belongs to one of three components,  C_L, C_0, and C_R,
such that γ _k=0 for k∈ C_0,  γ _k=-1 for k∈ C_L,
and γ _k=1 for k∈ C_R.
We also assume that u_kiid∼ N( μ ,σ ^2),
independently of γ _k.
Justification for both of these choices is given in
Section <ref>.
Thus, the inclusion of the kth predictor in the linear model is determined by
the value of γ _kiid∼ multinomial(-1,0,1;
  p_L,p_0,p_R).
Specification of the distribution of the responses is completed by
assuming that, conditional on the linear predictors, they are
independent draws from a particular exponential dispersion family,
with the most important special cases being the normal, binomial and
Poisson distributions.
Within this modeling framework the problem of
variable selection is cast as a classification problem, for which the main
interest lies in identifying which putative variables belong to C_L∪ C_R,
i.e., which latent variables, γ_k, are non-zero.

Let 𝐗 denote the N× J matrix with
jth column, 𝐱_j, containing the values of the jth
`locked in' predictor. Similarly, define 𝐙 to be the
N× K matrix with kth column 𝐳_k, the
corresponding vector for the kth putative predictor.  Then the
linear predictor model (<ref>) can be rewritten in matrix form
as

    η=𝐗β+𝐙Γ𝐮

where β =(β_1,...,β_J)', Γ≡ diag( γ _1,γ
_2,…,γ _K). Furthermore, the distributional assumptions
concerning {u_k} imply

    𝐙Γ𝐮 | Γ   ∼ N(
    𝐙Γμ,σ ^2𝐙Γ^2𝐙')

where μ=1_Kμ.

Mixture models provide simplicity and tractability and are very
popular in many applications. However, they are known to have
undesirable mathematical properties, such as unbounded likelihood and
lack of identifiability <cit.>. For example, in the
Gaussian case, if μ=0 then the product γ_ku_k in
model (<ref>) has marginal density of the form
h(x)=p_0· 0+(p_L+p_R)φ(x;0,σ), where φ(·;μ,σ)
denotes a normal density. Clearly, p_L and p_R cannot
be separately identified in such a model. However, when μ≠0,
the density is
h(x)=p_0· 0 +
p_Lφ(x;-μ,σ)+p_Rφ(x;μ,σ),
which is identifiable if p_L≠ p_R. Moreover, even if
p_L=p_R, the identifiability issue only concerns the sign of
μ which, in this case, does not affect the marginal distribution h.

Our hierarchical mixture model is similar to the well-known Spike and Slab
model of <cit.> which is implemented, for example, in the R
package BoomSpikeSlab <cit.> and in EMVS <cit.>.
The main difference is our choice of a three-way mixture model, in which
there are two non-null components, rather than one.
Compared with the spike and slab approach, our model  offers a couple of advantages.
To illustrate these advantages, it is helpful to plot the theoretical distributions
of γ_k u_k under the two models (Fig. <ref>). The two-component mixture
is depicted on the left, and our model on the right. In both cases, p_0=0.8.








The two-component mixture assumes that the non-null distribution
is symmetric, which implies a prior belief that the proportion of variables which are
positively correlated with the response is the same as the proportion
of predictors which are negatively correlated with the response. This may be an unreasonable
assumption, and the three-component mixture model is more flexible in
this regard.

In addition, the non-null component in the two-component mixture has much of its
mass around zero, which is counterintuitive because it is assumed that
variables in the non-null component have a non-zero effect. In
contrast, the three-component model assigns a very small probability
to non-null values near zero. Our mixture model also allows for the non-null components to
be highly concentrated, which may be especially useful in situation where there is
a single significant predictor.
The assumption of symmetry of the two nonnull components in our model may be relaxed, but it has
some important benefits. First, it implies that the model is invariant to the sign of columns
of 𝐙. Second, it allows us to borrow information across the two nonnull
components, which greatly contributes to computational stability, as well as power to detect true predictors,
and maintaining a low false positive rate. This is especially beneficial when at least one of the
nonnull components  consists of a small number of  predictors.

The three component prior can be viewed as similar to a non-local prior for which density functions are identically zero whenever a model parameter is equal to its null value.  Conversely, spike and slab priors are local priors where component densities are positive at null parameter values.  <cit.> demonstrate that model selection procedures based on non-local prior densities assign a posterior probability of one to the true model as the sample size n increases and certain regularity conditions on the design
matrix pertain.  Furthermore under the same conditions, they show that standard approaches based on local prior specifications result in the asymptotic assignment of a posterior probability of zero to the true model.














§ ESTIMATION AND VARIABLE SELECTION



 §.§ The Complete Data Likelihood

The marginal distribution of the responses according to the model
described in the previous section is determined by the
parameter vector θ={β, μ, σ^2, p, ϕ},
where ϕ is a dispersion parameter which, depending on the
GLM specification, may or may not be known.  Define
n_j=∑_k=1^KI(γ_k=j), for j=-1,0,1 corresponding to the left,
middle, and right components of the mixture model, respectively, and let
φ(·;μ,Σ) denote a multivariate normal
density with mean μ and covariance matrix Σ. Then assuming
a canonical link function, the complete data log likelihood

is

    ℓ(𝐲,  γ| θ)=log f_C(𝐲, γ | θ)
         =   log∫exp{∑_i=1^Nm_i/ϕ[η_iy_i-b(η_i)]
          }φ(𝐮;Γμ,σ^2𝐙Γ^2𝐙')d𝐮
           +∑_i=1^Nc(y_i,ϕ/m_i) +∑_j=-1,0,1n_jlog p_j ,

where the m_i's are known positive weights, b is the cumulant generator for the GLM satisfying
b'(η_i)=λ_i, and b”(η_i)= V(λ_i), where V is the GLM variance function.

In the Gaussian case the integral in (<ref>) has a closed form
and the complete data loglikelihood reduces to

    ℓ(𝐲,  γ| θ)
        =    -1/2
      (𝐲-𝐗β-𝐙Γμ)'
      (ϕ𝐖^-1+σ^2𝐙Γ^2𝐙')^-1
      (𝐲-𝐗β-𝐙Γμ)
                                           
            -1/2log|ϕ𝐖^-1+σ^2𝐙Γ^2𝐙'| +
       ∑_j=-1,0,1n_jlog p_j-Nlog(2π)/2,

where 𝐖=(w_i) and w_i≡ m_i.
For non-Gaussian GLMs an approximate complete data log-likelihood is
obtained by substituting the so-called 'working response' and
'iterative weight matrix' in place of 𝐲 and 𝐖 in (<ref>), with
components given by

    ỹ_i=g(λ̃_i)+g'(λ̃_i)(y_i-λ̃_i)
      w̃_i=m_i/g'(λ̃_i) .

The substitutions can be justified on
the basis of a Laplace approximation to the integral in
(<ref>), and are the basis of numerous algorithms
in the literature for fitting GLMs with random effects. See, for
example, <cit.>.



 §.§ A Generalized Alternating Maximization Algorithm

In principle, 
the ML estimate of the model parameter vector,
θ, can be obtained using the EM algorithm
<cit.> using the complete data log-likelihood given in
(<ref>),
with the Q-function given by
Q(θ;θ') = E_θ'{log
    f_c(𝐲,γ|θ) | 𝐲} where θ'
 denotes the current estimate of θ.
However, in this case the expectation is intractable, so we propose using the simple plug-in approximation

    E_θ'{log
        f_c(𝐲,γ|θ) | 𝐲}≈ E_θ'{log
        f_c(𝐲,γ'|θ) | 𝐲} ,

where γ' is obtained using one of the methods described later in
this section.

Since the expectation terms are approximated, the resulting iterative estimation procedure
will not fall within the EM framework, or even within the GEM framework <cit.>, and thus,
convergence results from these frameworks will not apply.
This type of EM variant, where the E-step is replaced with an approximation is an example of
the Generalized Alternating Minimization (GAM) framework of <cit.>, of which
EM and GEM are special cases. We will use the GAM theory to show the convergence of our algorithm,
which we describe in Algorithm <ref>.



This algorithm uses a 'likelihood-ratio' approach because a latent variable γ_k is changed from its
current value if that change increases the likelihood in a meaningful way
(holding all other latent indicators at their current values).
The set S consists of
all the variables which yield an improvement greater than  δ in the loglikelihood when their
current classification according to the three-component mixture model is
changed, while holding all other γ_ks at their current values.
If S is not empty, choosing k from this set in Line <ref>
is done according to one of the following methods:

 
  * Greedy: choose k∈ S for which d_j,k is largest.
 
  * Weighted probability: choose k∈ S with probability
 
    d_j,k/∑_(r,s)∈ Sd_r,s .




For the model in (<ref>) Algorithm <ref> converges to a stationary point in a finite number of steps.



A proof of the proposition is given in Appendix <ref>.



 §.§ Notes and Further Details

First, recall that in the binomial and Poisson models we work with the `working response and weights',
so ỹ_i and w̃_i have to be iteratively updated at each iteration of Algorithm <ref>.

Second, we emphasize that our algorithm modifies at most one coordinate in each iteration because,
changing more than one variable
may introduce multicollinearity. What this means is that even though each single variable
may increase the likelihood, changing a set of variables which are highly correlated
may cause the log-likelihood to decrease (and it may even approach -∞, because ℓ
involves the logarithm of the precision matrix.)
As a consequence, our algorithm automatically prevents selecting models with highly correlated
predictors. We do, however, keep track of predictors that are correlated with ones selected to be
in the model, since they are likely to be related to the outcome as well. We return to this point in the
Case Studies section.







Third, we may choose γ' in Line <ref> to be the variables selected by any other
method (FDR, SIS, EMVS, etc.)
Then, the log-likelihood of the final model selected by our algorithm will be greater than or equal to
the one obtained at the 0^th iteration.

Fourth, as a referee pointed out, the parameter space for (θ, γ) cannot assumed to be unimodal,
especially when some putative variables are correlated with each other or with columns in  𝐗. Therefore, it
is recommended to run Algorithm <ref> multiple times, using the
weighted probability approach for selecting the next putative variable to be updated.
This approach is feasible because our method is computationally efficient. When the posterior distribution
is multimodal there is no one `correct' model, and running the algorithm multiple times will allow users to obtain
different, but possibly equally relevant sets of significant predictors each time they fit the model.
This approach is demonstrated in Section <ref> (the riboflavin data example).



Finally, to obtain the maximum likelihood estimates θ' in Line <ref>,
the following formulas are used.
Let Σ=ϕ𝐖^-1+σ^2𝐙Γ^2𝐙'
and 𝐇=[𝐗,𝐙Γ1], the update
formula for β̃=(β',μ)' is given by

    β̃=(𝐇'Σ^-1𝐇)^-1𝐇' Σ^-1𝐲

so μ is simply the last element in β̃,
and the updates for the variance components are given by

    ϕ=τ_e/Nσ^2=τ_r/(𝐙Γ)

provided (𝐙Γ)>0 and μ=σ^2=0 otherwise,
where

    τ_e   =   [ϕ𝐈_N-ϕ^2Σ^-1]+
    ϕ^2(𝐲-𝐇β̃)'Σ^-2
    (𝐲-𝐇β̃)
    τ_r   =   [σ^2𝐈_K-σ^4Γ𝐙'
    Σ^-1𝐙Γ]+
    σ^4(𝐲-𝐇β̃)' Σ^-1𝐙Γ^2𝐙'
    Σ^-1(𝐲-𝐇β̃) .

(see Section 8.3.b in ).
Finally, maximizing (<ref>) with respect
p_L,p_0,p_R leads to the updates p_j=n_j/K for j=-1,0,1.

Small values of μ are counterintuitive because they suggest that the mean effects of
selected putative variables are close to zero. Moreover, as noted
earlier, small values of μ can lead to identifiability problems.
To prevent such problems one might consider adding a penalty term on
μ as proposed by <cit.>.
Another possible approach to preventing identifiability problems is to use a different
nonnull prior, such as a mixture of two log-normal distributions,
so that the probability that variables with effect size close to zero will be considered as
nonnull, will be practically 0.
However, as we shall see in the next subsection, one of the strengths of our method is
that it allows us to achieve significant dimension reduction via the Woodbury identity, which requires the
normality assumption. So, were we to use a different nonnull distribution, we would need to add an
extra step to the algorithm in order to normalize the nonnull components (via a Laplace approximation, for example).
In practice, however, we find that identifiability problems
are effectively avoided by choosing appropriate initial values
for our Generalized Alternating Maximization algorithm.
For further details regarding the initialization of the algorithm, see the
Supplementary Material.







































 §.§ Modifications for large N and K

The complete data
log-likelihood (<ref>) contains a large (N× N) matrix which has to be inverted
to compute the iterative approximate ML estimates.
However, using the Woodbury identity <cit.>,

    Σ^-1    =    (ϕ𝐖^-1+σ^2𝐙Γ^2𝐙')^-1
       =   1/ϕ[𝐖^1/2(𝐈_N+σ^2/ϕ𝐖^1/2𝐙Γ^2𝐙'𝐖^1/2)^-1𝐖^1/2] 
       =   1/ϕ[𝐖- σ^2/ϕ𝐖^1/2𝐙Γ'
    (𝐈_K +σ^2/ϕΓ'𝐙'𝐖𝐙Γ)^-1Γ'𝐙'𝐖^1/2].

This simplifies the computations considerably because the (k,l)th element of
Γ'𝐙'𝐖𝐙Γ is
proportional to γ_kγ_l. Specifically, suppose there are L variables for which
γ_k≠ 0. Define
Γ_L to be the L× L reduced matrix in
which rows and columns of Γ corresponding to excluded putative variables
have been eliminated. Similarly, define 𝐙_L by eliminating
the corresponding rows of 𝐙. Then

    Σ^-1=1/ϕ[𝐖-σ^2/ϕ𝐖^1/2𝐙_LΓ'_L
    (𝐈_L+σ^2/ϕΓ'_L𝐙'_L𝐖𝐙_LΓ_L)^-1Γ'_L𝐙'_L𝐖^1/2] .

Thus, inversion of the N× N matrix Σ is reduced to
inverting a much lower dimensional L× L matrix. Similar
simplifications, due to the exclusion of most putative variables, apply
to computation of the determinant term in (<ref>).

Further details about implementation considerations are provided in the Supplementary
Materials. Specifically, we discuss how to deal with correlation and
interactions among the putative variables, and we discuss computational challenges
stemming from the large number of variables, such as how to avoid loading
the entire matrix 𝐙 to memory, and how to parallelize the algorithm.



§ SIMULATIONS

We conducted a simulation study to verify that under the assumed model (<ref>)
the algorithm yields accurate parameter estimates (not shown here), and to compare the performance of our
algorithm in terms of power and accuracy with other methods.
Not surprisingly, as the sample size increases (even when N is still much smaller than
K), the parameter estimates become more accurate, the power to detect the non-null
variables increases, and the Type-I error rate decreases.

In this section we focus on the results of a simulation study that compares the
Type I and Type II errors when the data are not  necessarily generated according to model (<ref>).
We call the program that implements variable selection according to
model (<ref>) SEMMS (Scalable EMpirical Bayes Model Selection), and
 compare its performance with several well-known variable selection approaches:


  *  <cit.>, which implements three types of penalties:
LASSO <cit.>, SCAD <cit.>, and MCP
<cit.>.

  *  <cit.>, to fit GLM models with the lasso or elastic net regularization.

  *  (least angle regression) <cit.>.

  *  Sure Independence Screening <cit.>. This package combines a
screening process which reduces the number of variables to be considered in each
regularization step, while ensuring that only irrelevant predictors are eliminated. This package
allows to choose among three penalty types (LASSO, SCAD, or MCP).

  *  <cit.>, a Bayesian approach to variable selection with
a fast implementation via the EM algorithm.

  *  <cit.>, a Bayesian (spike and slab) variable selection method.

  *  <cit.> a Bayesian (MCMC) method based on <cit.> using moment
and inverse moment Bayes Factors.

  *  <cit.>, a Bayesian (MCMC) implementation of the
spike and slab model.

  * FDR <cit.> a simple one-predictor-at-a-time approach, controlling
the false discovery rate. In our simulations, we controlled the FDR at the 0.05 level.


The comparison with  and  is only done in the normal
response case, since these packages do not have an option to fit a GLM model. Similarly,
 only handles a normal linear regression model, but the package
 <cit.> includes the function  which can also
be used to analyze binomial response and Cox regression models (but not Poisson).

Generally, these software packages were used with their default values. With  and
, the number of selected variables was such that together they explained 90% of the variability or (null) deviance.
The packages  and   simply return the selected variables.
When using  we set the `bigp.smalln' option to TRUE, and
the selected variables were the non-zero generalized elastic net (gnet) coefficients.
With  we varied the spike variance parameter (twenty values
between 10^-10 and 10^-1 equally spaced on the logarithmic scale)
and the type of the prior distribution
was set to `betabinomial'. The `independent' parameter was set to FALSE,
since it yielded much better results than the default (independent=TRUE, which
resulted in zero selected predictors in most cases.)
Hence, rather than running  under the assumption that the regression
coefficients and the error variance are independent, a priori, a conjugate prior was used
<cit.>.

In the case of a normal response we show results from nine different scenarios with
varying dependence structures and number of significant predictors.
In each scenario K=1000 predictors were initially drawn independently from a Unif[-1,1]
distribution, but in scenarios 3, 4, and 5 they were modified in order to
induce correlations. Scenarios 1-5 are made under the assumption of our mixture model,
while 6-9 are not.
The true number of predictors in each scenario is denoted by L. The error terms, ϵ_i,
are generated as i.i.d. N(0,0.1) variates, except for N2 
where ϵ_i∼ N(0,0.25).



  * A single significant predictor (L=1) is related to the response:
Y_i=Z_1i+ϵ_i.

  * The response is the sum of eight i.i.d. predictors (L=8):
 Y_i=Z_1i+…+Z_8i+ϵ_i.

  * Y_i=Z_1i+…+Z_8i+ϵ_i (again, L=8)
but Z_2=Z_1+δ_i2,  Z_3=-2Z_1+δ_i3, Z_4=-Z_1+δ_i4, and
Z_6=-Z_5+δ_i6, where δ_ik∼ N(0,0.2), independently.
The correlation between Z_1 and each of Z_2, Z_3, Z_4 is approximately 0.95,
as is the correlation between Z_5 and Z_6.

  * Y_i=Z_1i+…+Z_14i+ϵ_i (thus, L=14)
and Z_2,…,Z_10 are drawn from a multivariate normal distribution with mean 0 and a
covariance matrix with a compound symmetry structure, 0.01· I_9+ 0.05· J_9
where I_9 is a 9× 9 identity matrix and J_9 is a 9× 9 matrix of 1's.
The average correlation between any pair from the set Z_2 …,Z_10 is 0.93.

  * Y_i=Z_1i+…+Z_20i+ϵ_i (thus, L=20)
where  Z_1,…, Z_20 are drawn from a multivariate normal distribution
with an autoregressive (AR1) structure, with ρ=0.95.

  * Y_i=∑_j=1^15β_jZ_ji+ϵ_i (thus, L=15) where β_j∼ N(0,1), i.i.d.
In this case, the coefficients of the non-null predictors are drawn from the `slab' component in
the spike-and-slab model.

  * Y_i=∑_j=1^15β_jZ_ji+ϵ_i (thus, L=15) where
β_j=5, 1, 2, 4, 9, 3, 4, 1, 3, 2, 4, 2, 3, 1, 7. In this case, the effect sizes
have dramatically different magnitudes.

  * Y_i=-∑_j=1^4β_jZ_ji+∑_j=4^10β_jZ_ji+ϵ_i (thus, L=10) where β_j=5, 7, 2, 4, 9, 3, 4, 1, 3, 2.
Again, the magnitudes of effect sizes are very different, but this time
four of the effects are negative, and six are positive.

  * Y_i=-∑_j=1^4β_jZ_ji+∑_j=4^10β_jZ_ji+ϵ_i (thus, L=10) where β_j=2, 2, 2, 2, 2, 2, 6, 6, 6, 6.
Similar setting as the previous scenario in the sense that the magnitudes of effect sizes
are very different effects have positive and negative signs, but the effect sizes are from
two non-symmetric point masses (unlike our model, which assumes two normal distributions with
symmetric means.)

In our simulations we used different sample sizes and different number of predictors. Results
shown here are for  K=1,000, and N=100  (Tables <ref> and <ref>) and N=50  (Table <ref>).
The median true positives and false positives were calculated from 30 replications of
each scenario. We observe several things:


  * For both values of N SEMMS achieves the best or nearly the best results. When N=100 it finds all the
true predictors in simulations N1-N5 and N7-N8, and no false positive ones in any scenario.
When N=50, simulations N2-N4 yield a small number of
false positives, but the results are still very good.

  * In simulation N6 (Table <ref>) SEMMS finds 7 of the 15 true predictors, but keep in mind that
under the
spike and slab model which was used to generate the data, a large proportion of the slab component
overlaps with the spike component, and it is expected that under this model some true predictors will
not be detectable. Note that   which uses the spike and slab model for finding significant
predictors yields the same result. In this scenario  with a SCAD or MCP penalty achieves
a slightly better result.

  * When the effect sizes vary dramatically and do not follow our mixture prior assumption,
but are not concentrated at non-symmetric mass points (N6-N8), then SEMMS performs very well.
When the effect sizes are  concentrated at non-symmetric mass points (N9) SEMMS tends to find
the largest effect. This is perhaps an unlikely scenario, since we can expect the effect sizes to
vary, in which case SEMMS performs very well. In any case, because SEMMS maintains a very low
false positive rate, it is possible to deal with
cases such  as N9 by running SEMMS sequentially, each time moving the detected covariates to
the set of `locked-in' variables until no additional effects are found.

  * The spike and slab methods seem very sensitive to deviations from their assumed
model. In simulations N7-N9 both  and  do not detect
all the true effects, and have a large number of false positives.

  * Perhaps a bit unexpected, in the N=100 case the one at a time approach (FDR)
performs quite well, especially in scenarios N1, N3, N4, and N5 where it is the second
or third best method in terms of
its overall error. However, FDR is quite conservative in simulations N6-N9, and
when N is small relative to K the power of this method is much lower
 in simulations N2-N4.

  * Between the  and  packages which offer similar regularization options
the former appears to be better at maintaining lower FP rate.

  * Between the two packages that offer a fast Bayesian (spike and slab) fitting,  achieves
a lower FP rate in simulations N1-N4 and N6, but higher in simulation N5.

  * Generally, both  and  have high FP rates.
















We obtained similar results in the binary response case.
In Table <ref> we show results from two scenarios, both with K=1000 and N=120.
We only use the methods which perform variable selection in the GLM framework,
namely, , ,  , and the one-at-a-time approach,
controlling the FDR. The two simulation scenarios are:


  * η_i=2Z_3i+2Z_6i+2Z_7i
but Z_2=Z_1+δ_i2,  Z_3=-2Z_1+δ_i3, Z_4=-Z_1+δ_i4, and
Z_6=-Z_5+δ_i6, where δ_ik∼ N(0,0.2), independently.
The correlation between Z_1 and each of Z_2, Z_3, Z_4 is approximately 0.95,
as is the correlation between Z_5 and Z_6. Thus, in this case L=7.

  * Autoregressive structure: we set L=10 and  η_i=2Z_1,i+2Z_101,i
but Z_1-Z_5 have an AR(1) structure with ρ=0.95 and so do  Z_101-Z_105.

  * Hub network - we create an N× K matrix so that the columns consist of
g non-overlapping hubs such that within each
hub all K/g nodes follow a multivariate normal distribution with compound symmetry correlation structure.
We pick one node in one of the hubs to be the response, and replace it with a vector of N 0/1 values, based on
a logistic model with the other L=(K/g)-1 nodes in that hub as predictors. In Table <ref>
we show results when L=9 (i.e. 100 hubs, each with 10 nodes), but we get very similar with L=4 and L=19.
For this configuration Table <ref> includes results with two sample sizes: N=120 and N=80.

As was the case with the normal response,
SEMMS appears to
have the best overall performance with the binary responses data. FDR has a slightly higher power in simulation B3 when N=120, but it
detects zero variables when N=80, whereas the performance of SEMMS remains approximately the same
when N is decreased from 120 to 80.


Similarly to the normal and binomial responses, our method has the lowest median false positive rate
and a high power to detect the true effects in the Poisson model. Table <ref> shows the results
of simulation (P1), where η=3+0.3Z_1 + 0.25Z_2 -0.22Z_3 -0.19Z_4 + 0.27Z_5 -0.17Z_6 -0.25Z_7
and N=120.
Simulation (P2) uses the same model, except that predictors Z_1-Z_5 have an AR(1) structure with ρ=0.95.
With P1, SEMMS has a median TP five predictors out of seven, and median of 0 false positives.
In  the P2 setting SEMMS detects all true predictors, while still having no false positives.
The competing methods , ,   all achieve good
results in terms of power in P1, detecting all seven predictors (and fewer in P2), but a higher false positive count, as compared with
SEMMS. Among , , and  the latter appears to yield lower false positive rates, yielding a median of five with all three penalty types in P1. In P2  has a median of 0.5 false positive,
but it detects only three of the seven true predictors.
Note that in the one-at-a-time approach, we used the  option, and not
  since the latter yielded a very high false positive rate.




In our simulation results we did not include the two MCMC-based methods
 ( and ) because they proved to be far too slow to be
 practical, given the existence of
fast alternatives, including Bayesian approaches like  and especially .
For example, when we ran simulations N1 and N5 with  but with a much smaller number of predictors
(K=100) scenario N1 took an average of 4 minutes to complete each replicate, and N5 took an
average of 1.6 hours
(with a sample size N=30 and 1,000 MCMC iterations in both scenarios.)
In the N1 case, the correct model was detected almost every time, but in the case of N5
the median number of true positives was 4 (when L is actually 20), and the median FP was 0.
It appears that  gets much slower as L increases, and less powerful
as compared with competing methods.

The  package appears to be somewhat faster, taking about 40 minutes to
complete 1,000 MCMC iterations when K=1000. However, it does not seem to have an advantage
over the best available methods in terms
of power or error rate. For example, in the case of simulation N5
the median number of predictors found by  was 1 (out of 20), with
0 false positives.







§ CASE STUDIES



 §.§ Normal Response – the Riboflavin Data


In a recent paper demonstrating modern approaches to high-dimensional statistics,
<cit.> analyzed a data set from <cit.> in which the response
variable is the logarithm of riboflavin (vitamin B12) production rate,
and there are normalized expression levels of 4,088 genes which are used
as explanatory variables. The sample size in the data set is N=71.
In addition to the fact that the number of putative variable greatly exceeds
the number of observations, many of the putative variables are highly correlated.
Out of 8,353,828 pairs of genes, there are 70,349 with correlation coefficient
greater than 0.8 (in absolute value).

<cit.> report that the Lasso with B=500 independent random subsamples of size ⌊ N/2 ⌋,
and with q=20 variables that enter the regularization step
first, yields three significant and stable genes: LYSC_at, YOAB_at, and YXLD_at.
The model with these three variables has an R^2 of 0.68, and AIC of 118.6.
Their multisample-split method yields one significant variable (YXLD_at), and the
projection  estimator, used with Ridge-type score yields no significant variables
at the FWER-adjusted 5% significance level.
The model with just YXLD_at has an R^2 of 0.37 and AIC of 162.

<cit.> also used this data set, to demonstrate their TREX model. Their final model
includes three genes: YXLE_at, YOAB_at, and YXLD_at. The R^2 of their
model is 0.62, and the AIC is 130.68. Two of the genes are highly correlated (YXLD_at and
YXLE_at) and yield a variance inflation factor of 23.7 each, so when fitting the final
linear regression model neither appears to be significant.

We ran SEMMS using both the greedy and the weighted probability methods as
described in Section <ref>.
The greedy method yielded six significant genes: LYSC_at, SPOIISA_at , XTRA_at,
YDDK_at, YURQ_at, and YXLD_at. The predicted level of the logarithm of
riboflavin production rate is given by the formula:

    Ŷ   =   -7.16 - 0.23· Z_LYSC_at + 0.19· Z_SPOIISA_at
    + 0.16· Z_XTRA_at
          
    - 0.23· Z_YDDK_at + 0.26· Z_YURQ_at
    -0.35· Z_YXLD_at


This model has AIC=70.9  and an R^2 of 0.85.
Figure <ref>  depicts the observed values
 (Y) vs. the fitted values from three models (<cit.>, <cit.>,
 and greedy SEMMS.)
In addition to having much smaller residuals than the two other methods,
SEMMS provides much better prediction for low values of riboflavin. The other two
methods seem to  over-estimate the riboflavin levels when the true (normalized) values
are small (less than -9).


We ran SEMMS using the weighted probability method 100 times.
The best model included six genes (CARB_at, SPOVAB_at, XHLA_at,
YCKE_at, YOAB_at, YXLD_at) and had an AIC of  53.9 and an
R^2 of 0.88.
In the 100 runs of the randomized SEMMS a total of 16 genes were selected, yielding
an AIC of 46.8 and an R^2 of 0.92.

We also used packages that performed well in our simulations.  found
eight predictors (AIC=68.2) and  with MCP penalty found four predictors
(AIC=86.3).

The notion of `the best selected model' may not always be appropriate, since
(i) the number of putative variables is large there is no way to evaluate all possible
models, and (ii) some selected predictors can be part of a network of highly correlated
variables.
To deal with (ii) SEMMS detects predictors that are highly correlated with ones
selected to be in the model, and reports them as well. This is illustrated in
Figure <ref>.
The variables selected by (greedy) SEMMS are shown as red diamonds and their
 coefficients in the fitted linear model as numbers next to the (dark blue) edges.
Variables that are highly correlated with a selected predictor are depicted as orange dots and
they are connected via a gray edge to another predictor if their correlation with that predictor
exceeds a user-defined threshold (in this case, we used 0.75).
With the riboflavin dataset it can be seen that the YURQ_at gene is co-expressed with a
large group of genes, and hence each one of those genes (or perhaps some
weighted average of their expression levels) could be considered as a relevant predictor
for the response. In this case study, with a correlation threshold of 0.75 between predictors
greedy SEMMS reports a total of 62 relevant variables.





 §.§ Binary Response – the BMI Data

<cit.> demonstrated an application of a LASSO-based variable selection
method for regression models with compositional covariates. The analysis aims
to identify a subset of 87 bacteria genera in the gut whose subcomposition
is associated with body-mass index (BMI). The data, which was introduced in
<cit.>, is compositional,
which using our previous notation means that ∑_j=1^87 z_ij=1 for each i.
The total number of samples is N=96.
To apply our method directly, without changing the model to account for the sum
constraint, we perform the log ratio transformation and replace the matrix 𝐙
with 𝐙^K=[log(z_ij/z_iK)]. The data contains many zero counts, so
<cit.> replace them with 0.5 before converting the data to be in
compositional form. Note that there may be preferable imputation methods, but we choose
to use the same method in <cit.> in order to have a meaningful comparison.
We use a subset of
45 bacteria which had non-zero counts in at least 10% of the samples (N=96).  The
omitted genera have minimal contribution to the overall distribution of
the proportions.

To demonstrate the application of SEMMS to the binary response case we create
an categorical BMI variable with 4 levels: underweight - less than 18.5 (n=5),
normal [18.5, 25) (n=30), overweight [25.5, 30) (n=25), and obese  ≥30 (n=10).
We used the binomial option in SEMMS
in order to find  which bacteria are associated with p_b=Pr(obese) and
which are associated with p_w=Pr(overweight), with the normal group as the baseline.
The underweight group is too small and is not included in the analysis.
We found six bacteria associated with obesity
(Acidaminococcus, Alistipes, Allisonella, Butyricimonas, Clostridium,
and Oxalobacter), and five with overweight
(Anaerofilum, Faecalibacterium, Oscillibacter, Turicibacter, and Veillonella).
We also ran SEMMS with BMI as a continuous response, using the normal model and
obtained the same four genera reported by <cit.>
(Acidaminococcus, Alistipes,  Allisonella, and Clostridium).
The fact that these four are a subset of the genera found to be associated with
p_b but have no overlap with the ones associated with p_w
suggest that different BMI levels are  associated with different
bacteria, and thus, a categorical analysis or perhaps quantile regression
may be more appropriate than the conditional mean models in regression analysis.




















 §.§ Survival Analysis – the NKI70 data

Our model can be used to deal with censored survival times.
We follow <cit.> who proposed using an
artificial Poisson model to fit Cox's proportional hazards (PH) regression model.
Recall that in the proportional hazard model, the hazard functions have the form
λ(t;𝐳)=λ_0(t)exp(𝐳β).
Suppose there are q deaths occurring after survival times t_1,…,t_q
(one death after each t_i.)
For each h=1,…,q let Y_h,j be independent Poisson random variables with
parameters μ_h,j, such that 
    μ_h,j=N_h,jexp(α_h+𝐳_h,jβ)

where N_h,j is the number of patients who are at risk after t_h
and whose explanatory vector is 𝐳_h,j(t_h).
Since the excess life of a patient beyond time t_h-δ t_h follows an exponential
distribution, the number of deaths in each group follows a Poisson distribution
with parameter μ_h,j.
In other words, rather than modeling the survival times,
we model the number of (instantaneous) deaths among the survivors at a given time
in a given group.
The mean number of instantaneous deaths among survivors in a group of
individuals sharing the same properties is associated with a linear combination of
the explanatory variables through the log link function.
This strategy allows for ties and censoring in the data.
Another approach to survival modeling is the  accelerated failure time (AFT) models.
More biostatisticians are using the AFT model because it is based on the linear model
and the estimated regression coefficients have a rather direct physical interpretation.
The AFT model can be formulated as variants of estimating equations that
correspond to generalized linear models via a Buckley-James transformed Gaussian estimating equations
<cit.>.

To demonstrate this approach, we used the NKI70 dataset from <cit.>, which is
available in the `' package <cit.>.
This dataset contains gene expression measurements for 70 genes,
obtained from lymph nodes of 144 breast cancer patients.
The 70 genes were determined as prognostic for metastasis-free survival
in earlier studies.  Large-scale microarray gene expression analyses have revealed
a signature set of genes that can predict breast-cancer prognosis
(<cit.>, <cit.>, <cit.>, <cit.>, <cit.>, <cit.>).
However, few genes overlapped in these assays, and only a few of
the breast-associated genes have been validated at the protein
level.

We analyzed which genes were associated with either death or recurrence
of metastasis in less than three months. We included the (scaled) age of
the patients as an explanatory variable, as well as the logarithm of the
time until death/recurrence or censoring where, for the purpose of
this analysis we considered any
individual who survived more than 3 months as censored.  SEMMS generated the following model


    η=logμ=-2.8-0.37log(t)-0.35· Z_SCUBE2 - 0.79· Z_KNTC2
    
     -0.56· Z_ZNF533+0.36· Z_IGFBP5.1+1.09· Z_PRC1

The residual deviance for this model is 39 (vs. 87.55 for an intercept-only model).
These selected genes have been  determined to be critical  as prognostic for metastasis-free survival
in an earlier studies.
In particular, a cross-platform comparison of breast-cancer
gene sets from these profiling studies revealed SCUBE2 as a common gene that has been validated at the protein
level  <cit.>.


















§ CONCLUSION

In 1996 Brad Efron stated that variable selection in regression is the most important
problem in statistics <cit.>. Since then many papers have been written
on the topic, as this continues to be a challenging problem in the age of high-throughput
sequencing in genomics, and as other types of `omics' data become available and more
affordable.

We have developed a model-based, empirical Bayes approach to variable
selection. We define a mixture model in which the putative
variables are modeled as random effects, and we demonstrate that this approach
results in high power to select the correct variables
while maintaining a low rate of false positive selections, in a variety
of situations.
Our algorithm is scalable and computationally efficient because of: (i) the
parsimony of the mixture model, since
the mixture model involves a very small number of parameters, and that number
remains constant regardless of the number
of putative variables; and (ii) the usage of the Generalized Alternating
Maximization algorithm to estimate the model parameters, as well as
an efficient dimension reduction trick via the Woodbury identity.
The Generalized Alternating Maximization algorithm not only converges
significantly faster than simulation-based methods, but also uses
memory more efficiently, since it has to keep only the posterior non/nonnull
probabilities for each predictor from only the two latest iterations,
as opposed to entire chains.
Finally, a simple modification to the algorithm prevents multicollinearity problems in
the fitted regression model.





SUPPLEMENTARY MATERIAL






  Title: The file SEMMS_a.pdf contains some implementation notes, instructions how
to obtain and install the package, and three examples with code and plots generated by
SEMMS. (pdf format)


  R-package: version 0.1.0 of  SEMMS is available online at <https://haim-bar.uconn.edu/software/>.
The GNU zipped tar file contains datasets used in the case studies, as well as one
simulated dataset and a file with ozone levels data for an additional case study of a well-studied
dataset.





§ APPENDIX

Proof of Proposition <ref>:
We denote the estimates obtained in step t of Algorithm <ref> by a superscript.
First, we show that the algorithm terminates after a finite number of steps.
Given γ^t, the estimate for θ from Line <ref> is unique, because the log-likelihood
is a convex function. Therefore, if at any t, S^t+1 is empty (that is, the loglikelihood will not improve
by changing γ) then
the algorithm terminates because the estimate for θ will not change in
the subsequent iteration.
So, we need to show that for all t≥ 0 and s>1 we cannot have
γ^t=γ^t+s (i.e., there are no loops in the sequence of γ^t).
Suppose that for some t there exists s>1 such that γ^t=γ^t+s. Then,
the convexity of the loglikelihood function implies that θ^t=θ^t+s. Together, this means that
ℓ^t(𝐲) = ℓ^t+s(𝐲), but in the algorithm we selected
the coordinates to modify in each step, k_t,…,k_t+s, so that in each
iteration the log-likelihood strictly increased, so we must have
ℓ^t+s(𝐲) - ℓ^t(𝐲) > 0.
Proposition 7 in <cit.> shows that minimization (and
consequently maximization) of a continuous function forms a closed
point-to-set map.  Their Lemma 9 shows the Cartesian product of
two closed point-to-set-maps is itself closed.  With these results in
hand one can apply Zangwill's generalized convergence theorem <cit.> to
the sequence (θ^t,  γ^t).  


chicago

