IGG (version 1.0)

igg: Inverse Gamma-Gamma Regression

Description

This function provides a fully Bayesian approach for obtaining a sparse estimate of the \(p \times 1\) vector, \(\beta\) in the univariate linear regression model,

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

where \(\epsilon \sim N_n (0, \sigma^2 I_n)\). This is achieved by placing the inverse gamma-gamma (IGG) prior on the coefficients of \(\beta\). In the case where \(p > n\), we utilize an efficient sampler from Bhattacharya et al. (2016) to reduce the computational cost of sampling from the full conditional density of \(\beta\) to \(O(n^2 p)\).

Usage

igg(X, y, a=NA, b=NA, sigma2=NA, max.steps=10000, burnin=5000)

Arguments

X

\(n \times p\) design matrix.

y

\(n \times 1\) response vector.

a

The parameter for \(IG(a,1)\). If not specified (a=NA), defaults to \(1/2+1/p\). User may specify a value for a between 0 and 1.

b

The parameter for \(G(b,1)\). If not specified (b=NA), defaults to \(1/p\). User may specify a value for b between 0 and 1.

sigma2

The variance parameter. If the user does not specify this (sigma2=NA), the Gibbs sampler will estimate this using Jeffreys prior. If \(\sigma^2\) is known or estimated separately (e.g. through empirical Bayes), the user may also specify it.

max.steps

The total number of iterations to run in the Gibbs sampler. Defaults to 10,000.

burnin

The number of burn-in iterations for the Gibbs sampler. Defaults to 5,000.

Value

The function returns a list containing the following components:

beta.hat

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

beta.med

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

beta.intervals

The lower and upper endpoints of the 95 percent credible intervals for all \(p\) estimates in \(\beta\).

igg.classifications

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

Details

The function performs sparse estimation of \(\beta\) in the standard linear regression model and variable selection from the \(p\) covariates. Variable selection is performed by assessing the 95 percent marginal posterior credible intervals. The lower and upper endpoints of the 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 | ( \lambda_i, \xi_i, \sigma^2) \sim N(0, \sigma^2 \lambda_i \xi_i), i = 1, ..., p,$$ $$\lambda_i \sim IG(a,1), i = 1, ..., p,$$ $$\xi_i \sim G(b,1), i = 1, ..., p,$$ $$\sigma^2 \propto 1/\sigma^2.$$

If \(\sigma^2\) is known or estimated separately, the Gibbs sampler does not sample from the full conditional for \(\sigma^2\).

References

Bai, R. and Ghosh, M. (2018). "The Inverse Gamma-Gamma Prior for Optimal Posterior Contraction and Multiple Hypothesis Testing." Submitted, arXiv:1711.07635.

Bhattacharya, A., Chakraborty, A., and Mallick, B.K. (2016). "Fast Sampling with Gaussian Scale Mixture Priors in High-Dimensional Regression." Biometrika, 69(2): 447-457.

Examples

# NOT RUN {
###################### 
# Load diabetes data #
######################
data(diabetes)
attach(diabetes)
X <- scale(diabetes$x)
y <- scale(diabetes$y)

################################
# Fit the IGG regression model #
################################
igg.model <- igg(X=X, y=y, max.steps=5000, burnin=2500)

##############################
# Posterior median estimates #
##############################
igg.model$beta.med

###########################################
# 95 percent posterior credible intervals #
###########################################
igg.model$beta.intervals

######################
# Variable selection #
######################
igg.model$igg.classifications

# }