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 through method: "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. This vignette introduces them through a worked example.

The package is installed from CRAN in the usual way.

install.packages("bayesQRsurvey")

library("bayesQRsurvey")
library("ggplot2")

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.02 sec
#> 
#> Coefficients (posterior means):
#>             tau=0.100 tau=0.500 tau=0.900
#> (Intercept)     4.249     6.167     7.948
#> age             3.833     3.536     3.351
#> I(age^2)       -0.377    -0.286    -0.204
#> sexGirls       -0.555    -0.586    -0.509
#> 
#> (sigma fixed at 1)
#> 
#> Acceptance rate by quantile:
#> tau=0.100 tau=0.500 tau=0.900 
#>        NA        NA        NA

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.167    5.818    6.503
#>          age  3.536    3.197    3.886
#>     I(age^2) -0.286   -0.360   -0.216
#>     sexGirls -0.586   -0.831   -0.378

Convergence diagnostics are computed for every level and stored in $diagnosis.

fit_ald$diagnosis[["tau=0.500"]]
#>      variable  rhat ess_bulk ess_tail
#> 1 (Intercept) 1.001      634      921
#> 2         age 1.000      526      902
#> 3    I(age^2) 1.000      518      814
#> 4    sexGirls 1.001      479      849

Graphical output

The plot() method takes a type argument. Trace plots and marginal posterior densities assess convergence; when which is omitted, every coefficient is displayed.

plot(fit_ald, type = "trace", tau = 0.5,
     color_palette = "grey", theme_style = "none")

plot(fit_ald, type = "density", tau = 0.5,
     color_palette = "grey", theme_style = "none")

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.305

Multiple-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 76 / 149.5 / 1252     1
#>  0.100     20/20  64 / 140.5 / 580     1
#>  0.150     20/20    49 / 163 / 345     1

The function plotQuantileRegion() intersects the half-spaces for a chosen covariate profile and draws the resulting contours over the observed responses. 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.

plotQuantileRegion(fit_mo, response = c("wgt", "hgt"),
                   datafile = Anthro, 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.