NormalBetaPrime (version 2.2)

nbp.VB: Variational EM Algorithm for Normal-Beta Prime Regression

Description

This function implements the variational EM approach for the normal-beta prime (NBP) model in the standard linear regression model,

$$y = X \beta + \epsilon,$$

where \(\epsilon \sim N_n (0, \sigma^2 I_n)\). This is achieved by placing the normal-beta prime (NBP) prior of Bai and Ghosh (2019) on the coefficients of \(\beta\). Mean field variational Bayes (MFVB) is used to approximate the posterior \(\pi(\beta | y)\) with an appropriate variational density \(q (\beta)\). The hyperparameters can be set deterministically by the user or they may be automatically selected by marginal maximum likelihood (MML).

It is recommended that the user use the Monte Carlo implementation, nbp, rather than the variational Bayes (VB) function, nbp.VB, for more accurate point esimates and posterior approximations.

Usage

nbp.VB(X, y, method.hyperparameters=c("fixed","mml"), a=0.5, b=0.5, c=1e-5, d=1e-5,  
    selection=c("dss", "intervals"), tol=0.001, n.iter=1000)

Arguments

X

\(n \times p\) design matrix. Should be centered.

y

\(n \times 1\) response vector. Should be centered.

method.hyperparameters

The method for estimating the shape parameters \((a,b)\). If "mml" is chosen, the function estimates the marginal maximum likelihood (MML) estimates of \((a,b)\) by the EM algorithm described in Bai and Ghosh (2019). If "fixed" is chosen, then \((a,b)\) are not estimated from the data.

a

Shape parameter for \(\beta'(a,b)\). The default is 0.5. The user may specify a different value for a (\(a>0\)). This is ignored if the method for estimating hyperparameters is "mml".

b

Shape parameter for \(\beta'(a,b)\). The default is 0.5. The user may specify a different value for b (\(b>0\)). This is ignored if the method for estimating hyperparameters is "mml".

c

The shape parameter for the \(IG(c,d)\) prior on unknown variance parameter, \(\sigma^2\). The default is \(10^{-5}\).

d

The rate parameter for the \(IG(c,d)\) prior on the unknown variance parameter, \(\sigma^2\). The default is \(10^{-5}\).

selection

The method of variable selection. "dss" implements the decoupled selection and shrinkage (DSS) method of Hahn and Carvalho (2015) to select variables. "intervals" performs variable selection by examining the 95 percent posterior credible intervals for each coefficient, \(\beta_i, i = 1, \ldots, p\).

tol

The convergence criterion. If the absolute value of the difference between the current ELBO and the previous ELBO falls below tol, then the variational EM algorithm terminates.

n.iter

The maximum number of coordinate ascent iterations to run. Defaults to 1000.

Value

The function returns a list containing the following components:

beta.hat

The posterior mean estimate of \(\beta\).

beta.var

The posterior variance estimates for each \(\beta_i, i=1, \ldots, p.\)

beta.intervals

The 95 percent credible intervals for all \(p\) estimates in \(\beta\).

nbp.classifications

A \(p\)-dimensional binary vector with "1" if the covariate is selected and "0" if it is deemed irrelevant.

sigma2.esimate

Estimate of unknown variance component \(\sigma^2\).

a.estimate

MML estimate of shape parameter \(a\). If \(a\) was fixed a priori, returns fixed \(a\).

b.estimate

MML estimate of shape parameter \(b\). If \(b\) was fixed a prior, returns fixed \(b\).

Details

The function implements the normal-beta prime (NBP) model of Bai and Ghosh (2019) using mean field variational Bayes (MFVB). The posterior variances and 95 percent credible intervals for each of the \(p\) covariates are also returned so that the user may assess uncertainty quantification. The full model is:

$$Y | (X, \beta ) \sim N_n(X \beta, \sigma^2 I_n),$$ $$\beta_i | \omega_i^2 \sim N(0, \sigma^2 \omega_i^2), i = 1, ..., p,$$ $$\omega_i^2 \sim \beta'(a,b), i = 1, ..., p,$$ $$\sigma^2 \sim IG(c,d),$$

where \(\beta'(a,b)\) denotes the beta prime density,

$$\pi(\omega_i^2) = \frac{\Gamma(a+b)}{\Gamma(a) \Gamma(b)} (\omega_i^2)^{a-1} (1+\omega_i^2)^{-a-b}. $$

References

Bai, R. and Ghosh, M. (2019). "On the beta prime prior for scale parameters in high-dimensional Bayesian Regression models." Pre-print, arXiv:1807.06539.

Blei, D. M., Kucukelbir, A., and McAuliffe, J. D. (2017). "Variational inference: A review for statisticians." Journal of the American Statistical Association, 112(518):859-877.

Hahn, P. R. and Carvalho, C. M. (2015). "Decoupling shrinkage and selection in Bayesian linear models: A posterior summary perspective." Journal of the American Statistical Association, 110(509):435-448.

Examples

# NOT RUN {
#################################
## Example on synethetic data. ## 
#################################
n <- 50
p <- 80
X <- matrix(rnorm(n*p,-3,3), nrow=n, ncol=p)
beta.true <- c(rep(2,5), rep(0,p-5)) # True beta has five entries of '2' and rest '0'

X <- scale(X) # Center and scale X
y <- crossprod(t(X), beta.true) + rnorm(n)

################################
# Fit the NBP regression model #
# using variational Bayes      #
################################
nbp.model <- nbp.VB(X=X, y=y, method.hyperparameters="mml", selection="dss")

nbp.model$beta.hat     # posterior mean estimates
nbp.model$beta.var     # posterior variance estimates
nbp.model$nbp.classifications    # Variables selected

# }