| Title: | Adaptive Bayesian Quantile Regression for Count Data |
|---|---|
| Description: | Implements Bayesian quantile regression for count data using the jittering technique for discrete data smoothing and an asymmetric Laplace distribution likelihood. Supports adaptive variable selection via a random-bridge penalty with a beta prior on the power parameter, as well as fixed-bridge and Lasso penalties. Utilizes Markov chain Monte Carlo with Gibbs sampling and adaptive Metropolis-Hastings algorithms for posterior inference, provides Gelman-Rubin convergence diagnostics, and predicts conditional quantiles for count responses. Methodology and applications are based on the following key references: Luo, Zhou, Hu, and Li (2026, Journal of Mathematics, 2026:1543166, <doi:10.1155/jom/1543166>), Koenker and Bassett (1978, Econometrica, 46, 33-50, <doi:10.2307/1913643>), Machado and Santos Silva (2005, Journal of the American Statistical Association, 100, 1226-1237, <doi:10.1198/016214505000000330>), Yu and Moyeed (2001, Statistics and Probability Letters, 54, 437-447, <doi:10.1016/S0167-7152(01)00124-9>), Polson, Scott, and Windle (2014, Journal of the Royal Statistical Society Series B, 76, 713-733, <doi:10.1111/rssb.12042>), Park and Casella (2008, Journal of the American Statistical Association, 103, 681-686, <doi:10.1198/016214508000000337>), and Roberts and Rosenthal (2009, Journal of Computational and Graphical Statistics, 18, 349-367, <doi:10.1198/jcgs.2009.06134>). |
| Authors: | Shikhar Tyagi [aut, cre] (ORCID: <https://orcid.org/0000-0003-1606-0844>), Arvind Pandey [aut], Bhupendra Singh [aut], Vrijesh Tripathi [aut] |
| Maintainer: | Shikhar Tyagi <[email protected]> |
| License: | GPL (>= 3) |
| Version: | 0.1.0 |
| Built: | 2026-10-05 10:49:45 UTC |
| Source: | https://github.com/cran/BayesQRCount |
Performs Bayesian quantile regression for count data using the jittering-based asymmetric Laplace distribution (ALD) likelihood with random-bridge, fixed-bridge, or Lasso penalty priors. Implements the methodology of Luo, Zhou, Hu, and Li (2026).
bqr_count( y, X, tau = 0.5, method = c("random_bridge", "fixed_bridge", "lasso"), xi_fixed = 0.5, N_iter = 20000, burn_in = 10000, thin = 10, n_chains = 3, lambda_prior = c(a = 1, b = 1), xi_prior = c(e = 1, f = 1), delta = 0.05, sigma0_sq = 0.01, seed = NULL )bqr_count( y, X, tau = 0.5, method = c("random_bridge", "fixed_bridge", "lasso"), xi_fixed = 0.5, N_iter = 20000, burn_in = 10000, thin = 10, n_chains = 3, lambda_prior = c(a = 1, b = 1), xi_prior = c(e = 1, f = 1), delta = 0.05, sigma0_sq = 0.01, seed = NULL )
y |
Integer vector of non-negative count response data. |
X |
Numeric matrix of covariates (n x p). An intercept column is added internally; do not include one. |
tau |
Numeric, quantile level in (0, 1). Default is 0.5. |
method |
Character string specifying the penalty method:
|
xi_fixed |
Numeric, power parameter for fixed-bridge penalty (default 0.5).
Used only when |
N_iter |
Integer, total number of MCMC iterations (default 20000). |
burn_in |
Integer, number of burn-in iterations to discard (default 10000). |
thin |
Integer, thinning interval (default 10). |
n_chains |
Integer, number of MCMC chains to run (default 3). |
lambda_prior |
Numeric vector of length 2, hyperparameters (a, b) for lambda ~ Gamma(a, b). Default is c(a = 1, b = 1). |
xi_prior |
Numeric vector of length 2, hyperparameters (e, f) for xi ~ Beta(e, f). Default is c(e = 1, f = 1). Used only for random_bridge. |
delta |
Numeric, probability for fixed (non-adaptive) component in the Metropolis-Hastings mixture proposal for xi. Default is 0.05. |
sigma0_sq |
Numeric, fixed proposal variance for the non-adaptive M-H component. Default is 0.01. |
seed |
Integer or NULL, random seed for reproducibility. |
An object of class "bqr_count", a list containing:
Array of posterior beta samples (n_saved x p x n_chains).
Matrix of posterior sigma samples (n_saved x n_chains).
Matrix of posterior lambda samples (n_saved x n_chains).
Matrix of posterior xi samples (n_saved x n_chains).
Data frame with posterior mean, median, SD, and 95 percent HPD intervals for all parameters.
Named numeric vector of posterior mean coefficients.
Numeric vector of fitted quantile values.
Logical vector indicating selected variables (HPD excludes 0).
The quantile level used.
The penalty method used.
Number of chains.
The response data.
The original covariate matrix (without intercept).
The full design matrix (with intercept).
The matched function call.
set.seed(42) n <- 100 X <- matrix(rnorm(n * 2), n, 2) y <- rpois(n, exp(0.5 + 0.3 * X[, 1] - 0.2 * X[, 2])) fit <- bqr_count(y, X, tau = 0.5, method = "random_bridge", N_iter = 500, burn_in = 200, thin = 2, n_chains = 1) print(fit)set.seed(42) n <- 100 X <- matrix(rnorm(n * 2), n, 2) y <- rpois(n, exp(0.5 + 0.3 * X[, 1] - 0.2 * X[, 2])) fit <- bqr_count(y, X, tau = 0.5, method = "random_bridge", N_iter = 500, burn_in = 200, thin = 2, n_chains = 1) print(fit)
Extract coefficients from bqr_count
## S3 method for class 'bqr_count' coef(object, ...)## S3 method for class 'bqr_count' coef(object, ...)
object |
A |
... |
Additional arguments. |
A numeric vector of coefficients.
set.seed(42) n <- 80 x1 <- rnorm(n) x2 <- rnorm(n) X <- cbind(x1, x2) y <- rpois(n, exp(0.5 + 0.8 * x1)) fit <- bqr_count(y, X, tau = 0.5, method = "lasso", N_iter = 500, burn_in = 200, thin = 2, n_chains = 1) coef(fit)set.seed(42) n <- 80 x1 <- rnorm(n) x2 <- rnorm(n) X <- cbind(x1, x2) y <- rpois(n, exp(0.5 + 0.8 * x1)) fit <- bqr_count(y, X, tau = 0.5, method = "lasso", N_iter = 500, burn_in = 200, thin = 2, n_chains = 1) coef(fit)
Compares two or more fitted "bqr_count" models by RMSE and MAE.
compare_models(...)compare_models(...)
... |
Two or more objects of class |
A data frame with columns: Model, Method, Tau, RMSE, MAE.
set.seed(42) n <- 80 x1 <- rnorm(n) x2 <- rnorm(n) X <- cbind(x1, x2) y <- rpois(n, exp(0.5 + 0.8 * x1)) fit1 <- bqr_count(y, X, tau = 0.5, method = "lasso", N_iter = 500, burn_in = 200, thin = 2, n_chains = 1) fit2 <- bqr_count(y, X, tau = 0.5, method = "random_bridge", N_iter = 500, burn_in = 200, thin = 2, n_chains = 1) compare_models(Lasso = fit1, RandomBridge = fit2)set.seed(42) n <- 80 x1 <- rnorm(n) x2 <- rnorm(n) X <- cbind(x1, x2) y <- rpois(n, exp(0.5 + 0.8 * x1)) fit1 <- bqr_count(y, X, tau = 0.5, method = "lasso", N_iter = 500, burn_in = 200, thin = 2, n_chains = 1) fit2 <- bqr_count(y, X, tau = 0.5, method = "random_bridge", N_iter = 500, burn_in = 200, thin = 2, n_chains = 1) compare_models(Lasso = fit1, RandomBridge = fit2)
Plots kernel density estimates of posterior distributions for selected parameters.
densplot(object, pars = NULL, ...)densplot(object, pars = NULL, ...)
object |
An object of class |
pars |
Character vector of parameter names to plot. Default plots all beta coefficients. |
... |
Additional graphical parameters. |
Invisible NULL. Called for side effects (plotting).
set.seed(42) n <- 80 x1 <- rnorm(n) x2 <- rnorm(n) X <- cbind(x1, x2) y <- rpois(n, exp(0.5 + 0.8 * x1)) fit <- bqr_count(y, X, tau = 0.5, method = "lasso", N_iter = 500, burn_in = 200, thin = 2, n_chains = 1) densplot(fit)set.seed(42) n <- 80 x1 <- rnorm(n) x2 <- rnorm(n) X <- cbind(x1, x2) y <- rpois(n, exp(0.5 + 0.8 * x1)) fit <- bqr_count(y, X, tau = 0.5, method = "lasso", N_iter = 500, burn_in = 200, thin = 2, n_chains = 1) densplot(fit)
Extract fitted values from bqr_count
## S3 method for class 'bqr_count' fitted(object, ...)## S3 method for class 'bqr_count' fitted(object, ...)
object |
A |
... |
Additional arguments. |
A numeric vector of fitted quantiles.
set.seed(42) n <- 80 x1 <- rnorm(n) x2 <- rnorm(n) X <- cbind(x1, x2) y <- rpois(n, exp(0.5 + 0.8 * x1)) fit <- bqr_count(y, X, tau = 0.5, method = "lasso", N_iter = 500, burn_in = 200, thin = 2, n_chains = 1) fitted(fit)set.seed(42) n <- 80 x1 <- rnorm(n) x2 <- rnorm(n) X <- cbind(x1, x2) y <- rpois(n, exp(0.5 + 0.8 * x1)) fit <- bqr_count(y, X, tau = 0.5, method = "lasso", N_iter = 500, burn_in = 200, thin = 2, n_chains = 1) fitted(fit)
Computes the potential scale reduction factor (R-hat) for all model parameters from a multi-chain Bayesian quantile regression fit.
gelman_rubin(object)gelman_rubin(object)
object |
An object of class |
A data frame with columns Parameter and Rhat.
set.seed(42) n <- 80 x1 <- rnorm(n) x2 <- rnorm(n) X <- cbind(x1, x2) y <- rpois(n, exp(0.5 + 0.8 * x1)) fit <- bqr_count(y, X, tau = 0.5, method = "lasso", N_iter = 500, burn_in = 200, thin = 2, n_chains = 2) gelman_rubin(fit)set.seed(42) n <- 80 x1 <- rnorm(n) x2 <- rnorm(n) X <- cbind(x1, x2) y <- rpois(n, exp(0.5 + 0.8 * x1)) fit <- bqr_count(y, X, tau = 0.5, method = "lasso", N_iter = 500, burn_in = 200, thin = 2, n_chains = 2) gelman_rubin(fit)
Plot coefficient effects across different quantiles
plot_effects(fits, var_names = NULL, ...)plot_effects(fits, var_names = NULL, ...)
fits |
A named list of |
var_names |
Optional character vector of variable names to plot. |
... |
Additional arguments passed to |
Invisible NULL.
set.seed(42) n <- 80 x1 <- rnorm(n) x2 <- rnorm(n) X <- cbind(x1, x2) y <- rpois(n, exp(0.5 + 0.8 * x1)) fit1 <- bqr_count(y, X, tau = 0.25, method = "lasso", N_iter = 500, burn_in = 200, thin = 2, n_chains = 1) fit2 <- bqr_count(y, X, tau = 0.5, method = "lasso", N_iter = 500, burn_in = 200, thin = 2, n_chains = 1) fit3 <- bqr_count(y, X, tau = 0.75, method = "lasso", N_iter = 500, burn_in = 200, thin = 2, n_chains = 1) plot_effects(list(Q25 = fit1, Q50 = fit2, Q75 = fit3))set.seed(42) n <- 80 x1 <- rnorm(n) x2 <- rnorm(n) X <- cbind(x1, x2) y <- rpois(n, exp(0.5 + 0.8 * x1)) fit1 <- bqr_count(y, X, tau = 0.25, method = "lasso", N_iter = 500, burn_in = 200, thin = 2, n_chains = 1) fit2 <- bqr_count(y, X, tau = 0.5, method = "lasso", N_iter = 500, burn_in = 200, thin = 2, n_chains = 1) fit3 <- bqr_count(y, X, tau = 0.75, method = "lasso", N_iter = 500, burn_in = 200, thin = 2, n_chains = 1) plot_effects(list(Q25 = fit1, Q50 = fit2, Q75 = fit3))
Plot diagnostics for bqr_count
## S3 method for class 'bqr_count' plot(x, which = 1, ...)## S3 method for class 'bqr_count' plot(x, which = 1, ...)
x |
A |
which |
Integer from 1 to 3 specifying the plot type: 1 = density, 2 = trace, 3 = residuals. |
... |
Additional arguments passed to specific plot functions. |
Invisible NULL.
set.seed(42) n <- 80 x1 <- rnorm(n) x2 <- rnorm(n) X <- cbind(x1, x2) y <- rpois(n, exp(0.5 + 0.8 * x1)) fit <- bqr_count(y, X, tau = 0.5, method = "lasso", N_iter = 500, burn_in = 200, thin = 2, n_chains = 1) plot(fit, which = 3)set.seed(42) n <- 80 x1 <- rnorm(n) x2 <- rnorm(n) X <- cbind(x1, x2) y <- rpois(n, exp(0.5 + 0.8 * x1)) fit <- bqr_count(y, X, tau = 0.5, method = "lasso", N_iter = 500, burn_in = 200, thin = 2, n_chains = 1) plot(fit, which = 3)
Computes predicted conditional quantiles for new or training count data using the fitted model posterior.
## S3 method for class 'bqr_count' predict( object, newdata = NULL, type = c("quantile", "mean"), interval = c("none", "credible"), level = 0.95, ... )## S3 method for class 'bqr_count' predict( object, newdata = NULL, type = c("quantile", "mean"), interval = c("none", "credible"), level = 0.95, ... )
object |
An object of class |
newdata |
Optional numeric matrix of new covariates (without intercept).
If |
type |
Character, type of prediction: |
interval |
Character, |
level |
Numeric, credible interval level (default 0.95). |
... |
Additional arguments (currently unused). |
A numeric vector of predicted values (if interval = "none"),
or a data frame with columns fit, lower, upper
(if interval = "credible").
set.seed(42) n <- 80 x1 <- rnorm(n) x2 <- rnorm(n) X <- cbind(x1, x2) y <- rpois(n, exp(0.5 + 0.8 * x1)) fit <- bqr_count(y, X, tau = 0.5, method = "lasso", N_iter = 500, burn_in = 200, thin = 2, n_chains = 1) predict(fit)set.seed(42) n <- 80 x1 <- rnorm(n) x2 <- rnorm(n) X <- cbind(x1, x2) y <- rpois(n, exp(0.5 + 0.8 * x1)) fit <- bqr_count(y, X, tau = 0.5, method = "lasso", N_iter = 500, burn_in = 200, thin = 2, n_chains = 1) predict(fit)
Print method for bqr_count objects
## S3 method for class 'bqr_count' print(x, ...)## S3 method for class 'bqr_count' print(x, ...)
x |
A |
... |
Additional arguments passed to |
The original object invisibly.
set.seed(42) n <- 80 x1 <- rnorm(n) x2 <- rnorm(n) X <- cbind(x1, x2) y <- rpois(n, exp(0.5 + 0.8 * x1)) fit <- bqr_count(y, X, tau = 0.5, method = "lasso", N_iter = 500, burn_in = 200, thin = 2, n_chains = 1) print(fit)set.seed(42) n <- 80 x1 <- rnorm(n) x2 <- rnorm(n) X <- cbind(x1, x2) y <- rpois(n, exp(0.5 + 0.8 * x1)) fit <- bqr_count(y, X, tau = 0.5, method = "lasso", N_iter = 500, burn_in = 200, thin = 2, n_chains = 1) print(fit)
Print method for summary.bqr_count
## S3 method for class 'summary.bqr_count' print(x, ...)## S3 method for class 'summary.bqr_count' print(x, ...)
x |
A |
... |
Additional arguments passed to |
The original object invisibly.
set.seed(42) n <- 80 x1 <- rnorm(n) x2 <- rnorm(n) X <- cbind(x1, x2) y <- rpois(n, exp(0.5 + 0.8 * x1)) fit <- bqr_count(y, X, tau = 0.5, method = "lasso", N_iter = 500, burn_in = 200, thin = 2, n_chains = 1) s <- summary(fit) print(s)set.seed(42) n <- 80 x1 <- rnorm(n) x2 <- rnorm(n) X <- cbind(x1, x2) y <- rpois(n, exp(0.5 + 0.8 * x1)) fit <- bqr_count(y, X, tau = 0.5, method = "lasso", N_iter = 500, burn_in = 200, thin = 2, n_chains = 1) s <- summary(fit) print(s)
Extract residuals from bqr_count
## S3 method for class 'bqr_count' residuals(object, ...)## S3 method for class 'bqr_count' residuals(object, ...)
object |
A |
... |
Additional arguments. |
A numeric vector of residuals.
set.seed(42) n <- 80 x1 <- rnorm(n) x2 <- rnorm(n) X <- cbind(x1, x2) y <- rpois(n, exp(0.5 + 0.8 * x1)) fit <- bqr_count(y, X, tau = 0.5, method = "lasso", N_iter = 500, burn_in = 200, thin = 2, n_chains = 1) residuals(fit)set.seed(42) n <- 80 x1 <- rnorm(n) x2 <- rnorm(n) X <- cbind(x1, x2) y <- rpois(n, exp(0.5 + 0.8 * x1)) fit <- bqr_count(y, X, tau = 0.5, method = "lasso", N_iter = 500, burn_in = 200, thin = 2, n_chains = 1) residuals(fit)
Computes F1-score, precision, recall, MSE, MAE, and average bias given the true regression coefficients. Useful for simulation studies.
selection_metrics(object, true_beta, level = 0.95)selection_metrics(object, true_beta, level = 0.95)
object |
An object of class |
true_beta |
Numeric vector of true regression coefficients (excluding intercept). Must match the number of covariates. |
level |
Numeric, HPD interval level (default 0.95). |
A data frame with columns: F1, Precision, Recall, MSE, MAE, Bias.
set.seed(42) n <- 80 x1 <- rnorm(n) x2 <- rnorm(n) X <- cbind(x1, x2) y <- rpois(n, exp(0.5 + 0.8 * x1)) fit <- bqr_count(y, X, tau = 0.5, method = "lasso", N_iter = 500, burn_in = 200, thin = 2, n_chains = 1) selection_metrics(fit, true_beta = c(0.8, 0))set.seed(42) n <- 80 x1 <- rnorm(n) x2 <- rnorm(n) X <- cbind(x1, x2) y <- rpois(n, exp(0.5 + 0.8 * x1)) fit <- bqr_count(y, X, tau = 0.5, method = "lasso", N_iter = 500, burn_in = 200, thin = 2, n_chains = 1) selection_metrics(fit, true_beta = c(0.8, 0))
A simulated dataset of 200 observations with 5 covariates and a Poisson count response, designed for demonstrating Bayesian quantile regression with variable selection. The true model has three non-zero coefficients (x1, x2, x5) and two zero coefficients (x3, x4).
data(sim_count_data)data(sim_count_data)
A data frame with 200 rows and 6 columns:
Integer, count response variable generated from Poisson distribution.
Numeric, covariate with true coefficient 0.8.
Numeric, covariate with true coefficient -0.5.
Numeric, noise covariate with true coefficient 0.
Numeric, noise covariate with true coefficient 0.
Numeric, covariate with true coefficient 0.3.
The data were generated using:
mu = exp(0.5 + 0.8*x1 - 0.5*x2 + 0*x3 + 0*x4 + 0.3*x5),
y ~ Poisson(mu), with covariates drawn independently from
standard normal distributions. The true intercept is 0.5. The
sparsity level is 0.4 (2 out of 5 covariates are zero).
Simulated data for package demonstration.
data(sim_count_data) head(sim_count_data) table(sim_count_data$y)data(sim_count_data) head(sim_count_data) table(sim_count_data$y)
Summary method for bqr_count objects
## S3 method for class 'bqr_count' summary(object, ...)## S3 method for class 'bqr_count' summary(object, ...)
object |
A |
... |
Additional arguments. |
A summary.bqr_count object.
set.seed(42) n <- 80 x1 <- rnorm(n) x2 <- rnorm(n) X <- cbind(x1, x2) y <- rpois(n, exp(0.5 + 0.8 * x1)) fit <- bqr_count(y, X, tau = 0.5, method = "lasso", N_iter = 500, burn_in = 200, thin = 2, n_chains = 1) summary(fit)set.seed(42) n <- 80 x1 <- rnorm(n) x2 <- rnorm(n) X <- cbind(x1, x2) y <- rpois(n, exp(0.5 + 0.8 * x1)) fit <- bqr_count(y, X, tau = 0.5, method = "lasso", N_iter = 500, burn_in = 200, thin = 2, n_chains = 1) summary(fit)
Produces trace plots of MCMC samples for selected parameters.
traceplot(object, pars = NULL, ...)traceplot(object, pars = NULL, ...)
object |
An object of class |
pars |
Character vector of parameter names to plot. Default plots all beta coefficients. |
... |
Additional graphical parameters passed to |
Invisible NULL. Called for side effects (plotting).
set.seed(42) n <- 80 x1 <- rnorm(n) x2 <- rnorm(n) X <- cbind(x1, x2) y <- rpois(n, exp(0.5 + 0.8 * x1)) fit <- bqr_count(y, X, tau = 0.5, method = "lasso", N_iter = 500, burn_in = 200, thin = 2, n_chains = 1) traceplot(fit)set.seed(42) n <- 80 x1 <- rnorm(n) x2 <- rnorm(n) X <- cbind(x1, x2) y <- rpois(n, exp(0.5 + 0.8 * x1)) fit <- bqr_count(y, X, tau = 0.5, method = "lasso", N_iter = 500, burn_in = 200, thin = 2, n_chains = 1) traceplot(fit)
Extracts variable selection results from a fitted model based on highest posterior density (HPD) intervals.
variable_select(object, level = 0.95)variable_select(object, level = 0.95)
object |
An object of class |
level |
Numeric, HPD interval level for selection (default 0.95). |
A list with components:
Integer vector of indices of selected covariates (relative to the original X matrix, excluding intercept).
Character vector of names of selected covariates.
Named numeric vector of posterior mean coefficients for selected covariates.
Integer, number of selected covariates.
set.seed(42) n <- 80 x1 <- rnorm(n) x2 <- rnorm(n) X <- cbind(x1, x2) y <- rpois(n, exp(0.5 + 0.8 * x1)) fit <- bqr_count(y, X, tau = 0.5, method = "lasso", N_iter = 500, burn_in = 200, thin = 2, n_chains = 1) variable_select(fit)set.seed(42) n <- 80 x1 <- rnorm(n) x2 <- rnorm(n) X <- cbind(x1, x2) y <- rpois(n, exp(0.5 + 0.8 * x1)) fit <- bqr_count(y, X, tau = 0.5, method = "lasso", N_iter = 500, burn_in = 200, thin = 2, n_chains = 1) variable_select(fit)