hsstan {hsstan}

Hierarchical shrinkage models

Description

Run the No-U-Turn Sampler (NUTS) as implemented in Stan to fit a hierarchical shrinkage model.

Usage

hsstan(
  x,
  covs.model,
  penalized = NULL,
  family = gaussian,
  iter = 2000,
  warmup = floor(iter/2),
  scale.u = 2,
  regularized = TRUE,
  nu = ifelse(regularized, 1, 3),
  par.ratio = 0.05,
  global.df = 1,
  slab.scale = 2,
  slab.df = 4,
  qr = TRUE,
  seed = 123,
  adapt.delta = NULL,
  keep.hs.pars = FALSE,
  ...
)

Arguments

x

Data frame containing outcome, covariates and penalized predictors. Continuous predictors and outcome variable should be standardized before fitting the models as priors assume them to have mean zero and unit variance.

covs.model

Formula containing the unpenalized covariates.

penalized

Names of the variables to be used as penalized predictors. Any variable that is already part of the covs.model formula will be penalized. If NULL or an empty vector, a model with only unpenalized covariates is fitted.

family

Type of model fitted: either gaussian() for linear regression (default) or binomial() for logistic regression.

iter

Total number of iterations in each chain, including warmup (2000 by default).

warmup

Number of warmup iterations per chain (by default, half the total number of iterations).

scale.u

Prior scale (standard deviation) for the unpenalized covariates.

regularized

If TRUE (default), the regularized horseshoe prior is used as opposed to the original horseshoe prior.

nu

Number of degrees of freedom of the half-Student-t prior on the local shrinkage parameters (by default, 1 if regularized=TRUE and 3 otherwise).

par.ratio

Expected ratio of non-zero to zero coefficients (ignored if regularized=FALSE). The scale of the global shrinkage parameter corresponds to par.ratio divided by the square root of the number of observations; for linear regression only, it's further multiplied by the residual standard deviation sigma.

global.df

Number of degrees of freedom for the global shrinkage parameter (ignored if regularized=FALSE). Larger values induce more shrinkage.

slab.scale

Scale of the regularization parameter (ignored if regularized=FALSE).

slab.df

Number of degrees of freedom of the regularization parameter (ignored if regularized=FALSE).

qr

Whether the thin QR decomposition should be used to decorrelate the predictors (TRUE by default). This is silently set to FALSE if there are more predictors than observations.

seed

Optional integer defining the seed for the pseudo-random number generator.

adapt.delta

Target average proposal acceptance probability for adaptation, a value between 0.8 and 1 (excluded). If unspecified, it's set to 0.99 for hierarchical shrinkage models and to 0.95 for base models.

keep.hs.pars

Whether the parameters for the horseshoe prior should be kept in the stanfit object returned (FALSE by default).

...

Further arguments passed to rstan::sampling(), such as chains (4 by default), cores (the value of options("mc.cores") by default), refresh (iter / 10 by default).

Value

An object of class hsstan containing the following fields:

stanfit

an object of class stanfit containing the output produced by Stan, including posterior samples and diagnostic summaries. It can be manipulated using methods from the rstan package.

betas

posterior means of the unpenalized and penalized regression parameters.

call

the matched call.

data

the dataset used in fitting the model.

model.terms

a list of names for the outcome variable, the unpenalized covariates and the penalized predictors.

family

the family object used.

hsstan.settings

the optional settings used in the model.

See Also

kfold() for cross-validating a fitted object.

Examples

data(diabetes)

# non-default settings for speed of the example
df <- diabetes[1:50, ]
hs.biom <- hsstan(df, Y ~ age + sex, penalized=colnames(df)[5:10],
                  chains=2, iter=250)