用于稀疏Probit回归的平均场变分贝叶斯方法

《Computational Statistics & Data Analysis》：Mean-field Variational Bayes for Sparse Probit Regression

　　奥古斯托·法萨诺|乔瓦尼·雷鲍多
都灵大学和卡洛·阿尔贝托学院，苏联联盟大道218/bis号，都灵，10134，意大利

**摘要**
我们考虑在probit链接下对二元结果进行贝叶斯变量选择，并对回归系数使用尖峰-平板先验。鉴于马尔可夫链蒙特卡洛（MCMC）采样器在高维情况下遇到的计算挑战，我们开发了一种平均场变分贝叶斯近似方法，其中所有变分因子都允许闭合形式更新，并且证据下界也可以以闭合形式获得。这反过来又允许开发出一种高效的坐标上升变分推断算法来找到变分参数的最优值。该方法产生了后验包含概率和参数估计，使得在单一框架内可以进行可解释的选择和预测。正如在模拟数据和真实数据应用中所展示的，所提出的方法能够成功识别出重要变量，并且比MCMC快几个数量级，同时保持相当的准确性。

**1. 引言**
从Albert和Chib（1993年）的开创性工作开始，贝叶斯二元回归在过去几年中受到了持续的关注，针对probit和logit链接出现了多种计算策略，包括采样（Durante, 2019; Holmes和Held, 2006; Polson等人, 2013; Zens等人, 2024）或变分近似（Consonni和Marin, 2007; Fasano等人, 2022; Jaakkola和Jordan, 2000; Zens和Steel, 2025）。另见Chopin和Ridgway（2017）以及Anceschi等人（2023）关于可能的计算方法的综述。

在高维问题中，一个核心的推断任务是变量选择，通常利用诱导稀疏性的先验。这种基于模型的方法在单一连贯的框架内产生后验包含概率和参数估计，但对后验计算的要求很高：因此，标准的马尔可夫链蒙特卡洛（MCMC）可能会遇到收敛问题，或者当p很大或p远大于n时速度变得非常慢。另见Chang等人（2016）；Chu等人（2025）关于贝叶斯probit模型中二元和多类别结果的MCMC策略及相关讨论。

在这项工作中，为了克服上述采样方法的局限性，我们开发了一种基于二元掩码的变分贝叶斯程序，用于稀疏probit回归，并采用尖峰-平板先验（Ormerod等人，2017）。为此，设yi∈{0, 1}为二元响应，xi∈Rp为单元i=1,…,n的p维协变量向量。对于j=1,…,p，定义变量γj∈{0, 1}表示预测变量j是否对输出有影响（γj=1）或没有影响（γj=0）。因此，利用probit模型关于潜在高斯变量部分观测值的特征（Albert和Chib，1993），我们将带有二元掩码的贝叶斯probit模型写为：
$$
y_i = 1 \{ z_i > 0 \}, \quad z_i \mid \beta, \gamma \sim \text{indN}(x_i^\top \Gamma \beta, 1), \quad i=1,\ldots,n, \quad \beta \sim N_p(0, \nu^2 I_p), \quad \gamma_j \sim \text{iid} \Bernoulli(\rho), \quad j=1,\ldots,p,
$$
其中1{A}表示事件A的指示函数，Ip是p维单位矩阵，β=(β_1,\ldots,β_p)^\top，γ=(γ_1,\ldots,γ_p)^\top，Γ=diag(γ)，而ν^2 > 0和ρ∈(0, 1)表示已知的超参数。因此，我们通过线性预测变量xi^\topΓβ来编码选择，同时保持β～N_p(0, ν^2 I_p)与γ独立。注意，如果我们设置β?=Γβ，则β?j=γ_jβ_j具有混合先验，其在0处的概率为1?ρ，在ν^2处的概率为ρ，这与β?的经典尖峰-平板先验表示一致。

我们的兴趣在于推断哪些变量真正影响输出以及这种影响的强度。设y=(y_1,\ldots,y_n)，因此我们研究后验包含概率（PIPs），Pr[γ_j=1|y]，以及γ_jβ_j的后验分布，j=1,\ldots,p。这些量没有闭合形式的表达式，通过MCMC方法获得可能计算成本很高。为了克服这些限制，在第2节中，我们提出了一种平均场变分贝叶斯方法，利用Ormerod等人（2017）为线性模型开发的方法。

**2. 稀疏probit的平均场变分贝叶斯**
在提出的方法中，我们用以下平均场变分贝叶斯（MFVB）近似来近似模型（1）的联合后验分布p(β, z, γ|y)：
$$
q(\beta, z, \gamma) = q(\beta) q(z) \prod_{j=1}^p q(\gamma_j).
$$
重要的是，平均场分解是在块（β, z, γ）之间以及γ块内部施加的，但在β内部不是。特别是，如下所示，q(β)是一个具有完整协方差矩阵的多变量高斯分布，因此允许回归系数之间的依赖性，从而减轻了通常与完全分解的平均场近似相关联的贝叶斯不确定性低估问题，例如在高维设置或存在共线性时。注意，在具有高斯先验的probit情况下，Fasano等人（2022）开发了一种部分分解的变分贝叶斯方法，旨在实现计算可扩展性，并在高维设置下提供理论保证。然而，他们方法的计算可行性依赖于对回归系数的高斯先验假设，并不直接适用于带有潜在包含指示器的尖峰-平板变量选择，这需要不同的算法解决方案。

最优因子的形式由平均场方程给出；参见例如Blei等人（2017）。根据这些方程，如果θ=(θ_1,\ldots,θ_K)是所有模型参数的向量，分为K个子组件θ_1,\ldots,θ_K，并且后验分布p(θ|y)由变分密度q(θ)=∏_{k=1}^K q(θ_k)近似，则最优变分近似q(θ_k)满足logq(θ_k)=Eq(θ?k)[logp(θ_k|θ?k, y)]+const，其中θ?k表示去除了第k个子组件θ_k的向量θ；这种符号约定在以下内容中也将使用。我们现在推导最优q(β)、q(z)和q(γ_j)，j=1,\ldots,p的形式。我们用?(t)和Φ(t)分别表示标准正态密度和累积分布函数。类似地，?_p(x; Σ)表示均值为0且协方差矩阵为Σ的p变量高斯随机变量的密度，在x处评估。

**2.1. 最优变分因子q(β)**
从最优平均场方程（Blei等人，2017）来看，最优q(β)必须满足：
$$
\log q(\beta) = Eq(z) q(\gamma) [logp(β|z, γ, y)] + const,
$$
其中完整的条件p(β|z, γ, y)=p(β|z, γ)与具有设计XΓ和单位误差方差的高斯回归中的后验一致。另见Ormerod等人（2017）。可以很容易地推导出来，因为一般来说，
$$
\log p(β|z, γ) = \log p(z|β, γ) + \log p(β) + const.
$$
设X=[x_1,\ldots, x_n]^\top \in R^{n\times p}，对于模型（1）有：
$$
\log p(z|β, γ) + \log p(β) = -\frac{1}{2}(z - X^\Gamma \beta)^\top (z - X^\Gamma \beta) - \frac{1}{2} \nu^2 \beta^\top \beta + const,
$$
其中二次项可以展开为：
$$
(z - X^\Gamma \beta)^\top (z - X^\Gamma \beta) = z^\top z - 2 \beta^\top \Gamma X^\top z + \beta^\top \Gamma X^\top X^\Gamma \beta = z^\top z - 2 \beta^\top \Gamma X^\top z + \beta^\top (\Gamma G^\Gamma) \beta,
$$
因此，
$$
\log p(β|z, γ) = -\frac{1}{2} \beta^\top (\nu^2 I_p + \Gamma G^\Gamma) \beta + \beta^\top \Gamma X^\top z + const,
$$
即，
$$
p(β|z, γ) = \phi_p(β - V^\Gamma X^\top z; V) \quad \text{其中} V = (\nu^2 I_p + \Gamma G^\Gamma)^{-1}.
$$
因此，我们有：
$$
\log q(β) = Eq(z) q(\gamma) [-\frac{1}{2} \beta^\top (\nu^2 I_p + \Gamma G^\Gamma) \beta + \beta^\top \Gamma X^\top z] + const = -\frac{1}{2} \beta^\top (\nu^2 I_p + Eq(\gamma)[\Gamma G^\Gamma]) \beta + \beta^\top Eq(\gamma)[\Gamma] X^\top Eq(z) + const,
$$
从中我们可以识别出多变量高斯核的对数。为了方便书写参数，我们定义z?=Eq(z)[z]和W=Eq(γ)[Γ]=diag(w)，其中w=Eq(γ)[γ]∈[0,1]^p且Ω=Eq(γ)[γγ^\top]=W(I_p?W)+ww^\top，使得Ω_{jj}=w_j=Eq(γ_j)[γ_j]且Ω_{jk}=w_jw_k对于j≠k。现在利用ΓG^\Gamma=(X^\Gamma)^\top(X^\Gamma)=[γ_1X:1\ldotsγ_pX:p]^\top[γ_1X:1\ldotsγ_pX:p]=[γ_jγ_kX:j^\topX:k]j,k=1\leq p=G⊙(γγ^\top)，其中⊙表示Hadamard（逐元素）乘积，X: j表示X的第j列，我们有：
$$
Eq(\gamma)[\Gamma G^\Gamma] = G⊙Ω.
$$
因此，最优变分因子q(β)满足：
$$
q(β) = \phi_p(β - \mu; \Sigma),
$$
其中Σ=(ν^2 I_p + G⊙Ω)^{-1}且μ=Σ WX^\top z?。注意，由于G=X^\top X?0且Ω=Eq(γ)[γγ^\top]?0，Schur乘积定理给出G⊙Ω?0。因此对于任何ν^2>0，ν^2 I_p+G⊙Ω是正定的。如果某些w_j退化为零或接近零的值，可以使用一些简化或近似，避免通过求逆p×p精度矩阵直接计算Σ所需的O(p^3)成本。更多细节见第3节和补充材料。

**2.2. 最优变分因子q(z)**
从z=(z_1,\ldots,z_n)^\top的完整条件开始，注意到：
$$
p(z|β, γ, y) \sim p(z|β, γ) p(y|z) \sim \prod_{i=1}^n \phi(z_i - x_i^\top \Gamma \beta) 1\{z_i \in A_yi},
$$
其中A_yi=(0,∞)如果yi=1，否则为(?∞,0]。也就是说，根据Albert和Chib（1993），在(β, γ, y)的条件下，zi作为具有均值xi^\topΓβ和单位方差的高斯随机变量独立分布，截断到集合A_yi。从最优性方程中，对于z∈A_y，其中A_y=A_y1×?×A_yn，得到：
$$
\log q(z) = \sum_{i=1}^n \log q(z_i),
$$
其中，对于z_i∈A_yi，
$$
\log q(z_i) = Eq(β) q(\gamma) [log\phi(z_i - x_i^\top \Gamma \beta)] + const = -\frac{1}{2} Eq(β) q(\gamma)[(z_i - x_i^\top \Gamma \beta)^2] + const = -\frac{1}{2} z_i^2 + z_i x_i x_i^\top Eq(β) q(\gamma)[\Gamma β] + const = -\frac{1}{2} z_i^2 + z_i x_i x_i^\top Eq(β)[\Gamma] Eq(β)[β] + const = -\frac{1}{2} z_i^2 + z_i x_i x_i^\top W^\mu + const = -\frac{1}{2} (z_i - x_i^\top W^\mu)^2 + const.
$$
因此，最优变分近似q(z)分解为q(z)=∏_{i=1}^n q(z_i)，其中q(z_i)∝?(z_i?mi)1{z_i∈A_yi}，其中mi=xi^\top W^\mu。因此，q(z_i)是截断到A_yi的高斯随机变量N(mi,1)的密度。设ki=2^i?1∈{?1,1}，截断高斯随机变量的标准公式给出：
$$
(3) z?_i = Eq(z_i)[z_i] = mi + ki \lambda(k_i m_i),
$$
(4) Eq(z_i)[z_i^2] = 1 + m_i z?_i,
$$
其中λ(t)=?(t)/Φ(t)表示逆Mills比率。我们只需要z?=(z?_1,\ldots,z?_n)^\top在(2)中；(4)中的第二时刻对于监控证据下界（ELBO）很有用，这在补充材料中有报告。

表达式(3)和(4)来自截断高斯分布的标准性质。实际上，设T～N(m,1)并在T>0的条件下，那么E[T]=m+λ(m)且Var(T)=1?mλ(m)?λ(m)^2。因此E[T^2]=Var(T)+E[T]^2=1?mλ(m)?λ(m)^2+(m+λ(m))^2=1+mE[T]。如果我们在T≤0的条件下，我们有E[T]=m?λ(?m)且Var(T)=1+mλ(?m)?λ(?m)^2。适应之前的计算，我们在这种情况下也有E[T^2]=1+mE[T]。方程(3)和(4)简洁地表示了这些结果。

**2.3. 最优变分因子q(γ_j)**
首先，注意到γ_j的完整条件形式为：
$$
p(\gamma_j|γ?j, z, β, y) \sim p(β)p(γ)p(z|β, γ)p(y|z) \sim p(γ)p(z|β, γ) \sim \rho \gamma_j (1?ρ) 1?γ_j \exp{-\frac{1}{2}(z - X^\Gamma \beta)^\top (z - X^\Gamma \beta)},
$$
因此：
$$
(5) \log p(\gamma_j|γ?j, z, β, y) = \gamma_j \logit(\rho) - \frac{1}{2} \beta^\top \Gamma X^\top X^\Gamma \beta + \beta^\top \Gamma X^\top z + const.
$$
在将γ_j从完整条件中分离出来并取变分期望以计算q(γ_j)之前，让我们回想一下，如果V～Bernoulli(p)，那么p(v)=p^v(1?p)^{1?v}∝(p^1?p)v∝exp{v·logit(p)}，因此logp(v)=v·logit(p)+const。现在让我们在(5)的第二项和第三项中分离出γ_j项。它满足XΓβ=γ_jβ_jX:j+∑_{j′≠j} γ_j′β_j′X:j′，因此β^\top \Gamma X^\top z=γ_jβ_jX:j^\top z+∑_{j′≠j} γ_j′β_j′X:j′^\top z=γ_jβ_jX:j^\top z+const，且β^\top \Gamma X^\top X^\Gamma β=(∑_{l=1}^p γ_l β_l X:l)^\top(∑_{k=1}^p γ_k β_k X:k)=γ_jβ_jX:j^\top ∑_{k=1}^p X:k γ_k β_k+(∑_{l≠j} γ_l β_l X:l^\top) γ_jβ_jX:j+const=γ_jβ_j^2X:j^\top X:j+2γ_jβ_jX:j^\top X:j+2γ_jβ_jX:j^\top X?jΓ?j,?jβ?j+const，
其中X?j对应于去除了第j列的X，Γ?j,?j是去除了第j行和第j列的矩阵Γ。综合起来，我们得到：
$$
\log p(\gamma_j|γ?j, z, β, y) = \gamma_j [logit(\rho) - \frac{1}{2} β_j^2 X:j^\top X:j + β_j X:j^\top (z - X?jΓ?j,?jβ?j)] + const.
$$
对所有其他变分分布取期望，我们得到最优q(γ_j)必须满足：
$$
\log q(γ_j) = Eq(β) q(γ?j) q(z) [logp(γ_j|γ?j, z, β, y)] + const = γ_j \cdot [logit(\rho) - \frac{1}{2}(Σ_jj + μ_j^2) G_jj + μ_j X:j^\top z??∑_{k≠j}(Σ_jk + μ_j μ_k) w_k G_jk] + const,
$$
从中我们得到最优变分因子q(γ_j)是Bernoulli(w_j)随机变量的概率质量函数，其中w_j=expit(η_j)，且
$$
(6) \eta_j = \logit(\rho) - \frac{1}{2}(Σ_jj + μ_j^2) G_jj + μ_j X:j^\top z??∑_{k≠j}(Σ_jk + μ_j μ_k) w_k G_jk.
$$

**3. 算法**
可以通过坐标上升变分推断（CAVI）算法获得变分参数的最优值，其中参数依次更新，直到它们的相对变化最大或ELBO的变化变得可以忽略不计。该程序在算法1中给出。

**下载：** 下载高分辨率图像（225KB）
**下载：** 下载全尺寸图像

**算法1. 稀疏probit回归的MFVB**
我们现在讨论一些关于该算法的计算和实际考虑。

**3.1. 大p的高效近似**
算法1使用精度矩阵Σ^?1=ν^2 I_p+G⊙Ω更新高斯因子q(β)，其中G和Ω如上所定义。直接计算Σ的成本为O(p^3)，对于大的p来说可能是禁止性的。实际上，当许多PIPs接近0或p远大于n时，可以使用以下简化（详细信息在补充材料中）。首先，当许多包含概率w_j接近零时，只有部分预测变量对线性预测变量有贡献；因此，可以将更新限制在活跃集Sw={j:w_j>ε}的大小为k，将w_j≤ε（例如，ε=10^?10）视为零，并保持相应的系数在其先验值。然后在对k维系统上执行q(β)的更新：如果k≤n，我们使用k×k的受限精度矩阵；如果k＞n，我们应用Woodbury恒等式，以便主要的线性代数涉及n×n矩阵。最后，在非常高的维度中，我们可以通过使用ηj的近似版本来近似q(γ_j)的更新，避免Σ的非对角项，从而节省额外的计算成本。注意，这样我们得到了一个介于我们的默认（具有完整协方差的q(β)）和完全分解的VB之间的变分算法，后者还假设多变量因子q(β)有更严格的独立结构。

**3.2. 初始化和停止规则**
迭代运行直到收敛，收敛定义为两次连续迭代之间变分系数均值和包含概率的最大相对变化，在欧几里得范数中测量，降至规定的容忍度以下。在我们的实验中，我们使用10^?4的容忍度。注意，在所有合成数据和真实数据的实验中，我们将变分参数初始化为μ_j=0和w_j=ρ对于所有j=1,\ldots,p（算法1，第3行）。将μ_j=0设置为中性起点，不会在潜在均值向量m=XWμ中引入初始方向偏差，因此与先验规范一致。将w_j初始化为ρ也是自然的，因为它与先验包含期望相匹配。在CAVI优化中，不同的初始化可能会导致不同的局部最优解。作为实际检查，可以从多个w的起始值运行算法，并保留最终ELBO最大的解。在我们的数值研究中，上述默认初始化提供了稳定的解，我们没有观察到这些初始化的小扰动对主要推断和预测摘要有显著变化。此外，当后验在多个几乎等价的变量选择配置上放置相似的质量（例如，在弱信号或高度相关的预测变量下），变分近似可能会集中在其中一个几乎等价的配置上，并产生过于自信的PIPs。在这种情况下，不同的初始化可能会选择不同的但实际上等价的子集，通常具有可比的ELBO值和预测性能，反映了实际的非可识别性而不是拟合模型的实质性不稳定性。因此，虽然这种行为可能会影响模型空间的贝叶斯不确定性量化，但在目标是预测或选择稀疏模型的点估计时，这不太成问题。超参数设置
为了运行算法，需要设置两个超参数 ν2（β 的先验方差）和 ρ（变量的先验包含概率）。在下面的模拟研究和应用中，这些超参数是根据以下数据驱动的程序设置的。首先，对于给定的 ρ 值，我们设置 ν2=ν02/(ρ·p)，以控制线性预测器的先验方差接近 ν02，假设预测器已经标准化（这在实践中总是建议的）或具有单位方差。ν02 的作用是控制线性预测器的先验分散，从而也控制隐含的先验预测概率。在标准化预测器和先验包含率 ρ 的情况下，重新缩放 ν2=ν02/(ρp) 可以使得 Var(xi?Γβ) 在先验上接近 ν02。接下来，我们将 ν02 设置为 25，相当于潜在 probit 规模上的标准差约为 5，这提供了一个信息量较弱的先验，允许较大的效应大小，同时避免了先验在零附近过于集中，遵循贝叶斯二元回归模型中的常见建议（Fasano 等人，2022 年；Gelman 等人，2008 年）。我们还在补充材料中报告了对一系列 ν02 值的敏感性分析。

然后，在可能将数据分为训练集和测试集之后，通过在训练集上进行 K 折交叉验证来选择 ρ 的最优值，其中观测值的分配按输出值 yi 分层，以使每个折中的 yi=1 的比例大致与整个训练数据集中的比例相同。对于可能的 ρ 值网格中的每个值，通常范围从 0.05 到 0.5，步长为 0.05，我们计算平均交叉验证偏差 (7)devCV=1K∑k=1Kdevk，其中 devk=?2∑i∈Ik[yilogp^i+(1?yi)log(1?p^i)]，Ik 是分配给第 k 折的观测值的索引集，p^i=Pr^[yi=1∣y?FOLD(i)] 是当模型在所有训练观测值集上拟合时，第 i 个训练观测值的成功估计概率，不包括属于第 i 个观测值的折中的观测值，其索引集表示为 FOLD(i)=Ikifi∈Ik。调用 W(?FOLD(i)) 和 μ(?FOLD(i)) 是从 MFVB 算法拟合 y?FOLD(i) 获得的 Γ 和 β 的变分均值估计，我们使用插值估计 p^i=Φ(xi?W(?FOLD(i))μ(?FOLD(i))) 来计算 p^i。

4. 模拟研究
我们研究了所提出算法在选择适当变量方面的性能，并将其性能与 MCMC 采样方法进行了比较，其中从 Gibbs 采样器获得了 10,000 个后验样本（在 1000 个燃烧期之后），其详细信息在补充材料中报告，该方法利用了完整条件表达式。我们考虑了两种不同的情况：一种是 p=200 和 n=1000，另一种更具挑战性的情况是 p=1000 和 n=500。在每种情况下，我们考虑了 50 个复制的数据集。在这两种情况下，每次复制中，非零系数的真实数量是 0.02·p，其中一半是分布在 ?3 和 ?1 之间的等间距点，另一半是分布在 1 和 3 之间的等间距点。协变量值是独立同分布的标准正态分布，结果输出值是根据模型 (1) 生成的。对于每个数据集，我们使用 MFVB 和 MCMC 计算参数的后验近似值，并比较结果。我们按照第 3 节中解释的方式固定 ν2 和 ρ，交叉验证程序中使用 K=5 折，ν02=25。然后我们比较了这些方法在运行时间、真正例率和假正例率以及测试集上的预测准确性方面的表现，测试集的大小与用于调整 ρ 的训练集不同。对于后验比较，如果估计的后验包含概率 Pr^[γj=1∣y] 超过 0.5，则认为变量 j 被包含；否则，该变量被排除。我们将真正例率计算为正确包含的变量数量除以活跃变量的总数，而假负例率是正确排除的变量数量除以非活跃变量的总数。对于具有协变量向量 xNEW 的新观测值 yNEW，使用插值估计 Pr^[yNEW=1∣y]=Φ(xNEW?Wμ) 来计算 MFVB 的预测概率，其中 W 是对角矩阵，包含 γ 的后验期望，μ 代表使用 MFVB 获得的 β 的后验期望。

4. 模拟研究
我们研究了所提出算法在选择适当变量方面的性能，并将其性能与 MCMC 采样方法进行了比较，其中从 Gibbs 采样器获得了 10,000 个后验样本（在 1000 个燃烧期之后），其详细信息在补充材料中报告，该方法利用了完整条件表达式。我们考虑了两种不同的情况：一种是 p=200 和 n=1000，另一种更具挑战性的情况是 p=1000 和 n=500。在每种情况下，我们考虑了 50 个复制的数据集。在这两种情况下，每次复制中，非零系数的真实数量是 0.02·p，其中一半是分布在 ?3 和 ?1 之间的等间距点，另一半是分布在 1 和 3 之间的等间距点。协变量值是独立同分布的标准正态分布，结果输出值是根据模型 (1) 生成的。对于每个数据集，我们使用 MFVB 和 MCMC 计算参数的后验近似值，并比较结果。我们按照第 3 节中解释的方式固定 ν2 和 ρ，交叉验证程序中使用 K=5 折，ν02=25。然后我们比较了这些方法在运行时间、真正例率和假负例率以及测试集上的预测准确性方面的表现，测试集的大小与用于调整 ρ 的训练集不同。对于后验比较，如果估计的后验包含概率 Pr^[γj=1∣y] 超过 0.5，则认为变量 j 被包含；否则，该变量被排除。我们将真正例率计算为正确包含的变量数量除以活跃变量的总数，而假负例率是正确排除的变量数量除以非活跃变量的总数。对于具有协变量向量 xNEW 的新观测值 yNEW，使用插值估计 Pr^[yNEW=1∣y]=Φ(xNEW?Wμ) 来计算 MFVB 的预测概率，其中 W 是对角矩阵，包含 γ 的后验期望，μ 代表使用 MFVB 获得的 β 的后验期望。对于 MCMC，预测概率是通过平均后验抽样来计算的。所有模拟都在配备 Apple M4 Pro 处理器的 12 核 Mac mini 上执行。

表 1 还报告了 βj 和 γj 的平均有效样本大小，作为 Gibbs 采样器的标准混合诊断指标。如表 1 所示，在 p < n 和 p > n 的设置中，MFVB 在第一种情况下快一个数量级以上，在第二种情况下快三个数量级以上。此外，在更简单的 p ? n 设置中，两种方法在变量选择和样本外偏差方面提供了几乎相同的结果，变量选择几乎完美。在更具挑战性的 p ? n 设置中，两种方法在真正例率和假负例率方面都表现良好，尽管两者之间存在一些小的差异。通过观察在重复模拟中作为数据生成回归参数 γj0βj0 的函数绘制的估计后验包含概率（以下简称 PIPs），可以更好地理解这些差异，对于第一种（图 1）和第二种（图 2）情况。总体而言，我们观察到 MFVB 的 PIPs 倾向于在边界 0 或 1 处退化，这与 Ormerod 等人（2017 年）对线性模型的理论结果一致，其中显示在适当的初始化下，它们收敛到真实的数据生成指标 γj0∈{0,1}。另一方面，MCMC 倾向于更不稀疏，可能包含比需要的更多变量：特别是对于 p ? n 的情况，这可能导致样本外偏差更高。

表 1. 按情况和方法的平均结果。括号中的标准差。
空单元格
p=200, n=1000
p=1000, n=500
指标
MFVB
MCMC
MFVB
MCMC
真正例率（%）
100 (0)
100 (0)
90.30 (7.17)
94.40 (4.70)
假负例率（%）
100 (0)
99.77 (0.35)
99.78 (0.23)
94.95 (3.60)
偏差
161.31 (15.44)
163.51 (15.92)
219.15 (25.59)
441.51 (29.47)
运行时间（秒）
0.93 (0.24)
48.54 (2.53)
27.34 (2.58)
609
58.84 (248)
54.62
βj 的有效样本大小—20
40.53 (542.30)
—
130
3.17 (187.32)
γj 的有效样本大小—3
394.44 (567.99)
—
339
9.19 (703.90)
下载：下载高分辨率图像（201KB）
下载：下载全尺寸图像

图 1. 对于 p=200, n=1000，通过 MCMC 和 MFVB 在 50 个模拟数据集中估计的真实参数值 γj0βj0 的后验包含概率（PIPs）。出于图形目的，0 和 1 之间的距离不是按比例的，我们将真实的 γj0βj0∈{±3,±1,0} 分布在它们实际值的附近。

图 2. 对于 p=1000, n=500，通过 MCMC 和 MFVB 在 50 个模拟数据集中估计的真实参数值 γj0βj0 的后验包含概率（PIPs）。出于图形目的，0 和 1 之间的距离不是按比例的，我们将真实的回归参数 γj0βj0 分布在它们实际值的附近。

现在更详细地考虑这两种情况，在第一种（p ? n）情况下，MFVB 的估计也与 MCMC 对非零 γj0βj0 值获得的 PIP 估计相匹配。如图 1 所示，在这种情况下，与问题维度相比，大量的观测值使得两种方法都能很好地学习哪些是真正活跃的变量。对于不活跃的变量，即数据生成参数 γj0 等于零的变量，MCMC 倾向于更加保守，估计的 PIPs 远离边界，而 MFVB 仍然显示出几乎确定性的分类，正如预期的那样。这种特殊属性对于高维设置中的 MFVB 似乎是有益的，特别是在 p ? n 的情况下，如第二种情况 p=1000, n=500 所示，结果在图 2 中报告。在这种情况下，两种方法对于真正例率和假负例率的表现都很好，尽管两者之间存在一些小的差异。通过观察在重复模拟中作为数据生成回归参数 γj0βj0 的函数绘制的估计后验包含概率（PIPs），可以更好地理解这些差异。一般来说，我们观察到 MFVB 的 PIPs 倾向于在边界 0 或 1 处退化，这与 Ormerod 等人（2017 年）对线性模型的理论结果一致，其中显示在适当的初始化下，它们收敛到真实的数据生成指标 γj0∈{0,1}。另一方面，MCMC 倾向于不太稀疏，可能包含比需要的更多变量：特别是对于 p ? n 的情况，这可能导致样本外偏差更高。

表 1. 按情况和方法的平均结果。括号中的标准差。
空单元格
p=200, n=1000
p=1000, n=500
指标
MFVB
MCMC
MFVB
MCMC
真正例率（%）
100 (0)
100 (0)
90.30 (7.17)
94.40 (4.70)
假负例率（%）
100 (0)
99.77 (0.35)
99.78 (0.23)
94.95 (3.60)
偏差
161.31 (15.44)
163.51 (15.92)
219.15 (25.59)
441.51 (29.47)
运行时间（秒）
0.93 (0.24)
48.54 (2.53)
27.34 (2.58)
609
58.84 (248)
54.62
βj 的有效样本大小—20
40.53 (542.30)
—
130
3.17 (187.32)
γj 的有效样本大小—3
394.44 (567.99)
—
339
9.19 (703.90)
下载：下载高分辨率图像（201KB）
下载：下载全尺寸图像

图 1. 对于 p=200, n=1000，通过 MCMC 和 MFVB 在 50 个模拟数据集中估计的真实参数值 γj0βj0 的后验包含概率（PIPs）。出于图形目的，0 和 1 之间的距离不是按比例的，我们将真实的 γj0βj0∈{±3,±1,0} 分布在它们实际值的附近。

现在更详细地考虑这两种情况，在第一种（p ? n）情况下，MFVB 的估计也与 MCMC 对非零 γj0βj0 值获得的 PIP 估计相匹配。如图 1 所示，在这种情况下，与问题维度相比，大量的观测值使得两种方法都能很好地学习哪些是真正活跃的变量。对于不活跃的变量，即数据生成参数 γj0 等于零的变量，MCMC 倾向于更加保守，估计的 PIPs 远离边界，而 MFVB 仍然显示出几乎确定性的分类，正如预期的那样。这种特殊属性对于高维设置中的 MFVB 似乎是有益的，特别是在 p ? n 的情况下，如第二种情况 p=1000, n=500 所示，结果在图 2 中报告。在这种情况下，两种方法对于绝对值较大的真正活跃变量的 PIPs 估计是一致的。对于真正活跃变量的 βj0 值不是足够远离零的变量，MCMC 的估计更加分散，尽管两种方法在应该包含哪些变量方面有一致的意见。最后，对于不活跃的变量，MCMC 倾向于更加保守，给出更高的 PIPs。实际上，我们观察到 MFVB 更倾向于简约的表示，而 MCMC 由于更加保守，倾向于更容易过拟合数据。因此，MFVB 观察到的 PIPs 的退化可能是一个有利的属性，可能在 p 远大于 n 的设置中很有用。

虽然接近退化的 PIPs 可以产生清晰的排名、简约的选择和高效的计算，但它也与平均场变分近似的已知限制一致，即低估了后验不确定性。特别是，假设的包含指标的分解可能会产生过于集中的变分后验，使得 PIPs 可能比 MCMC 更接近 0 或 1。这主要影响贝叶斯不确定性量化（例如，可信区间或包含指标的后验分散），并且当信号较弱或预测变量高度相关时可能导致过于自信的包含或排除决策。对于预测或选定模型的点估计，我们在实验中没有观察到不利影响，但我们提醒不要将 MFVB 的 PIPs 解释为完全校准的后验概率。请注意，MFVB 的独立性假设主要涉及包含指标，而回归系数的变分分布没有分解。由于 q(β) 保留了完整的协方差结构，近似可以表示与相关预测变量相关的系数的后验依赖性。

5. 实际数据应用
我们在两个实际数据应用上评估了所提出的方法：LSVT 语音康复数据集，最初由 Tsanas 等人（2013 年）分析，以及 Alzheimer 数据集，由 Craig-Schapiro 等人（2011 年）研究，并在 Fasano 等人（2022 年）中进一步考虑。前者数据集有 n=126 个观测值，经过清洗后有 p=309 个预测变量，结果代表帕金森病患者语音表现的可接受性。后者处理了一个大规模情况：它包含 n=333 个关于阿尔茨海默病早期发作的观测值，以及相关的 p=9036 个预测变量（通过在原始数据集中添加成对交互项获得），显示了所提出方法在 MCMC 计算上不可行的设置中的实际用途。在这两种应用中，我们应用了第 2 节中描述的 MFVB 程序，ρ 和 ν2 的设置如第 3 节中解释的那样，ν02=25，并使用整个数据集作为训练集。

5.1. LSVT 语音康复应用
我们首先在 LSVT 语音康复数据集上评估了所提出的 MFVB 方法，该数据集涉及将帕金森病（PD）患者的语音表现分类为“可接受”（yi=1）或“不可接受”（yi=0）。该数据集包含来自 14 名患者的 n=126 个发音，可在 UCI 仓库公开获取。最初由 Tsanas 等人（2013 年）分析，该数据集包含 310 个特征，其中两个（Data_length 和 Ea2）的变异性非常有限，因此从分析中移除。其余的协变量被标准化为零均值和单位方差。还添加了一个截距，得到 p=309 个预测变量。Tsanas 等人（2013 年）的原始目标是使用自动化预测变量复制专家的评分并识别有影响力的特征：我们的分析使用完全贝叶斯、稀疏且可解释的模型来实现相同的目标。我们将所提出的 MFVB 方法与上述超参数设置的 Gibbs 采样基线（MCMC）进行了比较，在相同的先验规范下。表 2 报告了 MFVB 和 MCMC 包含的变量，按后验包含概率（PIP）排序。如表所示，MFVB 包含了截距和七个额外的变量（即，对应变量的 wj > 0.5），而其余变量的包含概率可以忽略不计。在这八个选定的变量中，有五个是根据 MCMC 确定的五个最有影响力的变量，另外两个变量也出现在 MCMC 选定的变量列表中。因此，MFVB 确定的八个变量中有七个得到了采样器的验证，而 MFVB 仍然更加简约。实际上，MCMC 倾向于包含更多变量，将 PIP > 0.5 分配给 25 个变量，这与第 4 节中报告的模拟结果一致。

表 2. 对于 MFVB 和 MCMC，按后验包含概率（PIP）排序的 LSVT 语音康复应用中包含的特征：MFVB 与 MCMC。**粗体条目表示MFVB估计的PIP(wj)大于0.5的特征。**

| 特征 (MFVB) | PIP (wj) | wjμj | 特征 (MCMC) | PIPE[γjβj] |
|-----------|---------|--------|-----------|
| Intercept | 1.000 | -0.743 | Intercept | 1.000 | -1.280 |
| IMF->NSR_SEO | 1.000 | 0.693 | MFCC_0th coef | 0.742 | -0.429 |
| Shimmer->Ampl_abs | 1.000 | -0.599 | MFCC_1st coef | 0.681 | 0.351 |
| MFCC_0th coef | 1.000 | -0.576 | IMF->NSR_SEO | 0.679 | 0.358 |
| MFCC_1st coef | 1.000 | 0.501 | MFCC_7th coef | 0.670 | 0.294 |
| HNR->HNR_dB_Praat_std | 1.000 | -0.415 | 12th delta | 0.669 | -0.284 |
| MFCC_12th coef | 0.997 | 0.311 | MFCC_4th coef | 0.658 | 0.306 |
| MFCC_7th coef | 0.993 | 0.285 | MFCC_8th coef | 0.645 | 0.273 |
| DFA | 0.641 | 0.314 | IMF->NSR_entropy | 0.632 | 0.298 |
| Log energy | 0.605 | -0.270 | MFCC_9th coef | 0.594 | -0.233 |
| Jitter->pitch_TKEO_prc | 0.578 | 0.227 | entropy_log_5coef | 0.569 | -0.228 |
| HNR->HNR_dB_Praat_std | 0.562 | -0.209 | entropy_log_2coef | 0.549 | -0.204 |
| entropy_log_4coef | 0.546 | -0.206 | IMF->SNR_TKEO | 0.545 | 0.217 |
| GNE->std | 0.520 | -0.181 | entropy_log_7coef | 0.519 | 0.171 |
| Shimmer->Ampl_abs | 0.512 | -0.170 | VFER->NSR_SEO | 0.508 | 0.163 |
| entropy_log3_4coef | 0.508 | -0.160 | IMF->NSR_TKEO | 0.505 | 0.155 |
| Shimmer->Ampl_AM | 0.504 | -0.160 |

**考虑到样本外的预测准确性，我们通过5折交叉验证（使用与优化ρ不同的折叠方式）获得估计值。对于每个折叠，我们通过在剩余四个折叠上训练模型来计算该折叠中观测值的预测概率，并根据方程（7）获得折叠平均偏差。取折叠的平均值，结果得到可比的平均偏差值：MFVB为18.89，MCMC为19.29。因此，使用MFVB识别的配置预计具有可比的样本外性能，但它更为简洁，并且可能更好地突出最重要的变量。此外，MFVB运行时间为0.16秒，而MCMC采样器需要大约544秒，快三个数量级。在分类准确性方面，当新观测值的后验预测概率超过0.5时，预测y^NEW=1，MFVB的准确率为86.5%，MCMC为82.5%，这与Tsanas等人（2013年）报告的更复杂的机器学习技术的水平相当，且具有完全可解释的模型的优势。请注意，上述偏差是基于完整的预测概率的，而误分类率取决于阈值（此处为0.5）。总体而言，所提出的方法提供了一个连贯的、统一的框架，用于联合执行变量选择和预测，结合了可解释的结果和与更先进的黑盒程序相当的预测准确性。**

**5.2. 阿尔茨海默病应用**
第二个真实数据应用关注的是MCMC方法不可行的大规模设置。具体来说，该应用基于R包AppliedPredictiveModeling（Kuhn和Johnson，2018年）中提供的AlzheimerDisease数据集构建。该数据集首次由Craig-Schapiro等人（2011年）分析，涉及阿尔茨海默病的存在与否，分别编码为yi=1或yi=0，共有n=333名患者，作为人口统计数据和实验室结果的函数。根据Fasano等人（2022年）的方法，原始数据集通过添加所有预测变量之间的成对交互作用以及截距进行了修改，从而得到p=9036个预测变量。对于如此高的p值，MCMC是不可行的，如果天真地实现算法1，所提出的MFVB也可能面临计算问题。因此，我们通过设置协方差项Σjk=0（对于j≠k）在方程（6）中评估ηj，从而在γ更新时强制βj之间的独立性，从而实现了一些计算简化。详细信息在补充材料的第S2节中报告。这些简化将MFVB的运行时间缩短到了27秒。先验超参数ν02和ρ的设置与LSVT语音康复分析中的设置相同。

如表3所示，根据通常基于PIP超过0.5的包含标准，选择了五个预测变量以及截距，未选中的预测变量的估计PIP低于0.1。这五个预测变量是脑脊液（CSF）tau（tau）、CSF淀粉样蛋白-β42（Ab_42）、血管内皮生长因子（VEGF）、生长调节α蛋白（GRO_alpha）和胰腺多肽（Pancreatic_polypeptide）。这些结果进一步证实了科学文献中的发现，特别是Craig-Schapiro等人（2011年）的原始研究，其中这五个预测变量被所有考虑的机器学习方法包括在内，唯一的例外是Boosted Tree算法中的GRO_alpha。此外，估计效应μj的符号与原始论文的发现一致。有趣的是，与Craig-Schapiro等人（2011年）考虑的多种方法包含的变量相比，主要差异在于Log Matrix Metalloproteinase-10（MMP-10）预测变量，MFVB没有为其分配相关的PIP，而Craig-Schapiro等人（2011年）探索的方法却将其包括在内。这种差异的一个合理解释是MMP-10与tau之间的强相关性，即0.458（Craig-Schapiro，2011年）。在这里考虑的稀疏贝叶斯框架中，一旦tau被赋予较高的后验包含概率，MMP-10的额外贡献可能会被大幅吸收，使其边际包含概率可以忽略不计。这种行为与spike-and-slab先验倾向于通过丢弃那些预测信号已经被包含的变量来支持简洁表示的倾向一致。在Craig-Schapiro等人（2011年）的ROC分析中可以找到进一步的支持，其中，在相关设置中，将MMP-10添加到log tau/Ab_42比率中仅使曲线下面积（AUC）略有改善，从0.8443增加到0.8518，而所提出方法选择的其他预测变量，如VEGF，提供了相对较大的增益（AUC 0.8766）。关于ν02选择的平均样本外偏差和敏感性分析的进一步结果在补充材料中报告。

**表3. MFVB为阿尔茨海默病应用选择的特征，按后验包含概率排序。**

| 特征 (MFVB) | wj | wjμj | Intercept |
|-----------|------|--------|---------|
| tau | 1.000 | 0.600 | 1.000 |
| Ab_42 | 1.000 | -0.424 | 1.000 |
| VEGF | 1.000 | -0.400 | 0.999 |
| GRO_alpha | 0.999 | 0.340 | 0.997 |
| Pancreatic_polypeptide | 0.997 | 0.325 |

**6. 结论**
在本文中，我们提出了一种用于稀疏二元分类问题的平均场变分贝叶斯（MFVB）方法，该方法采用probit链接。通过使用对系数的spike-and-slab先验来引入稀疏性；另见Ormerod等人（2017年）。如模拟和真实数据应用所示，所提出的CAVI算法可以轻松实现，得益于封闭形式的更新，并且比MCMC程序快几个数量级，从而实现了巨大的计算优势。当观测数量n相对于p足够大时，MFVB和MCMC选择的变量相同，并且在参数估计和预测概率方面都得到可比的解决方案。在p>n的情景中，MFVB倾向于选择更简洁的解决方案，选择的变量较少，而MCMC则更为保守，包括更多的变量，尽管相关的后验包含概率（PIP）较小。进一步的研究方向可以考虑更结构化的近似，如Fasano等人（2022年）对具有高斯先验的贝叶斯probit模型所做的，或者更复杂的模型，包括例如混合效应模型，其中准确和可扩展的近似至关重要（Zhou等人，2024年）。用于重现结果的代码和数据可在GitHub仓库augustofasano/MFVB_Probit_VARIABLE_selection2中找到。

**关于手稿准备过程中生成式AI和AI辅助技术的声明**
在准备这项工作时，作者使用了ChatGPT来校对和润色语言。使用该工具/服务后，作者根据需要审查和编辑了内容，并对发表文章的内容负全责。