sparseGAM (version 1.0)

SSGL: Spike-and-Slab Group Lasso Regression

Description

This is a stand-alone function for group-regularized regression models in the exponential dispersion family with the spike-and-slab group lasso (SSGL) penalty. Let \(y_i\) denote the \(i\)th response and \(x_i\) denote a \(p\)-dimensional vector of covariates. We fit models of the form,

$$g(E(y_i)) = \beta_0 + x_i^T \beta, i = 1, ..., n,$$

where \(g\) is a monotone increasing link function. The identity link function is used for Gaussian regression, the logit link is used for binomial regression, and the log link is used for Poisson, negative binomial, and gamma regression.

If the covariates in each \(x_i\) are grouped according to known groups \(g=1, ..., G\), then this function may estimate some of the \(G\) groups of coefficients as all zero, depending on the amount of regularization.

Another implementation of the SSGL model for Gaussian regression models is available on Github at https://github.com/jantonelli111/SSGL. This package sparseGAM also implements the SSGL model for binomial, Poisson, negative binomial, and gamma regression.

Usage

SSGL(y, X, X.test, groups, 
     family=c("gaussian","binomial","poisson","negativebinomial","gamma"), 
     nb.size=1, gamma.shape=1, weights, nlambda0=20, lambda0, lambda1, a, b, 
     max.iter=100, tol = 1e-6, print.iter=TRUE)

Arguments

y

\(n \times 1\) vector of responses for training data.

X

\(n \times p\) design matrix for training data, where the \(j\)th column of X corresponds to the \(j\)th overall covariate.

X.test

\(n_{test} \times p\) design matrix for test data to calculate predictions. X.test must have the same number of columns as X, but not necessarily the same number of rows. If no test data is provided or if in-sample predictions are desired, then the function automatically sets X.test=X in order to calculate in-sample predictions.

groups

\(p\)-dimensional vector of group labels. The \(j\)th entry in groups should contain either the group number or the name of the factor level that the \(j\)th covariate belongs to. groups must be either a vector of integers or factors.

family

exponential dispersion family. Allows for "gaussian", "binomial", "poisson", "negativebinomial", and "gamma". Note that for "negativebinomial", the size parameter must be specified, while for "gamma", the shape parameter must be specified.

nb.size

known size parameter \(\alpha\) in \(NB(\alpha,\mu_i)\) distribution for the negative binomial responses. Default is nb.size=1. Ignored if family is not "negativebinomial".

gamma.shape

known shape parameter \(\nu\) in \(Gamma(\mu_i,\nu)\) distribution for gamma responses. Default is gamma.shape=1. Ignored if family is not "gamma".

weights

group-specific, nonnegative weights for the penalty. Default is to use the square roots of the group sizes.

nlambda0

number of spike hyperparameters \(L\). Default is nlambda0=20.

lambda0

grid of \(L\) spike hyperparameters \(\lambda_0\). The user may specify either a scalar or a vector. If the user does not provide this, the program chooses the grid automatically.

lambda1

slab hyperparameter \(\lambda_1\) in the SSGL prior. Default is lambda1=1.

a

shape hyperparameter for the \(Beta(a,b)\) prior on the mixing proportion in the SSGL prior. Default is a=1.

b

shape hyperparameter for the \(Beta(a,b)\) prior on the mixing proportion in the SSGL prior. Default is b=dim(X)[2].

max.iter

maximum number of iterations in the algorithm. Default is max.iter=100.

tol

convergence threshold for algorithm. Default is tol=1e-6.

print.iter

Boolean variable for whether or not to print the current nlambda0 in the algorithm. Default is print.iter=TRUE.

Value

The function returns a list containing the following components:

lambda0

\(L \times 1\) vector of spike hyperpameters lambda0 used to fit the model. lambda0 is displayed in descending order.

beta0

\(L \times 1\) vector of estimated intercepts. The \(k\)th entry in beta0 corresponds to the \(k\)th spike hyperparameter in lambda0.

beta

\(p \times L\) matrix of estimated regression coefficients. The \(k\)th column in beta corresponds to the \(k\)th spike hyperparameter in lambda0.

mu.pred

\(n_{test} \times L\) matrix of predicted mean response values \(\mu_{test} = E(Y_{test})\) based on the test data in X.test (or training data X if no argument was specified forX.test). The \(k\)th column in mu.pred corresponds to the predictions for the \(k\)th spike hyperparameter in lambda0.

classifications

\(G \times L\) matrix of classifications, where \(G\) is the number of groups. An entry of "1" indicates that the group was classified as nonzero, and an entry of "0" indicates that the group was classified as zero. The \(k\)th column of classifications corresponds to the \(k\)th spike hyperparameter in lambda0.

loss

vector of either the residual sum of squares ("gaussian") or the negative log-likelihood ("binomial", "poisson", "negativebinomial", "gamma") of the fitted model. The \(k\)th entry in loss corresponds to the \(k\)th spike hyperparameter in lambda0.

References

Bai R. (2021). "Spike-and-slab group lasso for consistent Bayesian estimation and variable selection in non-Gaussian generalized additive models." arXiv pre-print arXiv:2007.07021.

Bai, R., Moran, G. E., Antonelli, J. L., Chen, Y., and Boland, M.R. (2021). "Spike-and-slab group lassos for grouped regression and sparse generalized additive models." Journal of the American Statistical Association, in press.

Examples

# NOT RUN {
## Generate data
set.seed(12345)
X = matrix(runif(100*10), nrow=100)
n = dim(X)[1]
groups = c("A","A","A","B","B","B","C","C","D","D")
groups = as.factor(groups)
true.beta = c(-2.5,1.5,1.5,0,0,0,2,-2,0,0)

## Generate responses from Gaussian distribution
y = crossprod(t(X),true.beta) + rnorm(n)

## Generate test data
n.test = 50
X.test = matrix(runif(n.test*10), nrow=n.test)

## Fit SSGL model with 10 spike hyperparameters
## Note that if user does not specify lambda0, the SSGL function chooses a grid automatically.

SSGL.mod = SSGL(y, X, X.test, groups, family="gaussian", lambda0=seq(from=50,to=5,by=-5))

## Regression coefficient estimates
SSGL.mod$beta

# Predicted n.test-dimensional vectors mu=E(Y.test) based on test data, X.test. 
# The kth column of 'mu.pred' corresponds to the kth entry in 'lambda.'
SSGL.mod$mu.pred 

# Classifications of the 8 groups. The kth column of 'classifications'
# corresponds to the kth entry in 'lambda.'
SSGL.mod$classifications

# }
# NOT RUN {
## Example with binomial regression

## Generate binary responses
eta = crossprod(t(X), true.beta)
y = rbinom(n, size=1, prob=1/(1+exp(-eta)))

## Fit SSGL model with 10 spike hyperparameters
## Note that if user does not specify lambda0, the SSGL function chooses a grid automatically.

SSGL.mod = SSGL(y, X, X.test, groups, family="binomial", 
		lambda0=seq(from=10,to=1,by=-1))

## Predicted probabilities of success mu=E(Y.test) based on test data, X.test
SSGL.mod$mu.pred

## Classifications of the 8 groups. 
SSGL.mod$classifications
# }