Package {bqmm}

bqmm: Bayesian Multilevel Quantile Regression

Description

Fits Bayesian mixed-effects (multilevel) quantile regression models using the asymmetric Laplace working likelihood and Stan, with an lme4-style formula interface, one or several quantiles per call, optional LKJ-correlated random effects, post-hoc non-crossing rearrangement, and the Yang, Wang and He (2016) posterior-variance correction for valid fixed-effect inference.

Author(s)

Maintainer: Kailas Venkitasubramanian kailasv@gmail.com [copyright holder]

References

Yu, K. and Moyeed, R. A. (2001). Bayesian quantile regression. Statistics & Probability Letters, 54(3), 437-447.

Geraci, M. and Bottai, M. (2014). Linear quantile mixed models. Statistics and Computing, 24(3), 461-479.

Yang, Y., Wang, H. J. and He, X. (2016). Posterior inference in Bayesian quantile regression with asymmetric Laplace likelihood. International Statistical Review, 84(3), 327-344.

See Also

Useful links:

• https://github.com/kvenkita/bqmm

• https://kvenkita.github.io/bqmm/

• Report bugs at https://github.com/kvenkita/bqmm/issues

Random-effect standard deviations and correlations

Description

Random-effect standard deviations and correlations

Usage

## S3 method for class 'bqmm'
VarCorr(x, sigma = 1, ...)

Arguments

Value

A named numeric vector of posterior-median random-effect standard deviations (with a posterior-median correlation matrix attached as the "correlation" attribute for unstructured models), or NULL if the model has no random effects.

Examples

fit <- bqmm(distance ~ age + (1 | Subject), data = nlme::Orthodont,
            tau = 0.5, chains = 1, iter = 300, refresh = 0, seed = 1)
VarCorr(fit)

The asymmetric Laplace family for quantile regression

Description

A lightweight family object describing the asymmetric Laplace distribution (ALD) working likelihood used by bqmm(). It mirrors the role of a stats::family object but is intentionally minimal in this release.

Usage

ald(link = "identity")

Arguments

Value

An object of class "bqmm_family".

Examples

ald()

Coerce a bqmm fit to a matrix of posterior draws

Description

Coerce a bqmm fit to a matrix of posterior draws

Usage

## S3 method for class 'bqmm'
as.matrix(x, ...)

Arguments

Value

A matrix of posterior draws (draws in rows, parameters in columns), using the raw Stan parameter names.

Examples

fit <- bqmm(distance ~ age + (1 | Subject), data = nlme::Orthodont,
            tau = 0.5, chains = 1, iter = 300, refresh = 0, seed = 1)
dim(as.matrix(fit))

Convert a bqmm fit to a posterior draws object

Description

Convert a bqmm fit to a posterior draws object

Usage

## S3 method for class 'bqmm'
as_draws(x, ...)

Arguments

Value

A draws_array (from the posterior package) with tidy variable names: ⁠b_<name>⁠ for fixed effects, ⁠sd_<component>⁠ for random-effect SDs, and sigma.

Examples

fit <- bqmm(distance ~ age + (1 | Subject), data = nlme::Orthodont,
            tau = 0.5, chains = 1, iter = 300, refresh = 0, seed = 1)
as_draws(fit)

Bayesian multilevel quantile regression

Description

Fits a Bayesian mixed-effects quantile regression model using the asymmetric Laplace working likelihood and Stan. The interface follows lme4: random effects are written inline in the formula, e.g. y ~ x + (1 + x | group), and nested or crossed grouping factors are both supported.

Usage

bqmm(
  formula,
  data,
  tau = 0.5,
  family = ald(),
  prior = NULL,
  cov = c("diagonal", "unstructured"),
  adjust = TRUE,
  prior_only = FALSE,
  chains = 4,
  iter = 2000,
  warmup = floor(iter/2),
  cores = getOption("mc.cores", 1L),
  seed = NULL,
  control = list(adapt_delta = 0.95),
  ...
)

Arguments

Details

One or several quantiles may be requested through tau. A scalar returns a single bqmm fit; a vector fits each quantile independently and returns a bqmm_multi container.

Value

A bqmm object (single tau) or a bqmm_multi object (vector tau).

Examples

# A minimal fit; raise chains/iter for real analyses.
fit <- bqmm(distance ~ age + (1 | Subject), data = nlme::Orthodont,
            tau = 0.5, chains = 1, iter = 300, refresh = 0, seed = 1)
summary(fit)

Assemble the Stan data list for the correlated-RE model

Description

Builds the data for inst/stan/bqmm_corr.stan (a single grouping factor with M correlated coefficients). Requires the parsed formula to have exactly one random-effects term.

Usage

bqmm_corr_standata(parsed, tau, prior, prior_only = FALSE)

Arguments

Value

A named list for rstan::sampling() with the correlated-RE model.

Resolve data-scaled default priors

Description

Fills in any NULL elements of a bqmm_prior() using simple, robust scales derived from the response and design matrix. The intent is to keep sigma identified and the sampler away from divergences, not to be informative.

Usage

bqmm_default_priors(prior, y, K)

Arguments

Value

A fully-specified bqmm_prior with numeric hyperparameters.

Parse an lme4-style mixed-model formula

Description

Thin wrapper around lme4::glFormula() that extracts everything bqmm needs: the fixed-effect design matrix X, the random-effect design matrix Z (dense), the response y, and a mapping from each column of Z to a variance component. Reusing lme4's parser means nested and crossed random effects are handled for free.

Usage

bqmm_parse_formula(formula, data, na.action = stats::na.omit, contrasts = NULL)

Arguments

Value

A list with elements:

y

numeric response vector.

X

fixed-effect design matrix (N x K).

Z

random-effect design matrix (N x Q), dense.

sd_map

integer vector (length Q) mapping each Z column to a variance component in 1:G.

re_components

character labels for the G variance components.

re_terms

per-term metadata: grouping factor, coefficient names, number of levels.

cnms, flist, Gp

the raw lme4 random-effect structures.

fixed_names

column names of X.

Priors for a Bayesian quantile mixed model

Description

Builds the list of prior hyperparameters passed to Stan. Defaults are weakly informative and scaled to the data (see bqmm_default_priors()); any element supplied here overrides the default.

Usage

bqmm_prior(
  beta_mean = 0,
  beta_sd = NULL,
  sigma_scale = NULL,
  re_scale = NULL,
  lkj = 2
)

Arguments

Value

An object of class "bqmm_prior".

Examples

bqmm_prior(beta_sd = 5)

Posterior-median random effects as a Q-vector aligned with the columns of Z

Description

Handles both the diagonal model (b drawn as an S x Q matrix) and the correlated model (b drawn as an S x M x L array). For the correlated model the column order (level - 1) * M + coef matches the columns of Z.

Usage

bqmm_ranef_vector(object)

Residuals at the posterior-median fit (fixed + random part)

Description

Residuals at the posterior-median fit (fixed + random part)

Usage

bqmm_residuals_median(object)

Draw posterior samples for one quantile

Description

Thin wrapper around rstan::sampling() that targets the compiled bqmm Stan model (stanmodels$bqmm, generated by rstantools; see SETUP.md) and surfaces basic convergence checks.

Usage

bqmm_sample(
  standata,
  model = "bqmm",
  chains = 4,
  iter = 2000,
  warmup = floor(iter/2),
  cores = 1L,
  seed = NULL,
  control = list(adapt_delta = 0.95),
  ...
)

Arguments

Value

A stanfit object.

Assemble the Stan data list for one quantile

Description

Turns a parsed formula (from bqmm_parse_formula()) and resolved priors into the named list consumed by inst/stan/bqmm.stan.

Usage

bqmm_standata(parsed, tau, prior, prior_only = FALSE)

Arguments

Value

A named list suitable for rstan::sampling().

Map random-effect columns to variance components

Description

lme4 orders the columns of t(Zt) block by term (delimited by Gp), and within a term block the coefficients vary fastest within each level. This function turns that layout into an integer sd_map assigning each column to a ⁠(term, coefficient)⁠ variance component, plus human-readable labels.

Usage

build_sd_map(cnms, Gp)

Arguments

Value

List with sd_map, components (labels) and terms (metadata).

Warn on convergence problems

Description

Warn on convergence problems

Usage

check_convergence(fit)

Arguments

Extract model coefficients

Description

Alias for fixef(); returns the posterior-median fixed effects.

Usage

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

Arguments

Value

A named numeric vector of posterior-median fixed-effect coefficients.

Examples

fit <- bqmm(distance ~ age + (1 | Subject), data = nlme::Orthodont,
            tau = 0.5, chains = 1, iter = 300, refresh = 0, seed = 1)
coef(fit)

Coefficient-versus-tau matrix for a bqmm_multi fit

Description

Coefficient-versus-tau matrix for a bqmm_multi fit

Usage

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

Arguments

Value

A tau-by-coefficient matrix of posterior-median fixed effects, with one row per quantile.

Examples

fit <- bqmm(distance ~ age + (1 | Subject), data = nlme::Orthodont,
            tau = c(0.25, 0.75), chains = 1, iter = 250,
            refresh = 0, seed = 1)
coef(fit)

Numeric core of the Infinitesimal Jackknife variance

Description

Pure-numeric and Stan-free, so it can be unit tested directly.

Usage

compute_ij(beta_draws, loglik_draws, groups = NULL)

Arguments

Value

A symmetric K x K covariance matrix.

Yang-Wang-He multiplicative variance correction

Description

Computes V = Sigma_post %*% G %*% Sigma_post, the mixed-model form of the YWH correction (see ywh_adjust()). The ALD working-likelihood score for observation i is ⁠s_i = (1/sigma) x_i (tau - 1{resid_i < 0})⁠, so the meat is ⁠G = (1/sigma^2) sum_i x_i x_i' psi_i^2⁠ (independence) or ⁠G = (1/sigma^2) sum_g (sum_{i in g} x_i psi_i)(...)'⁠ (cluster-robust), with ⁠psi_i = tau - 1{resid_i < 0}⁠.

Usage

compute_ywh_multiplicative(Sigma_post, X, resid, sigma, tau, groups = NULL)

Arguments

Value

A symmetric K x K covariance matrix.

Numeric core of the Yang-Wang-He sandwich

Description

Pure-numeric and Stan-free, so it can be unit tested directly.

Usage

compute_ywh_sandwich(
  X,
  resid,
  tau,
  Sigma_post = NULL,
  groups = NULL,
  bandwidth = NULL
)

Arguments

Value

List with vcov, D0, D1, bandwidth (the bandwidth actually used, which may have been grown to keep the bread full rank).

Confidence (credible) intervals for the fixed effects

Description

Wald-type intervals built from the posterior-median estimates and the (optionally misspecification-corrected) fixed-effect covariance.

Usage

## S3 method for class 'bqmm'
confint(
  object,
  parm,
  level = 0.95,
  adjusted = TRUE,
  method = c("ywh", "ij"),
  cluster = TRUE,
  ...
)

Arguments

Value

A matrix with one row per coefficient and lower/upper interval columns.

Examples

fit <- bqmm(distance ~ age + (1 | Subject), data = nlme::Orthodont,
            tau = 0.5, chains = 1, iter = 300, refresh = 0, seed = 1)
confint(fit)

Density of the asymmetric Laplace distribution

Description

Density of the asymmetric Laplace distribution

Usage

dald(x, mu = 0, sigma = 1, tau = 0.5, log = FALSE)

Arguments

Value

Numeric vector of (log) density values.

Linear predictor (conditional tau-quantile) at the posterior median

Description

Linear predictor (conditional tau-quantile) at the posterior median

Usage

## S3 method for class 'bqmm'
fitted(object, ...)

Arguments

Value

Numeric vector of fitted conditional quantiles.

Examples

fit <- bqmm(distance ~ age + (1 | Subject), data = nlme::Orthodont,
            tau = 0.5, chains = 1, iter = 300, refresh = 0, seed = 1)
head(fitted(fit))

Posterior-median fixed effects

Description

Posterior-median fixed effects

Usage

## S3 method for class 'bqmm'
fixef(object, ...)

Arguments

Value

A named numeric vector of posterior-median fixed-effect coefficients.

Examples

fit <- bqmm(distance ~ age + (1 | Subject), data = nlme::Orthodont,
            tau = 0.5, chains = 1, iter = 300, refresh = 0, seed = 1)
fixef(fit)

Extract fixed-effect posterior draws

Description

Extract fixed-effect posterior draws

Usage

get_fixef_draws(object)

Hall-Sheather bandwidth for quantile sparsity estimation

Description

Hall-Sheather bandwidth for quantile sparsity estimation

Usage

hall_sheather_bandwidth(n, tau, alpha = 0.05)

Arguments

Value

A positive bandwidth.

Cluster index for the cluster-robust IJ

Description

Cluster index for the cluster-robust IJ

Usage

ij_group_index(object)

Infinitesimal Jackknife (IJ) standard errors

Description

Computes Infinitesimal Jackknife variance estimates for the fixed effects (Giordano & Broderick, 2023; Ji, Lee & Rabe-Hesketh, 2024) from a single MCMC run. The IJ influence of observation i is I_i = n * cov_post(beta, loglik_i) — the posterior covariance between the fixed-effect draws and the per-observation (conditional) log-likelihood draws — and the variance estimate is

Usage

ij_vcov(object, cluster = TRUE)

Arguments

Details

V_IJ = (1 / (n (n - 1))) sum_i (I_i - Ibar)(I_i - Ibar)',

with a cluster-robust version that aggregates influences within cluster j as ⁠I_j = (J / n) sum_{i in j} I_i⁠ and replaces n by the number of clusters J.

Value

A K x K covariance matrix for the fixed effects.

Caveat for hierarchical models

The influences use the conditional per-observation log-likelihood (given the random effects), whereas the IJ of Ji, Lee & Rabe-Hesketh (2024) is derived for a marginal model. For coefficients identified by within-cluster variation (e.g. slopes) the conditional IJ agrees well with the Yang-Wang-He sandwich. For coefficients identified by between-cluster variation (the intercept of a random-intercept model) it can under-estimate the variance — up-weighting a cluster mostly shifts that cluster's random effect rather than the fixed effect — and it is noisier; how much depends on the random-effect-to-noise ratio. For valid fixed-effect inference in the hierarchical model fitted by bqmm, prefer the default Yang-Wang-He sandwich (ywh_adjust()); method = "ij" is provided for benchmarking.

Check whether quantile predictions are monotone in tau

Description

Check whether quantile predictions are monotone in tau

Usage

is_noncrossing(preds)

Arguments

Value

Logical scalar: TRUE if every row is non-decreasing.

Pointwise log-likelihood draws

Description

Pointwise log-likelihood draws

Usage

## S3 method for class 'bqmm'
log_lik(object, ...)

Arguments

Value

An S x N matrix of pointwise log-likelihood values, suitable for use with the loo package.

Examples

fit <- bqmm(distance ~ age + (1 | Subject), data = nlme::Orthodont,
            tau = 0.5, chains = 1, iter = 300, refresh = 0, seed = 1)
dim(log_lik(fit))

Construct a bqmm object

Description

Construct a bqmm object

Usage

new_bqmm(stanfit, parsed, tau, prior, family, cov = "diagonal", call, formula)

Construct a bqmm_multi container

Description

Holds a list of independent bqmm() fits, one per quantile, and presents them jointly through S3 methods (e.g. coefficient-versus-tau paths).

Usage

new_bqmm_multi(fits, parsed, formula, call)

Arguments

Value

A bqmm_multi object.

Number of observations used in the fit

Description

Number of observations used in the fit

Usage

## S3 method for class 'bqmm'
nobs(object, ...)

Arguments

Value

An integer, the number of observations.

Examples

fit <- bqmm(distance ~ age + (1 | Subject), data = nlme::Orthodont,
            tau = 0.5, chains = 1, iter = 300, refresh = 0, seed = 1)
nobs(fit)

Plot a bqmm fit

Description

Default plot shows fixed-effect posterior intervals. With bayesplot installed, richer MCMC plots are available via as_draws().

Usage

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

Arguments

Value

Invisibly, x.

Examples

fit <- bqmm(distance ~ age + (1 | Subject), data = nlme::Orthodont,
            tau = 0.5, chains = 1, iter = 300, refresh = 0, seed = 1)
plot(fit)

Plot coefficient-versus-tau paths for a bqmm_multi fit

Description

Plot coefficient-versus-tau paths for a bqmm_multi fit

Usage

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

Arguments

Value

Invisibly, x.

Examples

fit <- bqmm(distance ~ age + (1 | Subject), data = nlme::Orthodont,
            tau = c(0.25, 0.75), chains = 1, iter = 250,
            refresh = 0, seed = 1)
plot(fit)

Draws of the expected response (conditional tau-quantile)

Description

Draws of the expected response (conditional tau-quantile)

Usage

## S3 method for class 'bqmm'
posterior_epred(object, ...)

Arguments

Value

An S x N matrix of posterior draws of the linear predictor, where S is the number of posterior draws and N the number of observations.

Examples

fit <- bqmm(distance ~ age + (1 | Subject), data = nlme::Orthodont,
            tau = 0.5, chains = 1, iter = 300, refresh = 0, seed = 1)
dim(posterior_epred(fit))

Draws from the posterior predictive distribution

Description

Draws from the posterior predictive distribution

Usage

## S3 method for class 'bqmm'
posterior_predict(object, ...)

Arguments

Value

An S x N matrix of posterior predictive draws of the response, where S is the number of posterior draws and N the number of observations.

Examples

fit <- bqmm(distance ~ age + (1 | Subject), data = nlme::Orthodont,
            tau = 0.5, chains = 1, iter = 300, refresh = 0, seed = 1)
dim(posterior_predict(fit))

Predictions from a bqmm fit

Description

Predictions from a bqmm fit

Usage

## S3 method for class 'bqmm'
predict(
  object,
  newdata = NULL,
  re.form = NULL,
  noncrossing = c("none", "rearrange"),
  ...
)

Arguments

Value

Numeric vector of predicted conditional quantiles.

Examples

fit <- bqmm(distance ~ age + (1 | Subject), data = nlme::Orthodont,
            tau = 0.5, chains = 1, iter = 300, refresh = 0, seed = 1)
head(predict(fit, re.form = NA))

Random generation from the asymmetric Laplace distribution

Description

Uses the normal-exponential location-scale mixture representation of the ALD (Kozumi and Kobayashi, 2011).

Usage

rald(n, mu = 0, sigma = 1, tau = 0.5)

Arguments

Value

Numeric vector of length n.

Posterior-median random effects

Description

Posterior-median random effects

Usage

## S3 method for class 'bqmm'
ranef(object, ...)

Arguments

Value

A numeric vector of posterior-median random effects aligned with the columns of the random-effects design matrix Z, or NULL if the model has no random effects.

Examples

fit <- bqmm(distance ~ age + (1 | Subject), data = nlme::Orthodont,
            tau = 0.5, chains = 1, iter = 300, refresh = 0, seed = 1)
ranef(fit)

Rearrange fitted quantiles to remove crossing

Description

Post-hoc monotonisation of estimated quantile curves (Chernozhukov, Fernandez-Val and Galichon, 2010): at any covariate point, sorting the fitted values across quantile levels into increasing order yields a valid, non-crossing set of quantiles and never increases estimation error. This is the v0.1 non-crossing strategy; joint constrained estimation is deferred.

Usage

rearrange_quantiles(preds)

Arguments

Value

A matrix of the same shape with each row sorted increasingly.

Examples

m <- rbind(c(1, 0.5, 2), c(0, 1, 0.8))   # some crossings
rearrange_quantiles(m)

Objects exported from other packages

Description

These objects are imported from other packages. Follow the links below to see their documentation.

lme4

fixef, ranef, VarCorr

posterior

as_draws

rstantools

log_lik, posterior_epred, posterior_predict

Residuals from a bqmm fit

Description

Residuals from a bqmm fit

Usage

## S3 method for class 'bqmm'
residuals(object, ...)

Arguments

Value

Numeric vector of response-minus-fitted residuals.

Examples

fit <- bqmm(distance ~ age + (1 | Subject), data = nlme::Orthodont,
            tau = 0.5, chains = 1, iter = 300, refresh = 0, seed = 1)
head(residuals(fit))

Asymmetric Laplace check (pinball) loss

Description

Asymmetric Laplace check (pinball) loss

Usage

rho_tau(u, tau)

Arguments

Value

Numeric vector of loss values u * (tau - (u < 0)).

Summarize a bqmm fit

Description

Produces a fixed-effect coefficient table (estimate, standard error and interval) together with random-effect standard deviations.

Usage

## S3 method for class 'bqmm'
summary(
  object,
  level = 0.95,
  adjusted = TRUE,
  method = c("ywh", "ij"),
  cluster = TRUE,
  ...
)

Arguments

Value

An object of class summary.bqmm.

Examples

fit <- bqmm(distance ~ age + (1 | Subject), data = nlme::Orthodont,
            tau = 0.5, chains = 1, iter = 300, refresh = 0, seed = 1)
summary(fit)

Summarize a bqmm_multi fit

Description

Summarize a bqmm_multi fit

Usage

## S3 method for class 'bqmm_multi'
summary(object, ...)

Arguments

Value

A list of summary.bqmm objects, one per quantile.

Examples

fit <- bqmm(distance ~ age + (1 | Subject), data = nlme::Orthodont,
            tau = c(0.25, 0.75), chains = 1, iter = 250,
            refresh = 0, seed = 1)
summary(fit)

Map raw Stan parameter names to tidy bqmm names

Description

Map raw Stan parameter names to tidy bqmm names

Usage

tidy_param_names(nms, object)

Variance-covariance of the fixed effects

Description

Variance-covariance of the fixed effects

Usage

## S3 method for class 'bqmm'
vcov(object, adjusted = TRUE, method = c("ywh", "ij"), cluster = TRUE, ...)

Arguments

Value

A K x K covariance matrix for the fixed effects.

Examples

fit <- bqmm(distance ~ age + (1 | Subject), data = nlme::Orthodont,
            tau = 0.5, chains = 1, iter = 300, refresh = 0, seed = 1)
vcov(fit)

Yang-Wang-He posterior-variance correction (mixed-model form)

Description

The asymmetric Laplace likelihood is a working (misspecified) likelihood, so the naive MCMC posterior covariance of the fixed effects is not the asymptotic variance of the quantile-regression estimator. Yang, Wang and He (2016) correct this with a multiplicative sandwich that re-uses the posterior covariance as the "bread":

Usage

ywh_adjust(object, meat = c("cluster", "independence"))

Arguments

Details

V_adj = Sigma_post %% G %% Sigma_post,

where Sigma_post is the posterior covariance of the fixed effects and G is the meat (variance of the ALD working-likelihood score). For the mixed model this is the right object to correct: Sigma_post already encodes the multilevel pooling, so the correction retains the random-effect contribution to fixed-effect uncertainty while fixing the misspecified ALD scale. Under correct specification G approximately equals Sigma_post^{-1} and the correction reduces to ~Sigma_post. See compute_ywh_multiplicative().

The pure Koenker-Powell sandwich (⁠[compute_ywh_sandwich()]⁠, valid for fixed-effect quantile regression) was found by simulation to under-cover the fixed effects of a mixed model, because it is computed on residuals with the random effects removed and therefore drops the between-cluster variance. The multiplicative form here covers at or slightly above the nominal level across homoscedastic and heteroscedastic designs (see tools/bakeoff.R). Validity is claimed for the fixed-effect block only; variance components keep their model-based posterior.

Value

A list with the adjusted fixed-effect covariance vcov, the naive posterior covariance vcov_naive, the meat G, and sigma.

Cluster index for the cluster-robust meat

Description

Uses the first (outermost) grouping factor when one is available.

Usage

ywh_group_index(object)