bayesreg_gl: MCMC Sampler for Linear Regression with Global-Local Priors

In drkowal/dsp: Dynamic Shrinkage Processes

MCMC Sampler for Linear Regression with Global-Local Priors

Description

Run the MCMC for linear regression with global-local priors on the regression coefficients. The priors include:

• the dynamic horseshoe prior ('DHS');

• the static horseshoe prior ('HS');

• the Bayesian lasso ('BL');

• the normal stochastic volatility model ('SV');

• the normal-inverse-gamma prior ('NIG').

In each case, the prior is a scale mixture of Gaussians. Sampling is accomplished with a (parameter-expanded) Gibbs sampler, mostly relying on a dynamic linear model representation.

Usage

bayesreg_gl(
  y,
  X,
  prior = "DHS",
  marginalSigma = TRUE,
  nsave = 1000,
  nburn = 1000,
  nskip = 4,
  mcmc_params = list("mu", "yhat", "beta", "evol_sigma_t2", "obs_sigma_t2", "dhs_phi",
    "dhs_mean"),
  computeDIC = TRUE,
  verbose = TRUE
)

Arguments

y	the n x 1 vector of response variables
X	the n x p matrix of predictor variables
prior	the prior distribution; must be one of 'DHS' (dynamic horseshoe prior), 'HS' (horseshoe prior), 'BL' (Bayesian lasso), or 'NIG' (normal-inverse-gamma prior)
marginalSigma	logical; if TRUE, marginalize over beta to sample sigma
nsave	number of MCMC iterations to record
nburn	number of MCMC iterations to discard (burin-in)
nskip	number of MCMC iterations to skip between saving iterations, i.e., save every (nskip + 1)th draw
mcmc_params	named list of parameters for which we store the MCMC output; must be one or more of: "mu" (conditional mean) "beta" (regression coefficients) "yhat" (posterior predictive distribution) "evol_sigma_t2" (evolution/prior error variance) "obs_sigma_t2" (observation error variance) "dhs_phi" (DHS AR(1) coefficient) "dhs_mean" (DHS AR(1) unconditional mean)
computeDIC	logical; if TRUE, compute the deviance information criterion DIC and the effective number of parameters p_d
verbose	logical; should R report extra information on progress?

Details

If p >= n, we use the fast sampling scheme of BHATTACHARYA et al (2016) for the regression coefficients, which is O(n^2*p). If p < n, we use the classical sampling scheme of RUE (2001), which is O(p^3).

Value

A named list of the nsave MCMC samples for the parameters named in mcmc_params

Note

The data y may contain NAs, which will be treated with a simple imputation scheme via an additional Gibbs sampling step. In general, rescaling y to have unit standard deviation is recommended to avoid numerical issues.

Examples

## Not run: 
# Case 1: p < n
n = 200; p = 50; RSNR = 10
X = matrix(rnorm(n*p), nrow = n)
beta_true = simUnivariate(signalName = 'levelshift', p, include_plot = FALSE)$y_true
#beta_true = c(rep(1, 10), rep(0, p - 10))
y_true = X%*%beta_true; sigma_true = sd(y_true)/RSNR
y = y_true + sigma_true*rnorm(n)
out = bayesreg_gl(y, X)

plot_fitted(lm(y ~ X - 1)$coef,
          mu = colMeans(out$beta),
          postY = out$beta,
          y_true = beta_true)

# Case 2: p > n
n = 200; p = 1000; RSNR = 10
X = matrix(rnorm(n*p), nrow = n)
#beta_true = simUnivariate(signalName = 'levelshift', p, include_plot = FALSE)$y_true
beta_true = c(rep(1, 20), rep(0, p - 20))
y_true = X%*%beta_true; sigma_true = sd(y_true)/RSNR
y = y_true + sigma_true*rnorm(n)
out = bayesreg_gl(y, X)

plot_fitted(rep(0, p),
          mu = colMeans(out$beta),
          postY = out$beta,
          y_true = beta_true)

## End(Not run)

drkowal/dsp documentation built on July 19, 2023, 11:42 a.m.