Package {bbqr}

bbqr: Bayesian Quantile Regression with Lasso and Adaptive Lasso

Description

Markov chain Monte Carlo samplers for Bayesian quantile regression, based on the asymmetric Laplace distribution and the location-scale mixture representation of Kozumi and Kobayashi (2011) doi:10.1080/00949655.2010.496117. A binary response and an observed continuous response are both supported, each with three penalty layers behind one interface: no penalty, following Benoit and Van den Poel (2012) doi:10.1002/jae.1216; the Bayesian lasso, following Benoit, Al-Hamzawi and Yu (2013) doi:10.1007/s00180-013-0439-0; and the Bayesian adaptive lasso of Rubio Garcia (2023) https://soar.wichita.edu/entities/publication/a2f86232-4704-4ec2-b685-751e7b04ec42. In the binary family each is available as published and in a corrected form, the default, in which every improper prior component is replaced by a proper one so that the posterior exists unconditionally; the continuous family ships the corrected form only. The continuous adaptive-lasso layer at its default reproduces the penalty of Alhamzawi, Yu and Benoit (2012) doi:10.1177/1471082X1101200304. A binary threshold model identifies the coefficient vector only up to a positive scale, so the binary samplers expose the identification anchor as an explicit argument, allowing fixing the scale of the error distribution, fixing a single coefficient, and constraining the norm of the coefficient vector to be compared directly; an observed response identifies the scale, so the continuous samplers have no anchor and draw it every sweep. The MCMC cores are written in Fortran and called from R.

Two response families

The package fits the same three penalty layers to two kinds of response, through two parallel entry points.

bbqr()

A binary response, observed only as 0/1 through a threshold. The latent response is drawn each sweep. Because the threshold identifies the coefficient vector only up to a positive scale, these fits take an anchor that says how the scale is fixed. bbqr_alasso(), bbqr_lasso() and bbqr_none() are the direct forms.

cbqr()

An observed continuous response. There is no latent response draw and no anchor: the data identify the scale, so sigma is sampled every sweep. The location-scale latents z and s are still drawn. These fits instead take q, the exponent on sigma in the local-scale prior rate, which the binary model cannot see because its sigma is held fixed. cbqr_alasso(), cbqr_lasso() and cbqr_none() are the direct forms.

The two families do not share a prior constructor. prior() builds a bbqr.prior object for the binary fitters; the continuous fitters take a plain named list, because their hierarchy has different components (no anchor, a live sigma, and the q exponent) and different reference defaults. Passing a prior() object to cbqr() is an error, and says so.

Method sources

The binary samplers implement Rubio Garcia (2023) for the adaptive lasso, and Benoit, Al-Hamzawi and Yu (2013) and Benoit and Van den Poel (2012) for the lasso and unpenalised layers. On the continuous side the published specification each layer reproduces is a different one: at its default q = 1 the adaptive lasso is the penalty of Alhamzawi, Yu and Benoit (2012), the lasso at q = 0 is Benoit, Al-Hamzawi and Yu (2013) carried over to an observed response, and the unpenalised layer is the sampler of Kozumi and Kobayashi (2011). See citation("bbqr").

In the binary family, each of the three penalties is shipped in two forms, selected by prior(penalty, model = ): "v3", the hierarchy as published, and "v5", the default, in which every improper prior component has been replaced by a proper one so that the posterior exists unconditionally. The published adaptive-lasso and lasso hierarchies have no posterior, for the reason given in prior(); they are kept so that earlier results can be reproduced, and the print method labels them. The continuous family ships one hierarchy, the proper one: there is no model argument and no published continuous specification being reproduced alongside it.

Getting started

bbqr() is the entry point for a 0/1 response; see its help page for the three penalties, the identification anchors and the derivation versions. cbqr() is the entry point for a continuous one. prior() sets hyperparameters and selects the hierarchy for the binary fitters only; the continuous fitters take a plain named list, documented at cbqr().

Relationship to bayesQR

The Fortran RNG wrapper (src/wrapper.c) and the general layout of the package follow Benoit's bayesQR, which is GPL (>= 2) licensed and is acknowledged in the Authors@R field of DESCRIPTION. The six MCMC kernels shipped here – three binary, three continuous – were written for this package.

Author(s)

Maintainer: Fernando Rubio Garcia j332v755@wichita.edu (Wichita State University)

Other contributors:

• Dries F. Benoit (Author of 'bayesQR', from which the Fortran RNG wrapper and package layout are derived) [contributor, copyright holder]

• Rahim Al-Hamzawi [contributor]

• Keming Yu [contributor]

• Dirk Van den Poel [contributor]

See Also

Useful links:

• https://github.com/fernandorubiogarcia/bbqr

• Report bugs at https://github.com/fernandorubiogarcia/bbqr/issues

The tau-calibrated reference intercept prior

Description

At the sigma = 1 anchor and with ⁠x'beta = 0⁠, the model reduces to Pr(y = 1) = G_tau(beta0) with G_tau the distribution function of -E for E ~ ALD(0, 1, tau). By the probability integral transform, G_tau(beta0) is exactly Uniform(0, 1) when beta0 ~ ALD(0, 1, 1 - tau), so the reference normal prior is that law's moment match:

Usage

.bbqr_b0_reference(quantile)

Arguments

Details

⁠mean = -theta = -(1 - 2 tau) / (tau (1 - tau))⁠ ⁠var = (1 - 2 tau + 2 tau^2) / (tau^2 (1 - tau)^2)⁠

The mean is not zero away from the median. Centring the intercept prior at the origin puts the prior-predictive mean of Pr(y = 1) at 0.73 rather than 0.5 at tau = 0.1, so it is a specification error rather than a neutral default. The variance is within 5% of the value that minimises the Kolmogorov-Smirnov distance to Uniform(0, 1) over scalings of itself, so no optimisation is performed: the closed form is used as it stands.

Value

A list with elements mean and var.

Bayesian binary quantile regression

Description

Fits a quantile regression model to a binary response by Markov chain Monte Carlo, using the asymmetric Laplace distribution (ALD) and the location-scale mixture representation of Kozumi and Kobayashi (2011). Three penalty layers are available through the penalty argument, the identification anchor is exposed through anchor, and the prior hierarchy is chosen by model on the prior() object. The defaults – model = "v5" at anchor = "sigma1" – sample a posterior that is proper unconditionally and identified by construction.

Usage

bbqr(
  formula,
  data = NULL,
  quantile = 0.5,
  penalty = c("alasso", "lasso", "none"),
  ndraw = 10000,
  burn = NULL,
  keep = 1,
  prior = NULL,
  anchor = NULL,
  sigma_fixed = 1,
  standardize = FALSE,
  use_boundaries = FALSE,
  clip = NULL,
  sampler = NULL,
  beta_init = NULL,
  beta0_init = 0
)

Arguments

Value

For a single quantile, an object of class "bbqr": a list with components

beta

Matrix of retained slope draws, one column per covariate.

beta0

Vector of retained intercept draws.

sigma

Vector of retained ALD inverse-scale draws. Constant at sigma_fixed under anchor = "sigma1"; under the norm anchors it varies, because the rescaling step absorbs the scale into it.

accept

Final Metropolis acceptance rate for the penalty hyperparameter, or NA when the sampler has no Metropolis step.

penalty, anchor, quantile

The settings used.

prior

The prior actually applied.

model

Which hierarchy was sampled: "v3", "v4" or "v5".

derived

TRUE when the chain targets a posterior one of the derivations writes down; FALSE when a projection anchor was combined with "v4" or "v5".

b0_prior

The intercept prior actually used, as mean, var and prec. Recorded because the "v5" reference depends on quantile, so it cannot be reconstructed from the prior object alone.

Penalty-specific draws are also returned: lambdasq, omega and delta for the adaptive lasso; lambdasq, tauhyper and delta for the lasso. For several quantiles, an object of class "bbqr.list" holding one such object per level.

Derivation version

Each penalty has a published specification and a corrected one, selected by prior(penalty, model = ). "v3" is the published hierarchy: for the adaptive lasso and the lasso its joint posterior does not exist, because the log-uniform hyperprior on the global shrinkage scale (\omega or \tau_h) is not integrable along the ray where the scale, the \lambda_j^2 and the latent s_j vanish together while the likelihood stays bounded below. "v5", the default, replaces every improper component with a proper one – \mathrm{Gamma}(2, 2) on the global scale, and a \tau-calibrated normal prior on the intercept – so that the posterior is proper with no condition on the design, on the number of covariates, or on the anchor. "v4" is the intermediate step for the adaptive lasso. The fitted object records which hierarchy was sampled in model, and whether the chain targets a derived posterior in derived. Under "v4" and "v5" the numerical clamps, argument floors and the historical sampler routine are refused, because the derivations say none of them is part of the target; "v3" keeps them reachable so that earlier results can be reproduced. See prior() for the hierarchies in full.

Choosing a penalty

The three penalties share the same latent-variable augmentation and differ only in what sits on top of the slopes:

"none"

A vague N(0, \sigma^2_\beta) prior on each slope and no shrinkage, following Benoit and Van den Poel (2012).

"lasso"

One shrinkage parameter shared by all slopes, following Benoit, Al-Hamzawi and Yu (2013).

"alasso"

A separate shrinkage parameter per slope, so that coefficients are penalised individually rather than uniformly, following Rubio Garcia (2023).

Identification

A binary threshold model identifies the coefficient vector only up to a positive scale: \mathrm{sign}(x'\beta + \varepsilon) is unchanged if \beta and the error scale are multiplied by the same positive constant. Some anchor must therefore be imposed, and the choice is not innocuous. The available anchors are:

"sigma1"

Fix the ALD inverse-scale at sigma_fixed (default 1) and never sample it. The default for every penalty. It closes the scale ray outright, so the posterior summaries of \beta are on a stated scale, and it is the convention under which the \tau-calibrated intercept prior of model = "v5" is exact. In the simulation study behind this package its point estimates matched "free" while its credible intervals were markedly narrower, because a sampled \sigma lets the unidentified scale inflate them.

"free"

Impose nothing beyond the prior: sample \sigma from its full conditional and let the prior alone settle the scale. This is the sampled-\sigma convention of the adaptive-lasso and lasso derivations, and it is a derived posterior under every model; it is simply not identified by the likelihood, so ⁠||beta||⁠ is a prior artefact. Not available for penalty = "none": with neither a penalty nor an anchor the likelihood is flat along the scale ray, so nothing identifies the model and only the vague prior stops the chain drifting.

"beta1"

Hold the first slope at 1 and never sample it. The first column of the design matrix must be a covariate with a genuinely non-zero coefficient for this to be sensible.

"norm1"

Rescale so that \|(\beta_0, \beta)\| = 1 after every sweep.

"normslopes"

Rescale so that \|\beta\| = 1 over the slopes alone, the convention used in the maximum-score literature.

The two norm anchors share one projection step and differ only in how the scaling constant is chosen. In both cases \beta_0 and the ALD inverse-scale \sigma are rescaled alongside the slopes. At the observed-data level the model is invariant under (\beta_0, \beta, \sigma) \mapsto (c\beta_0, c\beta, \sigma/c). These projections are identification experiments and are not steps in any of the derivations.

Note, however, that this invariance does not extend to the penalised posterior: the shrinkage layer is stated on a fixed scale, so repeatedly renormalising the coefficients moves them relative to the penalty. Anchors that constrain the norm of the coefficient vector should therefore be expected to behave differently from "sigma1", which leaves the coefficients where the sampler put them. Under model = "v4" or "v5" the three projection anchors ("beta1", "norm1", "normslopes") warn and set derived = FALSE on the fit: the projection is an exact reparameterisation of the likelihood but not of the prior, so the chain targets no written-down posterior. They remain available because the comparison is the point of exposing the anchor at all.

Scaling

The prior variances and the penalty are stated on the scale of the design matrix as supplied. If one covariate is measured in units a hundred times larger than another, a shared penalty shrinks them by very different relative amounts, and the recovered direction can be badly wrong. Setting standardize = TRUE centres and scales the covariates before sampling and maps the draws back afterwards; the back-transform preserves the linear predictor exactly, but the posterior itself differs because the penalty now acts on comparable coefficients. Unless the covariates are already on a common scale, prefer standardize = TRUE.

References

Benoit, D. F. and Van den Poel, D. (2012). Binary quantile regression: a Bayesian approach based on the asymmetric Laplace distribution. Journal of Applied Econometrics, 27(7), 1174–1188. doi:10.1002/jae.1216

Benoit, D. F., Al-Hamzawi, R. and Yu, K. (2013). Bayesian lasso binary quantile regression. Computational Statistics, 28(6), 2861–2873. doi:10.1007/s00180-013-0439-0

Kozumi, H. and Kobayashi, G. (2011). Gibbs sampling methods for Bayesian quantile regression. Journal of Statistical Computation and Simulation, 81(11), 1565–1578. doi:10.1080/00949655.2010.496117

Rubio Garcia, F. (2023). Bayesian adaptive lasso binary quantile regression with hybrid resampling for classification of imbalanced data. M.S. thesis, Wichita State University, Dept. of Mathematics, Statistics, and Physics. https://soar.wichita.edu/entities/publication/a2f86232-4704-4ec2-b685-751e7b04ec42

See Also

prior() for hyperparameters, summary.bbqr() for posterior summaries, plot.bbqr() for trace and density plots, and predict.bbqr() for fitted probabilities. cbqr() fits the same three penalties when the response is observed rather than thresholded.

Examples

set.seed(42)
n <- 200
X <- matrix(rnorm(n * 4), n, 4, dimnames = list(NULL, paste0("x", 1:4)))
y <- as.numeric(X %*% c(1.5, -1, 0, 0) + rnorm(n) > 0)
dat <- data.frame(y = y, X)

fit <- bbqr(y ~ x1 + x2 + x3 + x4, data = dat, quantile = 0.5,
            ndraw = 1000, burn = 200)
fit
summary(fit)

# Compare the three penalties at the same quantile
fits <- lapply(c("none", "lasso", "alasso"), function(p)
  bbqr(y ~ x1 + x2 + x3 + x4, data = dat, penalty = p,
       ndraw = 1000, burn = 200))
sapply(fits, coef)

# The published adaptive-lasso hierarchy, whose posterior does not exist,
# is still available for reproducing earlier results
old <- bbqr(y ~ x1 + x2 + x3 + x4, data = dat, ndraw = 1000, burn = 200,
            prior = prior("alasso", model = "v3"))
c(fit$model, old$model)

Adaptive-lasso binary quantile regression

Description

Low-level sampler for binary quantile regression with an adaptive-lasso penalty, in which each slope receives its own shrinkage parameter. Most users should call bbqr() with penalty = "alasso" instead; this function is exported for direct access to the sampler. With the defaults – model = "v5" on the prior, anchor = "sigma1", no clamps – every transition implements the corresponding full conditional in the V5 appendix: the ALD inverse-scale is held at sigma_fixed, the omega draw is ⁠Gamma(alpha_omega + k delta, beta_omega + sum(1 / lambdasq))⁠, the intercept carries its \tau-calibrated normal prior, and delta is updated by a random-walk Metropolis step on the log scale under a proper Gamma(alpha_delta, beta_delta) prior. That last prior is required rather than optional: see prior(). anchor = "free" samples the inverse-scale from its full conditional instead, which the appendix also covers.

Usage

bbqr_alasso(
  formula,
  data = NULL,
  quantile = 0.5,
  ndraw = 10000,
  burn = NULL,
  keep = 1,
  prior = NULL,
  anchor = c("sigma1", "free", "beta1", "norm1", "normslopes"),
  sigma_fixed = 1,
  standardize = FALSE,
  use_boundaries = FALSE,
  clip = NULL,
  sampler = c("gig_half", "v4_invgaus"),
  kappa_lam = 0,
  beta_init = NULL,
  beta0_init = 0
)

Arguments

Value

An object of class "bbqr", as documented in bbqr(), with two extra components: clip, the resolved logical(18) of switches, and clip_hits, the number of times each site's clamp actually bound. clip_hits is counted whether or not the switch was on, so a fit with every clamp off still reports how often each clamp would have fired. Denominators differ by site. ystar, z, v4_zarg and v4_zbound are out of ndraw * n; s, lambdasq, v4_sarg and v4_sbound out of ndraw * p; mh_sd out of burn; delta out of the accepted Metropolis proposals; omega and v4_rate_om out of ndraw. Multiply by the number of guards in a group where it holds more than one: v4_zarg has two, v4_sarg three, each rate group two. v4_sd is out of n alone, since v4 floored that standard deviation in the initialisation sweep only and left the main-loop assignment bare. sigma and v4_rate_sig stay at zero under anchor = "sigma1", where nothing ever reaches them.

Derivation version

Which of the three adaptive-lasso hierarchies is sampled is set by model on the prior() object, not by an argument here. "v4" and "v5" are derivation-faithful: their appendices state that no clamp, no sampler-argument floor, and no exp(-kappa_lam / lambda_j^2) tilt is part of the target, so a fit that names them while enabling one of those devices is refused rather than recorded. Concretely, under model = "v4" or "v5" this function requires clip = "none" (the default), kappa_lam = 0 and sampler = "gig_half". model = "v3" leaves all three switches reachable, because reproducing what the earlier builds actually ran needs them.

One further check is a warning rather than an error. Under model = "v4" with a free sigma, a flat intercept is admissible only if a > 1, and the historical default a = 1 sits exactly on the divergent boundary. That configuration is what Experiments v4 and v6 ran, so it stays available and warns instead of failing.

References

Rubio Garcia, F. (2023). Bayesian adaptive lasso binary quantile regression with hybrid resampling for classification of imbalanced data. M.S. thesis, Wichita State University, Dept. of Mathematics, Statistics, and Physics. https://soar.wichita.edu/entities/publication/a2f86232-4704-4ec2-b685-751e7b04ec42

See Also

bbqr(), and cbqr_alasso() for the same penalty layer fitted to an observed continuous response

Examples

set.seed(1)
n <- 150
X <- matrix(rnorm(n * 3), n, 3, dimnames = list(NULL, c("x1", "x2", "x3")))
y <- as.numeric(X %*% c(1, -1, 0) + rnorm(n) > 0)
d <- data.frame(y = y, X)
fit <- bbqr_alasso(y ~ x1 + x2 + x3, data = d, ndraw = 400, burn = 100)
coef(fit)

Lasso binary quantile regression

Description

Low-level sampler for binary quantile regression with a standard lasso penalty: one shrinkage parameter shared by all slopes. This implements the hierarchical sampler of Benoit, Al-Hamzawi and Yu (2013). Most users should call bbqr() with penalty = "lasso" instead; this function is exported for direct access to the sampler.

Usage

bbqr_lasso(
  formula,
  data = NULL,
  quantile = 0.5,
  ndraw = 10000,
  burn = NULL,
  keep = 1,
  prior = NULL,
  anchor = c("sigma1", "free", "beta1", "norm1", "normslopes"),
  sigma_fixed = 1,
  standardize = FALSE,
  use_boundaries = FALSE,
  beta_init = NULL,
  beta0_init = 0
)

Arguments

Details

The default anchor = "sigma1" holds the ALD inverse-scale at sigma_fixed; anchor = "free" samples it from its full conditional, which is what the published algorithm does. Which hierarchy is sampled – the published one, whose posterior does not exist, or the corrected "v5" default – is set by model on the prior() object.

Value

An object of class "bbqr", as documented in bbqr().

References

Benoit, D. F., Al-Hamzawi, R. and Yu, K. (2013). Bayesian lasso binary quantile regression. Computational Statistics, 28(6), 2861–2873. doi:10.1007/s00180-013-0439-0

See Also

bbqr(), and cbqr_lasso() for the same penalty layer fitted to an observed continuous response

Examples

set.seed(1)
n <- 150
X <- matrix(rnorm(n * 3), n, 3, dimnames = list(NULL, c("x1", "x2", "x3")))
y <- as.numeric(X %*% c(1, -1, 0) + rnorm(n) > 0)
d <- data.frame(y = y, X)
fit <- bbqr_lasso(y ~ x1 + x2 + x3, data = d, ndraw = 400, burn = 100)
coef(fit)

Unpenalised binary quantile regression

Description

Low-level sampler for binary quantile regression with a vague normal prior on the slopes and no shrinkage layer, implementing Benoit and Van den Poel (2012). Most users should call bbqr() with penalty = "none" instead.

Usage

bbqr_none(
  formula,
  data = NULL,
  quantile = 0.5,
  ndraw = 10000,
  burn = NULL,
  keep = 1,
  prior = NULL,
  anchor = c("sigma1", "beta1", "norm1", "normslopes"),
  standardize = FALSE,
  use_boundaries = FALSE,
  beta_init = NULL,
  beta0_init = 0
)

Arguments

Details

The default anchor = "sigma1" reproduces the published algorithm, which holds the ALD scale at 1. "free" is deliberately absent from the choices here, unlike the penalised fitters: with neither a penalty nor an anchor the likelihood is flat along the scale ray, so the model is not identified and only the vague prior stops the chain drifting.

Value

An object of class "bbqr", as documented in bbqr().

References

Benoit, D. F. and Van den Poel, D. (2012). Binary quantile regression: a Bayesian approach based on the asymmetric Laplace distribution. Journal of Applied Econometrics, 27(7), 1174–1188. doi:10.1002/jae.1216

See Also

bbqr(), and cbqr_none() for the same penalty layer fitted to an observed continuous response

Examples

set.seed(1)
n <- 150
X <- matrix(rnorm(n * 3), n, 3, dimnames = list(NULL, c("x1", "x2", "x3")))
y <- as.numeric(X %*% c(1, -1, 0) + rnorm(n) > 0)
d <- data.frame(y = y, X)
fit <- bbqr_none(y ~ x1 + x2 + x3, data = d, ndraw = 400, burn = 100)
coef(fit)

Bayesian quantile regression for a continuous response

Description

User-facing entry point for the three continuous penalty layers. The binary threshold model is fitted by bbqr() instead; the two share no arguments that refer to identification, because a continuous response identifies the scale and a binary one does not.

Usage

cbqr(
  formula,
  data = NULL,
  quantile = 0.5,
  penalty = c("alasso", "lasso", "none"),
  ndraw = 10000,
  burn = NULL,
  keep = 1,
  prior = NULL,
  q = NULL,
  standardize = TRUE,
  use_boundaries = FALSE,
  clip = NULL,
  sampler = NULL,
  beta_init = NULL,
  beta0_init = 0
)

Arguments

Value

A cbqr object, or a cbqr.list when quantile has length > 1. The components are:

family, penalty, q, quantile

what was fitted.

beta, beta0, sigma

the kept draws: a matrix of slopes, and vectors for the intercept and the ALD inverse-scale. Unlike a binary fit, sigma is always a live parameter here.

lambdasq, omega, delta

the shrinkage layer, "alasso" only: the per-coefficient scales, the global scale, and its shape.

lambdasq, tau_hyper, delta

the shrinkage layer, "lasso" only: one global scale, its hyperparameter, and the shape. Absent for "none", which has no shrinkage layer.

accept, acc

Metropolis acceptance rate for delta, and its running trace. NA when there is no Metropolis step.

ndraw, burn, keep, ndraw_kept

the MCMC controls as used.

n, nvar, names, terms, call

the data and the model frame the fit came from.

prior

the resolved prior list, including any defaults filled in.

standardized

the centring and scaling applied to X and y, or NULL. Draws are already mapped back, so this is a record rather than something to undo.

b0_prior

the intercept prior actually used.

clip, clip_hits, sampler

the clamp configuration, how often each site fired, and the latent-draw routine. "alasso" only.

References

Alhamzawi, R., Yu, K. and Benoit, D. F. (2012). Bayesian adaptive lasso quantile regression. Statistical Modelling, 12(3), 279–297. doi:10.1177/1471082X1101200304

Benoit, D. F., Al-Hamzawi, R. and Yu, K. (2013). Bayesian lasso binary quantile regression. Computational Statistics, 28(6), 2861–2873. doi:10.1007/s00180-013-0439-0

Kozumi, H. and Kobayashi, G. (2011). Gibbs sampling methods for Bayesian quantile regression. Journal of Statistical Computation and Simulation, 81(11), 1565–1578. doi:10.1080/00949655.2010.496117

See Also

bbqr() for the binary threshold model; cbqr_alasso(), cbqr_lasso() and cbqr_none() for the direct forms; predict.cbqr(), summary.cbqr() and plot.cbqr() for the fitted object. prior() does not configure these fitters.

Examples

set.seed(1)
d <- data.frame(x1 = rnorm(120), x2 = rnorm(120))
d$y <- 1 + 2 * d$x1 - d$x2 + rnorm(120)
fit <- cbqr(y ~ x1 + x2, d, ndraw = 400, burn = 100)
round(coef(fit), 2)

Continuous quantile regression with the adaptive lasso

Description

Fits the continuous-response adaptive-lasso hierarchy for an observed response. The response is observed, so there is no latent response draw and no identification anchor: sigma is identified by the data and is sampled every sweep.

Usage

cbqr_alasso(
  formula,
  data = NULL,
  quantile = 0.5,
  ndraw = 10000,
  burn = NULL,
  keep = 1,
  prior = NULL,
  q = NULL,
  standardize = TRUE,
  clip = NULL,
  sampler = NULL,
  beta_init = NULL,
  beta0_init = 0
)

Arguments

Value

An object of class cbqr.

References

Alhamzawi, R., Yu, K. and Benoit, D. F. (2012). Bayesian adaptive lasso quantile regression. Statistical Modelling, 12(3), 279–297. doi:10.1177/1471082X1101200304

Benoit, D. F., Al-Hamzawi, R. and Yu, K. (2013). Bayesian lasso binary quantile regression. Computational Statistics, 28(6), 2861–2873. doi:10.1007/s00180-013-0439-0

Kozumi, H. and Kobayashi, G. (2011). Gibbs sampling methods for Bayesian quantile regression. Journal of Statistical Computation and Simulation, 81(11), 1565–1578. doi:10.1080/00949655.2010.496117

See Also

cbqr() for the common interface, bbqr_alasso() for the binary counterpart. prior() configures that one, not this one.

Examples

set.seed(1)
d <- data.frame(x1 = rnorm(120), x2 = rnorm(120))
d$y <- 1 + 2 * d$x1 - d$x2 + rnorm(120)
fit <- cbqr_alasso(y ~ x1 + x2, d, ndraw = 400, burn = 100)
round(coef(fit), 2)

Continuous quantile regression with the lasso

Description

Fits the continuous-response lasso hierarchy for an observed response: a single shared lambda^2 carried as a rate with a Gamma prior and a tau_h hyperlayer. The prior rate here is sigma^q * lambda^2 / 2, with lambda^2 a single global scale multiplying; the adaptive-lasso rate is sigma^q / (2 * lambda_j^2), one scale per coefficient and dividing. The two forms are different by construction, not a typo.

Usage

cbqr_lasso(
  formula,
  data = NULL,
  quantile = 0.5,
  ndraw = 10000,
  burn = NULL,
  keep = 1,
  prior = NULL,
  q = NULL,
  standardize = TRUE,
  use_boundaries = FALSE,
  beta_init = NULL,
  beta0_init = 0
)

Arguments

Details

This is Benoit, Al-Hamzawi & Yu (2013)'s hierarchy, not the adaptive one with k = 1, and the two are not reachable from one another.

Value

An object of class cbqr.

References

Alhamzawi, R., Yu, K. and Benoit, D. F. (2012). Bayesian adaptive lasso quantile regression. Statistical Modelling, 12(3), 279–297. doi:10.1177/1471082X1101200304

Benoit, D. F., Al-Hamzawi, R. and Yu, K. (2013). Bayesian lasso binary quantile regression. Computational Statistics, 28(6), 2861–2873. doi:10.1007/s00180-013-0439-0

Kozumi, H. and Kobayashi, G. (2011). Gibbs sampling methods for Bayesian quantile regression. Journal of Statistical Computation and Simulation, 81(11), 1565–1578. doi:10.1080/00949655.2010.496117

See Also

cbqr() for the common interface, bbqr_lasso() for the binary counterpart.

Examples

set.seed(1)
d <- data.frame(x1 = rnorm(120), x2 = rnorm(120))
d$y <- 1 + 2 * d$x1 - d$x2 + rnorm(120)
fit <- cbqr_lasso(y ~ x1 + x2, d, ndraw = 400, burn = 100)
round(coef(fit), 2)

Continuous quantile regression, unpenalised

Description

Fits the unpenalised hierarchy for an observed response. Unlike its binary counterpart, which holds sigma at the identification anchor and never samples it, this kernel draws sigma every sweep from its Gamma full conditional given the latents, with shape ⁠a + 3n/2⁠.

Usage

cbqr_none(
  formula,
  data = NULL,
  quantile = 0.5,
  ndraw = 10000,
  burn = NULL,
  keep = 1,
  prior = NULL,
  standardize = TRUE,
  use_boundaries = FALSE,
  beta_init = NULL,
  beta0_init = 0
)

Arguments

Value

An object of class cbqr.

References

Alhamzawi, R., Yu, K. and Benoit, D. F. (2012). Bayesian adaptive lasso quantile regression. Statistical Modelling, 12(3), 279–297. doi:10.1177/1471082X1101200304

Benoit, D. F., Al-Hamzawi, R. and Yu, K. (2013). Bayesian lasso binary quantile regression. Computational Statistics, 28(6), 2861–2873. doi:10.1007/s00180-013-0439-0

Kozumi, H. and Kobayashi, G. (2011). Gibbs sampling methods for Bayesian quantile regression. Journal of Statistical Computation and Simulation, 81(11), 1565–1578. doi:10.1080/00949655.2010.496117

See Also

cbqr() for the common interface, bbqr_none() for the binary counterpart.

Examples

set.seed(1)
d <- data.frame(x1 = rnorm(120), x2 = rnorm(120))
d$y <- 1 + 2 * d$x1 - d$x2 + rnorm(120)
fit <- cbqr_none(y ~ x1 + x2, d, ndraw = 400, burn = 100)
round(coef(fit), 2)

Extract posterior mean coefficients

Description

Extract posterior mean coefficients

Usage

## S3 method for class 'bbqr'
coef(object, ...)

## S3 method for class 'bbqr.list'
coef(object, ...)

Arguments

Value

A named numeric vector of posterior means, intercept first. For a "bbqr.list", a matrix with one column per quantile.

Trace and density plots for a fitted binary quantile regression

Description

Trace and density plots for a fitted binary quantile regression

Usage

## S3 method for class 'bbqr'
plot(x, which = NULL, type = c("both", "trace", "density"), ...)

Arguments

Value

x, invisibly. Called for its side effect.

Examples

set.seed(1)
n <- 150
X <- matrix(rnorm(n * 3), n, 3, dimnames = list(NULL, c("x1", "x2", "x3")))
y <- as.numeric(X %*% c(1, -1, 0) + rnorm(n) > 0)
fit <- bbqr(y ~ x1 + x2 + x3, data = data.frame(y = y, X),
            ndraw = 500, burn = 100)
op <- graphics::par(no.readonly = TRUE)
plot(fit, which = "x1")
graphics::par(op)

Trace and density plots for a continuous quantile regression

Description

Same interface as plot.bbqr(), with one addition: sigma is a plottable parameter here. In the binary model it is constant under the default anchor = "sigma1", so there is usually nothing to look at; here it is always a live parameter.

Usage

## S3 method for class 'cbqr'
plot(x, which = NULL, type = c("both", "trace", "density"), ...)

Arguments

Value

x, invisibly. Called for the side effect of drawing.

See Also

plot.bbqr() for the binary counterpart.

Predicted probabilities and indices from a binary quantile regression

Description

Predicted probabilities and indices from a binary quantile regression

Usage

## S3 method for class 'bbqr'
predict(object, newdata = NULL, type = c("response", "link", "class"), ...)

Arguments

Details

Predictions use the posterior mean coefficients. Because the coefficient vector is identified only up to a positive scale, type = "response" depends on the identification anchor, whereas type = "class" does not.

Value

A numeric vector, one entry per row of newdata.

Examples

set.seed(1)
n <- 150
X <- matrix(rnorm(n * 3), n, 3, dimnames = list(NULL, c("x1", "x2", "x3")))
y <- as.numeric(X %*% c(1, -1, 0) + rnorm(n) > 0)
d <- data.frame(y = y, X)
fit <- bbqr(y ~ x1 + x2 + x3, data = d, ndraw = 500, burn = 100)
head(predict(fit, newdata = d, type = "response"))
mean(predict(fit, newdata = d, type = "class") == y)

Predictions from a continuous quantile regression

Description

Returns the fitted conditional quantile of the response, which is the linear predictor. That is a different quantity from what predict.bbqr() returns – P(y = 1 | x) from the ALD distribution function – which on a continuous fit saturates at 1 and is silently meaningless.

Usage

## S3 method for class 'cbqr'
predict(object, newdata = NULL, interval = FALSE, level = 0.95, ...)

Arguments

Value

A numeric vector, or a matrix with columns fit, lower, upper when interval = TRUE.

Prior and hyperparameter settings for the binary fitters

Description

Builds the list of prior hyperparameters and Metropolis tuning constants for one of the three penalty layers.

Usage

prior(
  penalty = c("alasso", "lasso", "none"),
  model = c("v5", "v4", "v3"),
  a = NULL,
  b = NULL,
  alpha_delta = NULL,
  beta_delta = NULL,
  alpha_omega = NULL,
  beta_omega = NULL,
  alpha_tau = NULL,
  beta_tau = NULL,
  b0_mean = NULL,
  b0_var = NULL,
  delta_init = 1,
  omega_init = 1,
  tau_init = 1,
  mh_sd0 = 0.05,
  target_acc = 0.3,
  beta_var = 100
)

Arguments

Details

This object is accepted by bbqr() and the three direct binary fitters, and by those only. The continuous fitters (cbqr() and friends) take a plain named list, not a bbqr.prior: their hierarchy has no anchor, a live \sigma, and the q exponent, and its reference defaults differ. Passing the result of this function to cbqr() fails with “Unknown prior component(s)”; cbqr() documents the components they do accept. The defaults are the "v5" hierarchy, in which every prior component is a probability distribution and the posterior is proper unconditionally; model = "v3" reproduces the specification of the source paper for each penalty instead. In most cases this function only needs to be called to change one or two values.

Value

A list of class "bbqr.prior".

References

Benoit, D. F. and Van den Poel, D. (2012). Binary quantile regression: a Bayesian approach based on the asymmetric Laplace distribution. Journal of Applied Econometrics, 27(7), 1174–1188. doi:10.1002/jae.1216

Benoit, D. F., Al-Hamzawi, R. and Yu, K. (2013). Bayesian lasso binary quantile regression. Computational Statistics, 28(6), 2861–2873. doi:10.1007/s00180-013-0439-0

See Also

bbqr() and the direct binary fitters, which accept this object; cbqr(), which does not.

Examples

# Default adaptive-lasso prior: the fully proper "v5" hierarchy
prior("alasso")

# The published specification, for reproducing earlier results. Its
# posterior does not exist, and the print method says so.
prior("alasso", model = "v3")

# A vaguer prior on the slopes for the unpenalised sampler
prior("none", beta_var = 1000)

Summarise a fitted binary quantile regression

Description

Posterior means, standard deviations and equal-tailed credible intervals for every coefficient. A coefficient whose interval excludes zero is flagged as selected, which is the variable-selection rule used with these samplers.

Usage

## S3 method for class 'bbqr'
summary(object, level = 0.95, ...)

Arguments

Value

An object of class "bbqr.summary".

Examples

set.seed(1)
n <- 150
X <- matrix(rnorm(n * 3), n, 3, dimnames = list(NULL, c("x1", "x2", "x3")))
y <- as.numeric(X %*% c(1, -1, 0) + rnorm(n) > 0)
fit <- bbqr(y ~ x1 + x2 + x3, data = data.frame(y = y, X),
            ndraw = 500, burn = 100)
summary(fit)

Summarise a fitted continuous quantile regression

Description

Posterior means, standard deviations and equal-tailed credible intervals for every coefficient, plus the ALD inverse-scale. A coefficient whose interval excludes zero is flagged.

Usage

## S3 method for class 'cbqr'
summary(object, level = 0.95, ...)

Arguments

Details

Unlike the binary summary, there is no identification caveat: a continuous response identifies the scale, so the coefficients are on the scale of y and are directly interpretable.

Value

An object of class "cbqr.summary".