Binary and continuous quantile regression, and the anchor you have to choose

The hardware and bandwidth for this mirror is donated by METANET, the Webhosting and Full Service-Cloud Provider.
If you wish to report a bug, or if you are interested in having us mirror your free-software or open-source project, please feel free to contact us at mirror[@]metanet.ch.

The package fits the same three penalties to two kinds of response: a binary one through bbqr(), and an observed continuous one through cbqr(). Most of this vignette is about the binary case, because that is where the hard problem is. If your response is observed, skip to When the response is observed — the anchor machinery below simply does not apply to you.

The problem this package exists to solve

Quantile regression on a binary outcome has an awkward property that continuous quantile regression does not. The model is

\[y_i = \mathbb{1}\{\beta_0 + x_i'\beta + \varepsilon_i > 0\},\]

and \(\mathrm{sign}(\beta_0 + x'\beta + \varepsilon)\) does not change if you multiply \(\beta_0\), \(\beta\) and the error scale by the same positive constant. The data therefore identify only the direction of \(\beta\), never its length.

Every method has to pin the scale down somehow, and existing packages each make that choice internally without exposing it. bbqr makes it an argument, so the choice can be compared rather than assumed.

There is a second problem, separate from the first and easy to confuse with it. The published hierarchies for the penalised samplers have no posterior: their prior on the global shrinkage scale is not integrable, and a Gibbs sampler on such a target runs without complaint while producing numbers that depend on chain length and seed. bbqr ships each hierarchy in two forms and defaults to the one that exists. That is the subject of the section Two hierarchies below; identification is fixed by the anchor, propriety by the prior, and neither substitutes for the other.

Fitting a model

The interface is one function. penalty picks the shrinkage layer, anchor picks the identification restriction, and model on the prior() object picks the hierarchy.

n <- 300
X <- matrix(rnorm(n * 5), n, 5, dimnames = list(NULL, paste0("x", 1:5)))
truth <- c(1.5, -1, 0.8, 0, 0)
y <- as.numeric(X %*% truth + rnorm(n) > 0)
dat <- data.frame(y = y, X)

fit <- bbqr(y ~ x1 + x2 + x3 + x4 + x5, data = dat,
            quantile = 0.5, penalty = "alasso",
            ndraw = 2000, burn = 500)
fit
#> Bayesian binary quantile regression
#> 
#> Call:
#> bbqr(formula = y ~ x1 + x2 + x3 + x4 + x5, data = dat, quantile = 0.5, 
#>     penalty = "alasso", ndraw = 2000, burn = 500)
#> 
#>   Quantile        0.5
#>   Penalty         adaptive lasso
#>   Hierarchy       model = "v5"
#>   Identification  sigma = 1
#>   Observations    300
#>   Draws           1500 kept of 2000 (burn 500, thin 1)
#> 
#> Posterior mean coefficients:
#> (Intercept)          x1          x2          x3          x4          x5 
#>     -0.0104      3.5833     -2.3567      2.4188     -0.1176      0.4388

summary() adds credible intervals and flags the coefficients whose interval excludes zero, which is the variable-selection rule used with these samplers.

summary(fit)
#> Bayesian binary quantile regression
#> 
#> Call:
#> bbqr(formula = y ~ x1 + x2 + x3 + x4 + x5, data = dat, quantile = 0.5, 
#>     penalty = "alasso", ndraw = 2000, burn = 500)
#> 
#>   Quantile        0.5
#>   Penalty         adaptive lasso
#>   Hierarchy       model = "v5"
#>   Identification  sigma = 1
#>   Observations    300
#>   Draws           1500 kept of 2000 (burn 500, thin 1)
#>   MH acceptance   33.1%
#> 
#> Posterior summary (95% credible intervals):
#>                Mean     SD   Lower   Upper   
#> (Intercept) -0.0104 0.2109 -0.4517  0.4220   
#> x1           3.5833 0.4443  2.7531  4.4569  *
#> x2          -2.3567 0.2980 -2.9102 -1.7082  *
#> x3           2.4188 0.3673  1.6543  3.0736  *
#> x4          -0.1176 0.2137 -0.5338  0.3324   
#> x5           0.4388 0.2348 -0.0507  0.8727   
#> 
#> * interval excludes zero
#> Note: coefficients are identified only up to a positive scale;
#>       see the Identification section of ?bbqr.

Note the reminder in the footer. The coefficients are identified only up to a positive scale, so their magnitudes are only interpretable relative to the anchor you chose. Ratios of coefficients, and the sign pattern, are what carry meaning across anchors.

Three penalties

The penalties share the latent-variable augmentation and differ only in what sits above the slopes: "none" puts a vague normal prior on each slope (Benoit and Van den Poel, 2012), "lasso" uses one shrinkage parameter for all slopes (Benoit, Al-Hamzawi and Yu, 2013), and "alasso" gives each slope its own (Rubio Garcia, 2023).

pens <- c("none", "lasso", "alasso")
est <- sapply(pens, function(p)
  coef(bbqr(y ~ x1 + x2 + x3 + x4 + x5, data = dat, penalty = p,
            ndraw = 2000, burn = 500)))
round(est, 3)
#>               none  lasso alasso
#> (Intercept)  0.010  0.052  0.005
#> x1           4.236  4.058  3.883
#> x2          -2.814 -2.658 -2.489
#> x3           3.013  2.828  2.638
#> x4          -0.189 -0.197 -0.090
#> x5           0.427  0.470  0.422

The columns are not on a common scale — each sampler lands wherever its own anchor puts it — so compare directions rather than levels:

cosine <- function(a, b) sum(a * b) / sqrt(sum(a^2) * sum(b^2))
round(apply(est[-1, ], 2, cosine, truth), 4)
#>   none  lasso alasso 
#> 0.9900 0.9899 0.9914

On data this clean the three penalties are nearly indistinguishable. That is worth knowing: the penalty is not usually where the leverage is.

Two hierarchies

Each penalty is available as published (model = "v3") and in a corrected form (model = "v5", the default). The print method of prior() says which is which:

prior("alasso")
#> bbqr prior for penalty = 'alasso' , model = 'v5'  [proper omega and intercept]
#> 
#>   a            2
#>   b            2
#>   alpha_delta  2
#>   beta_delta   2
#>   alpha_omega  2
#>   beta_omega   2
#>   b0_mean      <tau-calibrated>
#>   b0_var       <tau-calibrated>
#>   delta_init   1
#>   omega_init   1
#>   mh_sd0       0.05
#>   target_acc   0.3
prior("alasso", model = "v3")
#> bbqr prior for penalty = 'alasso' , model = 'v3'  [IMPROPER: p(omega) prop 1/omega]
#> 
#>   a            1
#>   b            1
#>   alpha_delta  2
#>   beta_delta   2
#>   alpha_omega  0
#>   beta_omega   0
#>   b0_mean      0
#>   b0_var       Inf
#>   delta_init   1
#>   omega_init   1
#>   mh_sd0       0.05
#>   target_acc   0.3

The defect in the published adaptive lasso is the hyperprior \(p(\omega) \propto 1/\omega\) on the global shrinkage scale. The hierarchy is a scale family in \(\omega\): on the ray where \(\omega\), every \(\lambda_j^2\) and every latent \(s_j\) shrink to zero together with \(\beta_j = O(\sqrt{\omega})\), the prior measure is \(d\omega/\omega\) while the likelihood stays bounded below, so the normalising integral diverges. That holds for every prior on \(\delta\) and at either anchor. The published lasso carries the same defect in \(\tau_h\). Nothing about the sampler is wrong — every full conditional is a proper distribution, which is exactly why the problem never announces itself.

"v5" puts a \(\mathrm{Gamma}(2, 2)\) prior on the global scale, a proper prior on \(\delta\) where the published one was flat, and a \(\tau\)-calibrated normal prior on the intercept. Every component is then a probability distribution, the observed likelihood is bounded by one, and the posterior exists with no condition on the design, on the number of covariates, or on the anchor.

The observable signature is where \(\omega\) spends its time. Under the published prior it piles up near zero; the \(\mathrm{Gamma}(2, 2)\) density vanishes linearly at the origin and stops it:

om <- sapply(c(v3 = "v3", v5 = "v5"), function(m)
  bbqr(y ~ x1 + x2 + x3 + x4 + x5, dat, ndraw = 4000, burn = 1000,
       prior = prior("alasso", model = m))$omega)
signif(apply(om, 2, quantile, c(0, 0.01, 0.1, 0.5)), 3)
#>           v3     v5
#> 0%  2.41e-07 0.0263
#> 1%  4.25e-03 0.1640
#> 10% 1.75e-01 0.4180
#> 50% 1.97e+00 1.0900

This is a diagnostic rather than a proof — a Gibbs sampler on an improper target is null-recurrent, so a short chain can miss the drift entirely — but it is the mechanism, and on longer chains the published hierarchy’s \(\omega\) drifts without bound while "v5"’s does not.

Under "v5" the package refuses every device its derivation excludes: the numerical clamps, the sampler-argument floors, the historical inverse-Gaussian routine for the latents and the \(\exp(-\kappa/\lambda_j^2)\) tilt. Those devices existed because the published hierarchy needed them, and "v3" keeps them reachable so that earlier runs can be reproduced. The fitted object records which hierarchy it sampled in model, and whether the chain targets a derived posterior in derived.

Five anchors

Anchors, and which penalties accept them.
anchor	restriction	available
sigma1	ALD inverse-scale held at 1	all
free	none beyond the prior	lasso, alasso
beta1	first slope held at 1	all
norm1	\|\|(beta0, beta)\|\| = 1	all
normslopes	\|\|beta\|\| = 1, slopes only	all

penalty = "none" with anchor = "free" is deliberately not offered. With neither a penalty nor an anchor the likelihood is flat along the scale ray, so nothing identifies the model and only the vague prior stops the chain drifting.

"sigma1" is the default for every penalty. It closes the scale ray outright, so the coefficient scale is stated rather than left to the prior, and it is the convention under which the \(\tau\)-calibrated intercept prior of "v5" is exact. "free" samples \(\sigma\) from its full conditional, as the adaptive-lasso and lasso derivations also allow; its point estimates are as good, but its credible intervals are markedly wider, because the likelihood never identifies \(\sigma\) and the prior alone settles where the chain sits along the ray.

The two norm anchors go through the same rescaling step and differ only in how the constant is chosen. Both carry \(\beta_0\) and the ALD inverse-scale \(\sigma\) along with the slopes. At the observed-data level the model is invariant under \((\beta_0, \beta, \sigma) \mapsto (c\beta_0, c\beta, \sigma/c)\). These projections are identification experiments rather than transitions of any derivation, and under "v5" the fit says so: the projection is an exact reparameterisation of the likelihood but not of the prior, so the chain does not target a derived posterior, and the fit warns and records derived = FALSE. You can see the constraint holding exactly, and \(\sigma\) absorbing the rescaling:

f_sig  <- bbqr(y ~ x1 + x2 + x3 + x4 + x5, dat, anchor = "sigma1",
               ndraw = 2000, burn = 500)
f_norm <- bbqr(y ~ x1 + x2 + x3 + x4 + x5, dat, anchor = "norm1",
               ndraw = 2000, burn = 500)
#> Warning: model = 'v5' with anchor = 'norm1': the post-draw projection is an
#> exact reparameterisation of the likelihood but not of the prior, so this chain
#> does not target a derived posterior. Reported as derived = FALSE.

c(sigma1_scale_fixed  = all(f_sig$sigma == 1),
  norm1_unit_norm     = all(abs(sqrt(f_norm$beta0^2 +
                                     rowSums(f_norm$beta^2)) - 1) < 1e-8),
  norm1_sigma_moves   = sd(f_norm$sigma) > 0,
  sigma1_derived      = f_sig$derived,
  norm1_derived       = f_norm$derived)
#> sigma1_scale_fixed    norm1_unit_norm  norm1_sigma_moves     sigma1_derived 
#>               TRUE               TRUE               TRUE               TRUE 
#>      norm1_derived 
#>              FALSE

When the anchor matters

On easy data the anchor is nearly irrelevant. (The projection-anchor warning shown above is suppressed in the comparisons that follow; it would otherwise appear once per projected fit.)

anch <- c("sigma1", "free", "beta1", "norm1", "normslopes")
easy <- sapply(anch, function(a)
  cosine(coef(suppressWarnings(
    bbqr(y ~ x1 + x2 + x3 + x4 + x5, dat, anchor = a,
         ndraw = 2000, burn = 500)))[-1], truth))
round(easy, 4)
#>     sigma1       free      beta1      norm1 normslopes 
#>     0.9915     0.9931     0.9882     0.9923     0.9915

Now make the problem hard — an extreme quantile, heavy-tailed errors, and a badly imbalanced outcome:

lat  <- as.vector(X %*% truth) + rt(n, df = 3)
yh   <- as.numeric(lat + quantile(-lat, 0.85) > 0)
hard <- data.frame(y = yh, X)

res <- sapply(anch, function(a)
  cosine(coef(suppressWarnings(
    bbqr(y ~ x1 + x2 + x3 + x4 + x5, hard, quantile = 0.05,
         anchor = a, ndraw = 3000, burn = 750)))[-1], truth))
round(res, 4)
#>     sigma1       free      beta1      norm1 normslopes 
#>     0.9788     0.9756     0.9839     0.9809     0.9848

The ordering here is a single draw and will move from seed to seed, but the pattern is the one that holds across a full factorial study: anchors that constrain the norm of the coefficient vector degrade on hard problems, while fixing the error scale or a single coordinate does not.

The mechanism is that the rescaling step is a symmetry of the likelihood but not of the penalised posterior. The shrinkage layer is stated on a fixed scale, so renormalising the coefficients every sweep moves them relative to the penalty. The prediction that follows — and it holds — is that the norm anchors are harmless when the prior on \(\beta\) is vague, and damaging when it is informative on the scale the coefficients occupy. With penalty = "none" at its default beta_var = 100 they cost nothing; tighten that prior far enough and they break in the same way.

Practical guidance

The defaults — model = "v5", anchor = "sigma1", no clamps — sample the derived V5 posterior: proper \(\mathrm{Gamma}(2, 2)\) priors on \(\omega\) and \(\delta\), the \(\tau\)-calibrated intercept prior, and \(\sigma\) held at 1. Use anchor = "free" when you want \(\sigma\) sampled from its full conditional as the appendix also allows, and expect wider intervals. Do not set beta_delta = 0: the flat prior leaves the adaptive-lasso posterior improper, so \(\delta\) drifts and the \(\lambda_j^2\) homogenise until the penalty is no longer adaptive.
Use model = "v3" only to reproduce published results. Its posterior does not exist, the print method says so, and the numbers it produces depend on chain length and seed.
Avoid "norm1" and "normslopes" with a penalty unless you specifically want to study them. They are offered because the comparison is the point of this package, not because they are recommended, and under "v5" they are flagged as derived = FALSE.
Use standardize = TRUE whenever the covariates are on different scales. The penalty is stated on the scale of the design matrix as supplied, so unscaled covariates are shrunk by very different relative amounts; and under "v5" there is no proper prior on \(\omega\) that would make the hierarchy equivariant in the predictor scale, so standardising is what makes the prior mean the same thing across data sets.
Compare directions, not magnitudes, across anchors.

Other outputs

p <- predict(fit, newdata = dat, type = "response")
round(quantile(p, c(0, .25, .5, .75, 1)), 3)
#>    0%   25%   50%   75%  100% 
#> 0.000 0.106 0.477 0.871 0.998

mean(predict(fit, newdata = dat, type = "class") == y)
#> [1] 0.8333333

Several quantiles at once returns one fit per level:

mfit <- bbqr(y ~ x1 + x2 + x3 + x4 + x5, dat, quantile = c(0.25, 0.5, 0.75),
             ndraw = 1500, burn = 400)
round(coef(mfit), 3)
#>             tau=0.25 tau=0.50 tau=0.75
#> (Intercept)   -2.353   -0.018    2.382
#> x1             5.625    3.921    5.045
#> x2            -3.577   -2.560   -3.497
#> x3             3.592    2.662    3.943
#> x4            -0.128   -0.101   -0.338
#> x5             0.527    0.425    0.459

When the response is observed

Everything above is about the binary case, where the threshold destroys the scale and you have to choose an anchor to get it back. If the response is observed rather than thresholded, that problem disappears, and cbqr() is the entry point. The three penalties carry over, and so does the corrected hierarchy. The prior() constructor and the model switch do not: continuous fits take a plain named list and ship a single hierarchy.

cd <- data.frame(x1 = rnorm(150), x2 = rnorm(150))
cd$y <- 2 + 1.5 * cd$x1 - cd$x2 + rnorm(150)

cfit <- cbqr(y ~ x1 + x2, cd, quantile = 0.5, penalty = "alasso",
             ndraw = 1500, burn = 400)
round(coef(cfit), 3)
#> (Intercept)          x1          x2 
#>       2.006       1.509      -0.932

Three things differ, and each is forced by the response being observed rather than chosen for convenience.

No anchor. The data identify the scale, so sigma is sampled every sweep; the anchors would be misspecifications rather than competing conventions. It is a reported parameter now:

round(mean(cfit$sigma), 3)
#> [1] 2.68

A new q, whose default differs by penalty. q is the exponent on sigma in the local-scale prior rate, \(\sigma^q / (2\lambda_j^2)\). The defaults are inherited from the published penalties, not picked: the lasso penalty \(\lambda|\beta_j|\) gives q = 0 and the adaptive penalty \((\sqrt{\sigma}/\lambda_j)|\beta_j|\) gives q = 1. Holding sigma fixed, as the binary model does, makes those two the same sampler — which is why the choice never had to be faced before. q = 2 is the response-scale equivariant option.

standardize defaults to TRUE here, and scales y as well as X. No proper hierarchy is exactly response-scale equivariant, so fixing the scale is part of the specification. Draws are mapped back to the original units, so coef() is on the scale of your data either way.

Finally, predict() on a cbqr fit returns the fitted conditional quantile rather than a probability:

round(head(predict(cfit, newdata = cd), 4), 3)
#> [1] 2.242 1.910 3.071 1.093

References

Rubio Garcia, F. (2023). Bayesian adaptive lasso binary quantile regression with hybrid resampling for classification of imbalanced data. M.S. thesis, Wichita State University. https://soar.wichita.edu/entities/publication/a2f86232-4704-4ec2-b685-751e7b04ec42

Benoit, D. F. and Van den Poel, D. (2012). Binary quantile regression: a Bayesian approach based on the asymmetric Laplace distribution. Journal of Applied Econometrics, 27(7), 1174–1188.

Benoit, D. F., Al-Hamzawi, R. and Yu, K. (2013). Bayesian lasso binary quantile regression. Computational Statistics, 28(6), 2861–2873.

Kozumi, H. and Kobayashi, G. (2011). Gibbs sampling methods for Bayesian quantile regression. Journal of Statistical Computation and Simulation, 81(11), 1565–1578.

These binaries (installable software) and packages are in development.
They may not be fully stable and should be used with caution. We make no claims about them.