Bayesian Quantile Regression using Stan

Exploring robust quantile estimation under uncertainty with Stan

Key Highlights

Asymmetric Laplace Distribution: The backbone for modeling conditional quantiles.
Robustness and Uncertainty: Handles outliers effectively while incorporating prior information.
Stan & brms Integration: Both custom Stan coding and the brms package provide flexible implementations.

Introduction to Bayesian Quantile Regression

Bayesian quantile regression is an advanced statistical modeling technique that goes beyond traditional mean regression by estimating various conditional quantiles of a response variable. Standard linear regression focuses on modeling the conditional mean of the dependent variable given predictors, whereas quantile regression estimates the conditional quantile functions.

When approached from a Bayesian perspective, this methodology makes use of prior distributions to incorporate existing knowledge and propagates uncertainty through probabilistic inference. With Bayesian quantile regression, practitioners obtain a fuller picture of the data by examining different parts of the response distribution, ensuring robustness especially in the presence of skewed distributions or outliers.

Core Concepts and Theoretical Framework

The Role of the Asymmetric Laplace Distribution

At the heart of Bayesian quantile regression using Stan is the asymmetric Laplace distribution (ALD). This distribution is chosen for its unique property: its maximum likelihood estimates correspond to the estimates obtained from classical quantile regression. For a given quantile level \( \tau \) (where \( 0 < \tau < 1 \)), the ALD is parameterized in such a way that its mode or location parameter aligns with the targeted conditional quantile.

The density function of the asymmetric Laplace distribution can be conceptualized in Bayesian settings as:

\( \displaystyle p(y_i|\mu_i, \sigma, \tau) \propto \frac{\tau(1-\tau)}{\sigma}\exp \left\{-\frac{\rho_\tau (y_i-\mu_i)}{\sigma}\right\} \)

where the check function \( \rho_\tau(u) \) is defined as:

\( \displaystyle \rho_\tau(u) = u(\tau - \mathbf{1}_{\{u < 0\}}) \)

This formulation makes the ALD particularly apt for modeling quantiles. In a Bayesian framework, the ALD is used as the likelihood function, while prior distributions are assigned to both the regression coefficients and the scale parameter.

Model Specification in Stan

Data Preparation

The first step is to prepare and format the data properly for Stan's inference engine. This includes defining the number of observations, the number of predictors, the matrix of predictors, and specifying the quantile of interest (denoted \( \tau \)). Stan requires data to be structured in a specific way, typically as vectors and matrices.

Parameter Declaration and Prior Setting

In a Bayesian setting, parameters such as the regression coefficients (\( \beta \)) and a scale parameter (\( \sigma \)) are endowed with prior distributions. Common choices include normal priors for coefficients to reflect initial beliefs of generally moderate effects and exponential priors for the scale to enforce positivity.

Likelihood Specification Using the ALD

The likelihood function in Stan is constructed by applying the asymmetric Laplace distribution to each observation. Although Stan does not have a built-in ALD function, it permits the definition of custom probability density functions. By defining a custom function (e.g., asymmetric_laplace_lpdf) in the functions block, one can implement the ALD manually.

Implementation Strategies

Implementing Bayesian quantile regression using Stan can follow two primary routes. The first involves writing custom Stan code to define the ALD and fit the model using Markov Chain Monte Carlo (MCMC) methods. A sample code structure might involve a functions block for the ALD, a data block for inputs, a parameters block, and the model block to tie together priors and the likelihood.

The second strategy leverages the brms package in R. This package provides an easier interface to Stan, allowing users to specify models using R formulas. With brms, the asymmetric Laplace family is designated through a specific family function (asym_laplace()), and then one calls the main model fitting function, benefitting from automated diagnostics and posterior analyses.

Advantages and Considerations

Advantages of Bayesian Quantile Regression

Bayesian quantile regression offers several advantages over classical techniques, especially in handling irregular data:

Robustness to Outliers: By modeling different parts of the response distribution, the procedure is less influenced by extreme values.
Incorporation of Prior Information: Prior distributions allow analysts to include existing knowledge or beliefs into the model.
Comprehensive Uncertainty Quantification: Bayesian methods provide full posterior distributions, which offer a more detailed uncertainty analysis compared to point estimates.
Flexibility: The approach can be adapted to various model specifications, including multilevel and non-linear models.

Considerations and Limitations

Despite its strengths, there are several challenges and considerations when using Bayesian quantile regression:

Computational Complexity: Running MCMC simulations, especially for large datasets or complex models, can be time-consuming.
Single Quantile Focus: Typically, the model is built to estimate one quantile at a time. If multiple quantiles are desired, separate model runs are needed for each.
Priors Sensitivity: The choice of priors can influence the posterior estimates, necessitating careful selection and sensitivity analyses.
Custom Implementation Requirements: When not using brms, users must define the ALD function manually in Stan, which requires a deeper understanding of both the theory and Stan’s syntax.

Practical Implementation Using Stan and brms

Using Custom Stan Code

Below is a simplified outline of custom Stan code for Bayesian quantile regression:


// Data block
data {
    int N; // number of observations
    int P; // number of predictors
    real<lower=0, upper=1> tau; // quantile of interest
    vector[N] y; // response variable
    matrix[N, P] x; // predictor matrix
}
// Parameters block
parameters {
    vector[P] beta; // regression coefficients
    real<lower=0> sigma; // scale parameter
}
// Transformed parameters
transformed parameters {
    vector[N] mu = x * beta; // linear predictor for each observation
}
// Functions block to define the asymmetric Laplace density
functions {
    real asymmetric_laplace_lpdf(real y, real mu, real sigma, real tau) {
        real rho = (fabs(y - mu) + (2 * tau - 1) * (y - mu)) / 2;
        return log(tau) + log1m(tau) - log(sigma) - rho / sigma;
    }
}
// Model block
model {
    // Priors
    beta ~ normal(0, 1);
    sigma ~ exponential(1);
    // Likelihood construction using the custom ALD pdf
    for (i in 1:N) {
        y[i] ~ asymmetric_laplace(mu[i], sigma, tau);
    }
}

Note: Since Stan does not provide a built-in asymmetric Laplace distribution, you must carefully define it within the functions block and then use it within the model block.

Leveraging the brms Package in R

The brms package offers a more streamlined approach for users who prefer working in R without diving deep into Stan code. In brms, you can specify a Bayesian quantile regression model by using the asym_laplace() family. Here is an example:

# R code using brms
library(brms)

# Define formula: response ~ predictors
formula <- bf(y ~ x1 + x2)

# Fit the Bayesian quantile regression model for the median (tau = 0.5)
fit <- brm(
    formula = formula,
    data = your_data,
    family = asym_laplace(),
    prior = c(
        set_prior("normal(0,1)", class = "b"),
        set_prior("exponential(1)", class = "sigma")
    ),
    iter = 2000,
    chains = 4
)

# Summary and diagnostics
summary(fit)

brms takes care of the underlying Stan code generation and sampling, making Bayesian quantile regression accessible even to those without extensive expertise in Stan's modeling language.

Comparison Table: Custom Stan vs. brms Implementation

Feature	Custom Stan Code	brms Package
Ease of Use	Requires direct Stan programming and defining custom functions	User-friendly R interface with automatic Stan code generation
Flexibility	Highly customizable with full control over model specification	Limited to model structures supported by brms, though highly versatile
Priors Management	Manually defined in the model block	Defined through a convenient R syntax
Computational Complexity	May require more detailed tuning of sampling parameters	Handles MCMC sampling efficiently with built-in defaults
Multiple Quantiles	Requires separate models for each quantile	Same limitation; individual models per quantile are typically required

Further Considerations and Advanced Topics

Handling Multiple Quantiles

When a practitioner is interested in understanding various parts of the response distribution, it is possible to fit separate Bayesian models for each quantile of interest. Although this means running multiple MCMC samplings, it provides granular insights into the relationship between predictors and the dependent variable across different segments of the outcome distribution.

Each model can be estimated independently using either the custom Stan implementation or brms. The outputs can later be aggregated and compared to give a comprehensive quantile profile.

Interpreting Results and Diagnostic Checks

As with any Bayesian model, thorough diagnostic checks play a vital role in ensuring the model's reliability. Key diagnostic tools include assessing convergence (using Rhat values and trace plots), reviewing posterior predictive checks, and performing sensitivity analyses with different prior specifications. The output from these checks informs whether the estimated quantiles are robust and reliable for decision-making.

This step is fundamental in validating that the Bayesian quantile regression model correctly captures the underlying data structure, thus lending greater confidence in any subsequent interpretations or decisions based on its results.

Practical Applications

Use Cases in Real-World Scenarios

Bayesian quantile regression has a myriad of applications across diverse fields:

Econometrics: Examining wage distributions by analyzing the effects of education or experience at different quantile levels.
Environmental Science: Modeling extreme weather events by focusing on upper quantiles of pollution or temperature data.
Medical Statistics: Assessing the impact of medical treatments on various quantiles of patient recovery times rather than simply relying on the mean.
Finance: Analyzing risk measures by modeling quantiles of returns, which is crucial for Value-at-Risk assessments.

These applications benefit from the robustness of quantile regression and the comprehensive uncertainty quantification inherent in Bayesian methods, making it a powerful tool in data analytics and decision-making under uncertainty.