The exact contents of citations.db main_text.text for this paper — one flattened LaTeX string, title through conclusion, appendix excluded, unmodified except for removing email addresses. This is what our citation measures are computed over.
128,817 characters
Debiased Bayesian Inference for High-dimensional Regression Models
\pdfbookmark[1]{Title}{title}
\title{Debiased Bayesian Inference for High-dimensional Regression Models}
\author{Qihui Chen\thanks{School of Management and Economics and Shenzhen Finance Institute, The Chinese University of Hong Kong, Shenzhen (CUHK-Shenzhen); [email removed]}\\CUHK-Shenzhen
\and
Zheng Fang\thanks{Department of Economics, Emory University; [email removed]} \\Emory University
\and
Ruixuan Liu\thanks{CUHK Business School, Chinese University of Hong Kong; [email removed]} \\CUHK}
\date{\today}
\maketitle
\begin{abstract}
There has been significant progress in Bayesian inference based on sparsity-inducing (e.g., spike-and-slab and horseshoe-type) priors for high-dimensional regression models. The resulting posteriors, however, in general do not possess desirable frequentist properties, and the credible sets thus cannot serve as valid confidence sets even asymptotically. We introduce a novel debiasing approach that corrects the bias for the entire Bayesian posterior distribution. We establish a new Bernstein-von Mises theorem that guarantees the frequentist validity of the debiased posterior. We demonstrate the practical performance of our proposal through Monte Carlo simulations and two empirical applications in economics.
\end{abstract}
\newpage
\section{Introduction}
Applied researchers now routinely work with regression models that feature a large number of covariates. A primary inferential goal in econometrics is to estimate the ceteris paribus effect of a specific variable while controlling for other variables \citep{BelloniChernoHansen2013HD, BelloniChernozhukovCherverikovWei2018Many}. The prevailing practice interprets the coefficient on a regressor as a causal effect, conditional on the included controls. As the plausibility of conditional unconfoundedness is often argued using a large set of covariates, practitioners have increasingly embraced high-dimensional regression models. This setting has been extensively studied, predominantly using frequentist methods.
Bayesian inference, on the other hand, has long been valued for its coherent framework for handling uncertainty in statistical analysis. As highlighted by \citet{Rubin1984Applied}, Bayesian methods provide direct answers to many empirical questions by quantifying uncertainty about unknown parameters conditional on the observed data.\footnote{In discussing the commonly-used rule-of-thumb confidence interval, \cite{Rubin1984Applied} states that ``the interval is-at least in my experience-nearly always interpreted Bayesianly, that is, as providing a fixed observed interval in which the unknown (parameter of interest) $\mu$ lies with 95\% probability.'' As advocated by \cite{Imbens2021Pvalue}, it may even be preferable to adopt Bayesian inference in cases where Bayesian and frequentist procedures lead to different conclusions.} This appeal has grown among applied researchers, who often seek probabilistic statements about particular parameters of interest given the specific dataset at hand. Bayesian posterior distributions conveniently encapsulate both sampling variation and parametric uncertainty, unifying estimation and inference in a way that aligns well with the needs of empirical work.
In addressing the challenges posed by high dimensionality, recent methodological advances have substantially expanded the Bayesian toolkit for regression models with many covariates. Notable progress includes the development of spike-and-slab priors \citep{MitchellBeauchamp_BayesianVariable_1988,GeorgeMcCulloch_VariableSelection_1993} and horseshoe priors \citep{CarvalhoPolsonScott2010Horseshoe}, as well as scalable approximate Bayesian inference techniques \citep{RaySzabo2022Variational}. These innovations have attracted a growing interest in economics, as illustrated by \citet{GiannoneLenzaPrimiceri2021Sparsity}. However, while much of this literature focuses on model selection and estimation, relatively less attention has been devoted to delivering valid inference for a particular parameter of interest in the presence of high-dimensional controls---a central need in many empirical applications. This paper aims to fill this gap.
Specifically, we develop a novel inferential procedure for high-dimensional linear regression that introduces a debiasing step for the Bayesian posterior distribution of the parameter of interest. Building on the concept of debiased point estimators, which is well established in the frequentist literature, we instead debias the entire posterior distribution, yielding a new debiased Bayesian approach. Our method is tailored to high-dimensional settings and constructs credible sets from the debiased posterior, while remaining faithful to Bayesian principles by conditioning on the observed data.
The core theoretical contribution of our work is the establishment of new Bernstein-von Mises (BvM) results for the debiased Bayesian procedure. Our framework allows the number of covariates $p$ to exceed the sample size $n$, with $p$ growing exponentially in $n$. These results formally justify the asymptotic normality of the debiased posterior and ensure that credible sets achieve correct frequentist coverage under repeated sampling, thereby helping to bridge the gap between Bayesian and frequentist inference. Importantly, our results show that, after debiasing, Bayesian credible sets can match the performance of debiased frequentist methods---such as double machine learning \citep{ChernozhukovChetverikovDemirerDufloHansenNeweyRobins2018Double}---in high-dimensional regimes. This effort also echoes an important point made by \citet{Rubin_BayesianBootstrap_1981} regarding the assessment of Bayesian procedures under repeated sampling. Our framework is notably general. It is agnostic to the choice of prior for the regression coefficients and permits flexible selection of pilot estimators for the precision matrix. The specific implementations we present in Section~\ref{Sec:42} illustrate how the method can be practically adapted to high-dimensional settings. This flexibility allows the procedure to be tailored to a wide range of empirical applications without sacrificing computational scalability.
Another key advantage of our approach is its compatibility with computationally efficient approximate Bayesian inference, such as variational Bayes \citep{RaySzabo2022Variational}. Traditional spike-and-slab priors, while theoretically appealing, often entail substantial computational cost due to the reliance on Markov chain Monte Carlo (MCMC) sampling. In contrast, variational Bayes delivers approximate posteriors within seconds, while still attaining optimal contraction rates in high-dimensional settings. However, without debiasing, credible sets constructed from such approximate posteriors are not guaranteed to achieve correct frequentist coverage. By incorporating our debiasing procedure, the variational Bayes posterior becomes a principled tool for valid uncertainty quantification.
Our approach parallels an important line of work in the frequentist literature, in particular the debiasing of LASSO-type estimators \citep{JavanmardMontanari2014CI,VandeGeer_OnAsymptotically_2014,ZhangZhang2014CI}. The LASSO estimator itself can be viewed as the posterior mode under a Laplace prior, yet it typically exhibits bias and therefore requires post-estimation correction. Whereas debiased frequentist methods adjust a point estimator and then rely on its asymptotic properties for inference, our Bayesian debiasing method extends this idea to the entire posterior distribution. Because the debiasing step is applied at the posterior level, we develop new technical arguments to establish its large-sample properties. Monte Carlo simulations show that our procedure delivers substantial improvements over standard (uncorrected) Bayesian methods and remains competitive with debiased frequentist procedures in terms of coverage and estimation accuracy. In particular, when regression coefficients are large, our method offers pronounced advantages over debiased frequentist approaches, providing practitioners with a powerful tool for empirical analysis in high-dimensional settings.
Our work contributes to the active literature on the frequentist validity of Bayesian inference in high-dimensional regression. Early contributions by \citet{Ghosal1999HDL} and \citet{Bontemps2011BvM} established asymptotic normality of the posterior in Gaussian linear regression when the number of covariates grows slowly with the sample size. However, these results do not accommodate sparsity and require the ambient dimension to remain smaller than the sample size (see \citet[p.~2563]{Bontemps2011BvM}). More recently, \citet{Castilloetal_BayesianSparse_2015} analyzed Bayesian posteriors under spike-and-slab priors, establishing mixed normality under strong signal conditions; see also \citet{WuNarisettyYang2023ConditionalBayes} for extensions. \citet{Yang2019HDL} introduced a debiasing approach based on reparameterization and targeted prior modification, assigning a data-dependent Gaussian prior to the parameter of interest. Building on this idea, \citet{Castilloetal2024VB} incorporated mean-field approximation strategies to distinguish between high-dimensional nuisance parameters and parameters of interest.
In contrast, our approach attains debiasing via a post-processing step applied to the posterior, preserving standard Bayesian modeling and computation and allowing practitioners to directly leverage existing Bayesian methods for high-dimensional regression.
Our paper also connects with broader developments in semiparametric Bayesian inference and correction methods for posterior or prior distributions \citep{RayVaart2020Semi,BreunigLiuYu2022DR,BreunigLiuYu2024DiD,BreunigLiuYu2025Bart,YiuFongHolmesRousseau2023Semiparametric}. These works depart from naive plug-in principles to achieve robust large-sample properties under weak conditions. We contribute to this literature by explicitly addressing high-dimensional regression with sparsity-inducing priors and by showing that debiasing can be performed on approximate posteriors. Theoretically, we also establish BvM results for a growing subset of coefficients, which is a novel contribution in this context. Recently, \citet{DiTragliaLiu2025BML} proposed a Bayesian analog of double machine learning for linear regression when the number of covariates is relatively small, in the sense that $p = o(n)$. Their approach relies on priors for reduced-form covariance matrices to remain consistent with the likelihood principle \citep{Walker_Parametrization_2023}. However, when $p$ is large relative to $n$, practical computation of posteriors for high-dimensional covariance matrices remains challenging.
The remainder of the paper is organized as follows. Section~\ref{Sec: Debiased Bayesian inferene} introduces our inferential procedure and discusses implementation details, while Section~\ref{Sec: Applications} presents two empirical applications. In Section~\ref{Sec: 4}, we develop our main theoretical results under high-level conditions, which are then verified under more primitive assumptions. Section~\ref{Sec: Monte Carlo} examines the finite-sample performance of our procedure through Monte Carlo simulations. Section~\ref{Sec: Conclusion} concludes. All proofs and supporting results are relegated to the appendices, which also contain additional simulation evidence.
\section{Debiased Bayesian Inference}\label{Sec: Debiased Bayesian inferene}
Consider the prototypical linear regression model:
\begin{align}\label{Eqn:Model}
Y_i=X_i^\intercal\beta_0+\varepsilon_i,\quad \mathbb{E}_0[X_{i}\varepsilon_i]=0, \quad i = 1,\ldots, n,
\end{align}
where $\{Y_i,X_i\}_{i=1}^n$ are independent and identically distributed (\textit{i.i.d.}) and $\mathbb{E}_0[\cdot]$ is the population expectation.
The covariate $X_i$ is a $p$-dimensional vector, with $p$ potentially exceeding the sample size $n$. The true regression coefficient vector $\beta_0$ is assumed to be sparse, in the sense that only a small subset of regressors have nonzero effects. We denote the population covariance matrix of $X_i$ by $\Omega_0:=\mathbb{E}_0[X_iX_i^\intercal]$, and its inverse, the precision matrix, by $\Theta_0:=\Omega_0^{-1}$.
Since the seminal work of \citet{Tibshirani_Lasso_1996}, the LASSO has become the default method for estimating the high-dimensional linear regression model, owing to its ability to perform model selection and coefficient estimation simultaneously. It is well known that the LASSO estimator, as a penalized least-squares estimator, can be interpreted as a Bayesian point estimator---specifically, the posterior mode under independent Laplace priors on the regression coefficients (see Remark \ref{Rem:DebiasPoint}). However, when one considers the entire posterior distribution rather than only its mode, the Laplace prior alone fails to induce sparsity: the resulting posterior remains non-sparse and does not contract toward the true parameter at the optimal rate, unlike its mode. As shown in Theorem 7 of \citet{Castilloetal_BayesianSparse_2015}, the full Bayesian posterior corresponding to the Laplace prior lacks desirable asymptotic properties.\footnote{\citet{Castilloetal_BayesianSparse_2015} remark that “the LASSO is essentially non-Bayesian, in the sense that the corresponding full posterior distribution is a useless object” (p.1988). In practice, sampled posterior distributions of regression coefficients under the Laplace prior are indeed non-sparse \citep{Baietal_ReviewSpikeSlab_2021}.}
Our method is inspired by the frequentist literature on debiased inference \citep{ZhangZhang2014CI}, which applies a debiasing step to the LASSO estimator. Let $\hat{\beta}^{\text{Pilot}}_n$ denote any pilot estimator of $\beta_0$, and let $\hat{\Theta}_n$ be a pilot estimator of the precision matrix $\Theta_0$. The celebrated debiased inference procedure is based on the following construction:
\begin{align}
\hat{\beta}_n^{\text{Debias}}=\hat{\beta}_n^{\text{Pilot}}+ \hat{\Theta}_n\left[\frac{1}{n}\sum_{i=1}^{n}X_i(Y_i-X_i^\intercal\hat{\beta}_n^{\text{Pilot}})\right].
\end{align}
It is well established that the individual coordinates of the debiased estimator $\hat{\beta}_n^{\text{Debias}}$ are asymptotically normal, enabling valid statistical inference.
Our debiased Bayesian inference procedure departs from the frequentist approach in several key respects. The Bayesian modeling perspective treats the unknown regression coefficient as the random vector $\beta$. We begin by assigning a sparsity-inducing prior, such as a spike-and-slab prior or a horseshoe-type prior, and computing the (approximate) initial posterior distribution of $\beta$, denoted by $\Pi_{\beta}(\,\cdot\,|\{Y_i,X_i\}_{i=1}^n)$. Under sparsity, the resulting posterior achieves (near-)optimal contraction rates \citep{Castilloetal_BayesianSparse_2015,GaoVaartZhou2020General,SongLiang_NearlyOptimal_2023}. We then debias the entire posterior distribution using a correction term analogous to the frequentist construction. In this step, however, we replace the uniform weights $1/n$ with Bayesian bootstrap weights. Specifically, the Bayesian bootstrap weights are defined by normalizing independent standard exponential random variables (independent of $\{Y_i,X_i\}_{i=1}^n$ and the prior of $\beta$):
\begin{align}
W_{ni}=\frac{\omega_i}{\sum_{j=1}^n \omega_j},\quad \{\omega_j\}_{j=1}^{n} \stackrel{i.i.d.}{\sim} \textup{Exp}(1)\quad, i=1,\ldots, n.
\end{align}
In summary, we study the debiased posterior of the form:
\begin{align}\label{Eqn:DebiasBayes}
\tilde{\beta}
= \beta + \hat{\Theta}_n \left[ \sum_{i=1}^n W_{ni} X_i( Y_i - X_i^\intercal \beta ) \right],
\end{align}
where $\beta\sim \Pi_{\beta}(\,\cdot\,|\{Y_i,X_i\}_{i=1}^n)$ and $\hat{\Theta}_n$ is as introduced above. The pseudo-code for generating draws from the debiased posterior is provided in Algorithm~\ref{Algorithm}.
\begin{algorithm}[H]
\caption{Debiased Bayesian Procedure}\label{Algorithm}
\begin{algorithmic}
\STATE \textbf{Input:} Data $\{Y_i,X_i\}_{i=1}^n$; number of posterior draws $B$.
\STATE \textbf{Prior Specification:} Select a sparsity-inducing prior for $\beta$.
\STATE \textbf{Pilot Estimation:} Choose a pilot estimator $\hat{\Theta}_n$ of the precision matrix $\Theta_0$.
\STATE \textbf{Posterior Computation:}
\FOR{$b=1,\ldots, B$}
\STATE (a) Draw $\beta^b$ from the initial posterior $\Pi_{\beta}(\,\cdot\,|\{Y_i,X_i\}_{i=1}^n)$.
\STATE (b) Generate Bayesian bootstrap weights:
\[W^b_{ni}=\frac{\omega^b_i}{\sum_{j=1}^n \omega^b_j},\quad \{\omega_j^b\}_{j=1}^{n} \stackrel{i.i.d.}{\sim} \textup{Exp}(1), \quad i=1,\ldots,n.\]
\STATE (c) Compute the debiased posterior draw:
\[\tilde{\beta}^b=\beta^b+ \hat{\Theta}_n\left[\sum_{i=1}^{n}W_{ni}^bX_i(Y_i-X_i^\intercal\beta^b)\right].\]
\ENDFOR
\STATE \textbf{Output: $\{\tilde{\beta}^{b}:b=1,\ldots,B\}$}
\end{algorithmic}
\end{algorithm}
Several remarks are in order. First, the weight replacement is essential. To see this, consider the case $p<n$, where one may take $\hat{\Theta}_n =(\sum_{i=1}^{n} X_{i}X_{i}^\intercal/n)^{-1}$. Without the weight replacement, the expression simplifies to $\tilde{\beta} = (\sum_{i=1}^{n}X_{i}X_{i}^{\intercal})^{-1}\sum_{i=1}^{n}X_{i}Y_{i}$, which exhibits no posterior uncertainty. Second, the debiasing step introduces no additional computational burden. Since the Bayesian bootstrap weights are independent of the posterior draws of $\beta$, steps (a) and (b) of Algorithm \ref{Algorithm} can be parallelized efficiently to save time. Moreover, the one-time computation of $\hat{\Theta}_n$ can be carried out simultaneously, which remains computationally efficient as known in the frequentist literature (see the additional discussion in Section \ref{Sec:43}). Third, step (c) of Algorithm~\ref{Algorithm} admits a natural interpretation as a
Bayesian version of the residual bootstrap. Specifically, consider the case $p<n$ with $\hat{\Theta}_n =(\sum_{i=1}^{n} X_{i}X_{i}^\intercal/n)^{-1}$. For a given posterior draw $\beta^b$, define bootstrap residuals $\varepsilon_i^{b\ast}$ as $\varepsilon_i^{b}=Y_i-X_i^\intercal \beta^b$ weighted by the Bayesian bootstrap
weights $W_{ni}^b$. Construct the pseudo responses $Y_{i}^{b\ast} = X_i^{\intercal}\beta^b + \varepsilon^{b\ast}$ and regress $Y_{i}^{b\ast}$ on $X_i$. The resulting estimator coincides with $\tilde{\beta}^{b}$. Thus, each debiased posterior draw $\tilde\beta^b$ can be viewed as the regression coefficient obtained from a ``Bayesian residual bootstrap,''
where residuals are reweighted via Bayesian bootstrap rather than
resampled as in the classical residual bootstrap. While other bootstrap weights can be employed, the Bayesian bootstrap weights provide a natural Bayesian interpretation (see Section \ref{Sec:43}).
Our debiased Bayesian procedure provides simultaneous point estimation and uncertainty quantification. For $j=1,\ldots,p$, let $\beta_{0,j}$ and $\tilde{\beta}^{b}_j$ be the $j$th coordinate of $\beta_{0}$ and $\tilde{\beta}^{b}$ respectively. For $0<\alpha<1$, a $100(1-\alpha)\%$ credible set for the regression coefficient $\beta_{0,j}$ is defined as:
\begin{align}\label{Eqn:CS}
\mathcal{C}_{n,j}(\alpha,B)
= \big[\,\hat{c}_{n,j}(\alpha/2,B),~~\hat{c}_{n,j}(1-\alpha/2,B)\,\big],
\end{align}
where $\hat{c}_{n,j}(\alpha,B)$ denotes the $\alpha$th quantile of the posterior draws $\{\tilde{\beta}_j^b : b=1,\ldots,B\}$. The Bayesian point estimator (posterior mean) is obtained by averaging the simulation draws: $\bar{\tilde \beta}_j = \sum_{b=1}^B \tilde{\beta}_j^b/B$ for each $j=1\ldots, p$.
\section{Insights from Empirical Applications}\label{Sec: Applications}
Our paper aims not only to present a comprehensive theoretical development, but also to demonstrate the empirical benefits through concrete empirical applications. Before presenting the theoretical properties of our debiased Bayesian inference method, we first demonstrate its practical utility through two empirical exmaples adapted from \citet{GiannoneLenzaPrimiceri2021Sparsity}. In both cases, our goal is not to determine the overall sparsity of the regressors, but rather to evaluate the significance of a particular regressor—an objective well suited to our inference procedure.
For each application, we begin by presenting standard Bayesian inference results for the coefficient of interest using both a spike-and-slab prior and a horseshoe-type prior. The prior specifications and hyperparameter settings follow those used in our Monte Carlo simulations. For the spike-and-slab prior, we obtain the posterior distribution via the variational Bayesian approximation described in \citet{RaySzabo2022Variational}, drawing 8,000 samples from the resulting approximate posterior. The posterior under the horseshoe prior is generated using the MCMC algorithm of \citet{KimLeeGupta2020BayesianSC}, with a total of 16,000 draws and the first 8,000 discarded as burn-in. In both applications, the spike-and-slab posterior collapses to a point mass at zero. Under the horseshoe prior, the posterior exhibits greater dispersion, though in the second application the density remains tightly concentrated near zero. In both cases, the 95\% credible intervals include zero, implying that the effects are insignificant at the $5\%$ level.
To implement the debiasing adjustment, we estimate the precision matrix using the nodewise LASSO regression method of \citet{VandeGeer_OnAsymptotically_2014}, adopting the default tuning parameters from the R package \texttt{hdi} \citep{DBMM2015hdi}. For the Bayesian bootstrap, we generate 8,000 additional sets of bootstrap weights in parallel with the posterior draws of the regression coefficients. After applying the debiasing step, the posterior densities become approximately Gaussian, providing empirical support for our theoretical results. Overall, the debiased inference reveals a significant effect in the first application but no significant effect in the second.
\subsection{Determinants of Economic Growth}
The influential work of \citet{Barro1991Growth} initiated a long-standing debate over what drives long-term economic growth across countries. Researchers have identified numerous potential predictors, many of which are included in the original data set compiled by \citet{Barro1991Growth}. In line with \citet{BelloniChernozhukovHansen2013InferHigh}, we use this data set to analyze average GDP growth between 1960 and 1985 for 90 countries. The data set contains 60 candidate predictors, covering pre-1960 measures of a wide range of socio-economic, institutional, and geographical factors. Using standard Bayesian inference methods, we find that all of the posterior results are statistically insignificant. However, applying our debiased Bayesian approach uncovers a significant negative association between the initial GDP level and subsequent growth. This suggests the presence of catch-up effect, everything else equal, which is in line with the neoclassical economic growth theory.
\begin{figure}[htbp]
\centering
\begin{tikzpicture}
\begin{axis}[
width=0.95\textwidth,
height=7cm,
xlabel={Coefficient value},
ylabel={Posterior density},
xlabel style={font=\small},
ylabel style={font=\small},
ticklabel style={font=\footnotesize},
axis lines=left,
xmin=-1.8, xmax=0.6,
xtick={-1.5, -1, -0.5, 0, 0.5,1},
ymin=0, ymax=2.8,
ytick={0,1,2},
legend style={
at={(0.1,0.99)},
anchor=north west,
draw=none,
fill=none,
font=\small,
align=left
},
legend cell align=left,
tick style={black},
legend columns=1,
every axis plot/.append style={thick},
]
\pgfplotstableread[col sep=comma]{DH1.csv}\DHone
\pgfplotstableread[col sep=comma]{DV1.csv}\DVone
\pgfplotstableread[col sep=comma]{HS1.csv}\HSone
\addplot[
dashed,
thick,
black
] coordinates {(0,0) (0,2.7)};
\addlegendentry{Spike-and-Slab}
\addplot[
color=green,
thick
]
table[
x=x,
y=y,
] {\HSone};
\addlegendentry{Horseshoe}
\addplot[
color=blue,
thick
]
table[
x=x,
y=y,
] {\DVone};
\addlegendentry{Debiased-Spike-and-Slab}
\addplot[
color=red,
thick
]
table[
x=x,
y=y,
] {\DHone};
\addlegendentry{Debiased-Horseshoe}
\addplot [
name path=E,
draw=none,
forget plot
]
table[
x=x,
y=y,
] {\HSone};
\addplot [
name path=F,
domain=-3:3,
draw=none,
forget plot
] {0};
\addplot [
green,
opacity=0.3
] fill between [of=E and F, soft clip={domain=-0.846580927:0.126281795}];
\addplot [
name path=A,
draw=none,
forget plot
]
table[
x=x,
y=y,
] {\DVone};
\addplot [
name path=B,
domain=-3:3,
draw=none,
forget plot
] {0};
\addplot [
blue,
opacity=0.3
] fill between [of=A and B, soft clip={domain=-1.252297236:-0.133349493}];
\addplot [
name path=C,
draw=none,
forget plot
]
table[
x=x,
y=y,
] {\DHone};
\addplot [
name path=D,
domain=-3:3,
draw=none,
forget plot
] {0};
\addplot [
red,
opacity=0.3
] fill between [of=C and D, soft clip={domain=-1.201603532:-0.1728713}];
\end{axis}
\end{tikzpicture}
\caption{Posterior distribution with 95\% credible set}
\end{figure}
\subsection{Decline in Crime Rates}
Using U.S. state-level data, \citet{DonohueLevitt2001Abortion} identified a strong relationship between the legalization of abortion following the 1973 Roe v.\ Wade decision and the subsequent decline in crime rates. In their analysis, the dependent variable is the change in log per capita murder rates between 1986 and 1997 across states. This outcome is regressed on a measure of the effective abortion rate—which is always included as a predictor alongside twelve year dummy variables—and a collection of controls capturing alternative determinants of crime, such as the number of police officers and prisoners per 1,000 residents, among others. \citet{BelloniChernozhukovHansen2014HighATE} extended the set of controls by incorporating these variables in various forms, including levels, differences, squared differences, cross-products, initial conditions, and interactions with linear and quadratic time trends, resulting in a database of 284 variables and 576 observations. When analyzing this expanded model, the coefficient on the abortion rate is consistently found to be insignificant across all methods. Note that the horseshoe induced posterior places almost all of its posterior mass at zero, which is close to the degenerate (at zero) posterior of the spike-and-slab. Therefore, these two debiased posteriors—whether based on the spike-and-slab or horseshoe priors—are largely indistinguishable in this application. In this example, standard Bayesian approaches place (nearly) all posterior mass at zero, while our debiased Bayesian procedure yields more nuanced posterior distributions that may better capture estimation uncertainty.
\begin{figure}[htbp]
\centering
\begin{tikzpicture}
\begin{axis}[
width=0.95\textwidth,
height=7cm,
xlabel={Coefficient value},
ylabel={Posterior density},
xlabel style={font=\small},
ylabel style={font=\small},
ticklabel style={font=\footnotesize},
axis lines=left,
xmin=-0.7, xmax=0.7,
xtick={-0.5, 0, 0.5},
ymin=0, ymax=10.5,
ytick={0,5,10},
legend style={
at={(0.1,0.99)},
anchor=north west,
draw=none,
fill=none,
font=\small,
align=left
},
legend cell align=left,
tick style={black},
legend columns=1,
every axis plot/.append style={thick},
]
\pgfplotstableread[col sep=comma]{DH2.csv}\DHone
\pgfplotstableread[col sep=comma]{DV2.csv}\DVone
\pgfplotstableread[col sep=comma]{HS2.csv}\HSone
\addplot[
dashed,
thick,
black
] coordinates {(0,0) (0,10.5)};
\addlegendentry{Spike-and-Slab}
\addplot[
color=green,
thick
]
table[
x=x,
y=y,
] {\HSone};
\addlegendentry{Horseshoe}
\addplot[
color=blue,
thick
]
table[
x=x,
y=y,
] {\DVone};
\addlegendentry{Debiased-Spike-and-Slab}
\addplot[
color=red,
thick
]
table[
x=x,
y=y,
] {\DHone};
\addlegendentry{Debiased-Horseshoe}
\addplot [
name path=E,
draw=none,
forget plot
]
table[
x=x,
y=y,
] {\HSone};
\addplot [
name path=F,
domain=-3:3,
draw=none,
forget plot
] {0};
\addplot [
green,
opacity=0.3
] fill between [of=E and F, soft clip={domain=-0.018099333:0.00054724}];
\addplot [
name path=A,
draw=none,
forget plot
]
table[
x=x,
y=y,
] {\DVone};
\addplot [
name path=B,
domain=-3:3,
draw=none,
forget plot
] {0};
\addplot [
blue,
opacity=0.3
] fill between [of=A and B, soft clip={domain=-0.43715268:0.384011518}];
\addplot [
name path=C,
draw=none,
forget plot
]
table[
x=x,
y=y,
] {\DHone};
\addplot [
name path=D,
domain=-3:3,
draw=none,
forget plot
] {0};
\addplot [
red,
opacity=0.3
] fill between [of=C and D, soft clip={domain=-0.434090556:0.382386431}];
\end{axis}
\end{tikzpicture}
\caption{Posterior distribution with 95\% credible set (blue and red overlap)}
\end{figure}
\section{Main Theoretical Results}\label{Sec: 4}
Various features of posterior distributions are used by empirical researchers for frequentist type inference. In particular, regions with high posterior probabilities, known as credible sets, are often interpreted as confidence sets. The frequentist validity of the debiased posterior distributions serves as a theoretical justification for the proposed Bayesian procedure. In this section, we present the main theoretical results, beginning with a generic BvM theorem formulated under a set of high-level conditions. The result applies broadly and is not restricted to specific prior choices.
For the reader’s convenience, we introduce notation that will be used throughout the paper. Let $I_p$ denote the $p\times p$ identity matrix, $e_j$ its $j$th column, and $E_J$ the matrix formed by any $J$ columns of $I_p$. For $1 \leq q < \infty$, we write $\|v\|_{q}$ for the standard $\ell_q$-norm of a vector $v$, that is, $\|v\|_{q}: = (\sum_j |v_j|^{q})^{1/q}$. We use $\|v\|_{0}$ to denote the number of nonzero entries of $v$, and $\|v\|_{\infty}$ for its sup-norm. For a matrix $A$, $\|A\|_{q}$ denotes the induced (operator) $\ell_q$-norm for $1\leq q\leq \infty$, and $\|A\|_{\max}$ the elementwise sup-norm. For two positive sequences $a_n$ and $b_n$, we write $a_n\lesssim b_n$ if $a_n\leq C b_n$ for some constant $C$ and all $n$, and $a_n\asymp b_n$ if $a_n\lesssim b_n$ and $b_n \lesssim a_n $. Let $Z^{(n)} := \{Y_i,X_i\}_{i=1}^n$ denote the data and $W^{(n)} := \{W_{ni}\}_{i=1}^n$ the Bayesian bootstrap weights. As previously, $\mathbb{E}_0[\cdot]$ indicates that the expectation is evalauted with respect to the distribution of $Z^{(n)}$ under $\beta=\beta_0$, while $\Pi_W(\,\cdot\,|Z^{(n)})$ denotes the conditional distribution of $W^{(n)}$ given $Z^{(n)}$. Since $W^{(n)}$ is independent of $Z^{(n)}$, this conditional distribution does not depend on $Z^{(n)}$. For a positive sequence $r_n$, the asymptotic symbols $o_{P_0}(r_n)$ and $O_{P_0}(r_n)$ are defined with respect to the same underlying probability measure $P_0$. The sub-Gaussian norm of a random variable $Z$, denoted by $\|Z\|_{\psi_2}$, is defined as $\|Z\|_{\psi_2}:=\inf\{t>0: \mathbb E_{0}[\exp(Z^2/t^2)]\leq 2\}$. For a random vector $Z$, its sub-Gaussian norm is defined as $\|Z\|_{\psi_2}:=\sup _{\|x\|_2=1}\| Z^{\intercal}x\|_{\psi_2}$.
\subsection{Bernstein-von Mises Theorem}
We begin by establishing the result under high-level assumptions that specify the posterior contraction rate of $\Pi_{\beta}(\,\cdot\,| Z^{(n)})$, the convergence rate of the precision matrix estimator $\hat{\Theta}_n$, a frequentist point estimator that centers the debiased posterior distribution, and suitable moment and regularity conditions.
\begin{ass}[Contraction rate of $\Pi_{\beta}(\,\cdot\,| Z^{(n)})$]\label{Assump:Beta}
There exist some constant $C>0$ and a positive sequence $\epsilon_n\to 0$ such that
\[\mathbb{E}_0\big[\Pi_{\beta}\big(\|\beta-\beta_0\|_1\geq C\epsilon_n \mid Z^{(n)}\big)\big]\to 0.\]
\end{ass}
\begin{ass}[Convergence rate of $\hat{\Theta}_n$]\label{Assump:Precision}
Let $\hat{\Omega}_n:= \sum_{i=1}^{n}X_iX_i^{\intercal}/n$. There exist positive sequences $\gamma_n\to 0$ and $\delta_n\to 0$ such that
\[\|\hat{\Theta}_n\hat{\Omega}_n - I_{p}\|_{\max}=O_{P_0}(\gamma_n) \text{ and } \|\hat{\Theta}_n-\Theta_0\|_{\infty}=O_{P_0}(\delta_n).\]
\end{ass}
\begin{ass}[Centering point estimator]\label{Assump:Frequentist}
There exists a point estimator $\hat\beta_n$ (depending only on $Z^{(n)}$) such that the following expansion holds
\[\hat\beta_n = \beta_0 + \Theta_0 \frac{1}{n} \sum_{i=1}^{n}X_i\varepsilon_i + \Delta_n, ~\text{ where } \|\Delta_n\|_{\infty} = o_{P_0}(n^{-1/2}).\]
\end{ass}
\begin{ass}[Moment and regularity conditions]\label{Assump:Moments}
(i) $\{Y_{i},X_{i}\}_{i=1}^{n}$ are i.i.d. and satisfy the model \eqref{Eqn:Model};
(ii) $X_{i}$'s are sub-Gaussian random vectors with $\|X_{i}\|_{\psi_2}\in(0,\infty)$;
(iii) $\varepsilon_{i}$'s are sub-Gaussian random variables with $\|\varepsilon_{i}\|_{\psi_2}\in(0,\infty)$;
(iv) $\mathbb{E}_0[|e_j^{\intercal}\Theta_0 X_{i}\varepsilon_i|^3]<\infty$ for each $j=1,\ldots, p$.
\end{ass}
Assumption \ref{Assump:Beta} imposes a contraction rate on the initial posterior distribution of $\beta$.
We later provide primitive conditions under which this assumption holds for both spike-and-slab priors and horseshoe-type priors \citep{Castilloetal_BayesianSparse_2015, SongLiang_NearlyOptimal_2023}.
Importantly, we allow for either the exact posterior or an approximate posterior for $\beta$.
The exact posterior arises directly from the Bayes’ rule, given a prior and the likelihood function.
In low-dimensional settings, standard MCMC algorithms can deliver fairly accurate approximations to this exact posterior. In contrast, in high-dimensional problems, it may no longer be feasible. For example, the point-mass spike-and-slab prior is often considered theoretically ideal for sparse Bayesian problems. However, exploring the full posterior over the entire model space using point-mass spike-and-slab priors can be computationally prohibitive, because of the combinatorial complexity of updating the discrete indicators whether to include each variable or not. Recently, the variational Bayesian approximation has become increasingly popular where one relies on an approximate posterior that minimizes the Kullback–Leibler divergence to the true posterior within a restricted family, such as the mean-field variational family. We provide sufficient conditions to ensure that Assumption \ref{Assump:Beta} holds for mean-field variational approximations under spike-and-slab priors \citep{RaySzabo2022Variational}.
Assumption \ref{Assump:Precision} specifies the convergence rate of the precision matrix estimator.
When the dimensionality $p$ is fixed, this condition is trivially satisfied by the plug-in estimator $\hat{\Theta}_n = \hat{\Omega}_n^{-1}$ with $\gamma_n = 0$ and $\delta_n = n^{-1/2}$ under mild moment assumptions.
In high-dimensional settings ($p \to \infty$), we will present primitive conditions under which the assumption holds when $\hat{\Theta}_n$ is obtained from the nodewise LASSO regressions of \citet{VandeGeer_OnAsymptotically_2014} or the CLIME approach of \cite{Caietal_ConstrainedSparse_2011}. Assumption \ref{Assump:Frequentist} requires the existence of a frequentist point estimator that serves as the center of the debiased posterior. It is satisfied by the OLS estimator when $p$ is fixed and by the debiased LASSO estimator in high-dimensional models \citep{VandeGeer_OnAsymptotically_2014} under the primitive conditions provided in Section \ref{Sec:42}. Besides LASSO, other types of pilot estimators have also been studied in the literature; see Remark \ref{Rem:DebiasPoint}. Finally, Assumption~\ref{Assump:Moments} imposes sub-Gaussianity on both the regressors and the error terms. These regularity conditions are standard in the high-dimensional inference literature and ensure well-behaved concentration of sample quantities.
The debiased posterior distribution is jointly determined by the initial posterior distribution of $\beta$ and the distribution of the Bayesian bootstrap weights. Specifically, the debiased posterior distribution is determined by
\begin{align}
\Pi(\,\cdot\,|Z^{(n)}) := \Pi_{\beta}(\,\cdot\,|Z^{(n)}) \times \Pi_W(\,\cdot\,|Z^{(n)}),
\end{align}
where $\Pi_{\beta}(\,\cdot\,|Z^{(n)})$ and $\Pi_W(\,\cdot\,|Z^{(n)})$ are, respectively, the initial posterior distribution of $\beta$ and the distribution of the Bayesian bootstrap weights conditional on the data.
For each $j=1,\ldots,p$, let $\mathcal{L}_{\Pi}\big(e_{j}^{\intercal}\sqrt{n}(\tilde\beta-\hat{\beta}_{n})\mid Z^{(n)}\big)$ denote the posterior law of $e_{j}^{\intercal}\sqrt{n}(\tilde\beta-\hat{\beta}_{n})$ given $Z^{(n)}$. We study the weak convergence of these posterior laws, which we measure using the bounded Lipschitz distance $d_{BL}$.
For probability measures $P$ and $Q$ on $\mathbf{R}^k$, the bounded Lipschitz distance is defined as
\begin{align}\label{Eqn:dBL}
d_{BL}(P,Q)
:= \sup\left\{
\left| \int f\, d(P-Q) \right|
: \|f\|_{BL}\le 1
\right\},
\end{align}
where the bounded Lipschitz norm of a measurable function $f:\mathbf{R}^k\to\mathbf{R}$ is given by
\begin{align}
\|f\|_{BL}
:= \sup_{x\in\mathbf{R}^k}|f(x)|
+ \sup_{x\neq y}\frac{|f(x)-f(y)|}{\|x-y\|_\infty}.
\end{align}
\begin{thm}\label{Thm:BvM1}
Suppose Assumptions \ref{Assump:Beta}-\ref{Assump:Moments} hold. If $\sqrt{n}\epsilon_n\gamma_{n}\to0$ and $\sqrt{n}(\epsilon_n\|\Theta_0\|_\infty+ \delta_n)(\sqrt{\log p/n} + \log^2 n\log^{2} p/n)\to 0$, then for each $j=1,\ldots, p$,
\[d_{BL}\left(\mathcal{L}_{\Pi}\big(e_{j}^{\intercal}\sqrt{n}(\tilde\beta-\hat{\beta}_{n})\mid Z^{(n)}\big), N(0,\sigma^2_{0,j})\right)\overset{P_{0}}{\to} 0,\]
where $\sigma^2_{0,j}: = e_{j}^{\intercal}\Theta_0 \mathbb{E}_0[X_{i}X_{i}^{\intercal}\varepsilon_{i}^{2}]\Theta_0 e_{j}$. That is, the posterior law of each coordinate of $\sqrt{n}(\tilde\beta-\hat{\beta}_{n})$ converges weakly to a normal distribution in probability.
\end{thm}
Theorem \ref{Thm:BvM1} establishes a coordinatewise BvM result for the debiased posterior.
In particular, the asymptotic variance $\sigma^{2}_{0,j}$ for the $j$th coordinate of $\sqrt{n}(\tilde{\beta}-\hat{\beta}_{n})$ coincides with the asymptotic variance for the centering point estimator in Assumption \ref{Assump:Frequentist} under general, possibly heteroskedastic, errors.
The rate conditions $\sqrt{n}\epsilon_n\gamma_{n}\to0$ and $\sqrt{n}(\epsilon_n\|\Theta_0\|_\infty+ \delta_n)(\sqrt{\log p/n} + \log^2 n\log^{2} p/n)\to 0$ ensure two key properties:
(i) the bias stemming from the initial posterior and the error from estimating the precision matrix are asymptotically negligible; and
(ii) the residual posterior uncertainty, after the debiasing correction, is asymptotically Gaussian. Importantly, the theorem accommodates a high-dimensional regime in which the dimensionality $p$ may grow exponentially with the sample size $n$.
As will be shown later, under spike-and-slab or horseshoe priors for $\beta$ and nodewise LASSO or CLIME estimation for $\Theta_0$, the quantities $\epsilon_n$, $\gamma_n$, and $\delta_n$ can all be of order $\sqrt{\log p/n}$ if the sparsity levels of $\beta_0$ and $\Theta_0$ are bounded (implying $\|\Theta_0\|_\infty$ is bounded). Hence, the result remains valid as long as $\log^{5/2} p$ grows at most linearly with $n$, up to logarithmic factors.
The BvM theorem enables the construction of marginal credible intervals for any coordinate of the regression coefficient from the debiased posterior.
For $0<\alpha<1$, let $q_{n,j}(\alpha)$ denote the $\alpha$th quantile of the debiased posterior distribution for the $j$th coefficient---that is, the $\alpha$th conditional quantile of $\tilde{\beta}_{j}$ (the $j$th coordinate of $\tilde{\beta}$) given $Z^{(n)}$. We study the asymptotic frequentist coverage of the $100(1-\alpha)\%$ credible intervals, defined for $j=1,\ldots,p$ as
\begin{align}\label{Eqn:CSinf}
\mathcal{C}_{n,j}(\alpha)
= \big[\,q_{n,j}(\alpha/2),~~q_{n,j}(1-\alpha/2)\,\big].
\end{align}
In practice, $\mathcal{C}_{n,j}(\alpha)$ is infeasible but can be approximated by its simulated counterpart $\mathcal{C}_{n,j}(\alpha,B)$ defined in \eqref{Eqn:CS} by taking a large $B$.
\begin{cor}\label{Cor:BvM1}
Under the conditions of Theorem~\ref{Thm:BvM1}, for each $j=1,\ldots,p$ and $\alpha\in(0,1)$, if in addition $\sigma_{0,j}^{2}$ is bounded away from zero, then
\[
P_0\big(\beta_{0,j}\in \mathcal{C}_{n,j}(\alpha)\big) \to 1-\alpha.
\]
\end{cor}
This corollary provides an important frequentist implication of Theorem \ref{Thm:BvM1}:
the Bayesian credible intervals derived from the debiased posterior achieve asymptotically correct frequentist coverage for each regression coefficient.
In other words, for large samples, the posterior-based uncertainty quantification coincides with the frequentist notion of confidence intervals, thereby unifying Bayesian and frequentist inference in high-dimensional settings.
We next consider simultaneous inference for a growing number of coefficients. Without loss of generality, we focus on the first $J$ coordinates. For this purpose, we strengthen Assumption \ref{Assump:Moments}(iv) to the following condition.
\begin{ass}[Moment condition for simultaneous inference]\label{Assump:MomentsSimutaneous}
(i) The eigenvalues of $\Sigma^2_{0,J}: = E_{J}^{\intercal}\Theta_0 \mathbb{E}_0[X_{i}X_{i}^{\intercal}\varepsilon_{i}^{2}]\Theta_0 E_{J}$ are bounded;
(ii) there exists a sequence $\nu_J>0$ such that $\{\mathbb E_0[\|E_{J}^{\intercal}\Theta_0X_{i}\varepsilon_i\|^{3}_{\infty}]\}^{1/3}=O(\nu_{J})$.
\end{ass}
Assumption \ref{Assump:MomentsSimutaneous}(i) is satisfied if, for example, $\mathbb{E}_0[\varepsilon_{i}^{2}|X_i]$ is constant and the eigenvalues of $\Omega_0$ are bounded away from zero.
Assumption \ref{Assump:MomentsSimutaneous} (ii) holds with $\nu_J = O(\log J)$ when the random vectors $\Theta_0X_{i}\varepsilon_i$ are sub-exponential.
\begin{thm}\label{Thm:BvM2}
Suppose Assumptions \ref{Assump:Beta}-\ref{Assump:Frequentist},
\ref{Assump:Moments}(i)-(iii), and \ref{Assump:MomentsSimutaneous} hold.
If
$\sqrt{n}\epsilon_n\gamma_{n}\to0$, $\sqrt{n}(\epsilon_n\|\Theta_0\|_\infty+ \delta_n)(\sqrt{\log p/n} + \log^2 n\log^{2} p/n)\to 0$,
and $\nu_{J}^{6}J^{3}\log^{9}(1+J)/n\to 0$, then
\[
d_{BL}\!\left(
\mathcal{L}_{\Pi}\big(E_{J}^{\intercal}\sqrt{n}(\tilde{\beta}-\hat{\beta}_{n})\mid Z^{(n)}\big),\,
N(0,\Sigma^2_{0,J})
\right)
\overset{P_{0}}{\longrightarrow} 0.
\]
\end{thm}
Theorem \ref{Thm:BvM2} extends the marginal BvM result in Theorem \ref{Thm:BvM1} to simultaneous inference on a growing subset of regression coefficients. Specifically, the result establishes joint asymptotic normality of the debiased posterior for any subvector of $\tilde{\beta}$ of size $J$, provided that $J$ increases with the sample size $n$ at a suitable rate. This rate is derived from the strong approximation result for the exchangeable bootstrap in \cite{FangSantosShaikhTorgovitsky2023LP}. The condition $\nu_J^{6}J^{3}\log^{9}(1+J)/n \to 0$ captures the complexity of the simultaneous inference problem, balancing the growth of the subset dimension against sample size and moment bounds. If we assume enough moment restrictions on $\|E_{J}^{\intercal}\Theta_0X_{i}\varepsilon_i\|_{\infty}$ (e.g., $\mathbb E_0[\|E_{J}^{\intercal}\Theta_0X_{i}\varepsilon_i\|^{4}_{\infty}]$ being bounded uniformly in $n$ and $p$), we can obtain Theorem \ref{Thm:BvM2} under $J/n^2\to\infty$ (up to logarithmic factors) by Theorem \ref{Thm:strong approx}.
\subsection{Primitive Conditions}\label{Sec:42}
We now verify the high-level conditions in Assumptions \ref{Assump:Beta}-\ref{Assump:Frequentist} by introducing a set of primitive conditions. For concreteness, we focus on a benchmark specification in which the regression coefficients are assigned a spike-and-slab prior, the precision matrix is estimated via the nodewise LASSO procedure, and the centering point estimator is constructed by the deibased LASSO. The general theoretical framework, however, extends to global–local shrinkage priors, including the horseshoe family, and other estimation approaches for the precision matrix, including the CLIME approach. Detailed results for these additional cases are provided in the Appendix.
We follow the construction of \cite{Castilloetal_BayesianSparse_2015}, who specify the spike-and-slab prior through the following hierarchical mixture for $\beta = (\beta_1,\ldots,\beta_p)^{\intercal}$,
\begin{align}\label{SSprior}
\pi(\beta|r,\lambda) = \prod_{j=1}^{p}\big[(1-r)\,\delta_0(\beta_j) + r\,\psi(\beta_j|\lambda)\big],
\end{align}
where $\delta_0$ denotes a point mass at zero, and $\psi(\beta_j|\lambda) = \frac{\lambda}{2}\exp(-\lambda|\beta_j|)$ is the Laplace density with hyperparameter $\lambda>0$. The mixing weight $r$ controls the proportion of nonzero coefficients\footnote{If $r$ were taken to be deterministic, it would represent the expected number of nonzero coefficients \textit{a priori}.} and thus governs the degree of model sparsity. Following \citet{Castilloetal_BayesianSparse_2015}, we assign a Beta hyper-prior
\begin{align}\label{SSpriorBeta}
r\sim \mathrm{Beta}(1,p^u),
\end{align}
with hyperparameter $u>1$, large values of which favor sparse models by placing most of the prior mass on small values of $r$.
The point mass in the spike-and-slab prior is designed for explicit variable selection. However, the posterior under this prior places mass over all $2^p$ candidate models, which is computationally prohibitive in high dimensions. We therefore adopt a variational Bayesian approximation scheme. The variational posterior is defined as the best approximation of the posterior within a given class. The crux is to choose this class both sufficiently rich to approximate the exact posterior well, while at the same time also simple enough so that the minimization problem can be efficiently solved. For this purpose, we consider the following mean-field family from \cite{RaySzabo2022Variational} for $\mu = (\mu_1,\ldots,\mu_p)^{\intercal}$, $\sigma^2 = (\sigma^2_1,\ldots,\sigma^2_p)^{\intercal}$, and $\gamma = (\gamma_1,\ldots,\gamma_p)^{\intercal}$,
\begin{align}\label{VBFamily}
\hspace{-0.2cm}\mathcal{P}_{MF}:=\left\{P_{\mu,\sigma^2,\gamma}=\bigotimes_{j=1}^p\big[\gamma_j\,N(\mu_j,\sigma_j^2)+(1-\gamma_j)\,\delta_0\big]:\mu_j\in\mathbf{R},\sigma^2_j>0,\gamma_j\in[0,1] \right\},
\end{align}
where $\gamma_j$ plays the role of a variational inclusion probability.
The resulting variational Bayesian approximate posterior is defined as the minimizer of the Kullback-Leibler divergence with respect to the exact posterior:
\begin{align}\label{VBPost}
\tilde{\Pi}_{\beta}(\,\cdot\,|Z^{(n)})= \underset{P_{\mu,\sigma^2,\gamma}\in \mathcal{P}_{MF}}{\arg \min}\mathrm{KL}\left(P_{\mu,\sigma^2,\gamma}~\|~\Pi_{\beta}(\,\cdot\,|Z^{(n)})\right),
\end{align}
where $\mathrm{KL}(\,P\,\|\,Q\,):=\int \log (dP/dQ) dP$ for two probability distributions $P$ and $Q$, and $\Pi_{\beta}(\,\cdot\,|Z^{(n)})$ is the exact posterior. Note that the above mean-field class enforces substantial independence in the approximated posterior. By doing so, it significantly reduces the model complexity, because there are only $p$ inclusion variables $\gamma_j$ to consider, rather than a total of $2^p$ models which the exact posterior puts mass on.
The mean-field class approximates the posterior distribution by a product-form distribution, thereby ignoring dependence among different coordinates. Nevertheless, it has been shown to achieve the desired rate of contraction to the true parameter \citep{RaySzabo2022Variational}, which is essential for our high-level conditions to ensure the asymptotic normality of the debiased posterior. Without the debiasing step, however, the variational posterior generally fails to deliver asymptotically correct coverage, even in low-dimensional parametric settings \citep{WangBlei2019VB}. This highlights the versatility of our debiasing proposal from a complementary perspective. The trade-off for relaxing the exactness of the posterior distribution lies in the substantial computational gains. Component-wise coordinate-ascent variational inference algorithms for $\tilde{\Pi}_{\beta}(\,\cdot\,|Z^{(n)})$ are given in \cite{RaySzabo2022Variational}. Beyond the mean-field family, our arguments also extend to other closely related variational classes that allow dependence among the nonzero coordinates (see Equation (9) in \citealp{RaySzabo2022Variational}), yielding analogous theoretical guarantee.
Turning to Assumption \ref{Assump:Precision}, we consider the nodewise LASSO regression following \citet{VandeGeer_OnAsymptotically_2014}. Let $X_{j,i}$ denote the $j$th component of $X_{i}$ and let ${X}_{-j,i}\in\mathbf{R}^{p-1}$ be the subvector of $X_{i}$ excluding its $j$th entry. For each $j=1,\ldots, p$, we define the following LASSO regression:
\begin{align}\label{Eqn: nodewise1}
\hat{\theta}_j:=\underset{\theta \in \mathbf{R}^{p-1}}{\arg \min} \frac{1}{n}\sum_{i=1}^{n} (X_{j,i}-{X}_{-j,i}^{\intercal} \theta)^2 +2 \lambda_j\|\theta\|_1,
\end{align}
with the corresponding residual variance given by
\begin{align}\label{Eqn: nodewise2}
\hat{\tau}_j^2:= \frac{1}{n}\sum_{i=1}^{n} (X_{j,i}-{X}_{-j,i} ^{\intercal}\hat\theta_{j})^2+\lambda_j\|\hat{\theta}_j\|_1,
\end{align}
where $\lambda_{j}>0$ is a tuning parameter. Writing $\hat{\theta}_j=(\hat{\theta}_{j, 1},\ldots, \hat{\theta}_{j, j-1}, \hat{\theta}_{j, j+1}, \ldots, \hat{\theta}_{j, p})$, we construct the estimated precision matrix $\hat{\Theta}_n$ as
\begin{align}\label{Eqn: nodewise3}
\hat{\Theta}_n = \left(\begin{array}{cccc}
1/\hat{\tau}_1^2 & -\hat{\theta}_{1,2}/\hat{\tau}_1^2 & \cdots & -\hat{\theta}_{1, p}/\hat{\tau}_1^2 \\
-\hat{\theta}_{2,1}/\hat{\tau}_2^2 & 1/\hat{\tau}_2^2 & \cdots & -\hat{\theta}_{2, p}/\hat{\tau}_2^2 \\
\vdots & \vdots & \ddots & \vdots \\
-\hat{\theta}_{p, 1}/\hat{\tau}_p^2 & -\hat{\theta}_{p, 2}/\hat{\tau}_p^2 & \cdots & 1/\hat{\tau}_p^2
\end{array}\right).
\end{align}
The estimator $\hat\Theta_n$ is motivated by noting that the restriction $\Omega_0^{-1}\Omega_0=I_p$ amounts to first order conditions of $p$ linear projections after parametrizing $\Omega_0^{-1}$ in the same structure as $\hat\Theta_n$.
We choose the centering point estimator as the debiased LASSO estimator. For $\hat{\Theta}_n$ given in \eqref{Eqn: nodewise1}-\eqref{Eqn: nodewise3}, the debiased LASSO estimator is given by
\begin{align}\label{DebiasedLASSO}
\hat{\beta}_n=\hat{\beta}_n^{\text{LASSO}}+ \hat{\Theta}_n\left[\frac{1}{n}\sum_{i=1}^{n}X_i(Y_i-X_i^\intercal\hat{\beta}_n^{\text{LASSO}})\right],
\end{align}
where $\hat{\beta}_n^{\text{LASSO}}$ is a LASSO estimator given by
\begin{align}\label{LASSO}
\hat{\beta}_n^{\text{LASSO}}:=\underset{\beta \in \mathbf{R}^{p}}{\arg \min} \frac{1}{n}\sum_{i=1}^{n} (Y_i-X_{i}^{\intercal}\beta)^2 +\rho\|\beta\|_1,
\end{align}
where $\rho>0$ is a tuning parameter.
We adopt the following standard notation in analyzing the high-dimensional regression. For a vector $\beta=(\beta_1,\cdots,\beta_p)^{\intercal}\in\mathbf{R}^p$ and a subset $S\subset\{1,\cdots,p \}$ of indices, let $\beta_S$ be the vector $(\beta_j)_{j\in S}\in \mathbf{R}^{|S|}$, where $|S|$ is the cardinality of $S$. Let $S_{\beta}:=\{ j\in \{1,\cdots,p \} :\beta_{j}\neq 0\}$ denote the index set of nonzero components of $\beta$, and $s_{\beta}:= |S_{\beta}|$ denote the sparsity level of $\beta$.
Recall that $\hat{\Omega}_n= \sum_{i=1}^{n}X_iX_i^\intercal/n$. Define $\nu(\hat{\Omega}_n):=\max_{1\leq j\leq p}\sqrt{e_{j}^{\intercal}\hat{\Omega}_n e_{j}}$ as the square root of the largest diagonal element of $\hat{\Omega}_n$.
\begin{ass}[Model Specification]\label{Assump:Error}
(i) $X_i$ and $\varepsilon_i$ are independent with $\mathbb E_{0}[X_i] =0$ and $\varepsilon_{i}\sim N(0,1)$; (ii) the eigenvalues of $\Omega_0$ are bounded away from zero and the diagonal entries of $\Omega_0$ are bounded, that is, $\max_{1\leq j \leq p}e_{j}^\intercal \Omega_0 e_j<\infty$.
\end{ass}
\begin{ass}[Hyper/Tuning Parameters]\label{Assump:Tuning}
(i) $\nu(\hat{\Omega}_n)\sqrt{n} /p\leq\lambda\leq 4\nu(\hat{\Omega}_n)\sqrt{n\log p}$ holds with probability one; (ii) $\max_{1\leq j\leq p} \lambda_j \asymp \sqrt{\log p / n})$; (iii) $\rho \asymp \sqrt{\log p / n}$.
\end{ass}
\begin{ass}[Sparsity and Dimensionality]\label{Assump:Sparsity}
The following conditions hold:
\[s_{\beta_0}\|\Theta_0\|_{\infty}\frac{\log p}{\sqrt{n}}\left(1 + \frac{\log^{3/2}p\log^2 n}{\sqrt{n}}\right)\to 0\]
and
\[\max_{1\leq j\leq p}s_{\Theta_0,j}\frac{\log p}{\sqrt{n}}\left(1 + \frac{\log^{3/2}p\log^2 n}{\sqrt{n}}\right)\to 0,\]
where $s_{\Theta_0,j}:=s_{\Theta_0e_j}$ denotes the sparsity level of the $j$th row/column of $\Theta_0$.
\end{ass}
All the assumptions are standard in the literature. Following \cite{Castilloetal_BayesianSparse_2015},
we assume $\varepsilon_i \sim N(0,1)$ for simplicity. In the case of an unknown variance, $\varepsilon_i \sim N(0,\sigma_0^2)$, one may rescale the data using an estimate of $\sigma_{0}^2$, thereby adopting an empirical Bayes approach. Alternatively, a fully Bayesian treatment can be employed by assigning a prior to $\sigma_{0}^2$, for instance, an inverse-Gamma prior. In contrast to \cite{Castilloetal_BayesianSparse_2015} and \cite{RaySzabo2022Variational}, we provide primitive conditions for their high-level conditions on compatibility numbers. Since $\|\Theta_0\|_{\infty
}=O(\max_{1\leq j\leq p}\sqrt{s_{\Theta_0,j}})$, Assumption \ref{Assump:Sparsity} involves the sparsity levels of $\beta_0$ and $\Theta_0$ and the dimension $p$. However, the requirements are slightly different from those in the frequentist inference literature. For example, \citet{VandeGeer_OnAsymptotically_2014} require $s_{\beta_0}\log p/\sqrt{n}\to 0$ and $\max_{1\leq j\leq p}s_{\Theta_0,j}\log p/\sqrt{n}\to 0$. When the sparsity levels of $\beta_0$ and $\Theta_0$ are bounded, Assumption \ref{Assump:Sparsity} requires $\log^{5/2}p = o(n/\log^2 n)$, which is slightly more retrictive than $\log p = o(\sqrt{n})$ required in \citet{VandeGeer_OnAsymptotically_2014}. As in the literature on frequentist inference, $p$ is allowed to grow exponentially with $n$.
\begin{pro}\label{Pro:Primitive}
Consider the prior in \eqref{SSprior}-\eqref{SSpriorBeta}, the precision matrix estimator in \eqref{Eqn: nodewise1}-\eqref{Eqn: nodewise3}, and the centering point estimator defined by \eqref{DebiasedLASSO}-\eqref{LASSO}. Suppose Assumptions \ref{Assump:Moments}-\ref{Assump:Sparsity} hold and $\nu_{J}^{6}J^{3}\log^{9}(1+J)/n\to 0$. The conclusions of Theorems \ref{Thm:BvM1} and \ref{Thm:BvM2} hold for the exact posterior. If $\lambda =O_{P_0}(\sqrt{n\log p}/s_{\beta_0})$, the conclusions also hold for the variational Bayesian approximate posterior in \eqref{VBFamily}-\eqref{VBPost}.
\end{pro}
\begin{rem}
\cite{Castilloetal_BayesianSparse_2015} provide thorough theoretical analysis of the exact posterior in high-dimensional settings. They establish a distributional approximation result, showing that the posterior can be approximated by a mixture of normal distributions. Under suitable conditions on the hyperparameter $\lambda$ and a standard \emph{beta-min} condition, the posterior credible sets achieve correct asymptotic coverage for nonzero coefficients, while assigning asymptotic point mass at zero for truly zero coefficients. These assumptions are substantially stronger. More importantly, the computational burden of evaluating the exact posterior makes it infeasible in real high-dimensional problems. As our Monte Carlo simulations demonstrate, when using the variational Bayes approximation, the method performs poorly in finite samples related to nonzero coefficients, highlighting the practical limitations of the existing spike-and-slab framework in large-scale problems. \qed
\end{rem}
\begin{rem}\label{Rem:DebiasPoint}
It is well known that the LASSO estimator corresponds to the posterior mode under the prior: for $\beta = (\beta_1,\ldots,\beta_p)^{\intercal}$,
\begin{align}
\pi(\beta|\lambda) = \prod_{j=1}^{p}\frac{\lambda}{2}\exp(-\lambda|\beta_j|),
\end{align}
where $\lambda>0$ corresponds to the tuning parameter in LASSO controlling the amount of shrinkage and sparsity. Consequently, the frequentist debiasing approach can be viewed as applying a correction to a Bayesian point estimator. Beyond using the LASSO as the pilot estimator, the literature has also explored other Bayesian point estimators. For instance, \citet{ZhangPolitis2022Ridge} focus on ridge regression, which corresponds to the posterior mode under Gaussian priors on the regression coefficients, while \citet{Baietal2022Group} study the posterior mode in the additive regression, arising from a continuous spike-and-slab prior \citep{Rockova_ContinuousSpikeSlab_2018}. Our general debiased Bayesian inference encompasses these and other classes of priors. \qed
\end{rem}
\subsection{Additional Discussions}\label{Sec:43}
Our proposal can be seen as performing one-step update based on the identity $\beta_0 = \chi(\beta_0,\Theta_0,F_0)$, where $F_0$ denotes the distribution of $(Y_i,X_i)$, and $\chi$ is a mapping from the space of $(\beta_0,\Theta_0,F_0)$ to that of $\beta_0$ defined by
\begin{align}\label{Eqn: one-step}
\chi(\beta,\Theta,F)=\beta+\Theta\int x\left [y-x^\intercal \beta\right]dF(y,x).
\end{align}
We consider the posterior law of $F$, denoted by $\Pi_{F}(\,\cdot\,|Z^{(n)})$, to be the Bayesian bootstrap law, which can be viewed as the posterior of $F$ in a nonparametric Bayesian model with a Dirichlet process prior having a degenerate (zero-mass) base measure \citep{Rubin_BayesianBootstrap_1981}. Then, \eqref{Eqn:DebiasBayes} can be expressed as
\begin{align}
\tilde{\beta} = \chi(\beta,\hat{\Theta}_n,F),
\end{align}
where $(\beta,F)|Z^{(n)} \sim \Pi_{\beta}(\,\cdot\,|Z^{(n)}) \times \Pi_{F}(\,\cdot\,|Z^{(n)})$.
Informallly speaking, the debiased posterior is obtained by plugging in the frequentist estimator $\hat{\Theta}_n$ for $\Theta_0$, the initial posterior for $\beta$, and the Bayesian bootstrap for $F$. Note that the evaluation of $\Pi_{F}(\,\cdot\,|Z^{(n)})$ boils down to $\Pi_{W}(\,\cdot\,|Z^{(n)})$, whose randomness involves the weighted exponential random variables. Since the two posteriors are conditionally independent, the debiased posterior distribution arises from the product of the initial posterior distribution of $\beta$ and the posterior distribution of $F$. Next, we discuss the comparison of our approach with a number of recent proposals in the literature.
Our debiased posterior aligns with the one-step posterior framework of \citet{YiuFongHolmesRousseau2023Semiparametric}, which studies scalar parameters of interest in general semiparametric models. Nonetheless, our proposal is distinct in several important ways.
First, we consider a high-dimensional regime where the number of covariates may exceed the sample size, and we allow for simultaneous inference for a growing number of coefficients. The one-step correction term in \citet{YiuFongHolmesRousseau2023Semiparametric} involves the influence function perturbed by the Bayesian bootstrap weights. In a linear regression model, the $j$th regression coefficient $\beta_{0,j}$ has the influence function $(y,x^{\intercal})^{\intercal}\mapsto e_{j}^{\intercal}\Theta_{0} x(y-x^\intercal \beta_0)$ at the truth.
This representation permits the application of their one-step posterior in fixed-dimensional settings. In contrast, our framework explicitly employs sparsity-inducing priors in high-dimensional regimes. Moreover, the debiasing construction applies to the entire vector of regression coefficients, allowing the number of covariates to exceed the sample size while both $n$ and $p$ diverge. The resulting general BvM theorems enable inference on subvectors of regression coefficients whose dimension can grow with $(n,p)$.
Second, we estimate $\Theta_0$ via a frequentist method rather than a Bayesian one, which substantially improves computational efficiency. This necessitates addressing non-trivial technical challenges and developing new theoretical results. Suppose the covariates $\{X_i\}_{i=1}^n$ are \textit{i.i.d.} draws from ${N}(0,\Theta_0^{-1})$. The likelihood of $\{X_i\}_{i=1}^n$ as a function of the precision matrix $\Theta_0$ is
\begin{align}
L(\{X_i\}_{i=1}^n|\Theta)
= \frac{\det(\Theta)^{n/2}}{(2\pi)^{np/2}}
\exp\left(-\frac{1}{2}\sum_{i=1}^{n} X_{i}^{\intercal}\Theta X_{i}\right).
\end{align}
A natural Bayesian approach would place a sparsity-inducing prior, such as a spike-and-slab prior or a horseshoe-type prior, on the entries of $\Theta_0$ \citep{Atchade2019Graph}. However, computing the posterior distribution of a high-dimensional precision matrix is prohibitive in practice. Our proposal instead adopts a more pragmatic strategy: we simply plug in a frequentist estimator. In particular, one may compute $\hat{\Theta}_n$ either by running $p$ nodewise LASSO regressions or by solving $p$ linear programming problems in CLIME, both of which are computationally efficient. This plug-in approach is feasible because our procedure only requires $\hat{\Theta}_n$ to satisfy certain convergence rate conditions. Since our primary inferential target is the regression coefficients $\beta_0$, and $\Theta_0$ only plays a supporting role, this frequentist plug-in approach is both theoretically valid and practically effective.
Third, we accommodate a broad class of bootstrap weights beyond the Bayesian bootstrap weights; while the Bayesian bootstrap weights provide a natural Bayesian interpretation, alternative weights can be incorporated. A closer look at proofs indicates that other common exchangeable weights may also be employed. The Bayesian bootstrap law corresponds to the posterior law when the empirical distribution of $\{Y_i,X_i\}_{i=1}^n$ is equipped with a Dirichlet process prior with degenerate (zero-mass) base measure. For more general base measures, however, the resulting posterior law becomes substantially more cumbersome both theoretically and computationally; see Equation (2.1) in \cite{RayVaart2021BvM} for a representation of the posterior in this case.
\section{Monte Carlo Simulations}\label{Sec: Monte Carlo}
In this section, we evaluate the finite-sample performance of the proposed debiased Bayesian inference method through a series of Monte Carlo experiments. Specifically, we examine the frequentist coverage of the Bayesian credible sets and assess estimation accuracy in terms of bias and root mean squared error (RMSE). Beyond the stylized linear model with standard Gaussian errors presented in Section \ref{Sec: 4}, we further demonstrate the robustness and effectiveness of our approach across a range of more complex data-generating processes. We highlight that current results are all based on the default choices of hyper-parameters and tuning parameters in the existing R packages, which do not require additional interventions from users.
\pgfplotstableread{
beta alpha Bayes DebiasBayes DebiasLasso
0 0.95 0.9956444 0.9379778 0.9522667
1 0.95 0.1120000 0.9060000 0.8870000
2 0.95 0.0750000 0.9050000 0.9100000
3 0.95 0.0500000 0.8940000 0.8820000
4 0.95 0.0670000 0.9070000 0.6290000
5 0.95 0.0820000 0.9080000 0.8400000
}\pfnoheter
\pgfplotstableread{
beta alpha Bayes DebiasBayes DebiasLasso
0 0.95 0.9947778 0.9392 0.9522222
1 0.95 0.0620000 0.9210 0.9480000
2 0.95 0.0530000 0.9090 0.9270000
3 0.95 0.0570000 0.8880 0.9110000
4 0.95 0.0600000 0.9020 0.6450000
5 0.95 0.0730000 0.9250 0.8710000
}\pfnohomochi
\pgfplotstableread{
beta alpha Bayes DebiasBayes DebiasLasso
0 0.95 0.9913333 0.9454444 0.9519111
1 0.95 0.3070000 0.8910000 0.9250000
2 0.95 0.4590000 0.8810000 0.8470000
3 0.95 0.4670000 0.8850000 0.7230000
4 0.95 0.5210000 0.8990000 0.6220000
5 0.95 0.5810000 0.9010000 0.7560000
}\pfnohomon
\pgfplotstableread{
beta alpha Bayes DebiasBayes DebiasLasso
0 0.95 0.9998632 0.9392316 0.9528421
1 0.95 0.0530000 0.9210000 0.9080000
2 0.95 0.0560000 0.9310000 0.9440000
3 0.95 0.0780000 0.8900000 0.8970000
4 0.95 0.0820000 0.9060000 0.6760000
5 0.95 0.0940000 0.9290000 0.8530000
}\ponoheter
\pgfplotstableread{
beta alpha Bayes DebiasBayes DebiasLasso
0 0.95 0.9999158 0.9371368 0.9528105
1 0.95 0.0390000 0.9400000 0.9590000
2 0.95 0.0460000 0.9130000 0.9290000
3 0.95 0.0470000 0.9020000 0.9310000
4 0.95 0.0400000 0.9020000 0.7100000
5 0.95 0.0530000 0.9250000 0.8700000
}\ponohomochi
\pgfplotstableread{
beta alpha Bayes DebiasBayes DebiasLasso
0 0.95 0.9995158 0.9480526 0.9601368
1 0.95 0.1520000 0.9210000 0.9330000
2 0.95 0.3130000 0.9000000 0.8700000
3 0.95 0.3330000 0.8840000 0.7330000
4 0.95 0.4080000 0.8970000 0.5900000
5 0.95 0.5460000 0.9140000 0.7280000
}\ponohomon
\pgfplotstableread{
beta alpha Bayes DebiasBayes DebiasLasso
0 0.95 0.9999333 0.9378872 0.953159
1 0.95 0.0390000 0.9270000 0.921000
2 0.95 0.0430000 0.9220000 0.944000
3 0.95 0.0570000 0.8740000 0.918000
4 0.95 0.0500000 0.9140000 0.728000
5 0.95 0.0530000 0.9290000 0.885000
}\ptonoheter
\pgfplotstableread{
beta alpha Bayes DebiasBayes DebiasLasso
0 0.95 0.999959 0.9368205 0.9530615
1 0.95 0.018000 0.9210000 0.9540000
2 0.95 0.029000 0.9190000 0.9410000
3 0.95 0.046000 0.8840000 0.9340000
4 0.95 0.045000 0.9170000 0.7540000
5 0.95 0.042000 0.9280000 0.8780000
}\ptonohomochi
\pgfplotstableread{
beta alpha Bayes DebiasBayes DebiasLasso
0 0.95 0.9997385 0.9457128 0.9647282
1 0.95 0.1190000 0.9230000 0.9440000
2 0.95 0.2240000 0.8780000 0.8640000
3 0.95 0.2200000 0.8500000 0.7250000
4 0.95 0.3010000 0.8790000 0.5640000
5 0.95 0.4730000 0.9240000 0.7400000
}\ptonohomon
\pgfplotstableread{
beta alpha Bayes DebiasBayes DebiasLasso
0 0.95 0.9998444 0.9439556 0.9522
1 0.95 0.1900000 0.9210000 0.8590
2 0.95 0.2890000 0.9180000 0.9380
3 0.95 0.3330000 0.9330000 0.9500
4 0.95 0.3780000 0.9250000 0.9410
5 0.95 0.5790000 0.9450000 0.9680
}\CovpfnoheterSigmaone
\pgfplotstableread{
beta alpha Bayes DebiasBayes DebiasLasso
0 0.95 0.9998889 0.9355111 0.9498444
1 0.95 0.0800000 0.9090000 0.9370000
2 0.95 0.1220000 0.9240000 0.9390000
3 0.95 0.1300000 0.9270000 0.9550000
4 0.95 0.1540000 0.9370000 0.9550000
5 0.95 0.2510000 0.9450000 0.9560000
}\CovpfnohomochiSigmaone
\pgfplotstableread{
beta alpha Bayes DebiasBayes DebiasLasso
0 0.95 0.9991778 0.9604222 0.9544222
1 0.95 0.6210000 0.9200000 0.9470000
2 0.95 0.8080000 0.9120000 0.9240000
3 0.95 0.8710000 0.9230000 0.9320000
4 0.95 0.9050000 0.9210000 0.9360000
5 0.95 0.9280000 0.9450000 0.9510000
}\CovpfnohomonSigmaone
\pgfplotstableread{
beta alpha Bayes DebiasBayes DebiasLasso
0 0.95 0.9998211 0.9430632 0.9519053
1 0.95 0.1470000 0.9070000 0.8520000
2 0.95 0.1970000 0.9230000 0.9500000
3 0.95 0.2510000 0.9220000 0.9500000
4 0.95 0.3170000 0.9270000 0.9470000
5 0.95 0.5290000 0.9490000 0.9590000
}\CovponoheterSigmaone
\pgfplotstableread{
beta alpha Bayes DebiasBayes DebiasLasso
0 0.95 0.9999579 0.9342526 0.9499684
1 0.95 0.0460000 0.9280000 0.9580000
2 0.95 0.0700000 0.9310000 0.9560000
3 0.95 0.0790000 0.9140000 0.9370000
4 0.95 0.1110000 0.9340000 0.9530000
5 0.95 0.2010000 0.9550000 0.9610000
}\CovponohomochiSigmaone
\pgfplotstableread{
beta alpha Bayes DebiasBayes DebiasLasso
0 0.95 0.9991053 0.9630316 0.9560211
1 0.95 0.4760000 0.8800000 0.9440000
2 0.95 0.6820000 0.9040000 0.9310000
3 0.95 0.8080000 0.9240000 0.9390000
4 0.95 0.8660000 0.9300000 0.9380000
5 0.95 0.9250000 0.9460000 0.9480000
}\CovponohomonSigmaone
\pgfplotstableread{
beta alpha Bayes DebiasBayes DebiasLasso
0 0.95 0.9994667 0.9422359 0.9518769
1 0.95 0.1050000 0.8900000 0.8770000
2 0.95 0.1460000 0.9030000 0.9400000
3 0.95 0.1870000 0.9250000 0.9410000
4 0.95 0.2320000 0.9280000 0.9520000
5 0.95 0.4480000 0.9430000 0.9550000
}\CovptnoheterSigmaone
\pgfplotstableread{
beta alpha Bayes DebiasBayes DebiasLasso
0 0.95 0.9998256 0.9352974 0.9508923
1 0.95 0.0280000 0.9030000 0.9450000
2 0.95 0.0430000 0.9160000 0.9470000
3 0.95 0.0590000 0.9200000 0.9530000
4 0.95 0.0680000 0.9360000 0.9580000
5 0.95 0.1570000 0.9540000 0.9650000
}\CovptnohomochiSigmaone
\pgfplotstableread{
beta alpha Bayes DebiasBayes DebiasLasso
0 0.95 0.9987333 0.9634718 0.9564872
1 0.95 0.3640000 0.8100000 0.9550000
2 0.95 0.6020000 0.8850000 0.9380000
3 0.95 0.7430000 0.9130000 0.9350000
4 0.95 0.8540000 0.9300000 0.9420000
5 0.95 0.9240000 0.9390000 0.9520000
}\CovptnohomonSigmaone
\pgfplotstableread{
beta Bayes DebiasBayes DebiasLasso
0 0.001405958 0.0004518953 0.001585543
1 0.188266654 0.0198290763 0.006154626
2 0.379296830 0.0079773669 0.006147921
3 0.551779900 0.0223169900 0.005791037
4 0.687083060 0.0496947506 0.029428829
5 0.798797334 0.0089568229 0.008365539
}\BiaspfnoheterSigmaone
\pgfplotstableread{
beta Bayes DebiasBayes DebiasLasso
0 0.0001724629 0.021572592 0.022899255
1 0.2298244455 0.021302251 0.005174676
2 0.4558012007 0.018473234 0.002541299
3 0.6727092236 0.017276890 0.005215583
4 0.8796584906 0.036355861 0.024209296
5 1.4444430016 0.005467414 0.001221988
}\BiaspfnohomochiSigmaone
\pgfplotstableread{
beta Bayes DebiasBayes DebiasLasso
0 0.001383397 0.005604515 0.003177875
1 0.123256509 0.019453508 0.011714049
2 0.084265179 0.003992994 0.008970469
3 0.064312110 0.008131654 0.019866299
4 0.066967318 0.014106107 0.030450883
5 0.055819965 0.006161213 0.022271943
}\BiaspfnohomonSigmaone
\pgfplotstableread{
beta Bayes DebiasBayes DebiasLasso
0 0.002200622 0.0005885165 0.0004309773
1 0.209249474 0.0309213007 0.0058489507
2 0.433571940 0.0500124761 0.0268647271
3 0.606398851 0.0267364013 0.0079332456
4 0.727287633 0.0412250312 0.0239651555
5 0.883873170 0.0031275531 0.0058530058
}\BiasponoheterSigmaone
\pgfplotstableread{
beta Bayes DebiasBayes DebiasLasso
0 0.0004517858 0.017746747 0.014082731
1 0.2440684313 0.015509021 0.011653418
2 0.4768083808 0.025150244 0.004387045
3 0.6944520531 0.036664305 0.026167961
4 0.8928678258 0.009734311 0.024369349
5 1.5364042550 0.005458272 0.002653375
}\BiasponohomochiSigmaone
\pgfplotstableread{
beta Bayes DebiasBayes DebiasLasso
0 0.0003041703 0.003508030 0.00230409
1 0.1652696245 0.034634330 0.01369987
2 0.1777309113 0.028141826 0.02851929
3 0.1271038110 0.010364185 0.02419351
4 0.1082425540 0.013614436 0.03484569
5 0.0426113377 0.009365502 0.01871877
}\BiasponohomonSigmaone
\pgfplotstableread{
beta Bayes DebiasBayes DebiasLasso
0 0.0001799846 0.009249999 0.007015938
1 0.2165986597 0.041992712 0.011857832
2 0.4451585143 0.041864287 0.010924753
3 0.6271699454 0.043673796 0.018468610
4 0.7847656468 0.027481907 0.008362704
5 1.0399637874 0.027770277 0.028437195
}\BiasptnoheterSigmaone
\pgfplotstableread{
beta Bayes DebiasBayes DebiasLasso
0 0.001088812 0.03254957 0.034389545
1 0.240153859 0.02675843 0.004916506
2 0.479685095 0.03712065 0.012825109
3 0.710580816 0.03511702 0.023680412
4 0.929041048 0.05068757 0.035217777
5 1.652809687 0.04848980 0.038397464
}\BiasptnohomochiSigmaone
\pgfplotstableread{
beta Bayes DebiasBayes DebiasLasso
0 0.001279871 0.005058094 0.004698813
1 0.182678905 0.046673645 0.016240946
2 0.215813474 0.031674831 0.023118812
3 0.183858651 0.023055634 0.029313589
4 0.140857600 0.016156247 0.031925811
5 0.065958831 0.012642770 0.038299084
}\BiasptnohomonSigmaone
\pgfplotstableread{
beta Bayes DebiasBayes DebiasLasso
0 0.01925207 0.01589084 0.09247683
1 0.22005008 0.05499430 0.04028210
2 0.41834090 0.08614162 0.12418467
3 0.64245754 0.14470681 0.16520591
4 0.79990089 0.12796052 0.32591465
5 0.40817257 0.08722685 0.18174314
}\BiaspfnoheterSigmatwo
\pgfplotstableread{
beta Bayes DebiasBayes DebiasLasso
0 0.01861526 0.02142069 0.10998463
1 0.23570318 0.05861051 0.03837186
2 0.44150229 0.07770014 0.11975520
3 0.65084910 0.13654482 0.14655078
4 0.83032609 0.12435710 0.33407180
5 0.43106131 0.09160899 0.18111704
}\BiaspfnohomochiSigmatwo
\pgfplotstableread{
beta Bayes DebiasBayes DebiasLasso
0 0.003968657 0.003084803 0.01835252
1 0.113172058 0.020778603 0.03020753
2 0.111124704 0.030065644 0.06907534
3 0.062742655 0.022823021 0.10995558
4 0.041286355 0.017870030 0.13628569
5 0.026227520 0.012266941 0.08989703
}\BiaspfnohomonSigmatwo
\pgfplotstableread{
beta Bayes DebiasBayes DebiasLasso
0 0.0003854495 0.01096429 0.09335733
1 0.2331656821 0.06218585 0.03936348
2 0.4632173206 0.08336527 0.10991463
3 0.6692672165 0.15808925 0.15691825
4 0.8253575326 0.13633374 0.31972620
5 0.4103119204 0.09383610 0.18014970
}\BiasponoheterSigmatwo
\pgfplotstableread{
beta Bayes DebiasBayes DebiasLasso
0 0.0000417 0.01837294 0.11257658
1 0.2421614 0.05504538 0.03360933
2 0.4672417 0.07278867 0.10560262
3 0.6924365 0.15432506 0.14150559
4 0.8859219 0.13005213 0.32324648
5 0.4421301 0.09837324 0.17942505
}\BiasponohomochiSigmatwo
\pgfplotstableread{
beta Bayes DebiasBayes DebiasLasso
0 0.003506978 0.000316633 0.01817540
1 0.184636804 0.037862667 0.03510348
2 0.248930128 0.055971252 0.07236356
3 0.212215880 0.047770395 0.12039333
4 0.123404768 0.030420800 0.15347303
5 0.041941850 0.016230821 0.10470971
}\BiasponohomonSigmatwo
\pgfplotstableread{
beta Bayes DebiasBayes DebiasLasso
0 0.0007001274 0.01035257 0.09325036
1 0.2377599729 0.05398724 0.02732927
2 0.4640292922 0.07979414 0.09458665
3 0.6941966000 0.16645638 0.14223467
4 0.8925797685 0.14864868 0.30480988
5 0.4387484356 0.09187280 0.15887714
}\BiasptnoheterSigmatwo
\pgfplotstableread{
beta Bayes DebiasBayes DebiasLasso
0 0.0010033 0.01765966 0.1096036
1 0.2443898 0.05233355 0.0271311
2 0.4785253 0.09218970 0.1101223
3 0.6974934 0.15740223 0.1260607
4 0.9078971 0.13117534 0.2983622
5 0.4587348 0.10312010 0.1683140
}\BiasptnohomochiSigmatwo
\pgfplotstableread{
beta Bayes DebiasBayes DebiasLasso
0 0.00540785 0.00046631 0.01902918
1 0.19907592 0.03883401 0.03311536
2 0.31594224 0.07267480 0.07249739
3 0.32860568 0.07351858 0.12647247
4 0.21528844 0.05249209 0.16446921
5 0.06759288 0.01927046 0.10592719
}\BiasptnohomonSigmatwo
\pgfplotstableread{
beta Bayes DebiasBayes DebiasLasso
0 0.08254837 0.5086305 0.5441416
1 0.25971430 0.2637163 0.2786553
2 0.46170014 0.3033568 0.3153846
3 0.65905194 0.3591929 0.3726307
4 0.84968007 0.4218454 0.4267199
5 1.29128758 0.4603910 0.4594298
}\RMSEpfnoheterSigmaone
\pgfplotstableread{
beta Bayes DebiasBayes DebiasLasso
0 0.06209489 0.6611826 0.6778838
1 0.25160765 0.2675794 0.2839598
2 0.48733017 0.3771388 0.3906908
3 0.72631279 0.4635989 0.4665322
4 0.95027629 0.5102706 0.5193839
5 1.73324747 0.5847290 0.5860584
}\RMSEpfnohomochiSigmaone
\pgfplotstableread{
beta Bayes DebiasBayes DebiasLasso
0 0.06815679 0.2451978 0.2804317
1 0.19176764 0.1098113 0.1133537
2 0.26271838 0.1660429 0.1713026
3 0.27672554 0.1957949 0.2049842
4 0.28410034 0.2180002 0.2322905
5 0.24572306 0.2394319 0.2590445
}\RMSEpfnohomonSigmaone
\pgfplotstableread{
beta Bayes DebiasBayes DebiasLasso
0 0.04539301 0.5078626 0.5396947
1 0.25982750 0.2650866 0.2876093
2 0.47776810 0.2918635 0.3088945
3 0.68868893 0.3547176 0.3620187
4 0.88052546 0.4314152 0.4370849
5 1.36814845 0.4550914 0.4526235
}\RMSEponoheterSigmaone
\pgfplotstableread{
beta Bayes DebiasBayes DebiasLasso
0 0.02811768 0.6416459 0.6581995
1 0.24828050 0.2525080 0.2779783
2 0.49103707 0.3650983 0.3790894
3 0.73751676 0.4876046 0.4917899
4 0.96280747 0.5483693 0.5549353
5 1.78887747 0.5720341 0.5695781
}\RMSEponohomochiSigmaone
\pgfplotstableread{
beta Bayes DebiasBayes DebiasLasso
0 0.0353110 0.2342876 0.2740181
1 0.2103085 0.1154194 0.1137970
2 0.3150441 0.1729335 0.1697419
3 0.3360473 0.2032149 0.2035512
4 0.3418442 0.2322440 0.2403927
5 0.2545867 0.2439282 0.2550367
}\RMSEponohomonSigmaone
\pgfplotstableread{
beta Bayes DebiasBayes DebiasLasso
0 0.06266688 0.5109558 0.5367663
1 0.25404252 0.2516392 0.2783975
2 0.48008979 0.2982676 0.3109914
3 0.70624598 0.3689816 0.3798191
4 0.91121555 0.4225943 0.4320719
5 1.47325863 0.4628739 0.4641672
}\RMSEptnoheterSigmaone
\pgfplotstableread{
beta Bayes DebiasBayes DebiasLasso
0 0.03904932 0.6827244 0.6942065
1 0.25422926 0.2530075 0.2775954
2 0.49637675 0.3846165 0.3996284
3 0.73710518 0.4641862 0.4659896
4 0.97600696 0.5270730 0.5407012
5 1.84099104 0.5693631 0.5718987
}\RMSEptnohomochiSigmaone
\pgfplotstableread{
beta Bayes DebiasBayes DebiasLasso
0 0.04626209 0.2380254 0.2807367
1 0.22161738 0.1220298 0.1164194
2 0.34942491 0.1761957 0.1717340
3 0.39417550 0.2128253 0.2157517
4 0.38773882 0.2362682 0.2429289
5 0.27160848 0.2506018 0.2717188
}\RMSEptnohomonSigmaone
\pgfplotstableread{
beta Bayes DebiasBayes DebiasLasso
0 0.09168509 0.2163082 0.2046474
1 0.27463131 0.3054994 0.3104504
2 0.47637559 0.2479545 0.2600198
3 0.70733907 0.2712974 0.2602138
4 0.90924111 0.2656003 0.3858627
5 0.49621814 0.2501143 0.2675998
}\RMSEpfnoheterSigmatwo
\pgfplotstableread{
beta Bayes DebiasBayes DebiasLasso
0 0.08779023 0.2299180 0.2245792
1 0.26106128 0.2656319 0.2720877
2 0.48299388 0.2550444 0.2657801
3 0.71286718 0.2826652 0.2627743
4 0.92871305 0.2820502 0.4032482
5 0.50143339 0.2504466 0.2700973
}\RMSEpfnohomochiSigmatwo
\pgfplotstableread{
beta Bayes DebiasBayes DebiasLasso
0 0.02958598 0.09069449 0.08294056
1 0.20510245 0.11858631 0.11920002
2 0.23966992 0.12577262 0.13439824
3 0.21783138 0.12206385 0.15986872
4 0.19774494 0.11978397 0.17797528
5 0.11721226 0.10809705 0.14070302
}\RMSEpfnohomonSigmatwo
\pgfplotstableread{
beta Bayes DebiasBayes DebiasLasso
0 0.01113178 0.2225581 0.2142471
1 0.26707937 0.2972665 0.3013392
2 0.48454461 0.2314504 0.2482455
3 0.71034970 0.2735636 0.2609650
4 0.91087610 0.2646515 0.3805860
5 0.47610880 0.2380804 0.2650259
}\RMSEponoheterSigmatwo
\pgfplotstableread{
beta Bayes DebiasBayes DebiasLasso
0 0.00019168 0.2336082 0.2315057
1 0.25160764 0.2581977 0.2692475
2 0.48637188 0.2539778 0.2669972
3 0.72393887 0.2779722 0.2557865
4 0.94569961 0.2771307 0.3924041
5 0.48940920 0.2476144 0.2717588
}\RMSEponohomochiSigmatwo
\pgfplotstableread{
beta Bayes DebiasBayes DebiasLasso
0 0.02842996 0.09851303 0.08989901
1 0.22966045 0.11707309 0.11997834
2 0.33936494 0.12602572 0.13516457
3 0.38408620 0.13598380 0.16585410
4 0.30282318 0.12907162 0.19242741
5 0.15056226 0.11501320 0.15130191
}\RMSEponohomonSigmatwo
\pgfplotstableread{
beta Bayes DebiasBayes DebiasLasso
0 0.0154824 0.2169399 0.2146040
1 0.2563912 0.2867962 0.2976811
2 0.4870332 0.2355932 0.2488787
3 0.7215580 0.2788700 0.2581164
4 0.9449017 0.2617744 0.3629847
5 0.4860447 0.2361386 0.2545549
}\RMSEptnoheterSigmatwo
\pgfplotstableread{
beta Bayes DebiasBayes DebiasLasso
0 0.02227146 0.2315107 0.2322819
1 0.25263231 0.2660762 0.2803563
2 0.49126995 0.2557743 0.2762522
3 0.72526834 0.2825386 0.2575012
4 0.95687267 0.2664669 0.3690936
5 0.49710271 0.2577885 0.2759601
}\RMSEptnohomochiSigmatwo
\pgfplotstableread{
beta Bayes DebiasBayes DebiasLasso
0 0.03176429 0.09775338 0.09127577
1 0.23455020 0.11603385 0.11992921
2 0.38856441 0.13797691 0.13812807
3 0.48110861 0.15255231 0.17382285
4 0.40918779 0.14570092 0.20269544
5 0.19329238 0.12088637 0.15347705
}\RMSEptnohomonSigmatwo
We compare our debiased Bayes method with both standard Bayesian inference based on the spike-and-slab prior and a frequentist benchmark---specifically, the debiased inference procedure using a LASSO pilot estimator. The methods under comparison are summarized as follows.
\begin{table}[H]
\centering
\caption{Methods Compared in the Simulation Study}
\label{tab:methods}
\begin{tabular}{@{}p{4cm}p{10cm}@{}}
\toprule
\textbf{Method} & \textbf{Description} \\
\midrule
\textbf{Bayes}
& Standard Bayesian inference using the spike-and-slab prior,
where the posterior is approximated via variational Bayes. \\[0.2cm]
\textbf{Debiased-Bayes}
& Our proposed debiased Bayesian inference procedure. \\[0.2cm]
\textbf{Debiased-LASSO}
& The debiased LASSO estimator of \citet{VandeGeer_OnAsymptotically_2014}. \\
\bottomrule
\end{tabular}
\end{table}
\noindent\textbf{Simulation design.} We fix the sample size at $n = 100$ and vary the ambient dimensionality of the covariates by setting $p = 50, 100,$ and $200$, respectively. The true regression coefficient vector $\beta_0$ is sparse with $5$ nonzero entries. Specifically, the first five elements of $\beta_0$ are set to $(0.25, 0.5, 0.75, 1, 2)$, while the remaining coefficients are zero.
\noindent
\textbf{Data-generating processes.}
We consider six simulation scenarios designed to examine the robustness of each method:
\begin{itemize}
\item[] \textbf{S1 (Homoskedastic Normal):} $\varepsilon_i \sim N(0, 1)$ for all $i = 1, \ldots, n$.
\item[] \textbf{S2 (Non-Normal Errors):} $\varepsilon_i \sim \chi^2(3) - 3$, i.e., centered chi-squared noise.
\item[] \textbf{S3 (Heteroskedastic Normal):} $\varepsilon_i \sim N(0, \sigma_{i}^2)$, where $\sigma_i = 1 + |X_{1,i}|$.
\end{itemize}
The heteroskedastic specification in S3 is similar to that in Section 4.2 of \citet{HouMaWang2023CompositeQuantile}. For S1–S3, the regressors are drawn independently as $X_i \sim N(0, \Theta_0^{-1})$, where $\Theta_0$ is a $p\times p$ diagonal matrix with entries $(1, 2, \ldots, p)$ on the diagonal.
\begin{itemize}
\item[] \textbf{S4–S6 (Correlated Covariates):}
The error terms are specified as in S1–S3, but the covariates $X_i$ are sampled from a correlated Gaussian distribution $N(0, \Theta_0^{-1})$, where $\Theta_0$ is a banded precision matrix with its $ij$th element defined by
\[
\theta_{0,ij} =
\begin{cases}
1, & \text{if } i = j, \\
0.5, & \text{if } |i - j| = 1, \\
0, & \text{otherwise}.
\end{cases}
\]
\end{itemize}
Each simulation scenario is replicated 1,000 times.
\noindent
\textbf{Implementation details.}
For the standard Bayesian method, we employ the variational Bayes algorithm of \citet{RaySzabo2022Variational} to approximate the posterior distribution. Each posterior sample consists of $B = 8{,}000$ draws from the variational Bayes approximated posterior. The spike-and-slab prior uses a Laplace slab with $\lambda = 1$ and a Beta$(1, p^u)$ hyperprior on the inclusion probability with $u = 1$. No additional tuning is performed. For the precision matrix estimator $\hat{\Theta}_n$, we adopt the nodewise LASSO regression of \citet{VandeGeer_OnAsymptotically_2014}, as implemented in the R package \texttt{hdi} \citep{DBMM2015hdi}. This estimator is used for both the Debiased Bayes and the Debiased LASSO methods to ensure comparability. Our numerical experiments suggest that the default tuning parameters in the package yield stable and robust performance. We anticipate that further optimization of hyper-parameters could potentially improve finite-sample performance, such exploration is beyond the scope of the present study.
Figures \ref{Fig: one} and \ref{Fig: two} illustrate the empirical coverage probabilities of the $95\%$ credible or confidence sets for the regression coefficients across the simulation scenarios. Specifically, Figure \ref{Fig: one} reports results under specifications S1 (first column), S2 (second column), and S3 (third column), with each row corresponding to a different dimensionality setting: $p = 50$, $p = 100$, and $p = 200$. On the horizontal axis, the label 0 represents the average coverage across all zero coefficients in $\beta_0$, whereas labels 1–5 correspond to the coverage probabilities for the nonzero coefficients $(0.25, 0.5, 0.75, 1, 2)$, respectively. Figure \ref{Fig: two} displays the analogous coverage results for specifications S4–S6, following the same layout and interpretation. We also compare the estimation accuracy of the posterior mean of the debiased posterior distribution against the standard Bayesian posterior mean and the frequentist debiased estimator. Figures \ref{Fig: three} and \ref{Fig: four} report the empirical bias of the three estimators, while Figures \ref{Fig: five} and \ref{Fig: six} present the corresponding RMSE.
\begin{figure}[!h]
\centering\scriptsize
\begin{tikzpicture}
\begin{groupplot}[group style={group name=myplots,group size=3 by 3,horizontal sep= 0.8cm,vertical sep=1.1cm},
grid = minor,
width = 0.375\textwidth,
xmax=5,xmin=0,
ymax=1,ymin=0,
every axis title/.style={below,at={(0.2,0.8)}},
xlabel=$\beta_0$,
x label style={at={(axis description cs:0.95,0.04)},anchor=south},
xtick={0,1,2,3,4,5},
ytick={0,0.5,0.95,1},
tick label style={/pgf/number format/fixed},
legend style={text=black,cells={align=center},row sep = 3pt,legend columns = -1, draw=none,fill=none},
cycle list={
{smooth,tension=0,color=black, mark=halfsquare*,every mark/.append style={rotate=270},mark size=1.5pt,line width=0.5pt},
{smooth,tension=0,color=blue, mark=halfsquare*,every mark/.append style={rotate=90},mark size=1.5pt,line width=0.5pt},
{smooth,tension=0,color=red, mark=10-pointed star,mark size=1.5pt,line width=0.5pt},
{smooth,tension=0,color=RoyalBlue1, mark=halfcircle*,every mark/.append style={rotate=90},mark size=1.5pt,line width=0.5pt},
{smooth,tension=0,color=RoyalBlue2, mark=halfcircle*,every mark/.append style={rotate=180},mark size=1.5pt,line width=0.5pt},
{smooth,tension=0,color=RoyalBlue3, mark=halfcircle*,every mark/.append style={rotate=270},mark size=1.5pt,line width=0.5pt},
{smooth,tension=0,color=RoyalBlue4, mark=halfcircle*,every mark/.append style={rotate=360},mark size=1.5pt,line width=0.5pt},
}
]
\nextgroupplot
\node[anchor=north] at (axis description cs: 0.25, 0.95) {\fontsize{5}{4}\selectfont \shortstack{
\\ \\
Homo-Normal\\
$p=50$
}};
\addplot[smooth,tension=0.5,color=NavyBlue, no markers,line width=0.25pt, densely dotted,forget plot] table[x = beta,y=alpha] from \CovpfnohomonSigmaone;
\addplot table[x = beta,y=Bayes] from \CovpfnohomonSigmaone;
\addplot table[x = beta,y=DebiasBayes] from \CovpfnohomonSigmaone;
\addplot table[x = beta,y=DebiasLasso] from \CovpfnohomonSigmaone;
\nextgroupplot
\node[anchor=north] at (axis description cs: 0.25, 0.95) {\fontsize{5}{4}\selectfont \shortstack{
\\ \\
Homo-Chi\\
$p=50$
}};
\addplot[smooth,tension=0.5,color=NavyBlue, no markers,line width=0.25pt, densely dotted,forget plot] table[x = beta,y=alpha] from \CovpfnohomochiSigmaone;
\addplot table[x = beta,y=Bayes] from \CovpfnohomochiSigmaone;
\addplot table[x = beta,y=DebiasBayes] from \CovpfnohomochiSigmaone;
\addplot table[x = beta,y=DebiasLasso] from \CovpfnohomochiSigmaone;
\nextgroupplot[legend style = {column sep = 7pt, legend to name = LegendMon1}]
\addplot[smooth,tension=0.5,color=NavyBlue, no markers,line width=0.25pt, densely dotted,forget plot] table[x = beta,y=alpha] from \CovpfnoheterSigmaone;
\addplot table[x = beta,y=Bayes] from \CovpfnoheterSigmaone;
\addplot table[x = beta,y=DebiasBayes] from \CovpfnoheterSigmaone;
\addplot table[x = beta,y=DebiasLasso] from \CovpfnoheterSigmaone;
\node[anchor=north] at (axis description cs: 0.25, 0.95) {\fontsize{5}{4}\selectfont \shortstack{
\\ \\
Hetero\\
$p=50$
}};
\nextgroupplot
\node[anchor=north] at (axis description cs: 0.25, 0.95) {\fontsize{5}{4}\selectfont \shortstack{
\\ \\
Homo-Normal\\
$p=100$
}};
\addplot[smooth,tension=0.5,color=NavyBlue, no markers,line width=0.25pt, densely dotted,forget plot] table[x = beta,y=alpha] from \CovponohomonSigmaone;
\addplot table[x = beta,y=Bayes] from \CovponohomonSigmaone;
\addplot table[x = beta,y=DebiasBayes] from \CovponohomonSigmaone;
\addplot table[x = beta,y=DebiasLasso] from \CovponohomonSigmaone;
\nextgroupplot
\node[anchor=north] at (axis description cs: 0.25, 0.95) {\fontsize{5}{4}\selectfont \shortstack{
\\ \\
Homo-Chi\\
$p=100$
}};
\addplot[smooth,tension=0.5,color=NavyBlue, no markers,line width=0.25pt, densely dotted,forget plot] table[x = beta,y=alpha] from \CovponohomochiSigmaone;
\addplot table[x = beta,y=Bayes] from \CovponohomochiSigmaone;
\addplot table[x = beta,y=DebiasBayes] from \CovponohomochiSigmaone;
\addplot table[x = beta,y=DebiasLasso] from \CovponohomochiSigmaone;
\nextgroupplot[legend style = {column sep = 3.5pt, legend to name = LegendMon2}]
\node[anchor=north] at (axis description cs: 0.25, 0.95) {\fontsize{5}{4}\selectfont \shortstack{
\\ \\
Hetero\\
$p=100$
}};
\addplot[smooth,tension=0.5,color=NavyBlue, no markers,line width=0.25pt, densely dotted,forget plot] table[x = beta,y=alpha] from \CovponoheterSigmaone;
\addplot table[x = beta,y=Bayes] from \CovponoheterSigmaone;
\addplot table[x = beta,y=DebiasBayes] from \CovponoheterSigmaone;
\addplot table[x = beta,y=DebiasLasso] from \CovponoheterSigmaone;
\nextgroupplot
\node[anchor=north] at (axis description cs: 0.25, 0.95) {\fontsize{5}{4}\selectfont \shortstack{
\\ \\
Homo-Normal\\
$p=200$
}};
\addplot[smooth,tension=0.5,color=NavyBlue, no markers,line width=0.25pt, densely dotted,forget plot] table[x = beta,y=alpha] from \ptonohomon;
\addplot table[x = beta,y=Bayes] from \CovptnohomonSigmaone;
\addplot table[x = beta,y=DebiasBayes] from \CovptnohomonSigmaone;
\addplot table[x = beta,y=DebiasLasso] from \CovptnohomonSigmaone;
\nextgroupplot
\node[anchor=north] at (axis description cs: 0.25, 0.95) {\fontsize{5}{4}\selectfont \shortstack{
\\ \\
Homo-Chi\\
$p=200$
}};
\addplot[smooth,tension=0.5,color=NavyBlue, no markers,line width=0.25pt, densely dotted,forget plot] table[x = beta,y=alpha] from \CovptnohomochiSigmaone;
\addplot table[x = beta,y=Bayes] from \CovptnohomochiSigmaone;
\addplot table[x = beta,y=DebiasBayes] from \CovptnohomochiSigmaone;
\addplot table[x = beta,y=DebiasLasso] from \CovptnohomochiSigmaone;
\nextgroupplot[legend style = {column sep = 3.5pt, legend to name = LegendMon3}]
\node[anchor=north] at (axis description cs: 0.25, 0.95) {\fontsize{5}{4}\selectfont \shortstack{
\\ \\
Hetero\\
$p=200$
}};
\addplot[smooth,tension=0.5,color=NavyBlue, no markers,line width=0.25pt, densely dotted,forget plot] table[x = beta,y=alpha] from \CovptnoheterSigmaone;
\addplot table[x = beta,y=Bayes] from \CovptnoheterSigmaone;
\addplot table[x = beta,y=DebiasBayes] from \CovptnoheterSigmaone;
\addplot table[x = beta,y=DebiasLasso] from \CovptnoheterSigmaone;
\addlegendentry{Bayes};
\addlegendentry{Debiased-Bayes};
\addlegendentry{Debiased-LASSO};
\end{groupplot}
\node at ($(myplots c2r1) + (0,-2.25cm)$) {\ref{LegendMon1}};
\node at ($(myplots c2r2) + (0,-2.25cm)$) {\ref{LegendMon2}};
\node at ($(myplots c2r3) + (0,-2.25cm)$) {\ref{LegendMon3}};
\end{tikzpicture}
\caption{Coverage rates corresponding to different values of $\beta_0$ under settings S1 (first column), S2 (second column), and S3 (third column). The dashed line represents the 95\% benchmark.} \label{Fig: one}
\end{figure}
\begin{figure}[!h]
\centering\scriptsize
\begin{tikzpicture}
\begin{groupplot}[group style={group name=myplots,group size=3 by 3,horizontal sep= 0.8cm,vertical sep=1.1cm},
grid = minor,
width = 0.375\textwidth,
xmax=5,xmin=0,
ymax=1,ymin=0,
every axis title/.style={below,at={(0.2,0.8)}},
xlabel=$\beta_0$,
x label style={at={(axis description cs:0.95,0.04)},anchor=south},
xtick={0,1,2,3,4,5},
ytick={0,0.5,0.95,1},
tick label style={/pgf/number format/fixed},
legend style={text=black,cells={align=center},row sep = 3pt,legend columns = -1, draw=none,fill=none},
cycle list={
{smooth,tension=0,color=black, mark=halfsquare*,every mark/.append style={rotate=270},mark size=1.5pt,line width=0.5pt},
{smooth,tension=0,color=blue, mark=halfsquare*,every mark/.append style={rotate=90},mark size=1.5pt,line width=0.5pt},
{smooth,tension=0,color=red, mark=10-pointed star,mark size=1.5pt,line width=0.5pt},
{smooth,tension=0,color=RoyalBlue1, mark=halfcircle*,every mark/.append style={rotate=90},mark size=1.5pt,line width=0.5pt},
{smooth,tension=0,color=RoyalBlue2, mark=halfcircle*,every mark/.append style={rotate=180},mark size=1.5pt,line width=0.5pt},
{smooth,tension=0,color=RoyalBlue3, mark=halfcircle*,every mark/.append style={rotate=270},mark size=1.5pt,line width=0.5pt},
{smooth,tension=0,color=RoyalBlue4, mark=halfcircle*,every mark/.append style={rotate=360},mark size=1.5pt,line width=0.5pt},
}
]
\nextgroupplot
\node[anchor=north] at (axis description cs: 0.25, 0.95) {\fontsize{5}{4}\selectfont \shortstack{
\\ \\
Homo-Normal\\
$p=50$
}};
\addplot[smooth,tension=0.5,color=NavyBlue, no markers,line width=0.25pt, densely dotted,forget plot] table[x = beta,y=alpha] from \pfnohomon;
\addplot table[x = beta,y=Bayes] from \pfnohomon;
\addplot table[x = beta,y=DebiasBayes] from \pfnohomon;
\addplot table[x = beta,y=DebiasLasso] from \pfnohomon;
\nextgroupplot
\node[anchor=north] at (axis description cs: 0.25, 0.95) {\fontsize{5}{4}\selectfont \shortstack{
\\ \\
Homo-Chi\\
$p=50$
}};
\addplot[smooth,tension=0.5,color=NavyBlue, no markers,line width=0.25pt, densely dotted,forget plot] table[x = beta,y=alpha] from \pfnohomochi;
\addplot table[x = beta,y=Bayes] from \pfnohomochi;
\addplot table[x = beta,y=DebiasBayes] from \pfnohomochi;
\addplot table[x = beta,y=DebiasLasso] from \pfnohomochi;
\nextgroupplot[legend style = {column sep = 7pt, legend to name = LegendMon12}]
\addplot[smooth,tension=0.5,color=NavyBlue, no markers,line width=0.25pt, densely dotted,forget plot] table[x = beta,y=alpha] from \pfnoheter;
\addplot table[x = beta,y=Bayes] from \pfnoheter;
\addplot table[x = beta,y=DebiasBayes] from \pfnoheter;
\addplot table[x = beta,y=DebiasLasso] from \pfnoheter;
\node[anchor=north] at (axis description cs: 0.25, 0.95) {\fontsize{5}{4}\selectfont \shortstack{
\\ \\
Hetero\\
$p=50$
}};
\nextgroupplot
\node[anchor=north] at (axis description cs: 0.25, 0.95) {\fontsize{5}{4}\selectfont \shortstack{
\\ \\
Homo-Normal\\
$p=100$
}};
\addplot[smooth,tension=0.5,color=NavyBlue, no markers,line width=0.25pt, densely dotted,forget plot] table[x = beta,y=alpha] from \ponohomon;
\addplot table[x = beta,y=Bayes] from \ponohomon;
\addplot table[x = beta,y=DebiasBayes] from \ponohomon;
\addplot table[x = beta,y=DebiasLasso] from \ponohomon;
\nextgroupplot
\node[anchor=north] at (axis description cs: 0.25, 0.95) {\fontsize{5}{4}\selectfont \shortstack{
\\ \\
Homo-Chi\\
$p=100$
}};
\addplot[smooth,tension=0.5,color=NavyBlue, no markers,line width=0.25pt, densely dotted,forget plot] table[x = beta,y=alpha] from \ponohomochi;
\addplot table[x = beta,y=Bayes] from \ponohomochi;
\addplot table[x = beta,y=DebiasBayes] from \ponohomochi;
\addplot table[x = beta,y=DebiasLasso] from \ponohomochi;
\nextgroupplot[legend style = {column sep = 3.5pt, legend to name = LegendMon22}]
\node[anchor=north] at (axis description cs: 0.25, 0.95) {\fontsize{5}{4}\selectfont \shortstack{
\\ \\
Hetero\\
$p=100$
}};
\addplot[smooth,tension=0.5,color=NavyBlue, no markers,line width=0.25pt, densely dotted,forget plot] table[x = beta,y=alpha] from \ponoheter;
\addplot table[x = beta,y=Bayes] from \ponoheter;
\addplot table[x = beta,y=DebiasBayes] from \ponoheter;
\addplot table[x = beta,y=DebiasLasso] from \ponoheter;
\nextgroupplot
\node[anchor=north] at (axis description cs: 0.25, 0.95) {\fontsize{5}{4}\selectfont \shortstack{
\\ \\
Homo-Normal\\
$p=200$
}};
\addplot[smooth,tension=0.5,color=NavyBlue, no markers,line width=0.25pt, densely dotted,forget plot] table[x = beta,y=alpha] from \ptonohomon;
\addplot table[x = beta,y=Bayes] from \ptonohomon;
\addplot table[x = beta,y=DebiasBayes] from \ptonohomon;
\addplot table[x = beta,y=DebiasLasso] from \ptonohomon;
\nextgroupplot
\node[anchor=north] at (axis description cs: 0.25, 0.95) {\fontsize{5}{4}\selectfont \shortstack{
\\ \\
Homo-Chi\\
$p=200$
}};
\addplot[smooth,tension=0.5,color=NavyBlue, no markers,line width=0.25pt, densely dotted,forget plot] table[x = beta,y=alpha] from \ptonohomochi;
\addplot table[x = beta,y=Bayes] from \ptonohomochi;
\addplot table[x = beta,y=DebiasBayes] from \ptonohomochi;
\addplot table[x = beta,y=DebiasLasso] from \ptonohomochi;
\nextgroupplot[legend style = {column sep = 3.5pt, legend to name = LegendMon32}]
\node[anchor=north] at (axis description cs: 0.25, 0.95) {\fontsize{5}{4}\selectfont \shortstack{
\\ \\
Hetero\\
$p=200$
}};
\addplot[smooth,tension=0.5,color=NavyBlue, no markers,line width=0.25pt, densely dotted,forget plot] table[x = beta,y=alpha] from \ptonoheter;
\addplot table[x = beta,y=Bayes] from \ptonoheter;
\addplot table[x = beta,y=DebiasBayes] from \ptonoheter;
\addplot table[x = beta,y=DebiasLasso] from \ptonoheter;
\addlegendentry{Bayes};
\addlegendentry{Debiased-Bayes};
\addlegendentry{Debiased-LASSO};
\end{groupplot}
\node at ($(myplots c2r1) + (0,-2.25cm)$) {\ref{LegendMon12}};
\node at ($(myplots c2r2) + (0,-2.25cm)$) {\ref{LegendMon22}};
\node at ($(myplots c2r3) + (0,-2.25cm)$) {\ref{LegendMon32}};
\end{tikzpicture}
\caption{Coverage rates corresponding to different values of $\beta_0$ under settings S4 (first column), S5 (second column), and S6 (third column). The dashed line represents the 95\% benchmark.} \label{Fig: two}
\end{figure}
We summarize our findings as follows. The standard Bayesian approach based on the spike-and-slab prior performs well in identifying the zero coefficients, as reflected by favorable average coverage, low bias, and small RMSE for these components. However, its performance deteriorates substantially for the nonzero coefficients,\footnote{By contrast, \citet{RaySzabo2022Variational} considered much larger nonzero entries of $\beta_0$ (as large as 10), whereas ours are more moderate.} and becomes more fragile when the error distribution departs from the baseline homoscedastic Gaussian specification. In contrast, our proposed debiased Bayes approach consistently outperforms the standard method, achieving markedly better estimation accuracy and more reliable uncertainty quantification across all designs. Relative to the frequentist debiased LASSO benchmark, our method delivers comparable performance under the independent-covariate settings (S1–S3). Notably, in the more challenging correlated-covariate scenarios (S4–S6), the debiased Bayes procedure tends to attain substantially improved coverage, lower bias, and smaller RMSE for relatively larger coefficients. Overall, these results demonstrate the robustness and accuracy of the debaised Bayes approach across diverse conditions and the importance of debiasing standard Bayesian procedures in high-dimensional settings. It provides clear improvements over conventional Bayesian inference and exhibits complementary strengths relative to the leading frequentist benchmark.
\begin{figure}[!h]
\centering\scriptsize
\begin{tikzpicture}
\begin{groupplot}[group style={group name=myplots,group size=3 by 3,horizontal sep= 0.8cm,vertical sep=1.1cm},
grid = minor,
width = 0.375\textwidth,
xmax=5,xmin=0,
ymax=1.66,ymin=0,
every axis title/.style={below,at={(0.2,0.8)}},
xlabel=$\beta_0$,
x label style={at={(axis description cs:0.95,0.04)},anchor=south},
xtick={0,1,2,3,4,5},
ytick={0,0.5,1,1.5},
tick label style={/pgf/number format/fixed},
legend style={text=black,cells={align=center},row sep = 3pt,legend columns = -1, draw=none,fill=none},
cycle list={
{smooth,tension=0,color=black, mark=halfsquare*,every mark/.append style={rotate=270},mark size=1.5pt,line width=0.5pt},
{smooth,tension=0,color=blue, mark=halfsquare*,every mark/.append style={rotate=90},mark size=1.5pt,line width=0.5pt},
{smooth,tension=0,color=red, mark=10-pointed star,mark size=1.5pt,line width=0.5pt},
{smooth,tension=0,color=RoyalBlue1, mark=halfcircle*,every mark/.append style={rotate=90},mark size=1.5pt,line width=0.5pt},
{smooth,tension=0,color=RoyalBlue2, mark=halfcircle*,every mark/.append style={rotate=180},mark size=1.5pt,line width=0.5pt},
{smooth,tension=0,color=RoyalBlue3, mark=halfcircle*,every mark/.append style={rotate=270},mark size=1.5pt,line width=0.5pt},
{smooth,tension=0,color=RoyalBlue4, mark=halfcircle*,every mark/.append style={rotate=360},mark size=1.5pt,line width=0.5pt},
}
]
\nextgroupplot
\node[anchor=north] at (axis description cs: 0.25, 0.95) {\fontsize{5}{4}\selectfont \shortstack{
\\ \\
Homo-Normal\\
$p=50$
}};
\addplot table[x = beta,y=Bayes] from \BiaspfnohomonSigmaone;
\addplot table[x = beta,y=DebiasBayes] from \BiaspfnohomonSigmaone;
\addplot table[x = beta,y=DebiasLasso] from \BiaspfnohomonSigmaone;
\nextgroupplot
\node[anchor=north] at (axis description cs: 0.25, 0.95) {\fontsize{5}{4}\selectfont \shortstack{
\\ \\
Homo-Chi\\
$p=50$
}};
\addplot table[x = beta,y=Bayes] from \BiaspfnohomochiSigmaone;
\addplot table[x = beta,y=DebiasBayes] from \BiaspfnohomochiSigmaone;
\addplot table[x = beta,y=DebiasLasso] from \BiaspfnohomochiSigmaone;
\nextgroupplot[legend style = {column sep = 7pt, legend to name = LegendMon13}]
\addplot table[x = beta,y=Bayes] from \BiaspfnoheterSigmaone;
\addplot table[x = beta,y=DebiasBayes] from \BiaspfnoheterSigmaone;
\addplot table[x = beta,y=DebiasLasso] from \BiaspfnoheterSigmaone;
\node[anchor=north] at (axis description cs: 0.25, 0.95) {\fontsize{5}{4}\selectfont \shortstack{
\\ \\
Hetero\\
$p=50$
}};
\nextgroupplot
\node[anchor=north] at (axis description cs: 0.25, 0.95) {\fontsize{5}{4}\selectfont \shortstack{
\\ \\
Homo-Normal\\
$p=100$
}};
\addplot table[x = beta,y=Bayes] from \BiasponohomonSigmaone;
\addplot table[x = beta,y=DebiasBayes] from \BiasponohomonSigmaone;
\addplot table[x = beta,y=DebiasLasso] from \BiasponohomonSigmaone;
\nextgroupplot
\node[anchor=north] at (axis description cs: 0.25, 0.95) {\fontsize{5}{4}\selectfont \shortstack{
\\ \\
Homo-Chi\\
$p=100$
}};
\addplot table[x = beta,y=Bayes] from \BiasponohomochiSigmaone;
\addplot table[x = beta,y=DebiasBayes] from \BiasponohomochiSigmaone;
\addplot table[x = beta,y=DebiasLasso] from \BiasponohomochiSigmaone;
\nextgroupplot[legend style = {column sep = 3.5pt, legend to name = LegendMon23}]
\node[anchor=north] at (axis description cs: 0.25, 0.95) {\fontsize{5}{4}\selectfont \shortstack{
\\ \\
Hetero\\
$p=100$
}};
\addplot table[x = beta,y=Bayes] from \BiasponoheterSigmaone;
\addplot table[x = beta,y=DebiasBayes] from \BiasponoheterSigmaone;
\addplot table[x = beta,y=DebiasLasso] from \BiasponoheterSigmaone;
\nextgroupplot
\node[anchor=north] at (axis description cs: 0.25, 0.95) {\fontsize{5}{4}\selectfont \shortstack{
\\ \\
Homo-Normal\\
$p=200$
}};
\addplot table[x = beta,y=Bayes] from \BiasptnohomonSigmaone;
\addplot table[x = beta,y=DebiasBayes] from \BiasptnohomonSigmaone;
\addplot table[x = beta,y=DebiasLasso] from \BiasptnohomonSigmaone;
\nextgroupplot
\node[anchor=north] at (axis description cs: 0.25, 0.95) {\fontsize{5}{4}\selectfont \shortstack{
\\ \\
Homo-Chi\\
$p=200$
}};
\addplot table[x = beta,y=Bayes] from \BiasptnohomochiSigmaone;
\addplot table[x = beta,y=DebiasBayes] from \BiasptnohomochiSigmaone;
\addplot table[x = beta,y=DebiasLasso] from \BiasptnohomochiSigmaone;
\nextgroupplot[legend style = {column sep = 3.5pt, legend to name = LegendMon33}]
\node[anchor=north] at (axis description cs: 0.25, 0.95) {\fontsize{5}{4}\selectfont \shortstack{
\\ \\
Heter\\
$p=200$
}};
\addplot table[x = beta,y=Bayes] from \BiasptnoheterSigmaone;
\addplot table[x = beta,y=DebiasBayes] from \BiasptnoheterSigmaone;
\addplot table[x = beta,y=DebiasLasso] from \BiasptnoheterSigmaone;
\addlegendentry{Bayes};
\addlegendentry{Debiased-Bayes};
\addlegendentry{Debiased-LASSO};
\end{groupplot}
\node at ($(myplots c2r1) + (0,-2.25cm)$) {\ref{LegendMon13}};
\node at ($(myplots c2r2) + (0,-2.25cm)$) {\ref{LegendMon23}};
\node at ($(myplots c2r3) + (0,-2.25cm)$) {\ref{LegendMon33}};
\end{tikzpicture}
\caption{Bias corresponding to different values of $\beta_0$ under settings S1 (first column), S2 (second column), and S3 (third column).} \label{Fig: three}
\end{figure}
\begin{figure}[!h]
\centering\scriptsize
\begin{tikzpicture}
\begin{groupplot}[group style={group name=myplots,group size=3 by 3,horizontal sep= 0.8cm,vertical sep=1.1cm},
grid = minor,
width = 0.375\textwidth,
xmax=5,xmin=0,
ymax=1,ymin=0,
every axis title/.style={below,at={(0.2,0.8)}},
xlabel=$\beta_0$,
x label style={at={(axis description cs:0.95,0.04)},anchor=south},
xtick={0,1,2,3,4,5},
ytick={0,0.5,1},
tick label style={/pgf/number format/fixed},
legend style={text=black,cells={align=center},row sep = 3pt,legend columns = -1, draw=none,fill=none},
cycle list={
{smooth,tension=0,color=black, mark=halfsquare*,every mark/.append style={rotate=270},mark size=1.5pt,line width=0.5pt},
{smooth,tension=0,color=blue, mark=halfsquare*,every mark/.append style={rotate=90},mark size=1.5pt,line width=0.5pt},
{smooth,tension=0,color=red, mark=10-pointed star,mark size=1.5pt,line width=0.5pt},
{smooth,tension=0,color=RoyalBlue1, mark=halfcircle*,every mark/.append style={rotate=90},mark size=1.5pt,line width=0.5pt},
{smooth,tension=0,color=RoyalBlue2, mark=halfcircle*,every mark/.append style={rotate=180},mark size=1.5pt,line width=0.5pt},
{smooth,tension=0,color=RoyalBlue3, mark=halfcircle*,every mark/.append style={rotate=270},mark size=1.5pt,line width=0.5pt},
{smooth,tension=0,color=RoyalBlue4, mark=halfcircle*,every mark/.append style={rotate=360},mark size=1.5pt,line width=0.5pt},
}
]
\nextgroupplot
\node[anchor=north] at (axis description cs: 0.25, 0.95) {\fontsize{5}{4}\selectfont \shortstack{
\\ \\
Homo-Normal\\
$p=50$
}};
\addplot table[x = beta,y=Bayes] from \BiaspfnohomonSigmatwo;
\addplot table[x = beta,y=DebiasBayes] from \BiaspfnohomonSigmatwo;
\addplot table[x = beta,y=DebiasLasso] from \BiaspfnohomonSigmatwo;
\nextgroupplot
\node[anchor=north] at (axis description cs: 0.25, 0.95) {\fontsize{5}{4}\selectfont \shortstack{
\\ \\
Homo-Chi\\
$p=50$
}};
\addplot table[x = beta,y=Bayes] from \BiaspfnohomochiSigmatwo;
\addplot table[x = beta,y=DebiasBayes] from \BiaspfnohomochiSigmatwo;
\addplot table[x = beta,y=DebiasLasso] from \BiaspfnohomochiSigmatwo;
\nextgroupplot[legend style = {column sep = 7pt, legend to name = LegendMon14}]
\addplot table[x = beta,y=Bayes] from \BiaspfnoheterSigmatwo;
\addplot table[x = beta,y=DebiasBayes] from \BiaspfnoheterSigmatwo;
\addplot table[x = beta,y=DebiasLasso] from \BiaspfnoheterSigmatwo;
\node[anchor=north] at (axis description cs: 0.25, 0.95) {\fontsize{5}{4}\selectfont \shortstack{
\\ \\
Hetero\\
$p=50$
}};
\nextgroupplot
\node[anchor=north] at (axis description cs: 0.25, 0.95) {\fontsize{5}{4}\selectfont \shortstack{
\\ \\
Homo-Normal\\
$p=100$
}};
\addplot table[x = beta,y=Bayes] from \BiasponohomonSigmatwo;
\addplot table[x = beta,y=DebiasBayes] from \BiasponohomonSigmatwo;
\addplot table[x = beta,y=DebiasLasso] from \BiasponohomonSigmatwo;
\nextgroupplot
\node[anchor=north] at (axis description cs: 0.25, 0.95) {\fontsize{5}{4}\selectfont \shortstack{
\\ \\
Homo-Chi\\
$p=100$
}};
\addplot table[x = beta,y=Bayes] from \BiasponohomochiSigmatwo;
\addplot table[x = beta,y=DebiasBayes] from \BiasponohomochiSigmatwo;
\addplot table[x = beta,y=DebiasLasso] from \BiasponohomochiSigmatwo;
\nextgroupplot[legend style = {column sep = 3.5pt, legend to name = LegendMon24}]
\node[anchor=north] at (axis description cs: 0.25, 0.95) {\fontsize{5}{4}\selectfont \shortstack{
\\ \\
Heter\\
$p=100$
}};
\addplot table[x = beta,y=Bayes] from \BiasponoheterSigmatwo;
\addplot table[x = beta,y=DebiasBayes] from \BiasponoheterSigmatwo;
\addplot table[x = beta,y=DebiasLasso] from \BiasponoheterSigmatwo;
\nextgroupplot
\node[anchor=north] at (axis description cs: 0.25, 0.95) {\fontsize{5}{4}\selectfont \shortstack{
\\ \\
Homo-Normal\\
$p=200$
}};
\addplot table[x = beta,y=Bayes] from \BiasptnohomonSigmatwo;
\addplot table[x = beta,y=DebiasBayes] from \BiasptnohomonSigmatwo;
\addplot table[x = beta,y=DebiasLasso] from \BiasptnohomonSigmatwo;
\nextgroupplot
\node[anchor=north] at (axis description cs: 0.25, 0.95) {\fontsize{5}{4}\selectfont \shortstack{
\\ \\
Homo-Chi\\
$p=200$
}};
\addplot table[x = beta,y=Bayes] from \BiasptnohomochiSigmatwo;
\addplot table[x = beta,y=DebiasBayes] from \BiasptnohomochiSigmatwo;
\addplot table[x = beta,y=DebiasLasso] from \BiasptnohomochiSigmatwo;
\nextgroupplot[legend style = {column sep = 3.5pt, legend to name = LegendMon34}]
\node[anchor=north] at (axis description cs: 0.25, 0.95) {\fontsize{5}{4}\selectfont \shortstack{
\\ \\
Hetero\\
$p=200$
}};
\addplot table[x = beta,y=Bayes] from \BiasptnoheterSigmatwo;
\addplot table[x = beta,y=DebiasBayes] from \BiasptnoheterSigmatwo;
\addplot table[x = beta,y=DebiasLasso] from \BiasptnoheterSigmatwo;
\addlegendentry{Bayes};
\addlegendentry{Debiased-Bayes};
\addlegendentry{Debiased-LASSO};
\end{groupplot}
\node at ($(myplots c2r1) + (0,-2.25cm)$) {\ref{LegendMon14}};
\node at ($(myplots c2r2) + (0,-2.25cm)$) {\ref{LegendMon24}};
\node at ($(myplots c2r3) + (0,-2.25cm)$) {\ref{LegendMon34}};
\end{tikzpicture}
\caption{Bias corresponding to different values of $\beta_0$ under settings S4 (first column), S5 (second column), and S6 (third column).} \label{Fig: four}
\end{figure}
\begin{figure}[!h]
\centering\scriptsize
\begin{tikzpicture}
\begin{groupplot}[group style={group name=myplots,group size=3 by 3,horizontal sep= 0.8cm,vertical sep=1.1cm},
grid = minor,
width = 0.375\textwidth,
xmax=5,xmin=0,
ymax=1.85,ymin=0,
every axis title/.style={below,at={(0.2,0.8)}},
xlabel=$\beta_0$,
x label style={at={(axis description cs:0.95,0.04)},anchor=south},
xtick={0,1,2,3,4,5},
ytick={0,0.5,1,1.5},
tick label style={/pgf/number format/fixed},
legend style={text=black,cells={align=center},row sep = 3pt,legend columns = -1, draw=none,fill=none},
cycle list={
{smooth,tension=0,color=black, mark=halfsquare*,every mark/.append style={rotate=270},mark size=1.5pt,line width=0.5pt},
{smooth,tension=0,color=blue, mark=halfsquare*,every mark/.append style={rotate=90},mark size=1.5pt,line width=0.5pt},
{smooth,tension=0,color=red, mark=10-pointed star,mark size=1.5pt,line width=0.5pt},
{smooth,tension=0,color=RoyalBlue1, mark=halfcircle*,every mark/.append style={rotate=90},mark size=1.5pt,line width=0.5pt},
{smooth,tension=0,color=RoyalBlue2, mark=halfcircle*,every mark/.append style={rotate=180},mark size=1.5pt,line width=0.5pt},
{smooth,tension=0,color=RoyalBlue3, mark=halfcircle*,every mark/.append style={rotate=270},mark size=1.5pt,line width=0.5pt},
{smooth,tension=0,color=RoyalBlue4, mark=halfcircle*,every mark/.append style={rotate=360},mark size=1.5pt,line width=0.5pt},
}
]
\nextgroupplot
\node[anchor=north] at (axis description cs: 0.25, 0.95) {\fontsize{5}{4}\selectfont \shortstack{
\\ \\
Homo-Normal\\
$p=50$
}};
\addplot table[x = beta,y=Bayes] from \RMSEpfnohomonSigmaone;
\addplot table[x = beta,y=DebiasBayes] from \RMSEpfnohomonSigmaone;
\addplot table[x = beta,y=DebiasLasso] from \RMSEpfnohomonSigmaone;
\nextgroupplot
\node[anchor=north] at (axis description cs: 0.25, 0.95) {\fontsize{5}{4}\selectfont \shortstack{
\\ \\
Homo-Chi\\
$p=50$
}};
\addplot table[x = beta,y=Bayes] from \RMSEpfnohomochiSigmaone;
\addplot table[x = beta,y=DebiasBayes] from \RMSEpfnohomochiSigmaone;
\addplot table[x = beta,y=DebiasLasso] from \RMSEpfnohomochiSigmaone;
\nextgroupplot[legend style = {column sep = 7pt, legend to name = LegendMon15}]
\addplot table[x = beta,y=Bayes] from \RMSEpfnoheterSigmaone;
\addplot table[x = beta,y=DebiasBayes] from \RMSEpfnoheterSigmaone;
\addplot table[x = beta,y=DebiasLasso] from \RMSEpfnoheterSigmaone;
\node[anchor=north] at (axis description cs: 0.25, 0.95) {\fontsize{5}{4}\selectfont \shortstack{
\\ \\
Hetero\\
$p=50$
}};
\nextgroupplot
\node[anchor=north] at (axis description cs: 0.25, 0.95) {\fontsize{5}{4}\selectfont \shortstack{
\\ \\
Homo-Normal\\
$p=100$
}};
\addplot table[x = beta,y=Bayes] from \RMSEponohomonSigmaone;
\addplot table[x = beta,y=DebiasBayes] from \RMSEponohomonSigmaone;
\addplot table[x = beta,y=DebiasLasso] from \RMSEponohomonSigmaone;
\nextgroupplot
\node[anchor=north] at (axis description cs: 0.25, 0.95) {\fontsize{5}{4}\selectfont \shortstack{
\\ \\
Homo-Chi\\
$p=100$
}};
\addplot table[x = beta,y=Bayes] from \RMSEponohomochiSigmaone;
\addplot table[x = beta,y=DebiasBayes] from \RMSEponohomochiSigmaone;
\addplot table[x = beta,y=DebiasLasso] from \RMSEponohomochiSigmaone;
\nextgroupplot[legend style = {column sep = 3.5pt, legend to name = LegendMon25}]
\node[anchor=north] at (axis description cs: 0.25, 0.95) {\fontsize{5}{4}\selectfont \shortstack{
\\ \\
Hetero\\
$p=100$
}};
\addplot table[x = beta,y=Bayes] from \RMSEponoheterSigmaone;
\addplot table[x = beta,y=DebiasBayes] from \RMSEponoheterSigmaone;
\addplot table[x = beta,y=DebiasLasso] from \RMSEponoheterSigmaone;
\nextgroupplot
\node[anchor=north] at (axis description cs: 0.25, 0.95) {\fontsize{5}{4}\selectfont \shortstack{
\\ \\
Homo-Normal\\
$p=200$
}};
\addplot table[x = beta,y=Bayes] from \RMSEptnohomonSigmaone;
\addplot table[x = beta,y=DebiasBayes] from \RMSEptnohomonSigmaone;
\addplot table[x = beta,y=DebiasLasso] from \RMSEptnohomonSigmaone;
\nextgroupplot
\node[anchor=north] at (axis description cs: 0.25, 0.95) {\fontsize{5}{4}\selectfont \shortstack{
\\ \\
Homo-Chi\\
$p=200$
}};
\addplot table[x = beta,y=Bayes] from \RMSEptnohomochiSigmaone;
\addplot table[x = beta,y=DebiasBayes] from \RMSEptnohomochiSigmaone;
\addplot table[x = beta,y=DebiasLasso] from \RMSEptnohomochiSigmaone;
\nextgroupplot[legend style = {column sep = 3.5pt, legend to name = LegendMon35}]
\node[anchor=north] at (axis description cs: 0.25, 0.95) {\fontsize{5}{4}\selectfont \shortstack{
\\ \\
Hetero\\
$p=200$
}};
\addplot table[x = beta,y=Bayes] from \RMSEptnoheterSigmaone;
\addplot table[x = beta,y=DebiasBayes] from \RMSEptnoheterSigmaone;
\addplot table[x = beta,y=DebiasLasso] from \RMSEptnoheterSigmaone;
\addlegendentry{Bayes};
\addlegendentry{Debiased-Bayes};
\addlegendentry{Debiased-LASSO};
\end{groupplot}
\node at ($(myplots c2r1) + (0,-2.25cm)$) {\ref{LegendMon15}};
\node at ($(myplots c2r2) + (0,-2.25cm)$) {\ref{LegendMon25}};
\node at ($(myplots c2r3) + (0,-2.25cm)$) {\ref{LegendMon35}};
\end{tikzpicture}
\caption{RMSE corresponding to different values of $\beta_0$ under settings S4 (first column), S5 (second column), and S6 (third column).} \label{Fig: five}
\end{figure}
\begin{figure}[!h]
\centering\scriptsize
\begin{tikzpicture}
\begin{groupplot}[group style={group name=myplots,group size=3 by 3,horizontal sep= 0.8cm,vertical sep=1.1cm},
grid = minor,
width = 0.375\textwidth,
xmax=5,xmin=0,
ymax=1,ymin=0,
every axis title/.style={below,at={(0.2,0.8)}},
xlabel=$\beta_0$,
x label style={at={(axis description cs:0.95,0.04)},anchor=south},
xtick={0,1,2,3,4,5},
ytick={0,0.5,1},
tick label style={/pgf/number format/fixed},
legend style={text=black,cells={align=center},row sep = 3pt,legend columns = -1, draw=none,fill=none},
cycle list={
{smooth,tension=0,color=black, mark=halfsquare*,every mark/.append style={rotate=270},mark size=1.5pt,line width=0.5pt},
{smooth,tension=0,color=blue, mark=halfsquare*,every mark/.append style={rotate=90},mark size=1.5pt,line width=0.5pt},
{smooth,tension=0,color=red, mark=10-pointed star,mark size=1.5pt,line width=0.5pt},
{smooth,tension=0,color=RoyalBlue1, mark=halfcircle*,every mark/.append style={rotate=90},mark size=1.5pt,line width=0.5pt},
{smooth,tension=0,color=RoyalBlue2, mark=halfcircle*,every mark/.append style={rotate=180},mark size=1.5pt,line width=0.5pt},
{smooth,tension=0,color=RoyalBlue3, mark=halfcircle*,every mark/.append style={rotate=270},mark size=1.5pt,line width=0.5pt},
{smooth,tension=0,color=RoyalBlue4, mark=halfcircle*,every mark/.append style={rotate=360},mark size=1.5pt,line width=0.5pt},
}
]
\nextgroupplot
\node[anchor=north] at (axis description cs: 0.25, 0.95) {\fontsize{5}{4}\selectfont \shortstack{
\\ \\
Homo-Normal\\
$p=50$
}};
\addplot table[x = beta,y=Bayes] from \RMSEpfnohomonSigmatwo;
\addplot table[x = beta,y=DebiasBayes] from \RMSEpfnohomonSigmatwo;
\addplot table[x = beta,y=DebiasLasso] from \RMSEpfnohomonSigmatwo;
\nextgroupplot
\node[anchor=north] at (axis description cs: 0.25, 0.95) {\fontsize{5}{4}\selectfont \shortstack{
\\ \\
Homo-Chi\\
$p=50$
}};
\addplot table[x = beta,y=Bayes] from \RMSEpfnohomochiSigmatwo;
\addplot table[x = beta,y=DebiasBayes] from \RMSEpfnohomochiSigmatwo;
\addplot table[x = beta,y=DebiasLasso] from \RMSEpfnohomochiSigmatwo;
\nextgroupplot[legend style = {column sep = 7pt, legend to name = LegendMon16}]
\addplot table[x = beta,y=Bayes] from \RMSEpfnoheterSigmatwo;
\addplot table[x = beta,y=DebiasBayes] from \RMSEpfnoheterSigmatwo;
\addplot table[x = beta,y=DebiasLasso] from \RMSEpfnoheterSigmatwo;
\node[anchor=north] at (axis description cs: 0.25, 0.95) {\fontsize{5}{4}\selectfont \shortstack{
\\ \\
Hetero\\
$p=50$
}};
\nextgroupplot
\node[anchor=north] at (axis description cs: 0.25, 0.95) {\fontsize{5}{4}\selectfont \shortstack{
\\ \\
Homo-Normal\\
$p=100$
}};
\addplot table[x = beta,y=Bayes] from \RMSEponohomonSigmatwo ;
\addplot table[x = beta,y=DebiasBayes] from \RMSEponohomonSigmatwo ;
\addplot table[x = beta,y=DebiasLasso] from \RMSEponohomonSigmatwo ;
\nextgroupplot
\node[anchor=north] at (axis description cs: 0.25, 0.95) {\fontsize{5}{4}\selectfont \shortstack{
\\ \\
Homo-Chi\\
$p=100$
}};
\addplot table[x = beta,y=Bayes] from \RMSEponohomochiSigmatwo;
\addplot table[x = beta,y=DebiasBayes] from \RMSEponohomochiSigmatwo;
\addplot table[x = beta,y=DebiasLasso] from \RMSEponohomochiSigmatwo;
\nextgroupplot[legend style = {column sep = 3.5pt, legend to name = LegendMon26}]
\node[anchor=north] at (axis description cs: 0.25, 0.95) {\fontsize{5}{4}\selectfont \shortstack{
\\ \\
Heter\\
$p=100$
}};
\addplot table[x = beta,y=Bayes] from \RMSEponoheterSigmatwo ;
\addplot table[x = beta,y=DebiasBayes] from \RMSEponoheterSigmatwo ;
\addplot table[x = beta,y=DebiasLasso] from \RMSEponoheterSigmatwo ;
\nextgroupplot
\node[anchor=north] at (axis description cs: 0.25, 0.95) {\fontsize{5}{4}\selectfont \shortstack{
\\ \\
Homo-Normal\\
$p=200$
}};
\addplot table[x = beta,y=Bayes] from \RMSEptnohomonSigmatwo;
\addplot table[x = beta,y=DebiasBayes] from \RMSEptnohomonSigmatwo;
\addplot table[x = beta,y=DebiasLasso] from \RMSEptnohomonSigmatwo;
\nextgroupplot
\node[anchor=north] at (axis description cs: 0.25, 0.95) {\fontsize{5}{4}\selectfont \shortstack{
\\ \\
Homo-Chi\\
$p=200$
}};
\addplot table[x = beta,y=Bayes] from \RMSEptnohomochiSigmatwo;
\addplot table[x = beta,y=DebiasBayes] from \RMSEptnohomochiSigmatwo;
\addplot table[x = beta,y=DebiasLasso] from \RMSEptnohomochiSigmatwo;
\nextgroupplot[legend style = {column sep = 3.5pt, legend to name = LegendMon36}]
\node[anchor=north] at (axis description cs: 0.25, 0.95) {\fontsize{5}{4}\selectfont \shortstack{
\\ \\
Hetero\\
$p=200$
}};
\addplot table[x = beta,y=Bayes] from \RMSEptnoheterSigmatwo;
\addplot table[x = beta,y=DebiasBayes] from \RMSEptnoheterSigmatwo;
\addplot table[x = beta,y=DebiasLasso] from \RMSEptnoheterSigmatwo;
\addlegendentry{Bayes};
\addlegendentry{Debiased-Bayes};
\addlegendentry{Debiased-LASSO};
\end{groupplot}
\node at ($(myplots c2r1) + (0,-2.25cm)$) {\ref{LegendMon16}};
\node at ($(myplots c2r2) + (0,-2.25cm)$) {\ref{LegendMon26}};
\node at ($(myplots c2r3) + (0,-2.25cm)$) {\ref{LegendMon36}};
\end{tikzpicture}
\caption{RMSE corresponding to different values of $\beta_0$ under settings S4 (first column), S5 (second column), and S6 (third column).} \label{Fig: six}
\end{figure}
\newpage
\section{Conclusion}\label{Sec: Conclusion}
In this paper, we develop a new debiased Bayesian inferential method for high-dimensional linear regression models. The construction resembles the frequentist debiasing step, whereas the key difference is that we correct the entire posterior distribution rather than the point estimator. Our approach is tailored to building confidence intervals based on sparsity-inducing priors such as spike-and-slab or horseshoe type priors of the regression coefficients. We establish the frequentist validity of our proposal in the general setup and also provide low-level conditions. It is straightforward to observe that our debiasing step easily extends to other parametric or semiparametric models. To demonstrate the versatility of our general methodology, we mention two possible extensions.
One may consider the generalized linear model for the conditional density function of the dependent variable given covariates as follows:
\begin{align*}
p(y|x)\propto\exp\left(y\cdot x^\intercal\beta_0-d(x^\intercal\beta_0)\right),
\end{align*}
where $d$ is a convex link function. In this case, our debiased Bayesian method builds on
\begin{align*}
\tilde{\beta}=\beta+\hat{\Theta}_n\left[\sum_{i=1}^{n}W_{ni}X_i(Y_i-d^\prime(X_i^\intercal\beta))\right],
\end{align*}
where $\hat{\Theta}_n$ stands for a regularized inverse of the Hessian matrix $n^{-1}\sum_{i=1}^nd^{\prime\prime} (X_{i}^{\intercal}\hat\beta_n)X_iX_i^\intercal$, given some pilot estimator $\hat{\beta}_n$.
One may also incorporate group structure which often occurs in additive models \citep{Baietal2022Group}, multivariate outcome variables regression, and regressors collected over mixed frequencies \citep{MoglianiSimoni2021BMidas}. For example, consider the linear regression specified by the following form:
\begin{equation}
Y=\sum_{g=1}^G X_g^{\intercal}\beta_{g,0}+\varepsilon,
\end{equation}
where $\beta_{g,0}$ is an $m_g\times 1$ vector of coefficients, and $X_g$ is an $m_g\times 1$ vector of covariates corresponding to group $g = 1, \cdots,G$. The number of groups $G$ is potentially larger than the sample size $n$. In this case, our debiasing procedure is based on
\begin{align*}
\tilde{\beta}_{g}=\beta_{g}+\mathcal{E}_g^\intercal \hat{\Theta}_n\left[\frac{1}{n}\sum_{i=1}^{n}X_i(Y_i-X_i^\intercal\beta)\right],
\end{align*}
where $\mathcal{E}_g$ here is the matrix formed by columns of $I_p$ such that $\mathcal{E}_g\beta=\beta_g$ with $\beta=(\beta_1^\intercal,\ldots,\beta_G^\intercal)^\intercal$ and $p=m_1+\cdots+m_G$, $X_i=(X_{1,i}^\intercal,\ldots,X_{G,i}^\intercal)^\intercal$, and the initial posterior for $\beta$ can be obtained by the group spike-and-slab prior from \cite{Baietal2022Group}.
For both extensions, we believe our debiased Bayesian inference offers interesting new insights. We will defer the detailed theoretical development, as well as practical performance to some future work.