BayesSIM (version 1.0.0)

gpFisher: Bayesian Single-Index Regression with Gaussian Process Link and von Mises-Fisher 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 with von Mises-Fisher prior, and the link \(f(\cdot)\) is represented with Gaussian process.

Usage

gpFisher(
  formula,
  data,
  prior = NULL,
  init = NULL,
  monitors = NULL,
  niter = 10000,
  nburnin = 1000,
  thin = 1,
  nchain = 1,
  setSeed = FALSE
)
gpFisher_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 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 All sampling parameters' samplers are 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.

1. Index vector: von Mises--Fisher prior for the projection vector \(\theta\) (index). index_direction is a unit direction vector of the von Mises--Fisher distribution. By default, the estimated vector from projection pursuit regression is assigned. index_dispersion is the positive concentration parameter. By default, 150 is assigned.

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): An 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 \(\theta\). By default, the vector is sampled from the von Mises--Fisher prior.

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

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

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.

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.

Examples

# \donttest{
set.seed(123)
N <- 60; p <- 2
x1 <- runif(N, -3, 5)
x2 <- runif(N, -3, 5)
beta1 <- 0.45; beta2 <- sqrt(1-beta1^2)
sigma <- sqrt(0.0036)
xlin <- x1*beta1 + x2*beta2
eta <- 0.1*xlin + sin(0.5*xlin)^2
y <- rnorm(N, eta, sigma)
x <- matrix(c(x1, x2), ncol = 2)
simdata <- data.frame(x = x, y = y)
colnames(simdata) <- c("X1", "X2", "y")

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

# Split version
models <- gpFisher_setup(y ~ ., data = simdata, nchain = 1)
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)
summary(fit2)
# }