cv.SFGAM {sparseGAM}

Sparse Frequentist Generalized Additive Models

Description

This function implements K-fold cross-validation for sparse frequentist generalized additive models (GAMs) with the group LASSO, group SCAD, and group MCP penalties. The identity link function is used for Gaussian GAMs, the logit link is used for binomial GAMs, and the log link is used for Poisson, negative binomial, and gamma GAMs.

Usage

cv.SFGAM(y, X, df=6, 
         family=c("gaussian","binomial", "poisson", "negativebinomial","gamma"),
         nb.size=1, gamma.shape=1, penalty=c("gLASSO","gMCP","gSCAD"), taper,
         nfolds=10, nlambda=100, lambda, max.iter=10000, tol=1e-4)

Arguments

y

n \times 1 vector of responses.

X

n \times p design matrix, where the jth column of X corresponds to the jth overall covariate.

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

penalty

group regularization method to use on the groups of basis coefficients. The options are "gLASSO", "gSCAD", and "gMCP". To implement sparse GAMs with the SSGL penalty, use the SBGAM function.

taper

tapering term \gamma in group SCAD and group MCP controlling how rapidly the penalty tapers off. Default is taper=4 for group SCAD and taper=3 for group MCP. Ignored if "gLASSO" is specified as the penalty.

nfolds

number of folds K to use in K-fold cross-validation. Default is nfolds=10.

nlambda

number of regularization parameters L. Default is nlambda=100.

lambda

grid of L regularization parameters. The user may specify either a scalar or a vector. If the user does not provide this, the program chooses the grid automatically.

max.iter

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

tol

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

Value

The function returns a list containing the following components:

lambda

L \times 1 vector of regularization parameters lambda used to fit the model. lambda is displayed in descending order.

cve

L \times 1 vector of mean cross-validation error across all K folds. The kth entry in cve corresponds to the kth regularization parameter in lambda.

cvse

L \times 1 vector of standard errors for cross-validation error across all K folds. The kth entry in cvse corresponds to the kth regularization parameter in lambda.

lambda.min

value of lambda that minimizes mean cross-validation error cve.

References

Breheny, P. and Huang, J. (2015). "Group descent algorithms for nonconvex penalized linear and logistic regression models with grouped predictors." Statistics and Computing, 25:173-187.

Wang, H. and Leng, C. (2007). "Unified LASSO estimation by least squares approximation." Journal of the American Statistical Association, 102:1039-1048.

Yuan, M. and Lin, Y. (2006). Model selection and estimation in regression with grouped variables. Journal of the Royal Statistical Society: Series B (Statistical Methodology), 68: 49-67.

Examples

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

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

## Fit sparse Gaussian generalized additive model to data with the MCP penalty
gam.mod = SFGAM(y, X, X.test, family="gaussian", penalty="gMCP")

## The model corresponding to the 75th tuning parameter
gam.mod$lambda[75] 
gam.mod$classifications[,75] ## The covariate index is listed first

## Plot first function f_1(x_1) in 75th model
x1 = X.test[,1] 
## Estimates of all 20 function evaluations on test data
f.hat = gam.mod$f.pred[[75]] 
## 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])))