ssnet: Bayesian Spike-and-Slab Elastic Net with Spatial Structure

In jmleach-bst/ssnet: Fit Spike-and-Slab Elastic Net GLM With Spatial Priors on Parameters

Bayesian Spike-and-Slab Elastic Net with Spatial Structure

Description

Fits generalized linear models with spike-and-slab priors whose (logit of) probability of inclusion in the model is either assigned an intrinsic autoregression as a prior to incorporate spatial information or is unstructured. The model is fit using an EM algorithm where the E-step is fit by glmnet() and the M-step is fit using the stan function optimizing (when IAR prior is employed).

Usage

ssnet(
  x,
  y,
  family = c("gaussian", "binomial", "multinomial", "poisson", "cox"),
  offset = NULL,
  epsilon = 1e-04,
  alpha = 0.95,
  type.multinomial = "grouped",
  maxit = 50,
  init = rep(0, ncol(x)),
  init.theta = 0.5,
  ss = c(0.04, 0.5),
  Warning = FALSE,
  group = NULL,
  iar.prior = FALSE,
  adjmat = NULL,
  iar.data = NULL,
  tau.prior = "none",
  tau.manual = NULL,
  stan_manual = NULL,
  opt.algorithm = "LBFGS",
  p.bound = c(0.01, 0.99),
  plot.pj = FALSE,
  im.res = NULL,
  verbose = FALSE,
  print.iter = FALSE
)

Arguments

x	Design, or input, matrix, of dimension nobs x nvars; each row is an observation vector. It is recommended that x have user-defined column names for ease of identifying variables. If missing, then colnames are internally assigned x1, x2, ... and so forth.
y	Scalar response variable. Quantitative for family = "gaussian", or family = "poisson" (non-negative counts). For family = "gaussian", y is always standardized. For family = "binomial", y should be either a factor with two levels, or a two-column matrix of counts or proportions (the second column is treated as the target class; for a factor, the last level in alphabetical order is the target class). For family="cox", y should be a two-column matrix with columns named 'time' and 'status'. The latter is a binary variable, with '1' indicating death, and '0' indicating right censored. The function Surv() in package survival produces such a matrix. When family = "multinomial", y follows the documentation for glmnet, but it is preferred that y is a factor with two or more levels.
family	Response type (see above).
offset	A vector of length nobs that is included in the linear predictor.
epsilon	A positive convergence tolerance; the iterations converge when |dev - dev_old|/(|dev| + 0.1) < e.
alpha	A scalar value between 0 and 1 determining the compromise between the Ridge and Lasso models. When alpha = 1 reduces to the Lasso, and when alpha = 0 reduces to Ridge.
type.multinomial	If "grouped" then a grouped lasso penalty is used on the multinomial coefficients for a variable. This ensures they are all in our out together. The default is "ungrouped"
maxit	An integer giving the maximal number of EM iterations.
init	A vector of initial values for all coefficients (not for intercept). If not given, it will be internally produced. If family = "multinomial" and the same initializations are desired for each response/outcome category then init can be a vector. If different initializations are desired, then init should be a list, each element of which contains a vector of initializations. The list should be named according the response/outcome category as they appear in y.
init.theta	A single value between 0 and 1 to initialize inclusion probabilities. When parameter groups is 2 or more, can be a vector of initial values for parameter inclusion probabilities. Default is 0.5 for all parameters. Currently, it is not supported to have parameter-specific initializations for inclusion probabilities. This may change in forthcoming updates.
ss	A vector of two positive scale values for the spike-and-slab mixture double-exponential prior, allowing for different scales for different predictors, leading to different amount of shrinkage. Smaller scale values give stronger shrinkage. While the smaller of the two input values will be treated as the spike scale, it is recommended to specify the spike scale as the first element of the vector.
Warning	Logical. If TRUE, shows the error messages of not convergence and identifiability.
group	A numeric vector, or an integer, or a list indicating the groups of predictors. If group = NULL, all the predictors form a single group. If group = K, the predictors are evenly divided into groups each with K predictors. If group is a numberic vector, it defines groups as follows: Group 1: (group[1]+1):group[2], Group 2: (group[2]+1):group[3], Group 3: (group[3]+1):group[4], ... If group is a list of variable names, group[[k]] includes variables in the k-th group. The mixture double-exponential prior is only used for grouped predictors. For ungrouped predictors, the prior is double-exponential with scale ss[2] and mean 0. Note that grouped predictors when family = "multinomial" is still experimental, so use with caution.
iar.prior	Logical. When TRUE, imposes intrinsic autoregressive prior on logit of the probabilities of inclusion. When FALSE, treats probabilities of inclusion as unstructured.
adjmat	A data.frame or matrix containing a "sparse" representation of the neighbor relationships. The first column should contain a numerical index for a given location. Each index will be repeated in this column for every neighbor it has. The indices for the location's neighbors are then specified in the second column. Any additional columns are ignored.
iar.data	A list of output from mungeCARdata4stan that contains the necessary inputs for the IAR prior. When unspecified, this is built internally assuming that neighbors are those variables directly above, below, left, and right of a given variable location. im.res must be specified when allowing this argument to be built internally. It is not recommended to use this argument directly, even when specifying a more complicated neighborhood stucture; this can be specified with the adjmat argument, and then internally converted to the correct format.
tau.prior	One of c("none", "manual", "cauchy"). This argument determines the precision parameter in the Conditional Autoregressive model for the (logit of) prior inclusion probabilities. When "none", the precision is set to 1; when "manual", the precision is manually entered by the user; when "cauchy", the inverse precision is assumed to follow a Cauchy distribution with mean 0 and scale 2.5. Note that at this stage of development, only the "none" option has been extensively tested, so the other options should be used with caution.
tau.manual	When tau.prior = "manual", use this argument to specify a common precision parameter.
stan_manual	A stan_model that is manually specified. Especially when fitting multiple models in succession, specifying the stan model outside this "loop" may avoid errors.
opt.algorithm	One of c("LBFGS", "BFGS", "Newton"). This argument determines which argument is used to optimize the term in the EM algorithm that estimates the probabilities of inclusion for each parameter. Optimization is performed by optimizing.
p.bound	A vector defining the lower and upper boundaries for the probabilities of inclusion in the model, respectively. Defaults to c(0.01, 0.99).
plot.pj	When TRUE, prints a series of 2D graphs of the prior probabilities of inclusion at each step of the algorithm. This should NOT be used for 3D data.
im.res	A 2-element vector where the first argument is the number of "rows" and the second argument is the number of "columns" in each subject's "image". Default is NULL.
verbose	Logical. If TRUE, prints out the number of iterations and computational time.
print.iter	Logical. When TRUE, prints a quite excessive amount intermediate output for every iteration. Default is (obviously) FALSE.

Value

An object of class c("elnet" "glmnet" "bmlasso" "GLM").

Note

If iar.data = NULL, i.e. is left unspecified, then provided that im.res is specified, the function proximity_builder() from the package sim2Dpredictr builds the appropriate list of data for optimization with stan. Currently, im.res can only handle 2D data. Future versions may allow images to be 3D. However, the function will work given any appropriately specified neighborhood matrix, whatever the original dimension.

While the type.multinomial option is included, it is only valid for traditional elastic net models. Thus far we have only extended the spike-and-slab models for grouped selection.

References

\insertRef

Banerjee:2015ssnet

\insertRef

Friedman:2007ssnet

\insertRef

Friedman:2010ssnet

\insertRef

Morris:2017ssnet

\insertRef

Morris:2019ssnet

\insertRef

Rockova+George:2018ssnet

\insertRef

Tang:2017ssnet

\insertRef

Tang:2018ssnet

Examples

library(sim2Dpredictr)
set.seed(223)

## sample size
n <- 30
## image dims
nr <- 4
nc <- 4

## generate data
cn <- paste0("x", seq_len(nr * nc))
tb <- rbinom(nr * nc, 1, 0.2)
tx <- matrix(rnorm(n * nr * nc), nrow = n, ncol = nr * nc,
             dimnames = list(seq_len(n), cn))
ty <- tx %*% tb + rnorm(n)

## build adjacency matrix
adjmat <- proximity_builder(im.res = c(nr, nc), type = "sparse")
## stan model information
model_info <- mungeCARdata4stan(adjmat$nb.index,
                                table(adjmat$location.index))

## fit model
ex_model <- ssnet(x = tx, y = ty, alpha = 0.5,
                  iar.prior = TRUE, iar.data = model_info,
                  family = "gaussian")

jmleach-bst/ssnet documentation built on March 4, 2024, 5:04 p.m.