| Title: | Bayesian Quantile Regression Models for Complex Survey Data Analysis |
|---|---|
| Description: | Provides Bayesian quantile regression models for complex survey data under informative sampling using survey-weighted estimators. Both single- and multiple-output models are supported. To accelerate computation, all algorithms are implemented in 'C++' using 'Rcpp', 'RcppArmadillo', and 'RcppEigen', and are called from 'R'. See Nascimento and Gonçalves (2024) <doi:10.1093/jssam/smae015> and Nascimento and Gonçalves (2026) <doi:10.1093/jssam/smaf040>. |
| Authors: | Tomás Rodríguez Taborda [aut, cre], Johnatan Cardona Jiménez [aut], Marcus L. Nascimento [aut], Kelly Cristina Mota Gonçalves [aut] |
| Maintainer: | Tomás Rodríguez Taborda <[email protected]> |
| License: | MIT + file LICENSE |
| Version: | 0.4.2 |
| Built: | 2026-10-28 06:04:40 UTC |
| Source: | https://github.com/torodriguezt/bayesqrsurvey |
The dataset comprises 985 observations and 8 variables on children aged 0-60 months from the Central-West region of Brazil. The data are a cleaned subset of the 2006 Brazilian National Demographic Health Survey of Women and Children (NDHS), a complex survey in which sampling units were selected in two stages: primary sampling units (PSUs) were census tracts, and secondary sampling units (SSUs) were households. Tract selection within each stratum was designed to ensure a minimum number of blood samples, given the prevalence of vitamin A deficiency in children. Outlier weights (> 30 kg) and rows with missing values have been removed.
data("Anthro")data("Anthro")
The data frame has the following components:
wgt : weight
hgt : height
wgt_ind : weight-for-age index
hgt_ind : height-for-age index
ruc : rural-urban classification (urban = 1 and rural = 2)
sex : factor with levels 0 (girl) and 1 (boy)
age : age in months
dweight : sampling weights
The 2006 Brazilian National Demographic Health Survey of Women and Children (NDHS) was published by the Brazilian Institute of Geography and Statistics.
bqr.svy implements Bayesian methods for estimating quantile regression models
for complex survey data analysis regarding single (univariate) outputs. To
improve computational efficiency, the Markov Chain Monte Carlo (MCMC) algorithms
are implemented in 'C++'.
bqr.svy( formula, weights = NULL, data = NULL, quantile = 0.5, method = c("ald", "score", "approximate"), prior = NULL, niter = 20000, burnin = 0, thin = 1, verbose = TRUE, estimate_sigma = FALSE, pi_matrix = NULL )bqr.svy( formula, weights = NULL, data = NULL, quantile = 0.5, method = c("ald", "score", "approximate"), prior = NULL, niter = 20000, burnin = 0, thin = 1, verbose = TRUE, estimate_sigma = FALSE, pi_matrix = NULL )
formula |
a symbolic description of the model to be fit. Offset terms are not supported. |
weights |
an optional numerical vector containing the survey weights. If |
data |
an optional data frame containing the variables in the model. |
quantile |
numerical scalar or vector containing quantile(s) of interest (default=0.5). |
method |
one of |
prior |
a |
niter |
number of MCMC draws. |
burnin |
number of initial MCMC draws to be discarded.(default = 0) |
thin |
thinning parameter, i.e., keep every keepth draw (default=1). |
verbose |
logical flag indicating whether to print progress messages (default=TRUE). |
estimate_sigma |
logical flag indicating whether to estimate the scale parameter
when method = "ald" (default=FALSE and |
pi_matrix |
an optional |
The bqr.svy function can estimate three types of models, where the quantile regression coefficients are defined at the super-population level, and their estimators are built upon the survey weights.
"ald" - The asymmetric Laplace distribution as working likelihood.
"score" - A score based likelihood function.
"approximate" - A pseudolikelihood function based on a Gaussian approximation.
An object of class "bqr.svy" with the components listed below.
The components that depend on the quantile (beta, draws,
diagnosis and accept_rate) have one entry per fitted quantile,
named "tau=0.100", "tau=0.500" and so on, also when a single
quantile was fitted.
beta |
Matrix of posterior mean estimates of the regression coefficients, with one row per coefficient and one column per quantile. |
draws |
Named list of posterior draw matrices, one per quantile. The
first columns are the regression coefficients, followed by the ALD scale
when applicable. Its column is named |
diagnosis |
Named list of convergence diagnostics, one per quantile. |
accept_rate |
Named vector of average acceptance rates, one per quantile
( |
quantile |
The quantile(s) fitted. |
weights, n
|
The survey weights supplied and the number of observations. |
warmup, thin, niter
|
MCMC control parameters used during sampling. |
prior |
Prior specification used. |
call, formula, terms, model
|
Model specification details. |
xlevels, contrasts
|
Factor levels and contrasts of the model frame,
used by |
runtime |
Elapsed runtime in seconds. |
method |
Estimation method |
estimate_sigma |
Logical flag indicating whether the scale parameter
|
To obtain the coefficients, fitted values, covariance matrices and related
quantities, use the generic accessor functions documented in
bqr.svy.methods, together with
posterior_interval for credible intervals and
diagnostics for convergence, instead of accessing these
components directly. All of them select quantiles numerically through a
tau argument.
Nascimento, M. L. & Gonçalves, K. C. M. (2024). Bayesian Quantile Regression Models for Complex Survey Data Under Informative Sampling. Journal of Survey Statistics and Methodology, 12(4), 1105-1130. doi:10.1093/jssam/smae015
# Generate population data set.seed(123) N <- 10000 x1_p <- runif(N, -1, 1) x2_p <- runif(N, -1, 1) y_p <- 2 + 1.5 * x1_p - 0.8 * x2_p + rnorm(N) # Generate sample data n <- 500 z_aux <- rnorm(N, mean = 1 + y_p, sd = .5) p_aux <- 1 / (1 + exp(2.5 - 0.5 * z_aux)) s_ind <- sample(1:N, n, replace = FALSE, prob = p_aux) y_s <- y_p[s_ind] x1_s <- x1_p[s_ind] x2_s <- x2_p[s_ind] w <- 1 / p_aux[s_ind] data <- data.frame(y = y_s, x1 = x1_s, x2 = x2_s, w = w) # Basic usage with default method ('ald') and priors (vague) fit1 <- bqr.svy(y ~ x1 + x2, weights = w, data = data) # Specify informative priors prior <- prior( beta_x_mean = c(2, 1.5, -0.8), beta_x_cov = diag(c(0.25, 0.25, 0.25)), sigma_shape = 1, sigma_rate = 1 ) fit2 <- bqr.svy(y ~ x1 + x2, weights = w, data = data, prior = prior) # Specify different methods fit_score <- bqr.svy(y ~ x1 + x2, weights = w, data = data, method = "score") fit_approx <- bqr.svy(y ~ x1 + x2, weights = w, data = data, method = "approximate")# Generate population data set.seed(123) N <- 10000 x1_p <- runif(N, -1, 1) x2_p <- runif(N, -1, 1) y_p <- 2 + 1.5 * x1_p - 0.8 * x2_p + rnorm(N) # Generate sample data n <- 500 z_aux <- rnorm(N, mean = 1 + y_p, sd = .5) p_aux <- 1 / (1 + exp(2.5 - 0.5 * z_aux)) s_ind <- sample(1:N, n, replace = FALSE, prob = p_aux) y_s <- y_p[s_ind] x1_s <- x1_p[s_ind] x2_s <- x2_p[s_ind] w <- 1 / p_aux[s_ind] data <- data.frame(y = y_s, x1 = x1_s, x2 = x2_s, w = w) # Basic usage with default method ('ald') and priors (vague) fit1 <- bqr.svy(y ~ x1 + x2, weights = w, data = data) # Specify informative priors prior <- prior( beta_x_mean = c(2, 1.5, -0.8), beta_x_cov = diag(c(0.25, 0.25, 0.25)), sigma_shape = 1, sigma_rate = 1 ) fit2 <- bqr.svy(y ~ x1 + x2, weights = w, data = data, prior = prior) # Specify different methods fit_score <- bqr.svy(y ~ x1 + x2, weights = w, data = data, method = "score") fit_approx <- bqr.svy(y ~ x1 + x2, weights = w, data = data, method = "approximate")
Generic accessor functions for objects of class "bqr.svy" returned by
bqr.svy, following the conventions of the standard regression
methods in stats.
## S3 method for class 'bqr.svy' coef(object, tau = NULL, ...) ## S3 method for class 'bqr.svy' fitted(object, tau = NULL, ...) ## S3 method for class 'bqr.svy' vcov(object, tau = NULL, ...) ## S3 method for class 'bqr.svy.vcov' print(x, digits = 4, ...) ## S3 method for class 'bqr.svy' nobs(object, ...) ## S3 method for class 'bqr.svy' as.matrix(x, tau = NULL, include_sigma = NULL, ...) ## S3 method for class 'bqr.svy' sigma(object, tau = NULL, ...) ## S3 method for class 'bqr.svy' formula(x, ...) ## S3 method for class 'bqr.svy' terms(x, ...) ## S3 method for class 'bqr.svy' model.matrix(object, ...) ## S3 method for class 'bqr.svy' weights(object, ...) ## S3 method for class 'bqr.svy' update(object, formula., ..., evaluate = TRUE)## S3 method for class 'bqr.svy' coef(object, tau = NULL, ...) ## S3 method for class 'bqr.svy' fitted(object, tau = NULL, ...) ## S3 method for class 'bqr.svy' vcov(object, tau = NULL, ...) ## S3 method for class 'bqr.svy.vcov' print(x, digits = 4, ...) ## S3 method for class 'bqr.svy' nobs(object, ...) ## S3 method for class 'bqr.svy' as.matrix(x, tau = NULL, include_sigma = NULL, ...) ## S3 method for class 'bqr.svy' sigma(object, tau = NULL, ...) ## S3 method for class 'bqr.svy' formula(x, ...) ## S3 method for class 'bqr.svy' terms(x, ...) ## S3 method for class 'bqr.svy' model.matrix(object, ...) ## S3 method for class 'bqr.svy' weights(object, ...) ## S3 method for class 'bqr.svy' update(object, formula., ..., evaluate = TRUE)
object, x
|
an object of class |
tau |
numeric vector of fitted quantiles to extract. |
... |
further arguments passed to or from other methods. |
digits |
integer; number of decimal places used when printing the
covariance matrix returned by |
include_sigma |
logical; whether to keep the scale parameter column in the returned draws. Defaults to whether the model estimated it. |
formula. |
a change to the model formula, in the form used by
|
evaluate |
logical; if |
coef returns the posterior means of the regression coefficients.
fitted evaluates the estimated conditional quantile function
at the observed
covariates.
vcov returns the posterior covariance matrix of the coefficients,
computed from the posterior draws. It is a numeric matrix, stored at full
precision; only its print method rounds, to digits decimal
places.
as.matrix returns the posterior draws themselves. They are the basis
for any further inference, such as posterior probabilities like
, contrasts between quantiles, or passing the draws to
packages such as coda or bayesplot.
sigma returns the scale parameter of the working likelihood. It is the
posterior mean when method = "ald" and estimate_sigma = TRUE,
exactly 1 when "ald" holds it fixed, and NA for the
"score" and "approximate" methods, which have no scale
parameter.
formula returns the formula as supplied, terms its terms
object, and model.matrix the design matrix the fit was built on.
weights returns the survey weights as supplied, without
normalisation. update refits the model with part of the call changed.
These are point summaries, evaluated at the posterior mean. The posterior
credible intervals of the coefficients are computed by
posterior_interval.
coef and fitted return a vector for a single quantile,
otherwise a matrix with one column per quantile. vcov returns a
matrix for a single quantile, otherwise a named list of matrices. Matrix
arithmetic on it behaves as on any numeric matrix; unclass drops the
class that governs printing.
as.matrix returns the draws by parameters matrix of posterior draws
for one quantile. sigma returns a scalar for a single quantile and a
named vector otherwise. nobs returns the number of observations and
weights the vector of survey weights.
bqr.svy, bqr.svy.posterior for the
posterior draws and credible intervals,
predict for the conditional quantile at
new covariate values, diagnostics,
summary
set.seed(123) n <- 300 d <- data.frame(x1 = runif(n, -1, 1), x2 = runif(n, -1, 1), w = runif(n, 1, 3)) d$y <- 2 + 1.5 * d$x1 - 0.8 * d$x2 + rnorm(n) fit <- bqr.svy(y ~ x1 + x2, weights = w, data = d, quantile = c(0.25, 0.5, 0.75), niter = 2000, burnin = 500) coef(fit) # one column per quantile coef(fit, tau = 0.5) # a single named vector vcov(fit, tau = 0.5) head(fitted(fit, tau = 0.5)) nobs(fit) head(weights(fit)) # the survey weights, as supplied # the posterior draws, one row per retained iteration draws <- as.matrix(fit, tau = 0.5) dim(draws) mean(draws[, "x1"] > 0) # P(beta_x1 > 0)set.seed(123) n <- 300 d <- data.frame(x1 = runif(n, -1, 1), x2 = runif(n, -1, 1), w = runif(n, 1, 3)) d$y <- 2 + 1.5 * d$x1 - 0.8 * d$x2 + rnorm(n) fit <- bqr.svy(y ~ x1 + x2, weights = w, data = d, quantile = c(0.25, 0.5, 0.75), niter = 2000, burnin = 500) coef(fit) # one column per quantile coef(fit, tau = 0.5) # a single named vector vcov(fit, tau = 0.5) head(fitted(fit, tau = 0.5)) nobs(fit) head(weights(fit)) # the survey weights, as supplied # the posterior draws, one row per retained iteration draws <- as.matrix(fit, tau = 0.5) dim(draws) mean(draws[, "x1"] > 0) # P(beta_x1 > 0)
Computes equal-tailed credible intervals for the regression coefficients of
a fitted "bqr.svy" model from its posterior draws.
## S3 method for class 'bqr.svy' posterior_interval(object, prob = 0.95, tau = NULL, ...)## S3 method for class 'bqr.svy' posterior_interval(object, prob = 0.95, tau = NULL, ...)
object |
an object of class |
prob |
the credible level. Default |
tau |
a single fitted quantile. Defaults to the only one in the fit, and is required when several were fitted. |
... |
further arguments passed to or from other methods. |
The values agree with those reported by
summary. They describe posterior uncertainty
about where the conditional quantile lies, and are not prediction intervals
for a new observation.
posterior_interval is the generic defined in rstantools and
used by brms and rstanarm, so the same call works across these
packages.
A matrix with one row per coefficient and columns holding the lower and upper bounds.
bqr.svy.methods, diagnostics,
summary
set.seed(123) n <- 300 d <- data.frame(x1 = runif(n, -1, 1), x2 = runif(n, -1, 1), w = runif(n, 1, 3)) d$y <- 2 + 1.5 * d$x1 - 0.8 * d$x2 + rnorm(n) fit <- bqr.svy(y ~ x1 + x2, weights = w, data = d, quantile = c(0.25, 0.5), niter = 2000, burnin = 500) posterior_interval(fit, tau = 0.5) posterior_interval(fit, tau = 0.5, prob = 0.9)set.seed(123) n <- 300 d <- data.frame(x1 = runif(n, -1, 1), x2 = runif(n, -1, 1), w = runif(n, 1, 3)) d$y <- 2 + 1.5 * d$x1 - 0.8 * d$x2 + rnorm(n) fit <- bqr.svy(y ~ x1 + x2, weights = w, data = d, quantile = c(0.25, 0.5), niter = 2000, burnin = 500) posterior_interval(fit, tau = 0.5) posterior_interval(fit, tau = 0.5, prob = 0.9)
Extracts the convergence diagnostics stored by bqr.svy and
mo.bqr.svy, selecting quantiles numerically.
diagnostics(object, ...) ## S3 method for class 'bqr.svy' diagnostics(object, tau = NULL, ...) ## S3 method for class 'mo.bqr.svy' diagnostics(object, tau = NULL, ...)diagnostics(object, ...) ## S3 method for class 'bqr.svy' diagnostics(object, tau = NULL, ...) ## S3 method for class 'mo.bqr.svy' diagnostics(object, tau = NULL, ...)
object |
a fitted object of class |
... |
further arguments passed to or from other methods. |
tau |
numeric vector of fitted quantiles. |
For "bqr.svy" fits the diagnostics are those of the MCMC sampler, namely
the potential scale reduction factor together with the bulk and
tail effective sample sizes, computed with the posterior package. As a rule of
thumb and effective sample sizes above 400 indicate that
the chain has converged and carries enough information for stable posterior
summaries.
For "mo.bqr.svy" fits estimation is by an EM algorithm rather than MCMC,
so the diagnostics report, for each quantile and direction, the number of
iterations taken and whether the tolerance was met before max_iter.
A data frame of diagnostics for a single quantile, or a named list of such data frames when several quantiles are selected.
bqr.svy, mo.bqr.svy,
bqr.svy.methods
set.seed(123) n <- 300 d <- data.frame(x1 = runif(n, -1, 1), x2 = runif(n, -1, 1), w = runif(n, 1, 3)) d$y <- 2 + 1.5 * d$x1 - 0.8 * d$x2 + rnorm(n) fit <- bqr.svy(y ~ x1 + x2, weights = w, data = d, quantile = c(0.25, 0.5, 0.75), niter = 2000, burnin = 500) diagnostics(fit, tau = 0.5) diagnostics(fit) # all fitted quantilesset.seed(123) n <- 300 d <- data.frame(x1 = runif(n, -1, 1), x2 = runif(n, -1, 1), w = runif(n, 1, 3)) d$y <- 2 + 1.5 * d$x1 - 0.8 * d$x2 + rnorm(n) fit <- bqr.svy(y ~ x1 + x2, weights = w, data = d, quantile = c(0.25, 0.5, 0.75), niter = 2000, burnin = 500) diagnostics(fit, tau = 0.5) diagnostics(fit) # all fitted quantiles
mo.bqr.svy implements a Bayesian approach to multiple-output quantile regression for complex survey data analysis. The method builds a quantile region based on a directional approach. To improve computational efficiency, an Expectation-Maximization (EM) algorithm is implemented instead of the usual Markov Chain Monte Carlo (MCMC).
mo.bqr.svy( formula, weights = NULL, data = NULL, quantile = 0.5, prior = NULL, U = NULL, gamma_U = NULL, n_dir = NULL, epsilon = 1e-06, max_iter = 1000, verbose = FALSE, estimate_sigma = FALSE, seed = NULL )mo.bqr.svy( formula, weights = NULL, data = NULL, quantile = 0.5, prior = NULL, U = NULL, gamma_U = NULL, n_dir = NULL, epsilon = 1e-06, max_iter = 1000, verbose = FALSE, estimate_sigma = FALSE, seed = NULL )
formula |
a symbolic description of the model to be fit. Offset terms are not supported. |
weights |
an optional numerical vector containing the survey weights. If |
data |
an optional data frame containing the variables in the model. Missing values must be removed or imputed before fitting. |
quantile |
numerical scalar or vector containing quantile(s) of interest (default=0.5). |
prior |
a |
U |
an optional |
gamma_U |
an optional list with length equal to |
n_dir |
numerical scalar corresponding to the number of directions (if |
epsilon |
numerical scalar indicating the convergence tolerance for the EM algorithm (default = 1e-6). |
max_iter |
numerical scalar indicating maximum number of EM iterations (default = 1000). |
verbose |
logical flag indicating whether to print progress messages (default=FALSE). |
estimate_sigma |
logical flag indicating whether to estimate the scale parameter
(default=FALSE, with |
seed |
optional integer seed passed to |
An object of class "mo.bqr.svy" with the components listed
below.
call |
The matched call |
formula |
The model formula |
terms |
The terms object |
xlevels, contrasts
|
Factor levels and contrasts of the model frame, kept so that the design matrix can be rebuilt consistently from new data |
quantile |
Vector of fitted quantiles |
prior |
List of priors used for each quantile |
fit |
List of fitted results for each quantile, each containing one sub-list per direction |
coefficients |
Coefficients organized by quantile |
sigma |
List of scale parameters by quantile and direction.
If |
algorithm, tolerance, max_iter
|
EM settings used during estimation |
n_dir |
Number of directions |
U |
Matrix of projection directions ( |
Gamma_list |
List of orthogonal complement bases, one per direction |
model |
The model frame |
weights |
The survey weights supplied |
n_obs |
Number of observations |
n_vars |
Number of covariates |
response_dim |
Dimension of the response |
estimate_sigma |
Logical flag indicating whether the scale parameter
|
To obtain the coefficients, scale estimates and related quantities, use the
generic accessor functions documented in mo.bqr.svy.methods,
together with diagnostics for EM convergence, instead of
accessing these components directly. Estimation is by EM, so a fit stores
posterior modes and carries no quantification of posterior uncertainty.
Nascimento, M. L. & Gonçalves, K. C. M. (2024). Bayesian Quantile Regression Models for Complex Survey Data Under Informative Sampling. Journal of Survey Statistics and Methodology, 12(4), 1105-1130. doi:10.1093/jssam/smae015
library(MASS) # Generate population data set.seed(123) N <- 10000 data <- mvrnorm(N, rep(0, 3), matrix(c(4, 0, 2, 0, 1, 1.5, 2, 1.5, 9), 3, 3)) x_p <- as.matrix(data[, 1]) y_p <- data[, 2:3] + cbind(rep(0, N), x_p) # Generate sample data n <- 500 z_aux <- rnorm(N, mean = 1 + y_p, sd = 0.5) p_aux <- 1 / (1 + exp(2.5 - 0.5 * z_aux)) s_ind <- sample(1:N, n, replace = FALSE, prob = p_aux) y_s <- y_p[s_ind, ] x_s <- x_p[s_ind, ] w <- 1 / p_aux[s_ind] data_s <- data.frame(y1 = y_s[, 1], y2 = y_s[, 2], x1 = x_s, w = w) # Basic usage with default priors when U and gamma_U are given fit1 <- mo.bqr.svy( cbind(y1, y2) ~ x1, weights = w, data = data_s, quantile = c(0.1, 0.2), U = matrix(c(0, 1, 1/sqrt(2), 1/sqrt(2)), 2), gamma_U = list(c(1, 0), c(1/sqrt(2), -1/sqrt(2))) ) # Basic usage with default priors when n_dir is given fit2 <- mo.bqr.svy( cbind(y1, y2) ~ x1, weights = w, data = data_s, quantile = c(0.1, 0.2), n_dir = 2 )library(MASS) # Generate population data set.seed(123) N <- 10000 data <- mvrnorm(N, rep(0, 3), matrix(c(4, 0, 2, 0, 1, 1.5, 2, 1.5, 9), 3, 3)) x_p <- as.matrix(data[, 1]) y_p <- data[, 2:3] + cbind(rep(0, N), x_p) # Generate sample data n <- 500 z_aux <- rnorm(N, mean = 1 + y_p, sd = 0.5) p_aux <- 1 / (1 + exp(2.5 - 0.5 * z_aux)) s_ind <- sample(1:N, n, replace = FALSE, prob = p_aux) y_s <- y_p[s_ind, ] x_s <- x_p[s_ind, ] w <- 1 / p_aux[s_ind] data_s <- data.frame(y1 = y_s[, 1], y2 = y_s[, 2], x1 = x_s, w = w) # Basic usage with default priors when U and gamma_U are given fit1 <- mo.bqr.svy( cbind(y1, y2) ~ x1, weights = w, data = data_s, quantile = c(0.1, 0.2), U = matrix(c(0, 1, 1/sqrt(2), 1/sqrt(2)), 2), gamma_U = list(c(1, 0), c(1/sqrt(2), -1/sqrt(2))) ) # Basic usage with default priors when n_dir is given fit2 <- mo.bqr.svy( cbind(y1, y2) ~ x1, weights = w, data = data_s, quantile = c(0.1, 0.2), n_dir = 2 )
Generic accessor functions for objects of class "mo.bqr.svy" returned
by mo.bqr.svy.
## S3 method for class 'mo.bqr.svy' coef(object, tau = NULL, direction = NULL, ...) ## S3 method for class 'mo.bqr.svy' sigma(object, tau = NULL, direction = NULL, ...) ## S3 method for class 'mo.bqr.svy' nobs(object, ...) ## S3 method for class 'mo.bqr.svy' vcov(object, ...) ## S3 method for class 'mo.bqr.svy' formula(x, ...) ## S3 method for class 'mo.bqr.svy' terms(x, ...) ## S3 method for class 'mo.bqr.svy' model.matrix(object, ...) ## S3 method for class 'mo.bqr.svy' weights(object, ...) ## S3 method for class 'mo.bqr.svy' update(object, formula., ..., evaluate = TRUE)## S3 method for class 'mo.bqr.svy' coef(object, tau = NULL, direction = NULL, ...) ## S3 method for class 'mo.bqr.svy' sigma(object, tau = NULL, direction = NULL, ...) ## S3 method for class 'mo.bqr.svy' nobs(object, ...) ## S3 method for class 'mo.bqr.svy' vcov(object, ...) ## S3 method for class 'mo.bqr.svy' formula(x, ...) ## S3 method for class 'mo.bqr.svy' terms(x, ...) ## S3 method for class 'mo.bqr.svy' model.matrix(object, ...) ## S3 method for class 'mo.bqr.svy' weights(object, ...) ## S3 method for class 'mo.bqr.svy' update(object, formula., ..., evaluate = TRUE)
object |
an object of class |
tau |
numeric vector of fitted quantiles to extract. |
direction |
optional integer vector selecting directions. |
... |
further arguments passed to or from other methods. |
x |
an object of class |
formula. |
a change to the model formula, in the form used by
|
evaluate |
logical; if |
Because the model is directional, results are indexed by quantile and by
direction. For a given quantile the coefficient matrix has one column per
direction, with rows for the covariate coefficients followed by the
coefficients on the projected responses. Quantiles are selected
numerically through tau and directions through direction.
Estimation uses an EM algorithm, so a fit stores the posterior mode of the
coefficients rather than a sample from the posterior. Unlike a fit produced
by bqr.svy, it holds no posterior draws, and a covariance
matrix cannot be computed from point estimates alone. vcov is
therefore defined only to fail with a message saying so, rather than being
left to dispatch to the default method and return something misleading. For
the same reason no posterior_interval method is provided. Use
diagnostics to check that the algorithm converged.
What the model estimates for a given covariate profile is a quantile region
of the response space, the intersection of the half-spaces fitted along each
direction, rather than a single value. That region is computed, drawn and
returned by plot through its xValue
argument.
formula, terms and model.matrix recover the model
specification, weights returns the survey weights as supplied, and
update refits the model with part of the call changed.
coef returns a matrix of coefficients, parameters by
directions, for a single quantile, and a named list of such matrices
otherwise. sigma returns the scale parameter for each direction,
which is 1 throughout when estimate_sigma = FALSE. nobs
returns the number of observations and weights the vector of survey
weights. vcov does not return a value; it throws an error explaining
that this class carries no posterior uncertainty. formula,
terms and model.matrix return the model specification, and
update returns the refitted model.
set.seed(123) n <- 250 d <- data.frame(x1 = runif(n, -1, 1), w = runif(n, 1, 3)) d$y1 <- 1 + 1.2 * d$x1 + rnorm(n) d$y2 <- 0.5 - 0.7 * d$x1 + rnorm(n) fit <- mo.bqr.svy(cbind(y1, y2) ~ x1, weights = w, data = d, quantile = c(0.25, 0.5), n_dir = 4) coef(fit, tau = 0.5) # parameters by directions coef(fit, tau = 0.5, direction = 1:2) # a subset of directions sigma(fit, tau = 0.5) nobs(fit) head(weights(fit)) # the survey weights, as supplied diagnostics(fit, tau = 0.5)set.seed(123) n <- 250 d <- data.frame(x1 = runif(n, -1, 1), w = runif(n, 1, 3)) d$y1 <- 1 + 1.2 * d$x1 + rnorm(n) d$y2 <- 0.5 - 0.7 * d$x1 + rnorm(n) fit <- mo.bqr.svy(cbind(y1, y2) ~ x1, weights = w, data = d, quantile = c(0.25, 0.5), n_dir = 4) coef(fit, tau = 0.5) # parameters by directions coef(fit, tau = 0.5, direction = 1:2) # a subset of directions sigma(fit, tau = 0.5) nobs(fit) head(weights(fit)) # the survey weights, as supplied diagnostics(fit, tau = 0.5)
Plot method for objects of class bqr.svy produced by bqr.svy().
It can display fitted quantile curves, coefficient-quantile profiles,
MCMC trace plots, and posterior densities.
## S3 method for class 'bqr.svy' plot( x, y = NULL, type = c("fit", "quantile", "trace", "density"), tau = NULL, which = NULL, add_points = TRUE, combine = TRUE, show_ci = FALSE, ci_probs = c(0.025, 0.975), at = NULL, grid_length = 200, points_alpha = 0.4, point_size = 2.5, line_size = 1.2, main = NULL, use_ggplot = TRUE, theme_style = c("minimal", "classic", "bw", "light", "none"), color_palette = c("viridis", "plasma", "set2", "dark2", "grey", "none"), add_h0 = FALSE, add_ols = FALSE, ols_fit = NULL, ols_weights = NULL, ... )## S3 method for class 'bqr.svy' plot( x, y = NULL, type = c("fit", "quantile", "trace", "density"), tau = NULL, which = NULL, add_points = TRUE, combine = TRUE, show_ci = FALSE, ci_probs = c(0.025, 0.975), at = NULL, grid_length = 200, points_alpha = 0.4, point_size = 2.5, line_size = 1.2, main = NULL, use_ggplot = TRUE, theme_style = c("minimal", "classic", "bw", "light", "none"), color_palette = c("viridis", "plasma", "set2", "dark2", "grey", "none"), add_h0 = FALSE, add_ols = FALSE, ols_fit = NULL, ols_weights = NULL, ... )
x |
Object of class |
y |
Ignored (S3 signature). |
type |
One of |
tau |
Quantile(s) to plot; must appear in |
which |
Variable name(s) or coefficient index(es) to display.
For |
add_points |
(fit) Logical; overlay observed data points. |
combine |
(fit) Logical; with several |
show_ci |
(fit) Logical; draw credible bands. |
ci_probs |
Length-2 numeric vector with the lower/upper probabilities of
the credible interval shown by the fit ribbon, the density plot bounds, and
the quantile-process band. Default |
at |
(fit) Named list of fixed values for non-plotted covariates (see Details). |
grid_length |
(fit) Integer; number of points in the predictor grid. |
points_alpha |
(fit) Point transparency in |
point_size |
(fit) Point size. |
line_size |
(fit/quantile) Line width for fitted/summary lines. |
main |
Optional main title. |
use_ggplot |
Logical; if |
theme_style |
(ggplot) One of |
color_palette |
(ggplot) One of |
add_h0 |
(quantile) Logical; add a horizontal reference at |
add_ols |
(quantile) Logical; add the OLS estimate (dotted line) for the selected coefficient. |
ols_fit |
(quantile) Optional precomputed |
ols_weights |
(quantile) Optional numeric vector of weights when fitting
OLS internally (length must match |
... |
Not used. |
The type argument selects one of four displays.
type = "fit" draws the fitted quantile curves versus a single
numeric predictor (selected via which). Optionally overlay
observed points and credible bands. Other covariates can be held
fixed via at.
type = "quantile" draws a single coefficient as a function
of the quantile . Optionally add a reference line at 0 and
the corresponding OLS estimate.
type = "trace" draws MCMC traces at a chosen . One
or several coefficients may be selected via which; when more
than one is shown the panels are arranged with facet_wrap.
type = "density" draws posterior densities at a chosen
, faceted over the selected coefficients as for
type = "trace".
Additional details on the arguments.
tau must be included in x$quantile. If NULL, all
available quantiles in the object are used.
For type = "fit", which must name a numeric column in
the original model. If NULL, the first numeric predictor
(different from the response) is chosen automatically.
For type = "fit", at is a named list
(list(var = value, ...)) used to fix other covariates while
plotting versus the selected predictor. Provide valid levels for
factors.
When use_ggplot = TRUE, a ggplot object is returned and the
appearance is controlled by theme_style and
color_palette. Otherwise, base graphics are used and the
function returns invisible(NULL).
invisible(NULL) for base R graphics, or a ggplot object if
use_ggplot = TRUE.
data(mtcars) fit <- bqr.svy(mpg ~ wt + hp + cyl, data = mtcars, quantile = c(0.5, 0.75), method = "ald", niter = 20000, burnin = 10000, thin = 5) plot(fit, type = "fit", which = "wt", show_ci = TRUE) plot(fit, type = "quantile", which = "wt", add_h0 = TRUE, add_ols = TRUE) plot(fit, type = "quantile", which = c("(Intercept)", "wt", "hp", "cyl"), add_h0 = TRUE, add_ols = TRUE) plot(fit, type = "trace", which = "wt", tau = 0.5) plot(fit, type = "density", which = "wt", tau = 0.5)data(mtcars) fit <- bqr.svy(mpg ~ wt + hp + cyl, data = mtcars, quantile = c(0.5, 0.75), method = "ald", niter = 20000, burnin = 10000, thin = 5) plot(fit, type = "fit", which = "wt", show_ci = TRUE) plot(fit, type = "quantile", which = "wt", add_h0 = TRUE, add_ols = TRUE) plot(fit, type = "quantile", which = c("(Intercept)", "wt", "hp", "cyl"), add_h0 = TRUE, add_ols = TRUE) plot(fit, type = "trace", which = "wt", tau = 0.5) plot(fit, type = "density", which = "wt", tau = 0.5)
Draws quantile regions from a mo.bqr.svy model fitted with exactly two
response variables. Models with other response dimensions cannot be plotted. The
method projects the data onto a grid, determines which points lie inside each
quantile region, and visualises the result using ggplot2.
The response names and the data are taken from the fitted object, so neither
has to be supplied again; response and datafile are only needed
to override them. A quantile region is conditional on the covariates, so
xValue must give the covariate profile at which to draw it, one value
per design-matrix column.
## S3 method for class 'mo.bqr.svy' plot( x, y = NULL, response = NULL, datafile = NULL, ngridpoints = 200, xValue = 1, paintedArea = TRUE, range_y = NULL, main = NULL, color_palette = c("default", "grey"), theme_style = c("minimal", "none"), point_alpha = 0.3, point_size = 1.2, line_size = 0.8, verbose = FALSE, ... )## S3 method for class 'mo.bqr.svy' plot( x, y = NULL, response = NULL, datafile = NULL, ngridpoints = 200, xValue = 1, paintedArea = TRUE, range_y = NULL, main = NULL, color_palette = c("default", "grey"), theme_style = c("minimal", "none"), point_alpha = 0.3, point_size = 1.2, line_size = 0.8, verbose = FALSE, ... )
x |
An object of class |
y |
Ignored (S3 signature). |
response |
Optional character vector of length 2 giving the names of the fitted response columns, in the desired axis order. Defaults to the order used to fit the model; the two names may be reversed. |
datafile |
Optional data frame containing at least the two columns named
in |
ngridpoints |
Integer; number of grid points per axis (default = 200). |
xValue |
Numeric vector of covariate values at which to evaluate the
regression, one per design-matrix column and in that order. For an
intercept-only model use |
paintedArea |
Logical; if |
range_y |
An optional |
main |
Optional character string for the plot title. When |
color_palette |
One of |
theme_style |
One of |
point_alpha |
Numeric in |
point_size |
Numeric; size of the data points (default = 1.2). |
line_size |
Numeric; line width for contour outlines (default = 0.8). |
verbose |
Logical; if |
... |
Not used. |
Two display modes are available.
paintedArea = TRUE (default) draws filled regions with
contour outlines, one per quantile magnitude, from the outermost
(smallest ) to the innermost.
paintedArea = FALSE draws contour lines only, distinguished
by colour and linetype.
In both modes the quantile levels are identified through a standard
ggplot2 legend (labelled = level), so its position and
styling follow the active theme.
The quantile regions are built from the directional approach described in
Nascimento & Gonçalves (2025). For each quantile
, the function evaluates every grid point against the
half-space constraints defined by the fitted directions and orthogonal
bases.
Invisibly, a list with the components listed below.
plot |
A |
data |
A data frame with columns |
mo.bqr.svy, mo.bqr.svy.methods,
prior
# Load the Anthro dataset (children anthropometric data, already cleaned) data(Anthro, package = "bayesQRsurvey") # --- Directions --- n_dir <- 12 angles <- (0:(n_dir - 1)) * 2 * pi / n_dir U_mat <- rbind(cos(angles), sin(angles)) gamma_list <- lapply(seq_len(n_dir), function(k) matrix(c(-sin(angles[k]), cos(angles[k])), ncol = 1)) # --- Fit --- fit <- mo.bqr.svy( cbind(wgt, hgt) ~ 1, data = Anthro, quantile = c(0.05, 0.10, 0.25), U = U_mat, gamma_U = gamma_list, max_iter = 2000, verbose = FALSE ) # --- Plot quantile regions (filled ribbons) --- # the responses and the data come from the fit; the model is intercept-only, # so the default xValue = 1 is the right covariate profile plot(fit) # Contour-only style plot(fit, paintedArea = FALSE)# Load the Anthro dataset (children anthropometric data, already cleaned) data(Anthro, package = "bayesQRsurvey") # --- Directions --- n_dir <- 12 angles <- (0:(n_dir - 1)) * 2 * pi / n_dir U_mat <- rbind(cos(angles), sin(angles)) gamma_list <- lapply(seq_len(n_dir), function(k) matrix(c(-sin(angles[k]), cos(angles[k])), ncol = 1)) # --- Fit --- fit <- mo.bqr.svy( cbind(wgt, hgt) ~ 1, data = Anthro, quantile = c(0.05, 0.10, 0.25), U = U_mat, gamma_U = gamma_list, max_iter = 2000, verbose = FALSE ) # --- Plot quantile regions (filled ribbons) --- # the responses and the data come from the fit; the model is intercept-only, # so the default xValue = 1 is the right covariate profile plot(fit) # Contour-only style plot(fit, paintedArea = FALSE)
Evaluates the estimated conditional quantile function
at new covariate
values, optionally with pointwise credible intervals.
## S3 method for class 'bqr.svy' predict( object, newdata, tau = NULL, interval = c("none", "credible"), level = 0.95, na.action = stats::na.pass, ... )## S3 method for class 'bqr.svy' predict( object, newdata, tau = NULL, interval = c("none", "credible"), level = 0.95, na.action = stats::na.pass, ... )
object |
an object of class |
newdata |
an optional data frame of covariate values. If omitted, the fitted values are returned. |
tau |
numeric vector of fitted quantiles to predict. |
interval |
|
level |
the credible level of each pointwise interval. |
na.action |
how to handle missing values in |
... |
further arguments passed to or from other methods. |
What is returned is an estimate of the conditional quantile. With
tau = 0.9, the fitted model estimates the value below which 90% of the
responses fall among units sharing those covariates. It summarizes the conditional
distribution, in the same sense that predict.lm returns
the conditional mean, and not a draw of a new response.
interval = "credible" adds equal-tailed, pointwise intervals obtained
from the pseudo-posterior draws for each of the three estimation methods.
Each draw of gives a draw of , and
lwr and upr are the percentiles of those values at the
requested level. These intervals describe uncertainty about the
conditional quantile under the working model. They are not simultaneous
bands and do not describe the variability of individual responses. Nominal
frequentist or sampling-design coverage is not guaranteed by the
pseudo-posterior construction.
There is no interval = "prediction". Lower and upper conditional
quantiles can be used to estimate a prediction interval for a new response,
but this method does not construct or calibrate such intervals and does not
generate posterior predictive draws. For method = "ald", the
asymmetric Laplace distribution is used as a working likelihood; its shape
away from is not assumed to describe the response distribution.
When several quantiles are requested, each was fitted separately and nothing constrains them to be ordered, so predicted quantiles may cross at some covariate values.
With interval = "none", a vector for a single quantile and a
matrix with one column per quantile otherwise. With
interval = "credible", a matrix with columns fit,
lwr and upr, as in predict.lm, for a
single quantile, and a named list of such matrices otherwise.
bqr.svy.methods,
posterior_interval,
predict.rq
set.seed(123) n <- 300 d <- data.frame(x1 = runif(n, -1, 1), x2 = runif(n, -1, 1), w = runif(n, 1, 3)) d$y <- 2 + 1.5 * d$x1 - 0.8 * d$x2 + rnorm(n) fit <- bqr.svy(y ~ x1 + x2, weights = w, data = d, quantile = c(0.1, 0.5, 0.9), niter = 2000, burnin = 500) nd <- data.frame(x1 = c(-0.5, 0, 0.5), x2 = 0) predict(fit, nd, tau = 0.9) predict(fit, nd, tau = 0.9, interval = "credible") predict(fit, nd) # one column per fitted quantileset.seed(123) n <- 300 d <- data.frame(x1 = runif(n, -1, 1), x2 = runif(n, -1, 1), w = runif(n, 1, 3)) d$y <- 2 + 1.5 * d$x1 - 0.8 * d$x2 + rnorm(n) fit <- bqr.svy(y ~ x1 + x2, weights = w, data = d, quantile = c(0.1, 0.5, 0.9), niter = 2000, burnin = 500) nd <- data.frame(x1 = c(-0.5, 0, 0.5), x2 = 0) predict(fit, nd, tau = 0.9) predict(fit, nd, tau = 0.9, interval = "credible") predict(fit, nd) # one column per fitted quantile
Print methods for the objects returned by bqr.svy and
mo.bqr.svy. They print an S3 object of class
bqr.svy or mo.bqr.svy to the console.
## S3 method for class 'bqr.svy' print(x, digits = 3, ...) ## S3 method for class 'mo.bqr.svy' print(x, ...)## S3 method for class 'bqr.svy' print(x, digits = 3, ...) ## S3 method for class 'mo.bqr.svy' print(x, ...)
x |
An object of class |
digits |
Integer specifying the number of decimal places to print. Defaults to |
... |
Additional arguments that are passed to the generic |
set.seed(123) N <- 10000 x1_p <- runif(N, -1, 1) x2_p <- runif(N, -1, 1) y_p <- 2 + 1.5 * x1_p - 0.8 * x2_p + rnorm(N) # Generate sample data n <- 500 z_aux <- rnorm(N, mean = 1 + y_p, sd = .5) p_aux <- 1 / (1 + exp(2.5 - 0.5 * z_aux)) s_ind <- sample(1:N, n, replace = FALSE, prob = p_aux) y_s <- y_p[s_ind] x1_s <- x1_p[s_ind] x2_s <- x2_p[s_ind] w <- 1 / p_aux[s_ind] data <- data.frame(y = y_s, x1 = x1_s, x2 = x2_s, w = w) # Fit a model fit1 <- bqr.svy(y ~ x1 + x2, weights = w, data = data, niter = 2000, burnin = 500, thin = 2) print(fit1)set.seed(123) N <- 10000 x1_p <- runif(N, -1, 1) x2_p <- runif(N, -1, 1) y_p <- 2 + 1.5 * x1_p - 0.8 * x2_p + rnorm(N) # Generate sample data n <- 500 z_aux <- rnorm(N, mean = 1 + y_p, sd = .5) p_aux <- 1 / (1 + exp(2.5 - 0.5 * z_aux)) s_ind <- sample(1:N, n, replace = FALSE, prob = p_aux) y_s <- y_p[s_ind] x1_s <- x1_p[s_ind] x2_s <- x2_p[s_ind] w <- 1 / p_aux[s_ind] data <- data.frame(y = y_s, x1 = x1_s, x2 = x2_s, w = w) # Fit a model fit1 <- bqr.svy(y ~ x1 + x2, weights = w, data = data, niter = 2000, burnin = 500, thin = 2) print(fit1)
prior creates prior distributions for both single (bqr.svy) and multiple-output
(mo.bqr.svy) Bayesian quantile regression models for complex survey data.
prior( beta_x_mean = NULL, beta_x_cov = NULL, sigma_shape = 0.001, sigma_rate = 0.001, beta_y_mean = NULL, beta_y_cov = NULL )prior( beta_x_mean = NULL, beta_x_cov = NULL, sigma_shape = 0.001, sigma_rate = 0.001, beta_y_mean = NULL, beta_y_cov = NULL )
beta_x_mean |
vector of prior means for the regression coefficients. (default = NULL). |
beta_x_cov |
prior covariance matrix for the regression coefficients. (default = NULL). |
sigma_shape |
shape parameter for the inverse Gamma prior on the ALD
scale |
sigma_rate |
parameter multiplying |
beta_y_mean |
prior means for the coefficients related to the variables that emerge from the product between the orthogonal basis and the outputs (default = NULL). |
beta_y_cov |
prior covariance matrix for the coefficients related to the variables that emerge from the product between the orthogonal basis and the outputs. (default = NULL). |
The function prior builds prior distributions for the three methods implemented in the function
bqr.svy and for the multiple-output quantile regression implemented in the function mo.bqr.svy.
Every nonspecified prior parameter will get the default value.
The scale prior has density proportional to
, where is
sigma_shape and is sigma_rate.
method = "ald" in function bqr.svy allow the specification of hyperparameters
beta_x_mean, beta_x_cov, sigma_shape, and sigma_rate.
method = "score" in function bqr.svy allow the specification of hyperparameters
beta_x_mean and beta_x_cov.
method = "approximate" in function bqr.svy allow the specification of hyperparameters
beta_x_mean and beta_x_cov.
In function mo.bqr.svy, the specification of hyperparameters beta_x_mean,beta_x_cov,
sigma_shape, sigma_rate, beta_y_mean, and beta_y_cov are allowed.
An object of class "prior".
# Simulate data set.seed(123) n <- 200 x1 <- rnorm(n, 0, 1) x2 <- runif(n, -1, 1) w <- runif(n, 0.5, 2) # survey weights y1 <- 2 + 1.5*x1 - 0.8*x2 + rnorm(n, 0, 1) y2 <- 1 + 0.5*x1 - 0.2*x2 + rnorm(n, 0, 1) data <- data.frame(y1 = y1, y2 = y2, x1 = x1, x2 = x2, w = w) # Define a general informative prior prior_general <- prior( beta_x_mean = c(2, 1.5, -0.8), beta_x_cov = diag(c(0.25, 0.25, 0.25)), sigma_shape = 3, sigma_rate = 2, beta_y_mean = 1, beta_y_cov = 0.25 ) # Estimate the model parameters with informative prior fit_ald <- bqr.svy(y1 ~ x1 + x2, weights = w, data = data, prior = prior_general, method = "ald") fit_scr <- bqr.svy(y1 ~ x1 + x2, weights = w, data = data, prior = prior_general, method = "score") fit_apx <- bqr.svy(y1 ~ x1 + x2, weights = w, data = data, prior = prior_general, method = "approximate") # Multiple-output method fit_mo <- mo.bqr.svy(cbind(y1, y2) ~ x1 + x2, weights = w, data = data, prior = prior_general, n_dir = 10) plot(fit_ald, type = "trace", which = "x1", tau = 0.5) plot(fit_ald, type = "trace", which = "x2", tau = 0.5) print(fit_mo)# Simulate data set.seed(123) n <- 200 x1 <- rnorm(n, 0, 1) x2 <- runif(n, -1, 1) w <- runif(n, 0.5, 2) # survey weights y1 <- 2 + 1.5*x1 - 0.8*x2 + rnorm(n, 0, 1) y2 <- 1 + 0.5*x1 - 0.2*x2 + rnorm(n, 0, 1) data <- data.frame(y1 = y1, y2 = y2, x1 = x1, x2 = x2, w = w) # Define a general informative prior prior_general <- prior( beta_x_mean = c(2, 1.5, -0.8), beta_x_cov = diag(c(0.25, 0.25, 0.25)), sigma_shape = 3, sigma_rate = 2, beta_y_mean = 1, beta_y_cov = 0.25 ) # Estimate the model parameters with informative prior fit_ald <- bqr.svy(y1 ~ x1 + x2, weights = w, data = data, prior = prior_general, method = "ald") fit_scr <- bqr.svy(y1 ~ x1 + x2, weights = w, data = data, prior = prior_general, method = "score") fit_apx <- bqr.svy(y1 ~ x1 + x2, weights = w, data = data, prior = prior_general, method = "approximate") # Multiple-output method fit_mo <- mo.bqr.svy(cbind(y1, y2) ~ x1 + x2, weights = w, data = data, prior = prior_general, n_dir = 10) plot(fit_ald, type = "trace", which = "x1", tau = 0.5) plot(fit_ald, type = "trace", which = "x2", tau = 0.5) print(fit_mo)
Summary methods for the objects returned by bqr.svy and
mo.bqr.svy. For a "bqr.svy" fit the posterior mean,
the posterior credible interval and the convergence diagnostics are reported
for every quantile. For a "mo.bqr.svy" fit, which is estimated by EM,
the posterior mode, the scale parameter and the number of iterations are
reported for every quantile and direction.
## S3 method for class 'bqr.svy' summary(object, probs = c(0.025, 0.975), digits = 3, ...) ## S3 method for class 'mo.bqr.svy' summary(object, digits = 4, ...)## S3 method for class 'bqr.svy' summary(object, probs = c(0.025, 0.975), digits = 3, ...) ## S3 method for class 'mo.bqr.svy' summary(object, digits = 4, ...)
object |
An object of class |
probs |
Two-element numeric vector with credible interval probabilities.
Default |
digits |
Integer; number of decimals used by printing helpers. Default |
... |
Not used. |
An object of class summary.bqr.svy with one block per .
An object of class summary.mo.bqr.svy.