--- title: "Bayesian quantile regression for complex survey data" author: | Tomás Rodríguez Taborda, Johnatan Cardona Jiménez, Marcus L. Nascimento, Kelly C. M. Gonçalves output: rmarkdown::html_vignette vignette: > %\VignetteIndexEntry{Bayesian quantile regression for complex survey data} %\VignetteEngine{knitr::rmarkdown} %\VignetteEncoding{UTF-8} --- ```{r, include = FALSE} knitr::opts_chunk$set(collapse = TRUE, comment = "#>", fig.width = 7, fig.height = 4.5, fig.align = "center") ``` ## 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. ```{r, eval = FALSE} install.packages("bayesQRsurvey") ``` ```{r, message = FALSE} 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. ```{r} 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. ```{r} data("Anthro", package = "bayesQRsurvey") Anthro$age <- Anthro$age / 12 Anthro$sex <- factor(Anthro$sex, levels = c("1", "0"), labels = c("Boys", "Girls")) str(Anthro) ``` ## 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. ```{r} 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 ``` 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. ```{r} print(summary(fit_ald), tau = 0.5) ``` Convergence diagnostics are computed for every level and stored in `$diagnosis`. ```{r} fit_ald$diagnosis[["tau=0.500"]] ``` ### 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. ```{r, fig.height = 5} plot(fit_ald, type = "trace", tau = 0.5, color_palette = "grey", theme_style = "none") ``` ```{r, fig.height = 5} 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. ```{r} 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. ```{r, fig.height = 5} 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. ```{r} 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. ```{r} 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) ``` ## 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`. ```{r} 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 ``` 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. ```{r} print(summary(fit_mo), coefficients = FALSE) ``` 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. ```{r, fig.height = 5.5} 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. ## 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](https://doi.org/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](https://doi.org/10.1093/jssam/smaf040)