Quantile Regressions in Stan: Part I

Quantile regressions are relatively new in the Bayesian literature first appearing in Yu and Moyeed (2001) who described an asymmetric Laplace distribution to estimate the conditional quantiles. This post is going to show the asymmetric Laplace and the augmented version. In part II I will show the score likelihood method.

The asymmetric Laplace distribution is given as

\[ \operatorname{ALD}(Y \mid x, \beta) = \frac{\tau^n (1 - \tau)^n}{\sigma^n}\exp \bigg( -\frac{\sum_{i=1}^n \rho_\tau (y_i - x_i^{\top}\beta)}{\sigma} \bigg) \]

where

\[ \rho_\tau(u) = \frac{|u| + (2\tau - 1)u}{2} \]

We may write this in Stan as

functions{
real q_loss(real q, vector u){
  return 0.5 * sum(abs(u) + (2 * q - 1) * u);
}

real ald_lpdf(vector y, real q, real sigma, vector q_est){
  int N = num_elements(y);
  
  return N * (log(q) + log1m(q) - log(sigma)) - q_loss(q, y - q_est) / sigma;
}
}
data {
  int N;                   // Number of observation
  int P;                   // Number of predictors
  real<lower=0, upper=1> q;
  vector[N] y;             // Response variable sorted
  matrix[N, P] x;
}
parameters {
  vector[P] beta;
  real<lower=0> sigma;
}
model {
  beta ~ normal(0, 4);
  sigma ~ exponential(1);
  y ~ ald(q, sigma, x * beta);
}

Let’s compare the 5th percentile to the Brq R package

library(Brq)
data("ImmunogG")
out_asym_laplace <- asymetric_laplace$sample(
  data = list(N = nrow(ImmunogG),
              P = 3,
              q = 0.05,
              y = ImmunogG$IgG,
              x = as.matrix(data.frame(alpha = 1, x = ImmunogG$Age, xsq = ImmunogG$Age^2))),
  seed = 12123123,
  parallel_chains = 4,
  iter_warmup = 500,
  iter_sampling = 500,
  refresh = 0,
  show_messages = FALSE,
  show_exceptions = FALSE
)
out_asym_laplace$summary()

# A tibble: 5 × 10
  variable     mean   median     sd    mad       q5       q95  rhat ess_bulk
  <chr>       <dbl>    <dbl>  <dbl>  <dbl>    <dbl>     <dbl> <dbl>    <dbl>
1 lp__     -671.    -670.    1.54   1.28   -673.    -669.      1.00     606.
2 beta[1]     0.599    0.617 0.255  0.278     0.152    0.985   1.00     380.
3 beta[2]     1.31     1.30  0.217  0.207     0.970    1.67    1.01     360.
4 beta[3]    -0.155   -0.153 0.0369 0.0329   -0.217   -0.0993  1.01     391.
5 sigma       0.165    0.164 0.0100 0.0100    0.150    0.182   1.01     738.
# ℹ 1 more variable: ess_tail <dbl>

y = ImmunogG$IgG
x = ImmunogG$Age
X=cbind(x, x^2)
fit = Brq(y ~ X , tau= 0.05,runs= 2000, burn=1000)
print(summary(fit))

Call:
Brq.formula(formula = y ~ X, tau = 0.05, runs = 2000, burn = 1000)

tau:[1] 0.05

            Estimate L.CredIntv  U.CredIntv
Intercept  0.6201064  0.1113146  1.03475824
Xx         1.2991424  0.8962901  1.68272470
X         -0.1550870 -0.2317421 -0.08612944

Using the covariates for extra information

In Lee, Chakraborty, and Su (2022) they detail a Bayesian approach to the envelope quantile regression. In all honesty, I briefly skimmed the paper and gathered that by adding a regression on the covariates and using this information in the asymmetric Laplace can help identification. This isn’t a recreation of their method but a quick detour to see if a simplified version of what they do increases ESS.

The idea is that we use the variability in the covariance of the regressors to help identify the coefficients on the ALD regression. I now have to separate out the intercept into \(\alpha\).

functions{
real q_loss(real q, vector u){
  return 0.5 * sum(abs(u) + (2 * q - 1) * u);
}

real ald_lpdf(vector y, real q, real sigma, vector q_est){
  int N = num_elements(y);
  
  return N * (log(q) + log1m(q) - log(sigma)) - q_loss(q, y - q_est) / sigma;
}
}
data {
  int N;                   // Number of observation
  int P;                   // Number of predictors
  real<lower=0, upper=1> q;
  vector[N] y;             // Response variable sorted
  matrix[N, P] x;
}
transformed data {
  array[N] vector[P] X;
  
  for (i in 1:P) {
    X[:, i] = to_array_1d(x[:, i]);
  }
}
parameters {
  vector[P] eta;
  real<lower=0> sigma;
  vector[P] mu;
  cholesky_factor_corr[P] L;
  vector<lower=0>[P] gamma;
  real alpha;
}
transformed parameters {
  vector[P] beta = eta .* gamma;
}
model {
  eta ~ normal(0, 4);
  sigma ~ exponential(1);
  y ~ ald(q, sigma, alpha + x * beta);
  alpha ~ normal(0, 4);
  
  gamma ~ exponential(1);
  mu ~ normal(0, 4);

  X ~ multi_normal_cholesky(mu, diag_pre_multiply(gamma, L));
}

Let’s fit the model.

out_envelope <- asymetric_laplace_envelope$sample(
  data = list(N = nrow(ImmunogG),
              P = 2,
              q = 0.05,
              y = ImmunogG$IgG,
              x = as.matrix(data.frame(x = ImmunogG$Age, xsq = ImmunogG$Age^2))),
  seed = 12123123,
  parallel_chains = 4,
  iter_warmup = 500,
  iter_sampling = 500,
  refresh = 0,
  show_messages = FALSE,
  show_exceptions = FALSE
)

out_envelope$summary(c("alpha", "beta", "sigma"))

# A tibble: 4 × 10
  variable   mean median      sd     mad     q5    q95  rhat ess_bulk ess_tail
  <chr>     <dbl>  <dbl>   <dbl>   <dbl>  <dbl>  <dbl> <dbl>    <dbl>    <dbl>
1 alpha     0.614  0.637 0.244   0.252    0.204  0.981  1.00     913.    1012.
2 beta[1]   1.29   1.29  0.203   0.201    0.982  1.62   1.00     863.     966.
3 beta[2]  -0.152 -0.151 0.0345  0.0329  -0.213 -0.102  1.00     965.     946.
4 sigma     0.166  0.165 0.00976 0.00983  0.151  0.182  1.00    1211.    1075.

Slightly different estimates, higher ESS, slower to fit, but lower standard deviations! Neat.

Scale-mixture representation

The above may also be written as a mixture of exponential and normal distributions. Letting

\[ \begin{aligned} z_i &\sim \operatorname{Exp}(1) \\ \sigma &\sim \operatorname{Exp}(1) \\ y_i &\sim \mathcal{N}(x_i^{\top}\beta_p + \sigma \theta z_i,\, \tau \sigma \sqrt{z_i}) \end{aligned} \] where

\[ \begin{aligned} \theta &= \frac{1 - 2 q}{q (1 - q)} \\ \tau &= \sqrt{\frac{2}{q(1 -q)}} \end{aligned} \]

The following Stan model implements the augmented approach.

data {
  int N;                   // Number of observation
  int P;                   // Number of predictors
  real<lower=0, upper=1> q;
  vector[N] y;             // Response variable sorted
  matrix[N, P] x;
}
transformed data {
  real theta = (1 - 2 * q) / (q * (1 - q));
  real tau = sqrt(2 / (q * (1 - q)));
}
parameters {
  vector[P] beta;
  vector<lower=0>[N] z;
  real<lower=0> sigma;
}
model {
  beta ~ normal(0, 2);
  sigma ~ exponential(1);
  
  // Data Augmentation
  z ~ exponential(1);
  y ~ normal(x * beta + sigma * theta * z, tau * sqrt(z) * sigma);
}

Let’s fit the model

out_aug_laplace <- augmented_laplace$sample(
  data = list(N = nrow(ImmunogG),
              P = 3,
              q = 0.05,
              y = ImmunogG$IgG,
              x = as.matrix(data.frame(alpha = 1, x = ImmunogG$Age, xsq = ImmunogG$Age^2))),
  seed = 12123123,
  parallel_chains = 4,
  iter_warmup = 500,
  iter_sampling = 500,
  refresh = 0,
  show_messages = FALSE,
  show_exceptions = FALSE
)

Warning: 57 of 2000 (3.0%) transitions ended with a divergence.
See https://mc-stan.org/misc/warnings for details.

out_aug_laplace$summary(c("beta"))

# A tibble: 3 × 10
  variable   mean median     sd    mad     q5     q95  rhat ess_bulk ess_tail
  <chr>     <dbl>  <dbl>  <dbl>  <dbl>  <dbl>   <dbl> <dbl>    <dbl>    <dbl>
1 beta[1]   0.606  0.619 0.243  0.260   0.193  0.990   1.01     303.     547.
2 beta[2]   1.30   1.30  0.206  0.198   0.968  1.62    1.01     335.     492.
3 beta[3]  -0.154 -0.154 0.0350 0.0325 -0.211 -0.0979  1.01     390.     669.

Ouch! Divergences. Although the parameter estimates look comparable. Let’s inspect the pairs plots between the augmented version and the asymmetric Laplace.

library(bayesplot)

This is bayesplot version 1.15.0

- Online documentation and vignettes at mc-stan.org/bayesplot

- bayesplot theme set to bayesplot::theme_default()

   * Does _not_ affect other ggplot2 plots

   * See ?bayesplot_theme_set for details on theme setting

source("../../../R/theme_blog.R")
use_blog_bayesplot()
pairs_np_style <- bayesplot::pairs_style_np(
  div_color = blog_colors$red,
  td_color = blog_colors$gold
)

np_aug_laplace <- nuts_params(out_aug_laplace)
np_asym_laplace <- nuts_params(out_asym_laplace)
mcmc_pairs(out_aug_laplace$draws(c("beta", "sigma")), np = np_aug_laplace, np_style = pairs_np_style)

mcmc_pairs(out_asym_laplace$draws(c("beta", "sigma")), np = np_asym_laplace, np_style = pairs_np_style)

The pairs plots don’t appear to show the divergence originating be due to beta or sigma. It’s most likely the augmentation variables z. As I’m up for time on this, I’d be curious to test out alternative specifications or see if anyone has encountered these issues with the augmented version and found a solution.

Try Student T

Warning

The Student T won’t result in the same posterior and is likely not ALD. However, it approaches a normal distribution with increasing \(\nu\).

Let’s use the student_t with \(\nu = 30\)

data {
  int N;                   // Number of observation
  int P;                   // Number of predictors
  real<lower=0, upper=1> q;
  vector[N] y;             // Response variable sorted
  matrix[N, P] x;
}
transformed data {
  real theta = (1 - 2 * q) / (q * (1 - q));
  real tau = sqrt(2 / (q * (1 - q)));
}
parameters {
  vector[P] beta;
  vector<lower=0>[N] z;
  real<lower=0> sigma;
}
model {
  beta ~ normal(0, 2);
  sigma ~ exponential(1);
  
  // Data Augmentation
  z ~ exponential(1);
  y ~ student_t(30, x * beta + sigma * theta * z, tau * sqrt(z) * sigma);
}

Let’s fit the same model

out_student_aug_laplace <- student_augmented_laplace$sample(
  data = list(N = nrow(ImmunogG),
              P = 3,
              q = 0.05,
              y = ImmunogG$IgG,
              x = as.matrix(data.frame(alpha = 1, x = ImmunogG$Age, xsq = ImmunogG$Age^2))),
  seed = 12123123,
  parallel_chains = 4,
  iter_warmup = 500,
  iter_sampling = 500,
  refresh = 0,
  show_messages = FALSE,
  show_exceptions = FALSE
)

Warning: 1 of 2000 (0.0%) transitions ended with a divergence.
See https://mc-stan.org/misc/warnings for details.

out_student_aug_laplace$summary(c("beta"))

# A tibble: 3 × 10
  variable   mean median     sd    mad     q5    q95  rhat ess_bulk ess_tail
  <chr>     <dbl>  <dbl>  <dbl>  <dbl>  <dbl>  <dbl> <dbl>    <dbl>    <dbl>
1 beta[1]   0.604  0.647 0.229  0.214   0.178  0.927  1.01     358.     462.
2 beta[2]   1.34   1.33  0.164  0.153   1.10   1.63   1.01     371.     733.
3 beta[3]  -0.156 -0.154 0.0267 0.0252 -0.202 -0.116  1.01     388.     715.

Great no divergences!

Conclusion

For the next post, we’ll take a look at an approximate likelihood using the data which asymptotically approaches the true likelihood with more data.

References

Lee, Minji, Saptarshi Chakraborty, and Zhihua Su. 2022. “A Bayesian Approach to Envelope Quantile Regression.” Statistica Sinica, January. https://doi.org/10.5705/ss.202020.0109.

Yu, Keming, and Rana A. Moyeed. 2001. “Bayesian Quantile Regression.” Statistics & Probability Letters 54 (4): 437–47. https://doi.org/https://doi.org/10.1016/S0167-7152(01)00124-9.