SSGL {SSGL} | R Documentation |
Spike-and-Slab Group Lasso for Group-Regularized Generalized Linear Models (GLMs)
Description
The SSGL
function implements maximum a posteriori (MAP) estimation for group-regularized GLMs with the spike-and-slab group lasso (SSGL) penalty of Bai et al. (2022) and Bai (2023). 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 regression. If the covariates in each x_i
are grouped according to known groups g=1, ..., G
, then this function can estimate some of the G
groups of coefficients as all zero, depending on the amount of regularization.
This function only returns point estimates. Please refer to the SSGL_gibbs
function if uncertainty quantification of the model parameters is desired. In general, we recommend using SSGL
for estimation and variable selection and SSGL_gibbs
for uncertainty quantification.
The SSGL
function also has the option of returning the generalized information criterion (GIC) of Fan and Tang (2013) for each regularization parameter in the grid lambda0
. The GIC can be used for model selection and serves as a useful alternative to cross-validation. The formula for the GIC and a given \lambda_0
is
DIC(\lambda_0) = \frac{1}{n} Deviance_{\lambda_0} + a_n \times \nu),
where Deviance_{\lambda_0}
is the deviance computed with the estimate of beta
based on spike hyperparameter \lambda_0
, \nu_0
is the number of nonzero elements in the estimated beta
, and a_n
is a sequence that diverges at a suitable rate relative to n
. As recommended by Fan and Tang (2013), we set a_n = \{\log(\log(n))\}\log(p)
.
If cross-validation is preferred for tuning \lambda_0
, please refer to the SSGL_cv
function.
Usage
SSGL(Y, X, groups, family=c("gaussian","binomial","poisson"),
X_test, group_weights, n_lambda0=25,
lambda0, lambda1=1, a=1, b=length(unique(groups)),
max_iter=100, tol = 1e-6, return_GIC=TRUE, print_lambda0=TRUE)
Arguments
Y |
|
X |
|
groups |
|
family |
exponential dispersion family of the response variables. Allows for |
X_test |
|
group_weights |
group-specific, nonnegative weights for the penalty. Default is to use the square roots of the group sizes. |
n_lambda0 |
number of spike hyperparameters |
lambda0 |
grid of |
lambda1 |
slab hyperparameter |
a |
shape hyperparameter for the |
b |
shape hyperparameter for the |
max_iter |
maximum number of iterations in the algorithm. Default is |
tol |
convergence threshold for algorithm. Default is |
return_GIC |
Boolean variable for whether or not to return the GIC. Default is |
print_lambda0 |
Boolean variable for whether or not to print the current value in |
Value
The function returns a list containing the following components:
lambda0 |
|
beta |
|
beta0 |
|
classifications |
|
Y_pred |
|
GIC |
|
lambda0_GIC_min |
The value in |
min_GIC_index |
The index of |
References
Bai, R. (2023). "Bayesian group regularization in generalized linear models with a continuous spike-and-slab prior." arXiv pre-print arXiv:2007.07021.
Bai, R., Moran, G. E., Antonelli, J. L., Chen, Y., and Boland, M.R. (2022). "Spike-and-slab group lassos for grouped regression and sparse generalized additive models." Journal of the American Statistical Association, 117:184-197.
Fan, Y. and Tang, C. Y. (2013). "Tuning parameter selection in high dimensional penalized likelihood." Journal of the Royal Statistical Society: Series B (Statistical Methodology), 75:531-552.
Examples
## 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)
beta_true = c(-2.5,1.5,1.5,0,0,0,2,-2,0,0)
## Generate responses from Gaussian distribution
Y = crossprod(t(X), beta_true) + 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: If you do not specify lambda0, the program will automatically choose a suitable grid.
SSGL_mod = SSGL(Y, X, groups, family="gaussian", X_test, 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 'Y_pred' corresponds to the kth entry in 'lambda.'
SSGL_mod$Y_pred
## Classifications of the 8 groups. The kth column of 'classifications'
## corresponds to the kth entry in 'lambda.'
SSGL_mod$classifications
## Plot lambda vs. GIC
plot(SSGL_mod$lambda0, SSGL_mod$GIC, type='l')
## Model selection with the lambda that minimizes GIC
SSGL_mod$lambda0_GIC_min
SSGL_mod$min_GIC_index
SSGL_mod$classifications[, SSGL_mod$min_GIC_index]
SSGL_mod$beta[, SSGL_mod$min_GIC_index]
## Example with Poisson regression
## Generate data
set.seed(1234)
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)
beta_true = c(-2.5,1.5,1.5,0,0,0,2,-2,0,0)
## Generate count responses
eta = crossprod(t(X), beta_true)
Y = rpois(n, exp(eta))
## Generate test data
n_test = 50
X_test = matrix(runif(n_test*10), nrow=n_test)
## Fit SSGL model
SSGL_poisson_mod = SSGL(Y, X, groups, family="poisson")
## Regression coefficient estimates
SSGL_poisson_mod$beta
## Predicted n_test-dimensional vectors mu=E(Y.test) based on test data, X_test.
## The kth column of 'Y_pred' corresponds to the kth entry in 'lambda.'
SSGL_poisson_mod$Y_pred
## Classifications of the 8 groups. The kth column of 'classifications'
## corresponds to the kth entry in 'lambda.'
SSGL_poisson_mod$classifications
## Plot lambda vs. GIC
plot(SSGL_poisson_mod$lambda0, SSGL_poisson_mod$GIC, type='l')
## Model selection with the lambda that minimizes GIC
SSGL_poisson_mod$lambda0_GIC_min
SSGL_poisson_mod$min_GIC_index
SSGL_poisson_mod$classifications[, SSGL_mod$min_GIC_index]
SSGL_poisson_mod$beta[, SSGL_mod$min_GIC_index]