SSLASSO {SSLASSO}

The Spike-and-Slab LASSO

Description

Spike-and-Slab LASSO is a spike-and-slab refinement of the LASSO procedure, using a mixture of Laplace priors indexed by lambda0 (spike) and lambda1 (slab).

The SSLASSO procedure fits coefficients paths for Spike-and-Slab LASSO-penalized linear regression models over a grid of values for the regularization parameter lambda0. The code has been adapted from the ncvreg package (Breheny and Huang, 2011).

Usage

SSLASSO(X, y, penalty = c("adaptive", "separable"), variance = c("fixed", "unknown"), 
lambda1, lambda0, nlambda = 100, theta = 0.5, sigma, a = 1, b,  
eps = 0.001, max.iter = 500,  counter = 10, warn = FALSE)

Arguments

X

The design matrix (n x p), without an intercept. SSLASSO standardizes the data by default.

y

Vector of continuous responses (n x 1). The responses will be centered by default.

penalty

The penalty to be applied to the model. Either "separable" (with a fixed theta) or "adaptive" (with a random theta, where theta ~ B(a,p)). The default is "adaptive".

variance

Whether the error variance is also estimated. Either "fixed" (with a fixed sigma) or "unknown" (with a random sigma, where p(sigma) ~ 1/sigma). The default is "fixed".

lambda1

Slab variance parameter. Needs to be less than lambda0. The default is lambda0 = 1.

lambda0

Spike penalty parameters (L x 1). Either a numeric value for a single run (L=1) or a sequence of increasing values for dynamic posterior exploration. The default is lambda0 = seq(1, nrow(X), length.out = 100).

nlambda

The number of lambda0 values. Default is 100.

theta

Prior mixing proportion. For "separable" penalty, this value is fixed. For "adaptive" penalty, this value is used as a starting value.

sigma

Error variance. For "fixed" variance, this value is fixed. For "unknown" variance, this value is used as a starting value.

a

Hyperparameter of the beta prior B(a,b) for the adaptive penalty (default a = 1).

b

Hyperparameter of the beta prior B(a,b) for the adaptive penalty (default b = ncol(X)).

eps

Convergence criterion: converged when difference in regression coefficients is less than eps (default eps = 0.001).

max.iter

Maximum number of iterations. Default is 500.

counter

Applicable only for the adaptive penalty. Determines how often the parameter theta is updated throughout the cycles of coordinate ascent. Default is 10.

warn

TRUE if warnings should be printed; FALSE by default

Details

The sequence of models indexed by the regularization parameter lambda0 is fitted using a coordinate descent algorithm. The algorithm uses screening rules for discarding irrelevant predictors along the lines of Breheny (2011).

Value

An object with S3 class "SSLASSO" containing:

beta

The fitted matrix of coefficients (p x L). The number of rows is equal to the number of coefficients p, and the number of columns is equal to L (the length of lambda0).

intercept

A vector of length L containing the intercept for each value of lambda0. The intercept is intercept = mean(y) - crossprod(XX, beta), where XX is the centered design matrix.

iter

A vector of length L containing the number of iterations until convergence at each value of lambda0.

lambda0

The sequence of regularization parameter values in the path.

penalty

Same as above.

thetas

A vector of length L containing the hyper-parameter values theta (the same as theta for "separable" penalty).

sigmas

A vector of length L containing the values sigma (the same as the initial sigma for "known" variance).

select

A (p x L) binary matrix indicating which variables were selected along the solution path.

model

A single model chosen after the stabilization of the regularization path.

Author(s)

Veronika Rockova <[e-mail]>, Gemma Moran <[e-mail]>

References

Rockova, V. and George, E.I. (2018) The Spike-and-Slab LASSO. Journal of the American Statistical Association.

Moran, G., Rockova, V. and George, E.I. (2018) On variance estimation for Bayesian variable selection. <https://arxiv.org/abs/1801.03019>

See Also

plot.SSLASSO

Examples

## Linear regression, where p > n

library(SSLASSO)

p <- 1000
n <- 100

X <- matrix(rnorm(n*p), nrow = n, ncol = p)
beta <- c(1, 2, 3, rep(0, p-3))
y = X[,1] * beta[1] + X[,2] * beta[2] + X[,3] * beta[3] + rnorm(n)

# Oracle SSLASSO with known variance

result1 <- SSLASSO(X, y, penalty = "separable", theta = 3/p)
plot(result1)

# Adaptive SSLASSO with known variance

result2 <- SSLASSO(X, y)
plot(result2)

# Adaptive SSLASSO with unknown variance

result3 <- SSLASSO(X, y, variance = "unknown")
plot(result3)