bsSpike: Bayesian Single-Index Regression with B-Spline Link and Spike-and-Slab Prior

Description

Fits a single-index model $Y_i \sim \mathcal{N}(f(X_i'\theta), \sigma^2), i = 1,\cdots,n$ where the link $f(\cdot)$ is represented by B-spline link function and the index vector $\theta$ has spike-and-slab prior.

Usage

bsSpike(
  formula,
  data,
  prior = NULL,
  init = NULL,
  monitors = NULL,
  niter = 10000,
  nburnin = 1000,
  thin = 1,
  nchain = 1,
  setSeed = FALSE
)
bsSpike_setup(
  formula,
  data,
  prior = NULL,
  init = NULL,
  monitors = NULL,
  niter = 10000,
  nburnin = 1000,
  thin = 1,
  nchain = 1,
  setSeed = FALSE
)

Value

A list typically containing:

coefficients: Mean estimates of index vector. Return list of model_setup does not include it.
ses: Mean standard error of index vector. Return list of model_setup does not include it.
residuals: Residuals with mean estimates of fitted values. Return list of model_setup does not include it.
fitted.values: Mean estimates of fitted values. Return list of model_setup does not include it.
linear.predictors: Mean estimates of single-index values. Return list of model_setup does not include it.
gof: Goodness-of-fit. Return list of model_setup function does not include it.
samples: Posterior draws of variables for computing fitted values of the model, including $\theta$, $\sigma^2$. Return list of model_setup does not include it.
input: List of data used in modeling, formula, input values for prior and initial values, and computation time without compiling.
defModel: Nimble model object.
defSampler: Nimble sampler object.
modelName: Name of the model.

Arguments

formula: an object of class formula. Interaction term is not allowed for single-index model.
data: an data frame.
prior: Optional named list of prior settings. Further descriptions are in Details section.
init: Optional named list of initial values. If the values are not assigned, they are randomly sampled from prior or designated value. Further descriptions are in Details section.
monitors: Optional character vector of additional monitor nodes. To check the names of the nodes, fit the model_setup function and then inspect the variable names stored in the model object using getVarMonitor.
niter: Integer. Total MCMC iterations (default 10000).
nburnin: Integer. Burn-in iterations (default 1000).
thin: Integer. Thinning for monitors (default 1).
nchain: Integer. Number of MCMC chains (default 1).
setSeed: Logical or numeric argument. Further details are provided in runMCMC setSeed argument.

Details

Model The single–index model uses a $m$-order polynomial spline with $k$ interior knots as follows: $$f(t) = \sum_{j=1}^{m+k} B_j(t)\,\beta_j$$ on $[a, b]$ with $t_i = X_i' \theta, i = 1,\cdots, n$ and $\|\theta\|_2 = 1$. $\{\beta_j\}_{j=1}^{m+k}$ are spline coefficient and $a_\theta$ and $ b_\theta$ are boundary knots where $a_\theta = min(t_i, i = 1, \cdots, n) - \delta$, and $b_\theta = max(t_i, i = 1,\cdots, n) + \delta$. $\theta$ is a p-dimensional index vector subject to a spike-and-slab prior for variable selection with binary indicator variable $\nu$.

Priors

The variable selection indicator $\nu$ has Beta–Bernoulli hierarchy prior. Beta hyper-prior on the Bernoulli–inclusion probability $w$, inducing $p(\nu) \propto \mathrm{Beta}(r_1 + n_\nu, r_2 + p - n_\nu)$ where $n_\nu = \Sigma_{i=1}^{p}I(\nu_i = 1)$. $r_1, r_2$ are shape and rate parameter of beta distribution.
Slab coefficients $\theta$ have normal distribution with zero mean and $\sigma_\theta^2$ variance.
Conditioned on $\theta$ and $\sigma^2$, the link coefficients $\beta = (\beta_1,\ldots,\beta_{m+k})^\top$ follow normal distribution with estimated mean vector $\hat{\beta}_{\theta} = (X_{\theta}'X_{\theta})^{-1}X_{\theta}'Y$ and covariance $\sigma^2 (X_{\theta}^\top X_{\theta})^{-1}$, where $X_{\theta}$ is the B-spline basis design matrix.
Inverse gamma prior on $\sigma^2$ with shape parameter $a_\sigma$ and rate parameter $b_\sigma$.

Sampling Samplers are automatically assigned by nimble.

Prior hyper-parameters These are the prior hyper-parameters set in the function. You can define new values for each parameter in prior_param.

Index vector: index_nu_r1, index_nu_r2 gives the shape and rate parameter of beta-binomial prior, respectively. For slab prior, normal distribution with zero mean is assigned for selected variables $\theta$. index_sigma_theta is for variance of $\theta$, and it is assigned 0.25 by default.
Link function: B-spline basis and coefficient of B-spline setup.
- basis: For the basis of B-spline, link_basis_df is the number of basis functions (default 21), link_basis_degree is the spline degree (default 2) and link_basis_delta is a small jitter for boundary-knot spacing control (default 0.01).
- beta: For the coefficient of B-spline, multivariate normal prior is assigned with mean link_beta_mu, and covariance link_beta_cov. By default, $\mathcal{N}_p\!\big(0, \mathrm{I}_p\big)$
Error variance (sigma2): Inverse gamma prior is assigned to $\sigma^2$ where sigma2_shape is shape parameter and sigma2_rate is rate parameter of inverse gamma distribution. (default sigma2_shape = 0.001, sigma2_rate = 100)

Initial values These are the initial values set in the function. You can define new values for each initial value in init_param.

Index vector:
- index_pi Initial selecting variable probability. (default: 0.5)
- index_nu Initial vector of inclusion indicators . By default, each nu is randomly drawn by $Bernoulli(1/2)$
- index Initial vector of index. By default, each element of index vector, which is chosen by $\nu$, is proposed by normal distribution.
Link function: Initial spline coefficients (link_beta). By default, $\big(X_{\theta}^\top X_{\theta} + \rho\, \mathrm{I}\big)^{-1} X_{\theta}^\top Y$ is computed, where $X_{\theta}$ is the B-spline basis design matrix.
Error variance (sigma2): Initial scalar error variance (default 0.01).

References

Antoniadis, A., Grégoire, G., & McKeague, I. W. (2004). Bayesian estimation in single-index models. Statistica Sinica, 1147-1164.

Hornik, K., & Grün, B. (2014). movMF: an R package for fitting mixtures of von Mises-Fisher distributions. Journal of Statistical Software, 58, 1-31.

McGee, G., Wilson, A., Webster, T. F., & Coull, B. A. (2023). Bayesian multiple index models for environmental mixtures. Biometrics, 79(1), 462-474.

Examples

Run this code

# \donttest{
set.seed(123)
n <- 200; d <- 4
theta <- c(2, 1, 1, 1); theta <- theta / sqrt(sum(theta^2))
f <- function(u) u^2 * exp(u)
sigma <- 0.5
X <- matrix(runif(n * d, -1, 1), nrow = n)
index_vals <- as.vector(X %*% theta)
y <- f(index_vals) + rnorm(n, 0, sigma)
simdata <- data.frame(x = X, y = y)
colnames(simdata) <- c(paste0("X", 1:4), "y")

# One tool version
fit1 <- bsSpike(y ~ ., data = simdata,
                niter = 5000, nburnin = 1000, nchain = 1)

# Split version
models <- bsSpike_setup(y ~ ., data = simdata)
Ccompile <- compileModelAndMCMC(models)
nimSampler <- get_sampler(Ccompile)
initList <- getInit(models)
mcmc.out <- runMCMC(nimSampler, niter = 5000, nburnin = 1000, thin = 1,
                   nchains = 1, setSeed = TRUE, inits = initList,
                   summary = TRUE, samplesAsCodaMCMC = TRUE)
fit2 <- as_bsim(models, mcmc.out)
summary(fit2)
# }

Run the code above in your browser using DataLab