SBGAM {sparseGAM}

Sparse Bayesian Generalized Additive Models

Description

This function implements sparse Bayesian generalized additive models (GAMs) with the spike-and-slab group lasso (SSGL) penalty. Let y_i denote the ith response and x_i denote a p-dimensional vector of covariates. GAMs are of the form,

g(E(y_i)) = \beta_0 + \sum_{j=1}^{p} f_j (x_{ij}), 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. With the SSGL penalty, some of the univariate functions f_j(x_j) will be estimated as \hat{f}_j(x_j) = 0, depending on the size of the spike hyperparameter \lambda_0 in the SSGL prior. The functions f_j(x_j), j = 1, ..., p, are modeled using B-spline basis expansions.

There is another implementation of sparse Gaussian GAMs with the SSGL penalty available at https://github.com/jantonelli111/SSGL, which uses natural cubic splines as the basis functions. This package sparseGAM uses B-spline basis functions and also implements sparse GAMs with the SSGL penalty for binomial, Poisson, negative binomial, and gamma regression.

For implementation of sparse frequentist GAMs with the group LASSO, group SCAD, and group MCP penalties, use the SFGAM function.

Usage

SBGAM(y, X, X.test, df=6, 
      family=c("gaussian","binomial","poisson","negativebinomial","gamma"), 
      nb.size=1, gamma.shape=1, 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 jth column of X corresponds to the jth 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.

df

number of B-spline basis functions to use in each basis expansion. Default is df=6, but the user may specify degrees of freedom as any integer greater than or equal to 3.

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 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".

nlambda0

number of spike hyperparameter 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 hyperparameters lambda0 used to fit the model. lambda0 is displayed in descending order.

f.pred

List of L n_{test} \times p matrices, where the kth matrix in the list corresponds to the kth spike hyperparameter in lambda0. The jth column in each matrix in f.pred is the estimate of the jth function evaluated on the test data in X.test for the jth covariate (or training data X if X.test was not specified).

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 for X.test). The kth column in mu.pred corresponds to the predictions for the kth spike hyperparameter in lambda0.

classifications

p \times L matrix of classifications. An entry of "1" indicates that the corresponding function was classified as nonzero, and an entry of "0" indicates that the function was classified as zero. The kth column of classifications corresponds to the kth spike hyperparameter in lambda0.

beta0

L \times 1 vector of estimated intercepts. The kth entry in beta0 corresponds to the kth spike hyperparameter in lambda0.

beta

dp \times L matrix of estimated basis coefficients. The kth column in beta corresponds to the kth 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 kth entry in loss corresponds to the kth 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

## Generate data
set.seed(12345)
X = matrix(runif(100*5), nrow=100)
n = dim(X)[1]
y = 3*sin(2*pi*X[,1])-3*cos(2*pi*X[,2]) + rnorm(n)

## Test data with 30 observations
X.test = matrix(runif(30*5), nrow=30)

## Fit sparse Bayesian generalized additive model to data with the SSGL penalty
## and 5 spike hyperparameters
SBGAM.mod = SBGAM(y, X, X.test, family="gaussian", lambda0=seq(from=50,to=10,by=-10))

## The model corresponding to the 1st spike hyperparameter
SBGAM.mod$lambda[1] 
SBGAM.mod$classifications[,1] 

## Plot first function f_1(x_1) in 2nd model
x1 = X.test[,1] 
## Estimates of all 20 function evaluations on test data
f.hat = SBGAM.mod$f.pred[[1]] 
## Extract estimates of f_1 
f1.hat = f.hat[,1] 

## Plot X_1 against f_1(x_1)
plot(x1[order(x1)], f1.hat[order(x1)], xlab=expression(x[1]), 
     ylab=expression(f[1](x[1])))