EBglmnet Vignette

Acronyms and Concepts

A good prior distribution is a balance of the aforementioned properties as shown in the Figure 1. Generally, NE prior leads to more non-zero predictors with smaller absolute coefficient values, while NEG prior leads to less non-zero predictors with stronger signals in terms of both absolute coefficient values as well as significance level.

Prior distributions that penalise posterior distributions

The following names should not be confused with the lasso and elastic net method from the comparison package glmnet:

EBglmnet: package that implements EBlasso and EBEN methods.

EBlasso: empirical Bayesian method with lasso prior distribution, which includes two sets of prior distributions: NE and NEG.

EBEN: empirical Bayesian method with elastic net prior distribution.

lasso prior: the hierarchical prior distribution that is equivalent with lasso penalty term when the marginal probability distribution for the regression coefficients is considered.

elastic net prior: the hierarchical prior distribution that is equivalent with elastic net penalty term when the marginal probability distribution for the regression coefficients is considered.

EBlasso-NE: EBlasso method with NE prior.

EBlasso-NEG: EBlasso method with NEG prior.

Generalized Linear Models (GLMs)

In a GLM

η = μI + Xβ, where $\bf X$ is an $\it{n}\times\it{p}$ matrix containing $\it{p}$ variables for $\it{n}$ samples ($\it{p}$ can be $\gg \it{n}$). η is an $\it{n}\times 1$ linear predictor and is related to the response variable y through a link function $\it{g}$: E($\mathit{\mathbf{y}} |\bf X$)=$\it{g}^{-1}$(μI + Xβ), and β is a $\it{p}\times 1$ vector of regression coefficients. Depending on certain assumption of the data distribution on y, the GLM is generally inferred through finding the set of model parameters that maximize the model likelihood function p(y|μ, β, φ), where φ denotes the other model parameters of the data distribution. However, such Maximum Likelihood (ML) approach is no longer applicable when $\it{p}\gg \it{n}$. With Bayesian lasso and Bayesian elastic net (EN) prior distribution on β, EBglmnet solves the problem by inferring a sparse posterior distribution for $\hat{\boldsymbol{\beta}}$, which includes exactly zero regression coefficients for irrelevant variables and both posterior mean and variance for non-zero ones. Comparing with the glmnet package, not only does EBglmnet provide features including both sparse outcome and hypothesis testing, simulation study and real data analysis in the reference papers also demonstrates the better performance in terms of Power of Detection (PD), False Discovery Rate (FDR), as well as Power Detecting Group Effects when applicable. While mathematical details of the EBlasso and EBEN methods can be found in the reference papers, the principle of the methods and differences on the prior distributions will be briefly introduced here.

Lasso and its Bayesian Interpretation

Lasso applies a penalty term on the log likelihood function and solves for $\hat{\boldsymbol{\beta}}$ by maximizing the following penalized likelihood:

$$ \hat{\boldsymbol{\beta}} = \arg_{\boldsymbol{\beta}}\max\left[\log\mathit{p}(\mathit{\boldsymbol{y}}|\mu, \boldsymbol{\beta}, \varphi) -\lambda||\boldsymbol{\beta}||_1\right]. $$

The $\it{l_1}$ penalty term can be regarded as a mixture of hierarchical prior distribution:

$$ \beta_j \sim \mathit{N}(0,\sigma_j^2),\\ \sigma_j^2 \sim \exp(\lambda), j = 1, \dots, p, $$

and maximizing the penalized likelihood function is equivalent to maximizing the marginal posterior distribution of β :

$$ \hat{\boldsymbol{\beta}} = \arg_{\boldsymbol{\beta}}\max \log \mathit{p}(\boldsymbol{\beta}|\mathit{\boldsymbol{y}},\mathbf{X},\mu,\lambda, \varphi)\\ \approx\arg_{\boldsymbol{\beta}}\max \log\int\left[\mathit{p}(\mathit{\boldsymbol{y}}|\mu, \boldsymbol{\beta}, \varphi)\cdot (2\pi)^{-p/2}\lvert\mathbf{A}\rvert^{1/2}\exp\{-\frac{1}{2}\boldsymbol{\beta}^T\mathbf{A}\boldsymbol{\beta}\}\cdot \prod^p_{j=1}\lambda\exp\{-\lambda\sigma_j^2\}\right]d\boldsymbol{\sigma}^2, $$

where A is a diagonal matrix with σ⁻² on its diagonal. Of note, lasso integrates out the variance information σ² and estimates a posterior mode $\hat{\boldsymbol{\beta}}$. The $\it{l_1}$ penalty ensures that a sparse solution can be achieved. In Bayesian lasso (Park and Casella, 2008), the prior probability distribution is also conditional on the residual variance so that it has a unique mode.

Empirical Bayesian Lasso (EBlasso)

EBglmnet keeps the variance information, while still enjoying the sparse property by taking a different and slightly complicated approach as showed below using EBlasso-NE as an example:

In contrary to the marginalization on β, the first step in EBlasso-NE is to obtain a marginal posterior distribution for σ² :

$$ \mathit{p}(\boldsymbol{\sigma}^2|\mathit{\boldsymbol{y}},\mathbf{X},\mu,\lambda, \varphi) = \int c \left[\mathit{p}(\mathit{\boldsymbol{y}}|\mu, \boldsymbol{\beta}, \varphi)\cdot (2\pi)^{-p/2}\lvert\mathbf{A}\rvert^{1/2}\exp\{-\frac{1}{2}\boldsymbol{\beta}^T\mathbf{A}\boldsymbol{\beta}\}\cdot \prod^p_{j=1}\lambda\exp\{-\lambda\sigma_j^2\}\right]d\boldsymbol{\beta}, $$

where c is a constant. While the integral in lasso is achieved through the conjugated Normal + Exponential (NE) prior, the integral in EBlasso-NE is completed through mixture of two normal distributions: $\it{p}$(β|σ²) and p(y|μ, β, φ), and the latter one is typically approximated to a normal distribution through Laplace approximation if itself is not a normal PDF. Then the estimate $\hat{\boldsymbol{\sigma}}^2$ can be obtained by maximizing this marginal posterior distribution, which has the following form:

$$ \hat{\boldsymbol{\sigma}}^2 = \arg_{\boldsymbol{\sigma}^2}\max \log \mathit{p}(\boldsymbol{\sigma}^2|\mathit{\boldsymbol{y}},\mathbf{X},\mu,\lambda, \varphi)\\ = \arg_{\boldsymbol{\sigma}^2}\max \left[ \log\mathit{p}(\mathit{\boldsymbol{y}}|\mu, \boldsymbol{\sigma}^2, \varphi,\lambda)-\sum^p_{j=1}\lambda\sigma_j^2 + c\right]. $$

Given the constraint that σ² > 0, the above equation essentially maximizes the $\it{l_1}$ penalized marginal likelihood function of σ², which images the $\it{l_1}$ penalty in lasso with the beauty of producing a sparse solution for $\hat{\boldsymbol{\sigma}}^2$. Note that if σ̂_j² = 0, β̂_j will also be zero and variable x_j will be excluded from the model. Finally, with the sparse estimate of $\hat{\boldsymbol{\sigma}}^2$, the posterior estimate of $\hat{\boldsymbol{\beta}}$ and other nuance parameters can then be obtained accordingly.

Hierarchical Prior Distributions in EBglmnet

EBglmnet Implementation and Usage

The EBglmnet algorithms use greedy coordinate descent, which successively optimizes the objective function over each parameter with others fixed, and cycles repeatedly until convergence. Key algorithms are implemented in C/C++ with matrix computation using the BLAS/LAPACK packages. Due to closed form solutions for $\hat{\boldsymbol{\sigma}}^2$ in all prior setups and other algorithmic and programming techniques, the algorithms can compute the solutions very fast.

We recommend using EBlasso-NEG when there are a large number of candidate effects (eg.,  ≥ 10⁶ number of effects such as whole-genome epistasis analysis and GWAS), and using EBEN when groups of highly correlated variables are of interest.

The authors of EBglmnet are Anhui Huang and Dianting Liu. This vignette describes the principle and usage of EBglmnet in R. Users are referred to the papers in the Reference section for details of the algorithms.

Gaussian Model

EBglmnet assumes a Gaussian distribution on y by default, i.e., p(y|μ, β, φ) = N(μI + Xβ, σ₀²I), where φ = σ₀² is the residual variance. In the above example, both μ̂ and $\hat{\sigma_0}^2$ are listed in the output:

output$Intercept
output$residual

Binomial Model

If there are two possible outcomes in the response variable, a binomial distribution assumption on y is available in EBglmnet, which has p(y|μ, β, φ) following a binomial distribution and φ ∈ ∅. Same as the widely-used logistic regression model, the link function is $\eta_i = logit(p_i)=\log(\frac{Pr(y_i)=1}{1-Pr(y_i=1)}), i = 1, \dots, n$. To run EBglmnet with binomial models, users need to specify the parameter family as binomial:

yBinary = y>mean(y);
output = EBglmnet(x,yBinary,family="binomial", hyperparameters = c(0.1, 0.1))

For illustration purpose, the above codes created a binary variable yBinary by set the cutoff at the mean Life Expectancy value.

Back to Top

Gaussian Model with Main Effects

We will demontrate the application of EBglmnet in multiple QTL mapping using a simulated F2 population, which is available along with EBglmnet package. The genotype of the F2 population is simulated from cross of two inbred lines. A total of p = 481 gentic markers were simulated to be evenly spaced on a large chromosome of 2400 centi-Morgan (cM) with an interval of d = 5 cM. Theoretically, two adjacent markers have a correlation coefficient R = e^−2d = 0.9048 since the Haldane map function is assumed. The dummy variable for the three genotypes, AA, Aa and aa of individual i at marker j was defined as x_ij = 1, 0, −1, respectively.

data(BASIS)#this is the predictors matrix
N = nrow(BASIS)
p = ncol(BASIS)
j = sample((p-2),1)
cor(BASIS[,c(j,j+1,j+2)]) #Correlation structure among predictors 

##           m183      m184      m185
## m183 1.0000000 0.9178161 0.8342258
## m184 0.9178161 1.0000000 0.9168214
## m185 0.8342258 0.9168214 1.0000000

Let us first simulate a quantitative phenotype with population mean 100 and residual variance that contribute to 10% of population variance. Ten QTLs were simulated from the 481 markers, and effect sizes are randomly generated from [2,3]. We assume that QTLs were coincided with markers. If QTLs were not on markers, they may still be detected given the above correlation structure, although a slightly larger sample size may be needed to give the same PD:

set.seed(1);
Mu = 100; #population mean;
nTrue = 10; # we assume 10 out of the 481 predictors are true causal ones
trueLoc = sort(sample(p,nTrue));
trueEff = runif(nTrue,2,3); #effect size from 2-3
xbeta = BASIS[,trueLoc]%*%trueEff;
s2 =  var(xbeta)*0.1/0.9 #residual variance with 10% noise
residual = rnorm(N,mean=0,sd=sqrt(s2))
y = Mu + xbeta + residual;  

To demonstrate the performance of EBglmnet in analyzing dataset with p > n, we will analyze the simulated dataset using a smaller sample size n = 300:

n = 300;
index = sample(N,n); 
CV = cv.EBglmnet(x=BASIS[index,],y=y[index],family="gaussian",prior= "lassoNEG",nfold= 5)

With 5 fold CV, EBlasso-NEG identified the following effects:

CV$fit
trueLoc

Comparing with the true QTL locations, EBlasso-NEG successfully identified all QTLs. Of note, an identified marker that is <20cM from a true QTL is generally not considered as a false effect, given the small distance and high genotype correlation. Since there is limited prior information in this simulation, other methods in EBglmnet will yield similar results.

Binomial Model and Separation Problem

In many genetics and population analysis such as binary QTL mapping or GWAS, both genetic effects (eg., X takes values of 1, 0, -1 denoting three genotype values AA, Aa, aa) and response variable (eg., y takes values of 0 and 1 denoting the two phenotype classes) are discrete values. Separation problem (complete separation and semi-complete separation) occurs often in such data, especially when epistasis is considered due to the larger number of variables p′ = p(p + 1)/2. Of note, separation is a problem in logistic regression where there exist some coefficients β such that y_i = 1 whenever x_i^Tβ > 0, and y_i = 0 whenever x_i^Tβ < 0, i = 1, …, n. Unless the phenotype is a Mendelian trait, finding a genetic factor/set of factors that perfectly predict the phenotype outcome is too good to be true. While separation problem has been well documented in many statiscal textbook (eg., Ch9 of Altman M, Gill J, Mcdonald M P. (2005)), it is less studied in high dimensional sparse modeling methods. We next examine the problem using lasso as an example.

n = 300;
set.seed(1)
index = sample(nrow(BASIS),n)
p = ncol(BASIS);
m = p*(p+1)/2;
#1. simulate true causal effects
nMain = 10;
nEpis  = 10;
mainLoc = sample(p,nMain);
episLoc = sample(seq((p+1),m,1),nEpis);
trueLoc = sort(c(mainLoc,episLoc)); #a vector in [1,m]
nTrue = length(trueLoc);
trueLocs = ijIndex(trueLoc, p); #two columns denoting the pair (i,j)
#2. obtain true predictors
basis = matrix(0,n,nTrue);
for(i in 1:nTrue)
{
    if(trueLocs[i,1]==trueLocs[i,2])
    {
        basis[,i] = BASIS[index,trueLocs[i,1]]
    }else
    {
        basis[,i] = BASIS[index,trueLocs[i,1]]*BASIS[index,trueLocs[i,2]]
    }
}
#3. simulate true effect size   
trueEff  = runif(nTrue,2,3);
#4. simulate response variables
xbeta = basis%*%trueEff;
vary = var(xbeta);
Pr = 1/(1+ exp( -xbeta));
y = rbinom(n,1,Pr);

CV = cv.EBglmnet(x=BASIS[index,],y=y,family="binomial",prior="lasso",nfold=5)

## Empirical Bayes LASSO Logistic Model (Normal + Exponential prior) 5 fold cross-validation

ind = which(CV$fit[,5]<=0.1)#p-value cutoff
CV$fit[ind,]

##       predictor       beta posterior variance  t-value     p-value
##  [1,]        54  0.3574937         0.04064460 1.773238 0.077207261
##  [2,]        56  0.6570230         0.05495714 2.802648 0.005399570
##  [3,]        72  0.5370808         0.03918611 2.713148 0.007051139
##  [4,]        91  0.5854443         0.04522653 2.752890 0.006268404
##  [5,]       210  0.5296459         0.03667812 2.765554 0.006036039
##  [6,]       225  0.5406408         0.05317435 2.344540 0.019704733
##  [7,]       227  0.4852761         0.04927319 2.186168 0.029578514
##  [8,]       239  0.4313054         0.04712116 1.986904 0.047845428
##  [9,]       240  0.3445439         0.04333749 1.655056 0.098962334
## [10,]       270  0.5462582         0.04604777 2.545621 0.011410310
## [11,]       293 -0.3105339         0.02784121 1.861080 0.063714242
## [12,]       336  0.3962398         0.03394442 2.150671 0.032303559

X = matrix(0,n,m);
X[,1:p] = BASIS[index,];
kk = p + 1;
for(i in 1:(p-1))
{
    for(j in (i+1):p)
    {
        X[,kk] = BASIS[index,i] * BASIS[index,j];
        kk = kk + 1;
    }
}

library(glmnet);

## Loading required package: Matrix

## Loaded glmnet 4.1-8

alpha = 1
lambdaRatio = 1e-4; #same as in EBlasso
cv = cv.glmnet(X, y, alpha = alpha,family="binomial",nfolds = 5,lambda.min.ratio=lambdaRatio)
nLambda = length(cv$lambda)
nLambda

## [1] 81

nbeta = rep(0,nLambda);
fit0 = cv$glmnet.fit;
for(i in 1:nLambda)
{
  nbeta[i] = length(which(fit0$beta[,i]!=0))
}
plot(nbeta,xlab=expression(paste(lambda, " in lasso selection path(n=300,p=115,921)")),
     ylab="No. of nonzero effects",xaxt="n")#
ticks = seq(1,nLambda,10)
axis(side=1, at= ticks,labels=round(cv$lambda[ticks],5), las=1,cex.axis = 0.5)
title("Number of nonzero effects in lasso selection path")

lambda= cv$lambda.min
coefs = fit0$beta
ind = which(cv$lambda==cv$lambda.min)
beta = coefs[,ind]
betaij = which(beta!=0)
Xdata = X[,betaij];
colnames(Xdata) = betaij; 
refit = glm(y ~ Xdata, family = binomial(link = "logit"))#separation occurs

## Warning: glm.fit: algorithm did not converge

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

EBglmnet Algorithms:

Anhui Huang, Shizhong Xu, and Xiaodong Cai. (2015).
Empirical Bayesian elastic net for multiple quantitative trait locus mapping.
Heredity, Vol. 114(1), 107-115.

Anhui Huang, Shizhong Xu, and Xiaodong Cai. (2014a).
Whole-genome quantitative trait locus mapping reveals major role of epistasis on yield of rice.
PLoS ONE, Vol. 9(1) e87330.

Anhui Huang, Eden Martin, Jeffery Vance, and Xiaodong Cai (2014b).
Detecting genetic interactions in pathway-based genome-wide association studies.
Genetic Epidemiology, 38(4), 300-309.

Anhui Huang, Shizhong Xu, and Xiaodong Cai. (2013).
Empirical Bayesian LASSO-logistic regression for multiple binary trait locus mapping.
BMC Genetics, 14(1),5.

Xiaodong Cai, Anhui Huang, and Shizhong Xu (2011).
Fast empirical Bayesian LASSO for multiple quantitative trait locus mapping.
BMC Bioinformatics, 12(1),211.