GBASS: Generalized Bayesian Adaptive Smoothing Splines

Description

Bayesian nonlinear regression under a range of likelihood models using generalized Bayesian adaptive smoothing splines. Robust regression with Student's t likelihoods, quantile regression, and related latent-scale models are included as special cases.

Author(s)

Maintainer: Kellin Rumsey knrumsey@lanl.gov

Construct a prior specification for GBASS

Description

Construct a prior specification for GBASS

Usage

build_prior(
  type = c("GIG", "GBP"),
  p,
  a,
  b,
  lower_bound = NULL,
  prop_sigma = NULL,
  adapt = NULL,
  adapt_delay = NULL,
  adapt_thin = NULL
)

build_GIG(
  p,
  a,
  b,
  lower_bound = NULL,
  prop_sigma = NULL,
  adapt = NULL,
  adapt_delay = NULL,
  adapt_thin = NULL
)

build_GBP(
  p,
  a,
  b,
  lower_bound = NULL,
  prop_sigma = NULL,
  adapt = NULL,
  adapt_delay = NULL,
  adapt_thin = NULL
)

Arguments

Value

A named list containing the prior specification.

Generalized Bayesian MARS

Description

Fits a generalized Bayesian MARS (GBASS) model for nonlinear regression under flexible latent-scale likelihoods.

Usage

gbass(
  X,
  y,
  w_prior = list(type = "GIG", p = 0, a = 0, b = 0),
  v_prior = list(type = "GIG", p = -15, a = 0, b = 30),
  maxInt = 3,
  maxBasis = 1000,
  npart = NULL,
  nmcmc = 10000,
  nburn = 9000,
  thin = 1,
  moveProbs = rep(1/3, 3),
  a_tau = 1/2,
  b_tau = NULL,
  a_lambda = 1,
  b_lambda = 1,
  m_beta = 0,
  s_beta = 0,
  scale = 1,
  Iw0 = rep(1, maxInt),
  Zw0 = rep(1, ncol(X)),
  verbose = TRUE
)

Arguments

Details

The priors for w and v_i must belong to the generalized inverse Gaussian ("GIG") or generalized beta prime ("GBP") families.

Each prior list may contain:

1. type: either "GIG" or "GBP".

2. p, a, b: prior hyperparameters.

3. lower_bound: optional lower bound for the parameter support. For backward compatibility, lb is also accepted.

4. prop_sigma: proposal standard deviation on the log scale for Metropolis updates, when applicable.

5. adapt, adapt_delay, adapt_thin: optional controls for adaptive Metropolis updates when applicable.

For w, prop_sigma and the adaptive Metropolis fields are mainly relevant when w_prior$type = "GBP". In the "GIG" case with nonzero beta, w is sampled using the modified half-normal formulation rather than a random-walk Metropolis step.

Retained draws are taken at iterations nburn + 1, nburn + 1 + thin, .... Thus nburn is interpreted as the number of initial iterations discarded as burn-in.

Value

An object of class "gbass" containing posterior draws and fitted model information.

Note

Current implementation notes:

1. Basis function parameters are stored as lists.

2. Knot locations use continuous uniform proposals.

3. Basis coefficients use a ridge-type prior.

Examples

ff1 <- function(x) 10.391*((x[1]-0.4)*(x[2]-0.6) + 0.36)

n <- 100
p <- 4
X <- matrix(runif(n * p), nrow = n)
y <- apply(X, 1, ff1)
mod <- gbass(X, y, nmcmc = 1000, nburn = 900, thin = 2)

Convert GBASS object to BASS-like object

Description

Converts an object of class gbass to an object with class bass, so that downstream BASS utilities such as Sobol decomposition can be used.

Usage

gbass2bass(gm)

gm2bm(gm)

Arguments

Value

An object with class c("bass", "gbass").

Sobol decomposition for GBASS models

Description

Computes Sobol sensitivity indices for a fitted gbass model by first converting it to a compatible BASS object with gbass2bass(), and then calling BASS::sobol().

Usage

gsobol(object, ...)

Arguments

Details

This is a thin wrapper around BASS::sobol() for GBASS models. Users who want direct access to the converted BASS object can call gbass2bass() explicitly and then use BASS::sobol() themselves.

Value

An object of class "bassSob" returned by BASS::sobol().

See Also

gbass2bass, sobol

Examples

ff1 <- function(x) 10.391*((x[1]-0.4)*(x[2]-0.6) + 0.36)

n <- 100
p <- 4
X <- matrix(runif(n * p), nrow = n)
y <- apply(X, 1, ff1)

mod <- gbass(X, y, nmcmc = 1000, nburn = 900)

# Direct wrapper
sob <- gsobol(mod)

# Equivalent manual conversion
bm <- gbass2bass(mod)
sob2 <- BASS::sobol(bm)

HBASS - Bayesian MARS with Horseshoe or Strawderman-Berger likelihood

Description

Wrapper around gbass() using generalized beta prime latent-scale priors.

Usage

hbass(
  X,
  y,
  w_prior = list(type = "GBP", p = 1, a = 1/2, b = 1/2),
  likelihood = "h",
  ...
)

Arguments

Details

The Horseshoe and Strawderman-Berger likelihoods can be expressed using generalized beta prime latent-scale distributions. This function provides a convenient wrapper for these likelihoods. For additional flexibility, use gbass() directly.

Value

An object of class c("hbass", "gbass") containing posterior draws and fitted model information.

Method-of-moments estimation for the Normal-Wald distribution

Description

Estimates Normal-Wald parameters from sample moments or from a supplied vector of moments.

Usage

nw_est_mom(
  data = NULL,
  stats = NULL,
  mu = NA,
  delta = NA,
  beta = NA,
  alpha = NA,
  triangle = FALSE,
  ...
)

Arguments

Details

This function computes method-of-moments estimates for the Normal-Wald distribution. If data is supplied, sample mean, variance, skewness, and kurtosis are computed automatically. Otherwise, the user may provide these moments directly through stats.

If some parameters are fixed, only the remaining parameters are estimated. When triangle = TRUE, the returned values are the transformed summary quantities \text{steepness} = (1 + |\gamma|)^{-1/2} and a corresponding asymmetry measure.

Value

If triangle = FALSE, a named numeric vector with entries mu, delta, beta, and alpha. If triangle = TRUE, a named numeric vector with entries asymmetry and steepness.

Examples

n <- 500
y <- rgamma(n, 3, 1.5) + rlnorm(n, 1, 0.5)
z <- (y - mean(y)) / sd(y)
skew <- mean(z^3)
kurt <- mean(z^4)

nw_est_mom(stats = c(NA, NA, skew, kurt), mu = 0, delta = 1, triangle = TRUE)
nw_est_mom(stats = c(NA, var(y), skew, kurt), mu = 0, triangle = TRUE)

Calibrate a prior for the Normal-Wald gamma parameter

Description

Selects prior hyperparameters for gamma in the Normal-Wald model using probability statements about the steepness parameter \xi = (1 + \gamma)^{-1/2}.

Usage

nw_gamma_prior(
  q1 = 0.1,
  q2 = 0.9,
  p1 = 0.5,
  p2 = 0.05,
  par0 = NULL,
  lambda = 0
)

Arguments

Details

This function chooses hyperparameters for a prior on gamma by solving a simple calibration problem based on the steepness parameter (1 + \gamma)^{-1/2}. The optimization is carried out with optim using the Nelder-Mead method.

Value

A named numeric vector with entries m_gamma and s_gamma.

Examples

nw_gamma_prior()
nw_gamma_prior(q1 = 0.2, q2 = 0.8, p1 = 0.4, p2 = 0.1)

Plot Normal-Wald shape parameters on the asymmetry-steepness triangle

Description

Visualizes posterior draws of the Normal-Wald shape parameters in terms of asymmetry and steepness. The plotted coordinates are \text{steepness} = (1 + \gamma)^{-1/2} and \text{asymmetry} = \beta (\beta^2 + \gamma^2)^{-1/2} \times \text{steepness}.

Usage

nw_triangle(obj, add = FALSE, details = FALSE, ...)

Arguments

Details

This plot provides a simple geometric summary of the shape of the Normal-Wald likelihood. The upper boundary corresponds to symmetric models, while points away from zero on the horizontal axis indicate asymmetry.

Value

Invisibly returns a data frame with columns asymmetry and steepness.

Examples

# mod is a fitted nwbass model

# Simple example
n <- 200
X <- lhs::maximinLHS(n, 2)
f <- 20 * apply(X, 1, function(x) sin(pi * x[1]) + x[2]^2)
eps <- rgamma(n, 3, 1.5) - 2
y <- f + eps

mod <- nwbass(X, y,
              m_beta=0, s_beta=10,
              m_gamma = nw_gamma_prior()[1], s_gamma = nw_gamma_prior()[2])
nw_triangle(mod)

Generalized Bayesian MARS with a Normal-Wald likelihood

Description

Fits a generalized BMARS model under a Normal-Wald likelihood. This provides flexible nonlinear regression with a unimodal, potentially skewed error model.

Usage

nwbass(
  X,
  y,
  w_prior = list(type = "GIG", p = -0.1, a = 0, b = 0.1),
  maxInt = 3,
  maxBasis = 1000,
  npart = NULL,
  nmcmc = 10000,
  nburn = 9000,
  thin = 1,
  moveProbs = rep(1/3, 3),
  a_tau = 1/2,
  b_tau = NULL,
  a_lambda = 1,
  b_lambda = 1,
  m_beta = 0,
  s_beta = 0,
  lag_beta = 1001,
  m_gamma = 1,
  s_gamma = 0,
  scale = 1,
  Iw0 = rep(1, maxInt),
  Zw0 = rep(1, ncol(X)),
  verbose = TRUE
)

nwbass2(...)

Arguments

Details

The latent local-scale prior is fixed internally to the Normal-Wald form v_i \sim \mathrm{GIG}(-1/2, \gamma^2, 1). Unlike gbass(), nwbass() does not expose a user-specified v_prior.

The w_prior list should contain:

1. type: either "GIG" or "GBP".

2. p, a, b: hyperparameters for the prior.

3. lower_bound: optional lower bound for w. For backward compatibility, lb is also accepted.

4. prop_sigma: proposal standard deviation on the log scale for Metropolis updates of w. This is only used when w is updated by Metropolis-Hastings, such as the "GBP" case.

5. adapt, adapt_delay, adapt_thin: optional controls for adaptive Metropolis updates of w when applicable.

Retained draws are taken at iterations nburn + 1, nburn + 1 + thin, .... Thus nburn is interpreted as the number of initial iterations discarded as burn-in.

Value

An object of class c("nwbass", "gbass") containing posterior draws and fitted model information.

Examples

n <- 200
X <- lhs::maximinLHS(n, 2)
f <- 20 * apply(X, 1, function(x) sin(pi * x[1]) + x[2]^2)
eps <- rgamma(n, 3, 1.5) - 2
y <- f + eps

mod <- nwbass(X, y,
              m_beta=0, s_beta=10,
              m_gamma = nw_gamma_prior()[1], s_gamma = nw_gamma_prior()[2])

Plot method for gbass objects

Description

Plot method for gbass objects

Usage

## S3 method for class 'gbass'
plot(x, ...)

Arguments

Value

No return value, called for its side effects.

Predict method for GBASS objects

Description

Returns posterior draws of either: (i) the linear predictor / mean surface, or (ii) the full posterior predictive distribution.

Usage

## S3 method for class 'gbass'
predict(
  object,
  newdata = NULL,
  mcmc.use = NULL,
  predictive = TRUE,
  bias_correct = FALSE,
  samples = 1,
  ...
)

Arguments

Details

If predictive = FALSE and bias_correct = FALSE, this returns draws of B(x)a.

If predictive = FALSE and bias_correct = TRUE, this returns draws of B(x)a + E(error | posterior draw), i.e. the mean response under the fitted GBASS error model.

If predictive = TRUE, this returns posterior predictive draws by simulating a fresh latent local variance v_new and Gaussian draw for each posterior sample.

For qbass objects, bias_correct = TRUE is usually not what you want, because qbass is typically being used for quantile regression rather than mean regression.

Currently, posterior predictive draws are implemented for GIG-based models. If the fitted object uses a GBP prior for v, predictive = TRUE will stop.

Value

a matrix with rows corresponding to posterior draws and columns corresponding to rows of newdata.

Print a gbass object

Description

Print a gbass object

Usage

## S3 method for class 'gbass'
print(x, ...)

Arguments

Value

The input object, invisibly.

QBASS - Bayesian MARS with an asymmetric Laplace likelihood

Description

Wrapper around gbass() for quantile regression.

Usage

qbass(
  X,
  y,
  q = 0.5,
  w_prior = list(type = "GIG", p = -0.1, a = 0, b = 0.1),
  ...
)

Arguments

Details

Performs quantile regression for quantile q using the asymmetric Laplace representation. For many quantiles, fitting separate models in parallel may be convenient.

Value

An object of class c("qbass", "gbass") containing posterior draws and fitted model information.

Generalized Inverse Gaussian Generator

Description

This function generates samples from the GIG(p, a, b) distribution with density function f(x) = cx^(p-1)exp(-1/2*(a*x + b/x))

Usage

rgig2(p, a, b)

Arguments

Details

A uniformly bounded rejection sample based on Hörmann et. al. (2014). Special cases include Gamma (b=0, p>0), Inverse Gamma (a=0, p<0) and Inverse Gaussian (p=-1/2).

Value

A numeric vector of random draws from the generalized inverse Gaussian distribution.

References

Hörmann, Wolfgang, and Josef Leydold. "Generating generalized inverse Gaussian random variates." Statistics and Computing 24.4 (2014): 547-557.

Examples

x <- rep(NA, 1000)
for(i in 1:1000) x[i] <- rgig2(2, 3, 0.5)
hist(x)

Random generation from the modified half-normal distribution

Description

Generates random samples from the modified half-normal distribution with density proportional to x^{\alpha - 1} \exp(-\beta x^2 + \gamma x) for x > 0.

Usage

rmhn(n = 1, alpha, beta, gamma)

Arguments

Details

This distribution arises in latent-scale updates for the generalized Bayesian adaptive smoothing spline model. The function is intended mainly for internal model fitting, but it may also be useful on its own.

Value

A numeric vector of length n containing random draws from the modified half-normal distribution.

Examples

x <- rmhn(1000, alpha = 3, beta = 1, gamma = 0.5)
hist(x, breaks = 30, freq = FALSE)

TBASS - Bayesian MARS with a Student's t likelihood

Description

Wrapper around gbass() for a Student's t error model.

Usage

tbass(X, y, df = 5, ...)

Arguments

Details

Uses an inverse-gamma latent-scale representation through a GIG prior on the local variance components. If df > 2, the scale is set to (df - 2) / df so that the marginal variance matches w.

Value

An object of class c("tbass", "gbass") containing posterior draws and fitted model information.