BoomSpikeSlab (version 1.2.6)

shrinkage.regression: Shrinking Regression Coefficients

Description

Fits a Bayesian regression model with a shrinkage prior on the coefficient. The model is

$$% y_i \sim N(x_i \beta, \sigma^2) \\ % 1 / \sigma^2 \sim Gamma(df/2, ss/2) \\ % g_1(\beta) \sim N(b1, v1) \\ % g_2(\beta) \sim N(b2, v2) \\ % \dots $$

In this notation, \(g_k(\beta) \sim N(b_k, v_k)\) indicates that the subset of coefficients in group k are a priori independent draws from the specified normal distribution. In addition, each subset-level prior may include a hyperprior, in which case the subset-level prior parameters will be updated as part of the MCMC. The hyperprior has the form of independent priors on the mean and precision parameters: $$% b_i ~ N(prior.mean, prior.variance) \\ % 1 / v_i ~ Chisq(df, guess.at.sd). \\ % $$

Usage

ShrinkageRegression(response, predictors, coefficient.groups,
                    residual.precision.prior = NULL,
                    suf = NULL, niter, ping = niter / 10,
                    seed = NULL)
CoefficientGroup(indices, mean.hyperprior = NULL, sd.hyperprior = NULL,
                 prior = NULL)

Value

ShrinkageRegression returns a list containing MCMC draws from the posterior distribution of model parameters. Each of the following is a matrix, with rows corresponding to MCMC draws, and columsn to distinct parameters.

• coefficients: regression coefficients.

• residual.sd: the residual standard deviation from the regression model.

• group.means: The posterior distribution of the mean of each coefficient group. If no mean hyperprior was assigned to a particular group, then the value here will be a constant (the values supplied by the prior argument to CoefficientGroup for that group).

• group.sds: The posterior distribution of the standard deviation of each coefficient group. If no sd.hyperprior was assigned to a particular group, then the value here will be a constant (the values supplied by the prior argument to CoefficientGroup for that group).

CoefficientGroup is a configuration utility used to define which coefficients should be shrunk together. It returns an object (list) formatted in the manner expected by

ShrinkageRegression.

Arguments

response

The numeric vector of responses.

predictors

The matrix of predictors, including an intercept term, if desired.

coefficient.groups

A list of objects of type CoefficientGroup, defining the pattern in which the coefficients should be shrunk together. Each coefficient must belong to exactly one CoefficientGroup.

residual.precision.prior

An object of type SdPrior describing the prior distribution of the residual standard deviation.

suf

An object of class RegressionSuf containing the sufficient statistics for the regression model. If this is NULL then it will be computed from response and predictors. If it is supplied then response and predictors are not used and can be left missing.

niter

The desired number of MCMC iterations.

ping

The frequency with which to print status updates.

seed

The integer-valued seed (or NULL) to use for the C++ random number generator.

indices

A vector of integers giving the positions of the regression coefficients that should be viewed as exchangeable.

mean.hyperprior

A NormalPrior object describing the hyperprior distribution for the average coefficient.

sd.hyperprior

An SdPrior object describing the hyperprior distribution for the standard deviation of the coefficients.

prior

An object of type NormalPrior giving the initial value of the distribution describing the collection of coefficients in this group. If either hyperprior is NULL then the corresponding prior parameter will not be updated. If both hyperpriors are non-NULL then this parameter can be left unspecified.

Author

Steven L. Scott

Examples

b0 <- -1
b1 <- rnorm(20, 3, .2)
b2 <- rnorm(30, -4, 7)
nobs <- 10000
beta <- c(b0, b1, b2)

X <- cbind(1, matrix(rnorm(nobs * (length(beta) - 1)), nrow = nobs, ncol = length(beta) - 1))
y.hat <- X %*% beta
y <- rnorm(nobs, y.hat, .5)

groups <- list(intercept = CoefficientGroup(1, prior = NormalPrior(0, 100)),
               first = CoefficientGroup(2:21,
                                        mean.hyperprior = NormalPrior(0, 100),
                                        sd.hyperprior = SdPrior(.2, 1)),
               second = CoefficientGroup(22:51,
                                         mean.hyperprior = NormalPrior(0, 100),
                                         sd.hyperprior = SdPrior(7, 1)))

model <- ShrinkageRegression(y, X, groups,
                             residual.precision.prior = SdPrior(.5, 1),
                             niter = 1000)