📊 statistics

Robust high-dimensional Bayesian regression with non-Gaussian errors under global--local shrinkage priors

This paper introduces a robust Bayesian framework for high-dimensional multivariate regression that employs scale-location mixture errors and horseshoe+ priors to simultaneously achieve sparse coefficient estimation and residual dependence graph recovery, offering superior performance over Gaussian methods in the presence of heavy tails, outliers, and skewness while maintaining efficiency under normality.

Original authors: Mohammad Arashi

Published 2026-06-09

📖 5 min read🧠 Deep dive

CC BY 4.0

Original authors: Mohammad Arashi

Original paper licensed under CC BY 4.0 (http://creativecommons.org/licenses/by/4.0/). ✨ This is an AI-generated explanation of the paper below. It is not written or endorsed by the authors. For technical accuracy, refer to the original paper. Read full disclaimer

Imagine you are trying to predict the future of a complex system, like the stock market or the economy, using a massive amount of data. You have hundreds of different "predictors" (like interest rates, oil prices, or weather) trying to explain dozens of "outcomes" (like stock returns or employment numbers).

In the world of statistics, the standard tool for this job is a method called Multivariate Regression. Think of this tool as a highly skilled but somewhat fragile detective. It's great at finding patterns, but it has two major weaknesses that this paper fixes:

It assumes everything is "Normal": The standard detective assumes that data points usually cluster neatly around an average, like people's heights. But in the real world (especially with money or genes), data is often "heavy-tailed." This means there are wild, unpredictable outliers—like a sudden market crash or a massive gene mutation—that don't fit the neat curve. When these outliers show up, the standard detective gets confused and draws the wrong conclusions.
It treats clues and connections separately: The detective looks at which predictors matter (the clues) and how the outcomes are connected to each other (the relationships) as two separate jobs. But in reality, these are deeply linked. If you get the clues wrong, you'll get the relationships wrong, too.

The Paper's Solution: A "Super-Detective" with a Flexible Mindset

The authors propose a new, robust framework (a "Super-Detective") that solves both problems at once. Here is how it works, using simple analogies:

1. The "Flexible Lens" (Handling Outliers)

Instead of assuming data follows a perfect bell curve, this new model uses a Scale-Location Mixture.

The Analogy: Imagine the standard model is a camera with a fixed focus. If a sudden, bright flash (an outlier) happens, the whole photo gets ruined.
The New Model: This model has a "smart lens" (based on the Student-t distribution). When it sees a wild, crazy data point, it doesn't panic. Instead, it says, "Okay, this point is weird, so I'll just turn down the volume on it." It assigns a low "weight" to that outlier, effectively ignoring its noise while still listening to the rest of the data. This allows it to handle heavy tails and sudden crashes without breaking.

2. The "Double-Pruned Garden" (Handling Complexity)

The model deals with two types of "noise" that need to be cut away:

Noise in the Clues: Many of the hundreds of predictors are actually useless.
Noise in the Connections: Many of the relationships between the outcomes are fake or non-existent.

The authors use a special tool called the Horseshoe+ Prior.

The Analogy: Imagine you have a garden with thousands of plants (data points). You want to keep only the rare, valuable flowers (true signals) and cut down the weeds (noise).
The Standard Tool (Lasso): This is like a lawnmower that cuts everything a little bit. It tends to accidentally trim the tips off the valuable flowers, making them smaller than they really are.
The New Tool (Horseshoe+): This is a pair of magical shears. It aggressively cuts down the tiny weeds (shrinking noise to zero) but leaves the big, valuable flowers completely untouched. It's "sharper" and more precise than the old tools, ensuring that if a signal is real, the model doesn't accidentally shrink it away.

3. Doing Two Jobs at Once

The biggest innovation is that this model does the "clue selection" and the "connection mapping" simultaneously.

The Analogy: Instead of hiring one person to find the clues and another to map the relationships, this model is a single agent who does both at the same time. Because it knows the clues are being filtered, it draws a more accurate map of how the outcomes are connected.

What the Paper Found (The Results)

The authors tested this new detective in four different ways:

When things are normal: If the data is clean and follows a bell curve, the new model works just as well as the old standard. It doesn't lose speed or accuracy.
When things are messy (Heavy Tails): When the data has wild outliers (like a financial crash), the old models get confused and produce bad maps. The new model stays calm, ignores the noise, and produces a much more accurate picture.
When things are lopsided (Skewness): Sometimes data isn't just messy; it's lopsided (like a pile of sand leaning to one side). The new model can detect this shape and adjust, whereas the old models fail.
Speed: They built two versions of the detective. One is a "slow and thorough" version (Gibbs sampler) that is perfect for smaller problems. The other is a "fast and efficient" version (Variational Inference) that is 10 to 70 times faster, making it usable for huge datasets.

Real-World Tests

The authors tested their model on two real-world datasets:

Macroeconomics (FRED-MD): They looked at economic indicators over decades. The model automatically identified the 2008 financial crisis and the 2020 pandemic as "outliers" and down-weighted them, while still finding the true relationships between economic factors.
Stock Market (S&P 500): They analyzed daily stock returns. The model successfully filtered out the "noise" of daily volatility and found a very small, clear network of how specific energy stocks were connected to each other, ignoring the noise from the rest of the market.

The Bottom Line

This paper introduces a statistical tool that is tougher (it ignores wild outliers), sharper (it finds true signals without shrinking them), and smarter (it figures out clues and connections at the same time). It proves that you don't have to choose between being robust to errors and being precise; you can have both.

Technical Summary: Robust High-Dimensional Bayesian Regression with Non-Gaussian Errors under Global–Local Shrinkage Priors

1. Problem Statement

The paper addresses the challenges of multivariate linear regression in high-dimensional settings where the number of predictors ( $p$ ) and response variables ( $q$ ) are large relative to the sample size ( $n$ ). The classical approach relies on two restrictive assumptions that are frequently violated in empirical sciences (e.g., macroeconomics, finance, genomics):

Gaussian Error Law: Standard estimators assume errors follow a multivariate normal distribution. However, real-world data often exhibit heavy tails (leptokurtosis), abrupt outliers, and asymmetry (skewness). Under Gaussian assumptions, a few gross outliers can distort coefficient estimates and, critically, corrupt the estimation of the conditional-dependence graph (the precision matrix $\Omega$ ).
Separate Sparsity Handling: Existing methods typically address sparsity in the regression coefficients ( $B$ ) and sparsity in the inverse error covariance ( $\Omega$ ) in isolation. There is a lack of integrated frameworks that simultaneously perform variable selection and conditional-dependence graph recovery while maintaining robustness to non-Gaussian errors.

2. Methodology

The authors propose a unified Bayesian framework, termed RHS-MR (Robust Horseshoe+ Multivariate Regression), which integrates robust error modeling with joint shrinkage priors.

A. Robust Likelihood via Scale–Location Mixtures
Instead of a Gaussian error, the model employs a flexible scale–location mixture of normals. Conditional on latent variables, the error structure allows for:

Multivariate $t$ -distribution: Achieved by mixing a Gaussian with an inverse-gamma weight ( $w_i$ ), accommodating heavy tails.
Multivariate Skew-Normal and Skew- $t$ : Achieved by introducing a latent half-normal vector ( $t_i$ ) to shift the mean, accommodating both heavy tails and asymmetry.
This formulation renders the conditional posteriors conjugate, facilitating tractable computation.

B. Joint Horseshoe+ Shrinkage
The model places a horseshoe+ prior simultaneously on two distinct sparse objects:

Regression Coefficients ( $B$ ): Each entry $\beta_{jk}$ is assigned a local scale parameter, augmented by a second half-Cauchy factor ( $\eta_{jk}$ ) compared to the standard horseshoe. This creates a marginal prior with a taller spike at zero (aggressive noise shrinkage) and heavier tails (near-unbiased treatment of large signals).
Precision Matrix ( $\Omega$ ): The off-diagonal entries (representing conditional dependencies) receive a horseshoe+ prior, while diagonal entries receive a diffuse exponential prior. This couples the sparsity of the regression map with the sparsity of the conditional-dependence graph.

C. Posterior Computation
Two algorithms are developed for posterior inference:

Gibbs Sampler: A fully tractable sampler derived from the scale-mixture and parameter-expanded representations. It utilizes a column-wise block update for the precision matrix (inspired by the graphical horseshoe) and closed-form updates for all full conditionals.
Mean-Field Variational Inference: A faster approximation algorithm that factorizes across blocks, offering a speed-up of one to two orders of magnitude with comparable accuracy in well-identified regimes.

3. Key Contributions

The paper makes four primary contributions:

Model Formulation: It formulates a robust multivariate regression model where the Gaussian error is replaced by a scale–location mixture (covering $t$ , skew-normal, and skew- $t$ laws), all sharing a common precision matrix $\Omega$ .
Joint Shrinkage: It introduces the simultaneous application of the horseshoe+ prior to both $B$ and the off-diagonal entries of $\Omega$ , effectively coupling variable selection with graph structure recovery.
Computational Algorithms: It derives a closed-form Gibbs sampler and a scalable variational inference algorithm for these complex models.
Theoretical Guarantees:
- Joint Posterior Contraction: Establishes that the posterior distributions of $B$ and $\Omega$ contract to the truth at near-minimax rates, even under non-Gaussian errors.
- Selection Consistency: Proves consistency for both variable selection (in $B$ ) and edge selection (in $\Omega$ ).
- Risk Dominance: Demonstrates via a Kullback–Leibler risk bound that the horseshoe+ prior dominates the standard horseshoe in the sparse-signal regime due to its heavier tails reducing bias on large signals.
- Bounded Sensitivity: Proves that under the $t$ -error, the posterior mean of $B$ remains bounded even as a single outlier tends to infinity, whereas Gaussian estimators diverge linearly.

4. Empirical Results

The methodology is evaluated through extensive simulations and two real-world applications.

Simulation Study:

Error Regimes: The model was tested under Gaussian, $t_3$ , skew-normal, and skew- $t_4$ error distributions.
Performance: Under Gaussian errors, the robust model performs indistinguishably from its Gaussian counterpart. Under heavy-tailed and skewed errors, the robust model significantly outperforms Gaussian-based benchmarks (including penalized likelihood methods like MRCE and OLS). Specifically, under skew- $t_4$ errors, the robust model reduced coefficient estimation error by 43% compared to the Gaussian ablation.
Contraction: As sample size ( $n$ ) increased, the robust estimator's error contracted at the predicted rate ( $\sim n^{-1/2}$ ), while the Gaussian estimator's precision matrix error failed to contract under heavy tails.
Outlier Sensitivity: When a single observation was contaminated with increasing magnitude, the Gaussian estimator's error grew linearly, while the robust estimator's error remained flat, confirming the bounded sensitivity theorem.
Computation: The variational algorithm achieved speed-ups of 10x to 70x over the Gibbs sampler with minimal loss in accuracy.

Applications:

FRED-MD (Macroeconomic Data): Applied to 8 macroeconomic indicators and 30 predictors. The model automatically down-weighted observations during the 2008 financial crisis and the 2020–2021 pandemic, recovering a sparse residual dependence graph and achieving lower predictive error than the Gaussian fit.
S&P 500 (Equity Returns): Applied to daily log-returns of 12 stocks and 31 factors. The robust fit identified a sparse residual graph where 4 of 5 selected edges connected Energy sector stocks, reflecting common exposure to commodity prices. It also down-weighted volatile trading days, including the onset of the 2007–2008 turmoil.

5. Significance and Claims

The paper claims to provide the first integrated Bayesian framework that simultaneously addresses high-dimensional sparsity in both regression coefficients and conditional-dependence structures while maintaining robustness to heavy-tailed and asymmetric errors.

The authors emphasize that their approach:

Preserves Efficiency: It does not sacrifice performance when errors are Gaussian (robustness costs "essentially nothing").
Corrects Structural Bias: It prevents the corruption of the conditional-dependence graph ( $\Omega$ ) by outliers, a failure mode common in Gaussian procedures.
Theoretical Soundness: It offers rigorous guarantees on posterior contraction, selection consistency, and bounded sensitivity, bridging the gap between robust statistics and high-dimensional Bayesian shrinkage.
Practical Utility: The automatic down-weighting of outliers via posterior weights provides an interpretable diagnostic for data quality without requiring pre-specified crisis indicators.

The paper concludes that the integration of global–local shrinkage with scale-mixture error models offers a superior alternative for empirical sciences where data non-normality and high dimensionality coexist.