Learn R Programming

BayesSIM (version 1.0.0)

gpSphere: Bayesian Single-Index Regression with Gaussian Process Link and Unit Sphere Prior

Description

Fits a single–index model \(Y_i \sim \mathcal{N}(f(X_i'\theta), \sigma^2), i = 1,\cdots,n\), where the index \(\theta\) lies on the unit sphere, and the link \(f(\cdot)\) is represented with Gaussian process.

Usage

gpSphere(
  formula,
  data,
  prior = NULL,
  init = NULL,
  method = "FB",
  lowerB = NULL,
  upperB = NULL,
  monitors = NULL,
  niter = 10000,
  nburnin = 1000,
  thin = 1,
  nchain = 1,
  setSeed = FALSE
)

gpSphere_setup( formula, data, prior = NULL, init = NULL, method = "FB", lowerB = NULL, upperB = 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.

method

Character, gpSphere model has 3 different types of sampling method, fully Bayesian method ("FB"), empirical Bayes approach ("EB"), and empirical Gibbs sampler ("EG"). Assign one sampler method. Empirical sampling approach is recommended in high-dimensional data. By default, fully Bayesian approach is assigned.

lowerB

This parameter is only for gpSphere model. Numeric vector of element-wise lower bounds for the "L-BFGS-B" method. When the empirical Bayes or Gibbs sampler method is used, the marginal likelihood is optimized via optim(method = "L-BFGS-B"). The vector must be ordered as c(index vector, lengthscale, amp, sigma2). Note that sigma2 is included only for the empirical Bayes method (omit it for Gibbs). By default, the lower bounds are -1 for the index vector and -1e2 for logarithm of lengthscale, amp, and (if present) sigma2.

upperB

This parameter is only for gpSphere model. Numeric vector of element-wise upper bounds for the "L-BFGS-B" method. When the empirical Bayes or Gibbs sampler method is used, the marginal likelihood is optimized via optim(method = "L-BFGS-B"). The vector must be ordered as c(index vector, lengthscale, amp, sigma2). Note that sigma2 is included only for the empirical Bayes method (omit it for Gibbs). By default, the upper bounds are 1 for the index vector and 1e2 for logarithm of lengthscale, amp, and (if present) sigma2.

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 Gaussian process with zero mean and and covariance kernel \(\eta \cdot \text{exp}(-\frac{(t_i-t_j)^2}{l})\) as a link function, where \(t_i, t_j, j = 1, \ldots, n\) is index value. Index vector should be length 1.

Priors

  • von Mises–Fisher prior on the index \(\theta\) with direction and concentration.

  • Covariance kernel: Amplitude, \(\eta\), follows log normal distribution with mean \(a_\eta\) and variance \(b_\eta\). Length-scale parameter follows gamma distribution with shape parameter \(\alpha_l\) and rate parameter \(\beta_l\).

  • Inverse-Gamma prior on \(\sigma^2\) with shape parameter \(a_\sigma\) and rate parameter \(b_\sigma\).

Sampling In the fully Bayesian approach, \(\theta\), \(l\), and \(\eta\) are updated via the Metropolis–Hastings algorithm, while \(f\) and \(\sigma^2\) are sampled using Gibbs sampling.

In the empirical Bayes approach, \(\theta\), \(l\), \(\eta\), and \(\sigma^2\) are estimated by maximum a posteriori (MAP), and \(f\) is sampled from its full conditional posterior distribution.

In the empirical Gibbs sampler, \(\theta\), \(l\), and \(\eta\) are estimated by MAP, whereas \(f\) and \(\sigma^2\) are sampled via Gibbs sampling.

For estimation via MAP, effective sample size or potential scale reduction factor is meaningless.

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

  1. Index vector: Nothing to assign.

  2. Link function:

    • Length-scale:Gamma distribution is assigned for length-scale parameter, \(l\). link_lengthscale_shape is shape parameter (default 1/8) and link_lengthscale_rate is rate parameter of lengthscale. (default 1/8)

    • Amplitude: Log-normal distribution is assigned for amplitude parameter, \(\eta\). link_amp_a is mean (default -1), and link_amp_b is variance. (default 1)

  3. 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 = 1, sigma2_rate = 1)

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

  1. Index vector (index): Initial unit index vector. By default, vector is randomly drawn from normal distribution and standardized.

  2. Link function: link_lengthscale is initial scalar length-scale parameter. (default: 0.1) link_amp is initial scalar amplitude parameter. (default: 1)

  3. Error variance (sigma2): Initial scalar error variance. (default: 1)

References

Choi, T., Shi, J. Q., & Wang, B. (2011). A Gaussian process regression approach to a single-index model. Journal of Nonparametric Statistics, 23(1), 21-36.

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 <- gpSphere(y ~ ., method = "EB", data = simdata,
                 niter = 1000, nburnin = 100)

# Split version
models <- gpSphere_setup(y ~ ., method = "EB", data = simdata)
Ccompile <- compileModelAndMCMC(models)
nimSampler <- get_sampler(Ccompile)
initList <- getInit(models)
mcmc.out <- runMCMC(nimSampler, niter = 1000, nburnin = 100, thin = 1,
                    nchains = 1, setSeed = TRUE, inits = initList,
                    summary = TRUE, samplesAsCodaMCMC = TRUE)
fit2 <- as_bsim(models, mcmc.out)
# The estimates computed by MAP - standard error of the esitmate is meaningless.
summary(fit2)
# }

Run the code above in your browser using DataLab