Bayesian quantile regression for complex survey data
Overview
bayesQRsurvey fits Bayesian quantile regression models to data collected under informative sampling designs, where the probability of selection is related to the response and ignoring it biases the estimates. The survey weights enter the working likelihood directly, so that the sampling design is accounted for in the posterior rather than through a subsequent correction.
The package covers two model families:
bqr.svy()fits models with a single response by MCMC, following Nascimento and Gonçalves (2024). Three working likelihoods are available throughmethod:"ald", based on the asymmetric Laplace distribution, which is the default and the fastest; and"score"and"approximate", which build a design-based variance into the likelihood itself.mo.bqr.svy()fits models with a vector response by an EM algorithm, following the directional approach of Nascimento and Gonçalves (2026). Here a quantile is not a curve but a nested family of convex regions.
Both functions return objects with print(),
summary() and plot() methods, together with
the generic accessor functions used for regression fits in base R
(coef(), fitted(), vcov(),
nobs(), weights(), predict(), …);
see ?bqr.svy.methods and ?mo.bqr.svy.methods
for the full list. This vignette introduces them through a worked
example.
The package is installed from CRAN in the usual way.
Every plotting function accepts color_palette = "grey"
for a print-friendly greyscale rendering, and
theme_style = "none", which applies no theme of its own so
that the plot honours the one set by the user. The settings below are
used throughout.
theme_set(theme_classic(base_size = 14) +
theme(axis.text = element_text(colour = "black"),
legend.title = element_blank(),
legend.background = element_blank()))
scale_tau <- scale_colour_manual(
values = c("0.1" = "grey70", "0.5" = "grey40", "0.9" = "black"),
breaks = c("0.9", "0.5", "0.1"),
labels = c("0.9" = expression(tau == 0.9),
"0.5" = expression(tau == 0.5),
"0.1" = expression(tau == 0.1)))Data
The Anthro data set holds anthropometric measurements of
985 children, together with the survey weight in dweight.
Age is recorded in months, so it is converted to years, and
sex is labelled for readability.
data("Anthro", package = "bayesQRsurvey")
Anthro$age <- Anthro$age / 12
Anthro$sex <- factor(Anthro$sex, levels = c("1", "0"),
labels = c("Boys", "Girls"))
str(Anthro)
#> 'data.frame': 985 obs. of 8 variables:
#> $ wgt : num 11.3 8.1 11.5 10 6.6 17.1 10 7.5 14 9.4 ...
#> $ hgt : num 79.7 68.6 86.1 73.1 65.4 102 79.8 68.7 92.7 80.7 ...
#> $ wgt_ind: num -0.982 -0.77 0.346 1.372 -0.858 ...
#> $ hgt_ind: num -2.984 -1.312 0.546 0.753 -0.201 ...
#> $ ruc : num 1 1 1 1 1 1 1 1 1 1 ...
#> $ sex : Factor w/ 2 levels "Boys","Girls": 1 1 2 2 2 2 2 2 1 2 ...
#> $ age : num 2.167 0.667 1.75 0.75 0.5 ...
#> $ dweight: num 1259 1259 2686 2038 1784 ...Single-output models
Weight is modelled as a quadratic function of age, adjusting for sex,
at three quantile levels. The sampling weights are passed through
weights and the levels through quantile.
Setting verbose = FALSE suppresses the progress bar. The
chains used here are short so that the vignette builds quickly; longer
runs are advisable in practice.
set.seed(50)
fit_ald <- bqr.svy(wgt ~ age + I(age^2) + sex, weights = dweight,
data = Anthro, quantile = c(0.1, 0.5, 0.9),
niter = 6000, burnin = 3000, thin = 1, verbose = FALSE)
fit_ald
#> Bayesian Quantile Regression for survey data
#> Method : ald
#> Quantile : 0.100 0.500 0.900
#> Formula : wgt ~ age + I(age^2) + sex
#> Runtime : 1.527 sec
#>
#> Coefficients (posterior means):
#> tau=0.100 tau=0.500 tau=0.900
#> (Intercept) 4.249 6.152 7.949
#> age 3.833 3.555 3.358
#> I(age^2) -0.377 -0.290 -0.206
#> sexGirls -0.551 -0.585 -0.510
#>
#> (sigma fixed at 1)The summary() method reports posterior means together
with 95% credible intervals. Passing tau to its
print() method restricts the output to a single quantile
level.
print(summary(fit_ald), tau = 0.5)
#>
#> Method: ald
#> Quantiles: 0.500
#>
#> == tau=0.500 ==
#> Draws: 3000 | Warmup: 3000 | Thin: 1
#>
#> variable mean lower_ci upper_ci
#> (Intercept) 6.152 5.793 6.485
#> age 3.555 3.216 3.919
#> I(age^2) -0.290 -0.365 -0.218
#> sexGirls -0.585 -0.800 -0.368Convergence diagnostics are computed for every level and retrieved
with diagnostics(), selecting the level numerically.
diagnostics(fit_ald, tau = 0.5)
#> variable rhat ess_bulk ess_tail
#> 1 (Intercept) 1.003 563 1004
#> 2 age 1.002 477 755
#> 3 I(age^2) 1.001 451 811
#> 4 sexGirls 1.005 522 756The estimates are reached through the usual accessors rather than the
components of the object: coef() returns the posterior
means, vcov() the posterior covariance matrix, and
as.matrix() the draws themselves. predict()
evaluates the estimated conditional quantile at new covariate values —
the quantile itself, not the average response — with an optional
credible band from the draws.
coef(fit_ald, tau = 0.5)
#> (Intercept) age I(age^2) sexGirls
#> 6.1522699 3.5545301 -0.2898247 -0.5847088
nd <- data.frame(age = 2, sex = c("Boys", "Girls"))
predict(fit_ald, newdata = nd)
#> tau=0.100 tau=0.500 tau=0.900
#> 1 10.406326 12.10203 13.84022
#> 2 9.855356 11.51732 13.33007
predict(fit_ald, newdata = nd, tau = 0.9, interval = "credible")
#> fit lwr upr
#> 1 13.84022 13.53485 14.16076
#> 2 13.33007 13.01466 13.62658Graphical output
The plot() method takes a type argument.
Trace plots and marginal posterior densities assess convergence; when
which is omitted, every coefficient is displayed.
Setting type = "fit" draws the fitted quantile curves
against a chosen predictor. Here color_palette = "none"
leaves the colours to the manual scale defined above.
plot(fit_ald, type = "fit", which = "age", add_points = FALSE,
color_palette = "none", theme_style = "none") +
scale_tau + labs(x = "Age (years)", y = "Weight (kg)") +
theme(legend.position = "inside", legend.position.inside = c(0.85, 0.18))
Fitting over a grid of quantile levels and setting
type = "quantile" shows how each coefficient varies across
the conditional distribution. The argument add_ols = TRUE
overlays the least-squares estimate for comparison.
set.seed(50)
fit_grid <- bqr.svy(wgt ~ age + I(age^2) + sex, weights = dweight,
data = Anthro, quantile = seq(0.1, 0.9, by = 0.2),
niter = 6000, burnin = 3000, thin = 1, verbose = FALSE)
plot(fit_grid, type = "quantile", add_ols = TRUE,
color_palette = "grey", theme_style = "none") + labs(x = "quantile")
The quadratic term attenuates towards the upper tail, a pattern that a single summary of the conditional distribution, such as the least-squares fit, cannot reveal.
Priors and estimation methods
Prior distributions are specified with prior() and
supplied through the prior argument. The default is
vague.
myprior <- prior(beta_x_mean = rep(0, 4), beta_x_cov = 25)
set.seed(50)
fit_prior <- bqr.svy(wgt ~ age + I(age^2) + sex, weights = dweight,
data = Anthro, quantile = 0.5,
niter = 6000, burnin = 3000, thin = 1,
prior = myprior, verbose = FALSE)The remaining two methods are selected through method.
Both are adaptive Metropolis-Hastings samplers and require longer chains
than the default method.
set.seed(50)
fit_score <- bqr.svy(wgt ~ age + I(age^2) + sex, weights = dweight,
data = Anthro, method = "score", quantile = 0.5,
niter = 20000, burnin = 5000, thin = 1, verbose = FALSE)
summary(fit_score)
#>
#> Method: score
#> Quantiles: 0.500
#>
#> == tau=0.500 ==
#> Acceptance rate (avg): 0.213
#> Draws: 15000 | Warmup: 5000 | Thin: 1
#>
#> variable mean lower_ci upper_ci
#> (Intercept) 6.113 5.641 6.579
#> age 3.599 3.121 4.085
#> I(age^2) -0.298 -0.402 -0.201
#> sexGirls -0.608 -0.877 -0.305Multiple-output models
With two responses a quantile is a region rather than a curve. The
model is specified by binding the responses on the left-hand side of the
formula, and the directions used to build the region are generated
automatically according to n_dir.
set.seed(50)
fit_mo <- mo.bqr.svy(cbind(wgt, hgt) ~ age + I(age^2) + sex,
weights = dweight, data = Anthro,
quantile = c(0.05, 0.10, 0.15),
n_dir = 20, max_iter = 2000, verbose = FALSE)
fit_mo
#>
#> Multiple-Output Bayesian Quantile Regression
#> --------------------------------------------------------
#> Formula : cbind(wgt, hgt) ~ age + I(age^2) + sex
#> Quantiles : 0.050, 0.100, 0.150
#> Directions : 20
#> Sample : 985 obs, 2 responses
#> Scale : sigma fixed at 1
#> Coefficients: 5 per direction
#> --------------------------------------------------------
#> Use summary() for coefficients and convergence details.Each direction contributes one half-space, and the quantile region is their intersection, so it is the region as a whole rather than any individual direction that is of interest. The compact form of the summary reports only whether the directional fits converged.
print(summary(fit_mo), coefficients = FALSE)
#>
#> Multiple-Output Bayesian Quantile Regression (Summary)
#> --------------------------------------------------------
#> Quantiles : 0.050, 0.100, 0.150
#> Directions : 20
#> Sample : 985 obs, 2 responses
#> Estimation : EM (posterior mode / MAP)
#> --------------------------------------------------------
#>
#> EM convergence by magnitude (iterations: min / median / max)
#>
#> tau converged iterations sigma
#> 0.050 20/20 79 / 163 / 1234 1
#> 0.100 20/20 71 / 143.5 / 576 1
#> 0.150 20/20 71 / 161.5 / 382 1The plot() method intersects the half-spaces for a
chosen covariate profile and draws the resulting contours over the
observed responses. The response names and the data are taken from the
fitted object. The profile is supplied through xValue, in
the order of the design-matrix columns, here (Intercept),
age, I(age^2) and sexGirls,
evaluated at a two-year-old boy.
plot(fit_mo, xValue = c(1, 2, 4, 0),
ngridpoints = 200, paintedArea = FALSE,
color_palette = "grey", theme_style = "none")
The contours are nested and convex, and they are read from the outside in: the smallest quantile magnitude traces the outermost region, and each larger one a more central region. Because every region is tied to a covariate profile, evaluating the function at different profiles shows how the joint distribution of weight and height shifts with age and sex, which a separate analysis of either response cannot recover.
References
Nascimento ML, Gonçalves KCM (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
Nascimento ML, Gonçalves KCM (2026). “A Bayesian Approach to Multiple-Output Quantile Regression Analysis under Informative Sampling.” Journal of Survey Statistics and Methodology, smaf040. doi:10.1093/jssam/smaf040