青岛大学，山东 青岛

收稿日期：2022年12月28日；录用日期：2023年1月24日；发布日期：2023年1月31日

摘要

尽管Spike and Slab方法广泛应用于贝叶斯变量选择，但其惩罚似然估计的潜力在很大程度上被忽视了。通过在贝叶斯模态中引入惩罚化似然观点，本文提出了新的Spike and Slab Lasso逻辑回归模型，将两个拉普拉斯密度的混合先验置于单个坐标上，可以自适应地收缩系数，即弱收缩重要预测量，强收缩不相关预测量，从而可以得到准确的估计和预测。同时，我们使用了预测–校正算法来求解Spike and Slab Lasso逻辑回归模型，并将该算法扩展到不可分离惩罚的情况。该算法利用凸优化的预测校正算法，沿着整个正则化路径有效的计算解，方便了模型选择，避免了正则化参数值不同时的独立优化。最后，模拟学习和实证结果表明本文所提模型比Lasso逻辑回归模型具有更优的性能。

关键词

Spike and Slab Lasso，逻辑回归，惩罚似然，预测校正算法，贝叶斯先验

The Spike and Slab Lasso Logistic Regression Model Based on Prediction Correction Algorithm

Qi Qi, Qi Zhang^*

Qingdao University, Qingdao Shandong

Received: Dec. 28^th, 2022; accepted: Jan. 24^th, 2023; published: Jan. 31^st, 2023

ABSTRACT

Although the Spike and Slab method is widely used in Bayesian variable selection, its potential for penalized likelihood estimation is ignored. By introducing the penalized likelihood into the Bayesian case, this paper proposes a new Spike and Slab Lasso logistic regression model, which places the mixed priors of two Laplace densities on a single coordinate, and can adaptively adjust the shrinkage coefficient, that is, the weak shrinkage important predictor and the strong shrinkage irrelevant predictor, so that accurate estimation and prediction can be obtained. At the same time, we use the prediction-correction algorithm for the Spike and Slab Lasso logistic regression model, and extend the algorithm to the case of non-separable penalty. The algorithm uses the prediction correction algorithm of convex optimization to effectively calculate the solution along the entire regularization path, which facilitates model selection and avoids independent optimization when the regularization parameter values are different. Finally, the simulation learning and empirical results show that the proposed model has better performance than lasso logistic regression model.

Keywords:Spike and Slab Lasso, Logistic Regression, Penalized Likelihood, Predictor-Corrector Algorithm, Bayesian Prior

Copyright © 2023 by author(s) and Hans Publishers Inc.

This work is licensed under the Creative Commons Attribution International License (CC BY 4.0).

http://creativecommons.org/licenses/by/4.0/

1. 引言

随着现代科技的飞速发展，空前规模的大数据，如基因和蛋白质组学数据，经济学和金融中的面板数据等，在当代的许多应用中遇到。在这些应用中，维数p可以与样本容量n相当，甚至比n大得多，这使得惩罚稀疏模型很适合分析这些数据，自从Tibshirani提出lasso [1] 以来，高维环境中变量选择的惩罚回归已经吸引了相当多的学者 [2] [3]。lasso和它的延伸是变量选择中常用的方法，这些方法对系数施加 $l_1$ 惩罚，惩罚的似然可以通过极快的优化算法来解决，例如LARS和循环坐标下降算法 [4] [5]。然而，lasso对所有系数使用单一惩罚，因此可以包括许多无关的预测因子，或过度收缩大系数。理想的方法是在大效应时诱导弱收缩，在不相关效应时诱导强收缩。几乎所有的回归都容易受到两个问题的影响：过度拟合和多重共线性。前者，参数优化到统计噪声，而不是适合的兴趣关系；后者，协变量组相互关联，致其难以识别个别影响，这些问题在高维协变量集的分析中尤其突出。一种流行的策略是具有先验分布的贝叶斯方法，通过在回归系数上放置某种先验，将系数缩小到零，这提高了模型的稳定性，并防止过拟合。特别地，一种吸引人的选择是对每个回归系数使用拉普拉斯先验分布，该回归系数对应目标函数中的一个 $l_1$ 惩罚，因此具有拉普拉斯先验分布的贝叶斯回归模型的MAP估计等价于Lasso回归的估计。在缺乏有力证据的情况下，拉普拉斯先验分布产生了精确的 ${\beta }_j=0$ 和对应于重要变量的子集的非零 ${\beta }_j$ 的惩罚点估计。通过这种方式，使用 $l_1$ 惩罚自然允许变量选择，其中只有协变量的子集被选择为对结果变量具有实质性预测作用。

也有人提出了变量选择的Spike and Slab方法 [6]，该方法直接产生于概率的考虑，涉及到线性回归的参数和模型空间上的设计先验层次。Gibbs抽样 [7] 用于识别具有高后验概率的有希望的模型，先验的选择通常是棘手的，尽管经验贝叶斯方法可以用来处理这个问题 [8]。最早的关于连续Spike and Slab先验的理论分析之一是由Ishwaran和Rao在2005年进行的 [9]，作者提出并研究了一类连续的双峰先验(两点student-t混合)，建立了后验均值对变量选择和多组分类的拟oracle误分类性能，并创造了“选择性收缩”一词来表示后验均值的渐近行为。到了2011年，Ishwaran和Rao在非正交低维设计的假定下进一步建立了两点高斯混合下后验均值的oracle性质 [10]。在另一个发展中，Narisetty和He在2014年在协变量数目分散的更一般设计中建立了高斯混合先验下贝叶斯因子的模型选择一致性 [11]。后来，到了2017年，V. Rockova和E. I. George在单变量线性回归背景下引入了一种新的稀疏正态均值框架 [12]，弥补了常用的频率主义策略(Lasso)和常用的贝叶斯策略(Spike and Slab)之间的差距，介绍了拉普拉斯先验和点质量的Spike and Slab先验之间的连续统一体——Spike and Slab Lasso (SSL)先验。在高维回归的背景下释放了由连续SSL先验产生的惩罚函数的潜力，超越了独立先验的框架。到了2018年，V. Rockova又在另一篇文章中证明了SSL全局后验模式在平方误差损失下是(接近)最大最小速率最优的 [13]，与Lasso相似。同时还证明了整个后验与全局模式保持同步，并集中(近似)于极大极小率，这是一致的单一拉普拉斯先验不具备的性质，利用适当的一类独立积先验(针对已知的稀疏度)和依赖的混合先验(针对未知的稀疏度)，可以获得极小极大速率的运算性。2020年，Ray Bai等人在SSL的基础上考虑了协变量的组结构 [14]，引入了Spike and Slab Group Lasso (SSGL)用于分组变量线性回归中的贝叶斯估计和变量选择，并进一步将SSGL扩展到稀疏广义可加模型。然而，SSL、SSGL和以往的大多数方法都是基于正态线性模型发展起来的，不能直接应用于其他模型，因此，将混合拉普拉斯先验的高维方法扩展到常规线性模型之外的框架为方法学和应用研究提供了重要的新方向 [15] [16]。

本文章节安排如下：在第2节中，我们介绍了Spike and Slab Lasso (SSL)逻辑回归模型，展示了如何将SSL方法应用于逻辑模型；在第3节中，我们介绍了预测–校正算法，并利用该算法求解Spike and Slab Lasso逻辑回归模型；第4节将逻辑回归模型下的SSL方法扩展到不可分离情况，并相应的将预测–校正算法扩展到不可分离情况；在第5节中，我们使用仿真实验来验证所提出的Spike and Slab Lasso方法，并与常见的Lasso方法进行比较。在第6节中，我们利用了南非心脏病数据来说明SSL逻辑回归模型的使用，建立了南非西开普省地区缺血性心脏病危险因素的强度，通过报告预测效果的均方误差(MSE)，曲线下面积(AUC)和误判率(misclassification)，证实了该方法的可行性与有效性。最后在第7节中对本文进行了总结。

2. Spike and Slab Lasso

2.1. 经典线性Spike and Slab Lasso模型

对于以下线性模型

$Y=X{\beta }_0+\epsilon$ (1)

其中Y为n维响应变量， $X_{n\times p}=\left(X_1,\cdots ,X_p\right)$ 是固定的p个潜在预测的回归矩阵， ${\beta }_0={\left({\beta }_{01},\cdots ,{\beta }_{0P}\right)}^{\prime }$ 是未知的p维回归系数向量， $\epsilon ~N_n\left(0,{\sigma }^2{I}_P\right)$ 为噪声向量。一般的Spike and Slab Lasso先验是这样的形式：

$\pi \left(\beta |\gamma \right)={\displaystyle \underset{i=1}{\overset{P}{\prod }}\left[{\gamma }_i{\phi }_1\left({\beta }_i\right)+\left(1-{\gamma }_i\right){\phi }_0\left({\beta }_i\right)\right]},\gamma ~\pi \left(\gamma \right)$ (2)

其中 $\gamma ={\left({\gamma }_1,\cdots ,{\gamma }_P\right)}^{\prime }$，${\gamma }_i\in \left\{0,1\right\}$ 是二进制向量的中间向量，索引 $2^P$ 个可能的模型。这里， ${\phi }_0\left(\beta \right)$ 作为建模无关(零)系数的Spike分布， ${\phi }_1\left(\beta \right)$ 作为建模大效应的Slab分布。对于Spike and Slab Lasso (SSL)，在上面的特殊变式(2)中引入拉普拉斯密度，即 ${\phi }_0\left(\beta \right)=\frac{{\lambda }_0}{2}{\text{e}}^{-{\lambda }_0\left|\beta \right|}$，${\phi }_1\left(\beta \right)=\frac{{\lambda }_1}{2}{\text{e}}^{-{\lambda }_1\left|\beta \right|}$。则在稀疏正态均值的背景下，这些拉普拉斯分布的两点混合将被称为Spike and Slab Lasso (SSL)先验。模型空间 $\pi \left(\gamma \right)$ 的灵活性极大地扩大了这些先验的范围，对我们来说，我们集中注意力于可交换的模型空间先验形式：

$\pi \left(\gamma |\theta \right)={\displaystyle \underset{j=1}{\overset{P}{\prod }}{\theta }^{{\gamma }_j}{\left(1-\theta \right)}^{1-{\gamma }_j}},\theta ~\pi \left(\theta \right)$ (3)

其中 $\theta =P\left({\gamma }_i=1|\theta \right)$ 为较大的 ${\beta }_{0j}$ 的先验期望分数。在 $\theta$ 条件下，SSL先验(2)可归结为混合物的独立乘积：

$\pi \left(\beta |\theta \right)={\displaystyle \underset{i=1}{\overset{P}{\prod }}\left[\theta {\phi }_1\left({\beta }_i\right)+\left(1-\theta \right){\phi }_0\left({\beta }_i\right)\right]}$ (4)

其中 $\theta \in \left(0,1\right)$ 为混合比例。SSL先验是将两个拉普拉斯密度的混合先验置于单个坐标 ${\beta }_j$ 上。选择一个质点峰值 ${\phi }_0\left({\beta }_j\right)=\infty I\left({\beta }_j=0\right)$，(当 ${\lambda }_0\to \infty$ 时获得)， ${\phi }_1\left({\beta }_j\right)\propto C>0$ (当 ${\lambda }_1\to 0$ 时获得)， $\log \pi \left(\beta |\theta \right)$ 坍缩为 $l_0$ 惩罚；在另一端，选择 ${\phi }_1\left({\beta }_j\right)={\phi }_0\left({\beta }_j\right)$ 产生参数 ${\lambda }_1={\lambda }_0$ 的lasso惩罚。因此，SSL先验的一个特征是它们能够在这两种处理中间形成非凹连续体。

2.2. 正则化逻辑回归

当响应变量为二元时，通常采用线性逻辑回归模型，设数据集为 $D=\left\{\left(x_1,y_1\right),\cdots ,\left(x_i,y_i\right),\cdots ,\left(x_n,y_n\right)\right\}$，其中 $X\in R^{n\times p}$ 是数据矩阵，其中p是特征(参数或属性)的数量 [17]，对于每个 $x_i\in R^p$，y是二进制结果向量，结果是 $y_i=0$ 或 $y_i=1$，逻辑回归模型通过预测因子的线性函数来表示类条件的期望：

$E\left(y_i=1|x_i,\beta \right)=p_i=\frac{{\text{e}}^{{\beta }_0+{x^{\prime }}_i\beta }}{1+{\text{e}}^{{\beta }_0+{x^{\prime }}_i\beta }}=\frac{1}{1+{\text{e}}^{-\left({\beta }_0+{x^{\prime }}_i\beta \right)}},i=1,\cdots ,n$

Logit变换是响应为1的概率的对数，定义为：

${\eta }_i=g\left(p_i\right)=\log \frac{p_i}{1-p_i}={\beta }_0+{x^{\prime }}_i\beta$

在矩阵形式下，logit转换函数表示为：

$\eta =X\beta$

Logit转换函数很重要，因为它是线性的，故具有线性回归模型的许多特性，在逻辑回归中，函数 $g(·)$ 也称为正则链接函数。现在，假设观测值是独立的，则似然函数为：

$L\left(\beta \right)={\displaystyle \underset{i=1}{\overset{n}{\prod }}{\left(p_i\right)}^{y_i}}{\left(1-p_i\right)}^{1-y_i}={{\displaystyle \underset{i=1}{\overset{n}{\prod }}\left(\frac{{\text{e}}^{{\beta }_0+{x^{\prime }}_i\beta }}{1+{\text{e}}^{{\beta }_0+{x^{\prime }}_i\beta }}\right)}}^{y_i}{\left(\frac{1}{1+{\text{e}}^{{\beta }_0+{x^{\prime }}_i\beta }}\right)}^{1-y_i}$ (5)

因此，对数似然函数为

$l\left(\beta \right)={\displaystyle \underset{i=1}{\overset{n}{\sum }}\left(y_i\log \left(\frac{{\text{e}}^{{\beta }_0+{x^{\prime }}_i\beta }}{1+{\text{e}}^{{\beta }_0+{x^{\prime }}_i\beta }}\right)+\left(1-y_i\right)\log \left(\frac{1}{1+{\text{e}}^{{\beta }_0+{x^{\prime }}_i\beta }}\right)\right)}$ (6)

已经有学者证明了逻辑回归的极大似然估计量满足极大似然的期望性质 [18]，但不幸的是，关于 $\beta$ 最大化 $l\left(\beta \right)$ 没有封闭解。因此，逻辑回归的极大似然估计是通过数值优化方法得到的，在本文中，使用了预测–校正算法，该算法在预测步中使用交替方向法得到预测点，然后对得到的部分预测点进行校正，从而保证了算法的收敛性。

2.3. 逻辑回归的可分离SSL惩罚

我们的方案的一个关键因素是利用Spike and Slab Lasso模态估计(变量选择程序的基础)和广义Lasso估计之间的联系，一个重要的第一步是理解可分离的SSL惩罚的结构。假设混合比例 $\theta$ 是固定的，这简化了SSL惩罚的操作，因为它是可根据 $\theta$ 条件分离的，即此时 ${\beta }_j$ 具有条件独立的先验，这导致了可分离的SSL惩罚。

定义1：给定 $\theta \in \left(0,1\right)$，可分离SSL惩罚被定义为：

$pe{n}_s\left(\beta |\theta \right)=\log \left[\frac{\pi \left(\beta |\theta \right)}{\pi \left(0_p|\theta \right)}\right]={\displaystyle \underset{i=1}{\overset{p}{\sum }}\log \left[\frac{\theta {\phi }_1\left({\beta }_j\right)+\left(1-\theta \right){\phi }_0\left({\beta }_j\right)}{\theta {\phi }_1\left(0\right)+\left(1-\theta \right){\phi }_0\left(0\right)}\right]}$ (7)

为了便于惩罚操作，我们将其居中，使 $pe{n}_s\left(0_p|\theta \right)=0$。由于 $\beta$ 与 $\theta$ 的条件无关，惩罚函数由单次函数建立：

$p\left({\beta }_j|\theta \right)=-{\lambda }_1\left|{\beta }_j\right|+\log \left[\frac{p_{\theta }^{\ast }\left(0\right)}{p_{\theta }^{\ast }\left({\beta }_j\right)}\right]$ (8)

$pe{n}_s\left(\beta |\theta \right)={\displaystyle \underset{j=1}{\overset{p}{\sum }}p\left({\beta }_j|\theta \right)}=-{\lambda }_1\left|\beta \right|+{\displaystyle \underset{j=1}{\overset{p}{\sum }}\log \left[\frac{p_{\theta }^{\ast }\left(0\right)}{p_{\theta }^{\ast }\left({\beta }_j\right)}\right]}$ (9)

其中

$p_{\theta }^{\ast }\left({\beta }_j\right)=\frac{\theta {\phi }_1\left({\beta }_j\right)}{\theta {\phi }_1\left({\beta }_j\right)+\left(1-\theta \right){\phi }_0\left({\beta }_j\right)}$ (10)

描述(9)将可分离的SSL惩罚写成Lasso惩罚和非凹惩罚的自适应和，使其最终是非凹的。与在 $\left|{\beta }_j\right|$ 和 $\lambda$ 中都是线性的Lasso惩罚不同，SSL惩罚是 $\left|{\beta }_j\right|$ 和 $\left({\lambda }_0,{\lambda }_1,\theta \right)$ 的非线性函数，尽管有明显的差异，但两个惩罚之间具有一定的联系，求导后，这个联系就暴露出来了，导数对应于一个隐偏差项，在估计中起着至关重要的作用，因此我们有以下引理：

引理1：可分离的SSL惩罚满足：

$\frac{\partial pe{n}_s\left(\beta |\theta \right)}{\partial \left|{\beta }_j\right|}\equiv -{\lambda }_{\theta }^{\ast }\left({\beta }_j\right)$ (11)

其中，

${\lambda }_{\theta }^{\ast }\left({\beta }_j\right)={\lambda }_1{p}_{\theta }^{\ast }\left({\beta }_j\right)+{\lambda }_0\left[1-p_{\theta }^{\ast }\left({\beta }_j\right)\right]$ (12)

另一方面，拉普拉斯先验求导为： $\frac{\partial }{\partial \left|{\beta }_j\right|}\log \phi \left({\beta }_j|\lambda \right)=-\lambda$。因此SSL偏差项是两个Lasso偏差项的凸组合，重要的是，这种组合是有适应性的，因为 $p_{\theta }^{\ast }\left({\beta }_j\right)$ 依赖于 ${\beta }_j$，这是Spike and Slab惩罚的独特特征，它对每个系数分别加权Spike和Slab的贡献。与之形成鲜明对比的是，无论系数大小如何，Lasso惩罚对每个系数的偏差都是相同的，这通常是收缩和偏差之间冲突的根源， ${\lambda }_{\theta }^{\ast }\left({\beta }_j\right)$ 的额外灵活性极大地缓解了这种冲突。

${\lambda }_{\theta }^{\ast }\left({\beta }_j\right)$ 是一个自适应的线性组合，由 $p_{\theta }^{\ast }\left({\beta }_j\right)$ 加权，影响较大的系数 $p_{\theta }^{\ast }\left({\beta }_j\right)$ 接近1并且收缩的更小，这是因为 ${\lambda }_{\theta }^{\ast }\left({\beta }_j\right)$ 主要由 ${\lambda }_1$ 驱动，为了避免过度收缩， ${\lambda }_1$ 被设置的很小。影响较小的系数则相反，它们具有较小的包含概率，因此 ${\lambda }_{\theta }^{\ast }\left({\beta }_j\right)$ 被较大的惩罚 ${\lambda }_0$ 所接管，也就是说，当 $\left|{\beta }_j\right|$ 较大时，倾向于收缩少量，当 $\left|{\beta }_j\right|$ 较小时，倾向于收缩大量，这是SSL估计器背后选择性收缩特性的一种表现。混合比例 $p_{\theta }^{\ast }\left({\beta }_j\right)$ 可以看做条件包含概率 $p\left({\gamma }_i=1|{\beta }_i,\theta \right)$，从而有：

$p_{\theta }^{\ast }\left({\beta }_j\right)=p\left({\gamma }_i=1|{\beta }_i,\theta \right)=\frac{1}{1+\frac{1-\theta }{\theta }\frac{{\lambda }_0}{{\lambda }_1}\exp \left[-\left|{\beta }_i\right|\left({\lambda }_0-{\lambda }_1\right)\right]}$

3. 预测校正算法

预测校正算法是实现数值延拓的基本策略之一 [19]。在许多方法中，预测校正算法通过实用初始条件(参数一个极值的解)显式的求出一系列解，并根据当前解继续求出相邻解，我们详细说明了如何使用预测校正算法跟踪曲线 $H\left(\beta ,{\lambda }_{\theta }^{\ast }\right)=0$ 到我们的问题设置。为了找到 $\beta ={\left({\beta }_0,{\beta }^{\prime }\right)}^{\prime }$，我们的准则使用 ${\lambda }_{\theta }^{\ast }$ 从一个固定的值开始减少，这等价于最小化以下问题：

$l\left(\beta ,{\lambda }_{\theta }^{\ast }\right)=\left\{-{\displaystyle \underset{i=1}{\overset{n}{\sum }}\left[y_i\log \left(p_i\right)+\left(1-y_i\right)\log \left(1-p_i\right)\right]}-pe{n}_s\left(\beta |\theta \right)\right\}$ (13)

假设 $\beta$ 的分量全都不为0，且 $l\left(\beta ,{\lambda }_{\theta }^{\ast }\right)$ 关于 $\beta$ 可导，定义函数H为：

$H\left(\beta ,{\lambda }_{\theta }^{\ast }\right)=\frac{\partial l}{\partial \beta }=-X^{\text{T}}\left(y-p\right)+{\lambda }_{\theta }^{\ast }\left(\beta \right)Sgn\left(\begin{array}{c}0\\ \beta \end{array}\right)$ (14)

其中假设 $\beta$ 的所有分量都不为零，X是 $n\times \left(p+1\right)$ 阶矩阵，包括列为1。即使已经假设了 $\beta$ 的所有分量都不为零，但 $\beta$ 的非零分量集随 ${\lambda }_{\theta }^{\ast }$ 的变化而变化， ${\lambda }_{\theta }^{\ast }$ 中， ${\lambda }_1\in \left(0,\infty \right)$，${\lambda }_2$ 是很小的，固定的正常数， $\theta \in \left(0,1\right)$ 是给定的。令 $\hat{\beta }$ 表示可分离SSL惩罚 $pe{n}_s\left(\beta |\theta \right)$ 下的全局后验模式，即

$\hat{\beta }\left({\lambda }_{\theta }^{\ast }\right)=\underset{\beta \in R^p}{\arg \min }\left\{-{\displaystyle \underset{i=1}{\overset{n}{\sum }}\left[y_i\log p_i+\left(1-y_i\right)\log \left(1-p_i\right)\right]}-pe{n}_s\left(\beta |\theta \right)\right\}$ (15)

引理2： $\hat{\beta }={\left({\hat{\beta }}_1,\cdots ,{\hat{\beta }}_p\right)}^{\prime }$ 为 $pe{n}_s\left(\beta |\theta \right)$ 下的全局后验模式(15)的KKT条件为：

${X^{\prime }}_j\left(y-\hat{p}\right)={\lambda }_{\theta }^{\ast }\left({\hat{\beta }}_j\right)sign\left({\hat{\beta }}_j\right)$ for ${\hat{\beta }}_j\ne 0$ $j=1,\cdots ,p$ (16)

$\left|{X^{\prime }}_j\left(y-\hat{p}\right)\right|\le {\lambda }_{\theta }^{\ast }\left({\hat{\beta }}_j\right)$ for ${\hat{\beta }}_j=0$ $j=1,\cdots ,p$ (17)

对于任意 $j=1,\cdots ,p$，当 ${\hat{\beta }}_j=0$ 时，KKT条件暗示了 $1^{\prime }\left(y-\hat{p}\right)=0$，即

$\hat{p}=\bar{y}1=g^{-1}\left({\hat{\beta }}_0\right)1$ (18)

引理3：当 ${\lambda }_{\theta }^{\ast }$ 超过某个阈值时，截距是唯一的非零系数，即

${\tilde{\beta }}_0=g\left(\bar{y}\right)$，且 $H\left({\left({\hat{\beta }}_0,0,\cdots ,0\right)}^{\prime },{\lambda }_{\theta }^{\ast }\right)=0$ 对于 ${\lambda }_{\theta }^{\ast }>\underset{j\in \left\{1,\cdots ,p\right\}}{\max }\left|{x^{\prime }}_j\left(y-\bar{y}1\right)\right|$。

我们把 ${\lambda }_{\theta }^{\ast }$ 的这个阈值表示为 ${\lambda }_{\theta ,\max }^{\ast }$，从变量 $j_0=\underset{j}{\arg \max }\left|{x^{\prime }}_j\left(y-\bar{y}1\right)\right|$ 开始，随着 ${\lambda }_{\theta }^{\ast }$ 的进一步减少，其他变量加入活动集。从 ${\lambda }_{\theta ,\max }^{\ast }$ 开始减少 ${\lambda }_{\theta }^{\ast }$，我们在预测步和校正步之间交替，第k次迭代的步骤如下：

1) 步长：确定下降量 ${\lambda }_{\theta }^{\ast }$，给定 ${\lambda }_{\theta ,k}^{\ast }$，我们近似第二大 ${\lambda }_{\theta }^{\ast }$，其中活动集发生变化，命名为 ${\lambda }_{\theta ,k+1}^{\ast }$。

2) 预测步：用 ${\lambda }_{\theta }^{\ast }$ 的下降量线性近似 $\beta$ 的相应改变，称为 ${\hat{\beta }}^{k+}$。

3) 校正步：找到 $\beta$ 中和 ${\lambda }_{\theta ,k+1}^{\ast }$ 配对的精确解(即 $\beta \left({\lambda }_{\theta ,k+1}^{\ast }\right)$ )，用作为初始值，称为 ${\hat{\beta }}^{k+1}$。

4) 活动集：测试以查看当前活动集是否必须是改进的，如果是，用更新后的活动集重复校正步。

3.1. 预测步

定义 $f\left({\lambda }_{\theta }^{\ast }\right)=H\left(\beta \left({\lambda }_{\theta }^{\ast }\right),{\lambda }_{\theta }^{\ast }\right)$，则在第k个预测步中，在生成的当前活动集的范围内， $f\left({\lambda }_{\theta }^{\ast }\right)$ 对于所有的 ${\lambda }_{\theta }^{\ast }$ 都是0，通过对 $f\left({\lambda }_{\theta }^{\ast }\right)$ 关于 ${\lambda }_{\theta }^{\ast }$ 求导有 $f^{\prime }\left({\lambda }_{\theta }^{\ast }\right)=\frac{\partial H}{\partial {\lambda }_{\theta }^{\ast }}+\frac{\partial H}{\partial \beta }\frac{\partial \beta }{\partial {\lambda }_{\theta }^{\ast }}=0$，从上式我们计算 $\frac{\partial \beta }{\partial {\lambda }_{\theta }^{\ast }}$，有

$\frac{\partial \beta }{\partial {\lambda }_{\theta }^{\ast }}=-{\left(\frac{\partial H}{\partial \beta }\right)}^{-1}\frac{\partial H}{\partial {\lambda }_{\theta }^{\ast }}=-{\left({X^{\prime }}_A{W}_k{X}_A\right)}^{-1}Sgn\left(\begin{array}{c}0\\ {\hat{\beta }}^k\end{array}\right)$ (19)

其中W是n阶对角阵，元素为 $p_i\left(1-p_i\right),i=1,\cdots ,n$，$W_k$ 和 $X_A$ 表示当前活动集中X的列。从而在第k个预测步中， $\beta \left({\lambda }_{\theta ,k+1}^{\ast }\right)$ 被近似通过

$\begin{array}{c}{\hat{\beta }}^{k+}={\hat{\beta }}^k+\left({\lambda }_{\theta ,k+1}^{\ast }-{\lambda }_{\theta ,k}^{\ast }\right)\frac{\partial \beta }{\partial {\lambda }_{\theta }^{\ast }}\\ ={\hat{\beta }}^k-\left({\lambda }_{\theta ,k+1}^{\ast }-{\lambda }_{\theta ,k}^{\ast }\right){\left({X^{\prime }}_A{W}_k{X}_A\right)}^{-1}Sgn\left(\begin{array}{c}0\\ {\hat{\beta }}^k\end{array}\right)\end{array}$ (20)

上式中 $\beta$ 仅由当前非零系数组成，这种线性化相当于对对数似然进行二次逼近，并通过加权Lasso步骤来扩展当前解 ${\hat{\beta }}^k$。下列定理说明，预测步近似可以通过使 $\left({\lambda }_{\theta ,k}^{\ast }-{\lambda }_{\theta ,k\text{+}1}^{\ast }\right)$ 小来任意接近真实解。

定理1：令 $h_{\theta ,k}^{\ast }={\lambda }_{\theta ,k}^{\ast }-{\lambda }_{\theta ,k\text{+}1}^{\ast }$，假设 $h_{\theta ,k}^{\ast }$ 足够小，活动集在 ${\lambda }_{\theta }^{\ast }={\lambda }_{\theta ,k}^{\ast }$ 和 ${\lambda }_{\theta }^{\ast }={\lambda }_{\theta ,k+1}^{\ast }$ 是相同的，则近似解 ${\hat{\beta }}^{k+}$ 和真实解 ${\hat{\beta }}^{k+1}$ 有 $O\left(h_{\theta ,k}^2\right)$ 的区别。

证明：因为 $\frac{\partial \beta }{\partial {\lambda }_{\theta }^{\ast }}=-{\left({X^{\prime }}_A{W}_k{X}_A\right)}^{-1}Sgn\left(\begin{array}{c}0\\ {\hat{\beta }}^k\end{array}\right)$ 是连续可微的对于 ${\lambda }_{\theta }^{\ast }\in \left({\lambda }_{\theta ,k+1}^{\ast },{\lambda }_{\theta ,k}^{\ast }\right]$

${\hat{\beta }}^{k+1}={\hat{\beta }}^k-h_{\theta ,k}^{\ast }{\frac{\partial \beta }{\partial {\lambda }_{\theta }^{\ast }}|}_{{\lambda }_{\theta ,k}^{\ast }}+O\left(h_{\theta ,k}^2\right)={\hat{\beta }}^{k+}+O\left(h_{\theta ,k}^2\right)$

3.2. 校正步

在第k次校正步中，我们使用先前的预测步中的近似作为初始值来计算精确解 $\beta \left({\lambda }_{\theta ,k+1}^{\ast }\right)$。

定理2：如果将在 ${\lambda }_{\theta ,k}^{\ast }$ 和 ${\lambda }_{\theta ,k+1}^{\ast }={\lambda }_{\theta ,k}^{\ast }-h_{\theta ,k}^{\ast }$ 处的解记为 ${\hat{\beta }}^k$ 和 ${\hat{\beta }}^{k+1}$，他们被连接，使得对某个 $\alpha \in \left[0,1\right]$，在 ${\lambda }_{\theta }^{\ast }={\lambda }_{\theta ,k}^{\ast }-\alpha h_{\theta ,k}^{\ast }$ 处，我们的估计为：

$\hat{\beta }\left({\lambda }_{\theta }^{\ast }-\alpha h_{\theta ,k}^{\ast }\right)={\hat{\beta }}^k+\alpha \left({\hat{\beta }}^{k+1}-{\hat{\beta }}^k\right)$ (21)

则 $\hat{\beta }\left({\lambda }_{\theta }^{\ast }-\alpha h_{\theta ,k}^{\ast }\right)$ 和真实解 $\beta \left({\lambda }_{\theta }^{\ast }-\alpha h_{\theta ,k}^{\ast }\right)$ 有 $O\left(h_{\theta ,k}^{\ast }\right)$ 的区别。

证明：因为 $\frac{\partial \beta }{\partial {\lambda }_{\theta }^{\ast }}$ 是连续可微的对于 ${\lambda }_{\theta }^{\ast }\in \left({\lambda }_{\theta ,k+1}^{\ast },{\lambda }_{\theta ,k}^{\ast }\right]$，则下式成立：

$\hat{\beta }\left({\lambda }_{\theta }^{\ast }-\alpha h_{\theta ,k}^{\ast }\right)={\hat{\beta }}^k-\alpha h_{\theta ,k}^{\ast }\left(\frac{{\hat{\beta }}^{k+1}-{\hat{\beta }}^k}{-h_{\theta ,k}^{\ast }}\right)={\hat{\beta }}^k-\alpha h_{\theta ,k}^{\ast }{\frac{\partial \beta }{\partial {\lambda }_{\theta }^{\ast }}|}_{{\lambda }_{\theta ，k}^{\ast }}+O\left(h_{\theta ,k}^2\right)$

类似的，在 ${\lambda }_{\theta }^{\ast }={\lambda }_{\theta ,k}^{\ast }-\alpha h_{\theta ,k}^{\ast }$ 处的真实解为：

$\beta \left({\lambda }_{\theta }^{\ast }-\alpha h_{\theta ,k}^{\ast }\right)={\hat{\beta }}^k-\alpha h_{\theta ,k}^{\ast }{\frac{\partial \beta }{\partial {\lambda }_{\theta }^{\ast }}|}_{{\lambda }_{\theta ,k}^{\ast }}+O\left(h_{\theta ,k}^2\right)$

3.3. 活动集

活动集A从引理3中的截距开始，在每个校正步之后，应用下面的检查程序去检查A是否应该被扩充：

$\left|X_j^{\text{T}}\left(y-p\right)\right|>{\lambda }_{\theta }^{\ast }\left(\beta \right)$，对于任意 $j\in A$ $\Rightarrow A\leftarrow A\cup \left\{j\right\}$

我们用改进后的活动集重复校正步直到活动集不再扩充。然后我们从活动集中删除带有零系数的变量，即 $\left|{\hat{\beta }}_j\right|=0$ 对于任意 $j\in A$ $\Rightarrow A\leftarrow A\backslash \left\{j\right\}$。

KKT最优条件(16)~(17)暗示了 ${\hat{\beta }}_j\ne 0\Rightarrow \left|{X^{\prime }}_j\left(Y-\hat{P}\right)\right|={\lambda }_{\theta }^{\ast }\left({\hat{\beta }}_j\right)$，结合(16)式和(18)式，我们有以下引理：

引理4：令 $\hat{p}$ 是y的来自一个校正步的估计， $\hat{c}$ 表示因子和当前残差相关权重：

$\hat{c}=X^{\prime }\left(y-\hat{p}\right)$ (22)

A中各因子(除截距外)的绝对相关系数为 ${\lambda }_{\theta }^{\ast }$，而 $A^c$ 中各因子的值小于 ${\lambda }_{\theta }^{\ast }$。也就是说，在每个校正步之后，如果对于任意 $l\in A^c$ 有 $\left|{\hat{c}}_l\right|>{\lambda }_{\theta }^{\ast }$，我们通过增加 $x_l$ 来扩大活动集。直到活动集不再进一步扩大，校正步停止重复。如果对于任意 $l\in A$，${\hat{\beta }}_l=0$，则我们从活动集中估计 $x_l$。

3.4. 步长

如果 ${\lambda }_{\theta }^{\ast }=0$，算法停止，如果 ${\lambda }_{\theta }^{\ast }>0$，我们在 ${\lambda }_{\theta }^{\ast }$ 中近似最小的减量，活动集将会被改进。随着 ${\lambda }_{\theta }^{\ast }$ 减少 $h_{\theta }^{\ast }$，变化量表示为a，在当前校正步的近似变化量如下：

$\begin{array}{c}c\left(h_{\theta }^{\ast }\right)=\hat{c}-h_{\theta }^{\ast }\left(a\right)\\ =\hat{c}-h_{\theta }^{\ast }{X}^{\prime }\hat{W}{X}_A{\left({X^{\prime }}_A\hat{W}{X}_A\right)}^{-1}Sgn\left(\begin{array}{c}0\\ \hat{\beta }\end{array}\right)\end{array}$ (23)

其中 $h_{\theta }^{\ast }>0$ 是 ${\lambda }_{\theta }^{\ast }$ 中给定的减少量，对于A中的因子，a的值为 $Sgn\left(\begin{array}{c}0\\ \hat{\beta }\end{array}\right)$ 的值，求 $A^c$ 的任意因子与A中任意因子的绝对相关系数相同的 $h_{\theta }^{\ast }$，我们解以下方程：

$\left|c_j\left(h_{\theta }^{\ast }\right)\right|=\left|{\hat{c}}_j-h_{\theta }^{\ast }{a}_j\right|={\lambda }_{\theta }^{\ast }-h_{\theta }^{\ast },\forall j\in A^c$ (24)

方程给出了 ${\lambda }_{\theta }^{\ast }$ 中的步长估计为

$h_{\theta }^{\ast }=\underset{j\in A^c}{{\min }^{+}}\left\{\frac{{\lambda }_{\theta }^{\ast }-{\hat{c}}_j}{1-a_j},\frac{{\lambda }_{\theta }^{\ast }+{\hat{c}}_j}{1+a_j}\right\}$ (25)

另外，检查活动集中是否有变量在 ${\lambda }_{\theta }^{\ast }$ 减少 $h_{\theta }^{\ast }$ 之前达到0，我们解出方程：

${\beta }_j\left({\tilde{h}}_{\theta }^{\ast }\right)={\hat{\beta }}_j+{\tilde{h}}_{\theta }^{\ast }{\left({X^{\prime }}_A\hat{W}{X}_A\right)}^{-1}Sgn\left(\begin{array}{c}0\\ \hat{\beta }\end{array}\right)=0,\forall j\in A$ (26)

若 $\forall j\in A$，$0<{\tilde{h}}_{\theta }^{\ast }<h_{\theta }^{\ast }$，我们期望相应的变量在任何其他变量加入活动集之前被消去，因此 ${\tilde{h}}_{\theta }^{\ast }$ 而不是 $h_{\theta }^{\ast }$ 作为下一个步长。

4. 不可分离情况

如前所述，惩罚由参数 $\left({\lambda }_0,{\lambda }_1,\theta \right)$ 的三元组索引，它们串联工作以产生理想的性质。在整篇论文中，我们假设 ${\lambda }_1$ 被设为一个很小的值，因此不需要调谐，两个参数 $\left({\lambda }_0,\theta \right)$ 将驱动惩罚的性能，它们的调谐将是最重要的。在可分离情况中，我们假设 $\theta \in \left(0,1\right)$ 是给定的，然而 $\theta$ 控制着大系数的期望比例，在缺乏真稀疏度水平q的先验信息的情况下，任意预先指定 $\theta$ 可能会无意识的高估\低估真稀疏度分数q/p，从而影响性能。

本节采用贝叶斯策略，赋予 $\theta$ 一个合适的先验 $\theta ~\pi \left(\theta \right)$，使惩罚可以在不需要设置 $\theta$ 接近q/p的情况下，实现一定的自适应和升压性能。假设 $\pi \left(\theta \right)$ 是一个一般的先验， $\beta$ 中的坐标是根据 $\pi \left(\theta \right)$ 的边际相依分布的。

$\begin{array}{c}\pi \left(\beta \right)={\displaystyle {\int }_0^1{\displaystyle \underset{j=1}{\overset{p}{\prod }}\left[\theta {\phi }_1\left({\beta }_j\right)+\left(1-\theta \right){\phi }_0\left({\beta }_j\right)\right]\pi \left(\theta \right)\text{d}\theta }}\\ ={\left(\frac{{\lambda }_1}{2}\right)}^p{\text{e}}^{-{\lambda }_1\left|\beta \right|}{\displaystyle {\int }_0^1\frac{{\theta }^p}{{\displaystyle {\prod }_{j=1}^p{p}_{\theta }^{*}\left({\beta }_j\right)}}\pi \left(\theta \right)\text{d}\theta }\end{array}$ (27)

重写(27)作为惩罚函数，我们得到以下不可分离的SSL惩罚变量。

4.1. 不可分离SSL惩罚

定义2：不可分离的Spike and Slab Lasso (NSSL)惩罚在 $\theta ~\pi \left(\theta \right)$ 下定义为：

$pe{n}_{NS}\left(\beta \right)=\log \left[\frac{\pi \left(\beta \right)}{\pi \left(0_p\right)}\right]=-{\lambda }_1\left|\beta \right|+\log \left[\frac{{\displaystyle \int \frac{{\theta }^p}{{\displaystyle \underset{j=1}{\overset{p}{\prod }}{p}_{\theta }^{*}\left({\beta }_j\right)}}\pi \left(\theta \right)\text{d}\theta }}{{\displaystyle \int \frac{{\theta }^p}{{\displaystyle \underset{j=1}{\overset{p}{\prod }}{p}_{\theta }^{*}\left(0\right)}}\pi \left(\theta \right)\text{d}\theta }}\right]$ (28)

再一次，我们把惩罚中心化，使 $pe{n}_{NS}\left(0\right)=0$。将(28)与(11)相比，NSSL惩罚仍然是(可分离的) Lasso部分和非凹部分的相加组成，但现在非凹部分将是不可分的。一般来说，(28)中的积分没有一个封闭形式的解，这似乎使惩罚的可处理型变得复杂化了，然而，当意识到先验(内含偏差项)的得分函数可以写成简单且非常直观的形式后，操作就变得非常简单了，这种形式出现在下面引理5的不可分离类比中。令 ${\beta }_{\backslash j}$ 表示去掉第j个分量之外的其余分量所组成的 $\beta$ 的子向量。

引理5：不可分离的Spike and Slab Lasso惩罚(28)的导数满足：

$\frac{\partial pe{n}_{NS}\left(\beta \right)}{\partial \left|{\beta }_j\right|}\equiv -{\lambda }^{*}\left({\beta }_j;{\beta }_{\backslash j}\right)$ (29)

其中，

${\lambda }^{*}\left({\beta }_j;{\beta }_{\backslash j}\right)=p^{*}\left({\beta }_j;{\beta }_{\backslash j}\right){\lambda }_1+\left[1-p^{*}\left({\beta }_j;{\beta }_{\backslash j}\right)\right]{\lambda }_0$ (30)

$p^{*}\left({\beta }_j;{\beta }_{\backslash j}\right)\equiv {\displaystyle {\int }_0^1{p}_{\theta }^{*}\left({\beta }_j\right)\pi \left(\theta |\beta \right)\text{d}\theta }$ (31)

我们现在稍微停顿一下，来看(30)和(12)的区别，不可分离惩罚的混合概率为 $p_{\theta }^{\ast }(·)$ 在 $\pi \left(\theta |\beta \right)$ 上平均得到的聚合概率 $p^{*}\left({\beta }_j;{\beta }_{\backslash j}\right)$。正是通过这种平均，惩罚才有机会了解 $\beta$ 的稀疏水平。对不可分离惩罚的初步认识表明，它的自适应机制通过条件分布在概率域内运行，这一点在可分离惩罚中完全缺失了。由于 $p_{\theta }^{*}\left({\beta }_j\right)$ 是 $\theta$ 的非线性函数，因此在(17)中撇去 $\theta$ 对平均混合权重 $p_{\theta }^{*}\left({\beta }_j;{\beta }_{\backslash j}\right)$ 的影响尚不明显。这个谜题在以下引理中展开，它为NSSL惩罚的实施和理论研究提供了简化。

引理6：给定 $\beta \in {ℝ}^p$ 以及先验 $\pi \left(\theta \right)$，则有：

$p^{*}\left({\beta }_j;{\beta }_{\backslash j}\right)=p_{{\theta }_j}^{*}\left({\beta }_j\right)$ (32)

其中，

${\theta }_j=E\left[\theta |{\beta }_{\backslash j}\right]$ (33)

上述引理的意义在于，通过简单的代入 $\theta =E\left[\theta |{\beta }_{\backslash j}\right]$，我们可以将对可分离情形的认识转化为不可分离情形的认识，上述两个式子也说明了完全贝叶斯公式如何将概率意义赋予NSSL惩罚元素。

4.2. 自适应权重

我们现在考虑条件均值 $E\left[\theta |{\hat{\beta }}_{\backslash j}\right]$ 的形式，当p充分大时，后验期望 $E\left[\theta |{\hat{\beta }}_{\backslash j}\right]$ 将非常相似且接近于 $E\left[\theta |\hat{\beta }\right]$。对于 $\theta$ 的先验，我们将使用标准beta先验 $\theta ~B\left(a,b\right)$。现在我们检验条件分布 $\pi \left(\theta |\hat{\beta }\right)$，这种条件分布将受到非零系数 $\hat{q}={\Vert \beta \Vert }_0$ 的数量及其大小的影响。假设 $\hat{\beta }$ 中的前 $\hat{q}$ 项是非零的，这个分布的密度为：

$\pi \left(\theta |\hat{\beta }\right)\propto {\theta }^{a-1}{\left(1-\theta \right)}^{b-1}{\left(1-\theta z\right)}^{p-\hat{q}}{\displaystyle \underset{j=1}{\overset{\hat{q}}{\prod }}\left(1-\theta x_j\right)}$ (34)

其中， $z=1-\frac{{\lambda }_1}{{\lambda }_0}$，$x_j=\left(1-\frac{{\lambda }_1}{{\lambda }_0}{\text{e}}^{\left|{\hat{\beta }}_j\right|\left({\lambda }_0-{\lambda }_1\right)}\right)$。这个分布是高斯超几何分布的一种推广，归一化常数写为多变量超几何函数的欧拉积分表示，因此它的期望可以写成：

$E\left[\theta |\hat{\beta }\right]=\frac{{\displaystyle {\int }_0^1{\theta }^a{\left(1-\theta \right)}^{b-1}{\left(1-\theta z\right)}^{p-\hat{q}}{\displaystyle \underset{j=1}{\overset{\hat{q}}{\prod }}\left(1-\theta x_j\right)\text{d}\theta }}}{{\displaystyle {\int }_0^1{\theta }^{a-1}{\left(1-\theta \right)}^{b-1}{\left(1-\theta z\right)}^{p-\hat{q}}{\displaystyle \underset{j=1}{\overset{\hat{q}}{\prod }}\left(1-\theta x_j\right)\text{d}\theta }}}$ (35)

虽然上面的表达式(35)似乎很难计算，但当 ${\lambda }_0$ 非常大时，它可以有一个简单的多的形式，根据 [12]，我们在引理7中获得了这个简单得多的形式：

引理7：假设 $\pi \left(\theta |\hat{\beta }\right)$ 是根据 $\pi \left(\theta |\hat{\beta }\right)\propto {\theta }^{a-1}{\left(1-\theta \right)}^{b-1}{\left(1-\theta z\right)}^{p-\hat{q}}{\displaystyle \underset{j=1}{\overset{\hat{q}}{\prod }}\left(1-\theta x_j\right)}$ 来分布的。令 $\hat{q}={\Vert \hat{\beta }\Vert }_0$，当 ${\lambda }_0\to \infty$ 时，

$E\left[\theta |\hat{\beta }\right]=\frac{a+\hat{q}}{a+b+G}$ (36)

我们注意到这个表达式本质上是在beta先验下的 $\theta$ 的通常的后验均值。直观的说，当 ${\lambda }_0$ 发散时，权重 $p_{\theta }^{\ast }\left({\beta }_j\right)$ 集中在0和1处，得到熟悉的 $E\left(\theta |\hat{\beta }\right)$ 的形式。

4.3. 不可分离情况的预测校正算法

在本节中，我们将预测–校正算法扩展到不可分离惩罚的情况。具有不可分离惩罚的逻辑回归模型通过最小化下列损失函数来寻找系数：

$l\left(\beta ,{\lambda }_{{\theta }_j}^{*}\right)=\left\{-{\displaystyle \underset{i=1}{\overset{n}{\sum }}\left[y_i\log \left(p_i\right)+\left(1-y_i\right)\log \left(1-p_i\right)\right]}-pe{n}_{NS}\left(\beta \right)\right\}$ (37)

$l\left(\beta ,{\lambda }_{{\theta }_j}^{*}\right)$ 关于 $\beta$ 的一阶导和二阶导函数分别如下所示：

$H\left(\beta ,{\lambda }_{{\theta }_j}^{*}\right)=\frac{\partial l}{\partial \beta }=-X^{\text{T}}\left(y-p\right)+{\lambda }_{{\theta }_j}^{*}\left(\beta \right)Sgn\left(\begin{array}{c}0\\ \beta \end{array}\right)$ (38)

$\frac{\partial H}{\partial \beta }=\frac{{\partial }^{2}l}{\partial \beta \partial {\beta }^{\prime }}=X^{\prime }WX$ (39)

其中W是n阶对角阵，元素为 $W_{ii}=p_i\left(1-p_i\right),i=1,\cdots ,n$。在贝叶斯框架中，通过将 $\theta$ 作为一个附加的模型参数，可将预测–校正算法扩展到不可分离惩罚的情况，现在只需要通过算法来更新 $\theta$，而不是使用 $\theta$ 的固定值。因此，在第k次迭代中，使用 $\theta ={\theta }^{\left(k\right)}$ 来更新 ${\beta }^{\left(k+1\right)}$，其中是 $\theta$ 的第k次迭代值，是在引理7中得到的。

5. 模拟学习

在本节中，我们使用仿真实验来验证所提出的Spike and Slab Lasso方法，并与常见的Lasso方法进行比较。我们模拟了两个数据集，用第一个数据集作为训练数据来拟合模型，第二个数据集作为测试数据来评估预测值，对于每个模拟设置，我们重复模拟50次，并总结这些重复的结果。

我们使用所提出的SSL逻辑回归方法分析了每个模拟数据集，并与Lasso逻辑回归方法进行了对比，以说明我们的方法相比Lasso方法的优越性。对于Lasso方法，我们通过十折交叉验证选择λ的最优值，λ决定了最优的Lasso，我们还报告了最优的Lasso模型的结果(见图1)。

对于每个数据集，令 $n=100$，$p=1000$，我们从多元高斯分布中生成一个数据矩阵X，均值为 $0_p$，$\Sigma ={\left({\sigma }_{ij}\right)}_{i,j=1}^p$，其中当 $i\ne j$ 时， ${\sigma }_{ij}=0.6$，否则， ${\sigma }_{ii}=1$。真向量 ${\beta }_0$ 通过指定回归系数 $\frac{1}{\sqrt{3}}\left\{-2.5,-2,-1.5,-1,1,1.5,2,2.5\right\}$ 到 $q=8$ 个随机方向，其余系数设为0，响应由 $N\left({\eta }_i,{1.6}^2\right)$ 生成，其中 ${\eta }_i={\beta }_0+{\displaystyle {\sum }_{j=1}^m{x}_{ij}{\beta }_j}$。我们考虑以下三种设置，其均方误差(MSE)图分别见图2，图3，图4：

1) 非自适应选择 $\theta =0.5$，明显高估了真稀疏分数8/1000；

2) 非自适应oracle选择 $\theta =8/1000$ ；

3) 自适应选择 $\theta ~B\left(1,1000\right)$ ；

图1. Lasso系数的估计值。图中的垂直线对应最佳Lasso模型(交叉验证)

图2. 设置1情况下的MSE

图3. 设置2情况下的MSE

图4. 设置3情况下的MSE

表1. Lambda0的不同取值下，三种设置和Lasso情况的MSE对比

从表1可以看出，在n = 100，p = 1000的设置下，设置3是最优的，其次是设置2，因为设置2是系数的真实稀疏度。说明当我们设置权重θ自适应于数据时，我们的模型能更好的适应数据的稀疏性。上述表中结果也表明SSL逻辑回归方法与Lasso进行比较时表现良好。

6. 南非心脏病数据集

在这里，我们提出了一个二元数据的分析，以说明SSL逻辑回归模型的使用效果。图5中的数据是冠状动脉危险因素研究(CORIS)基线调查的一个子集，该调查在南非西开普省的三个农村地区进行。该研究的目的是建立该高发地区缺血性心脏病危险因素的强度。数据代表15~64岁的白人男性，响应变量为调查时是否存在心肌梗死(MI)，(该地区MI的总体患病率为5.1%)。在我们的数据集中有160个病例以及302个对照样本。

图5是南非心脏病数据的散点图矩阵。每个图显示一对风险因素，病例和对照用颜色编码(红色表示病例)。变量心脏病的家族史(famhist)是二元的(是或否)。使用疾病/非疾病变量，在不可分离惩罚的情况下，我们可以拟合一个SSL逻辑回归模型。系数的精确路径图在图6到图8中给出。

图5. 南非心脏病数据的散点图矩阵

图6. 系数的精确路径集，其中x轴为系数的l₁范数，垂直间断表示活动集的修改位置

图7. 系数的精确路径集，其中x轴为 ${\lambda }_{\theta }^{\ast }$，垂直间断表示活动集的修改位置

图8. 系数的精确路径集，其中x轴为步骤数，垂直间断表示活动集的修改位置

表2. 南非心脏病数据使用三种方法的对比

对于该数据集，我们分别建立Spike and Slab Lasso逻辑回归模型和Lasso逻辑回归模型，并且均采用预测校正算法来求解，其中 ${\lambda }_1$ 被设置为0.1， ${\lambda }_0$ 由十折交叉验证获得，结果如表2所示。可以看出，尽管Lasso逻辑回归模型的步骤数少于SSL逻辑回归，但无论θ是否固定，SSL逻辑回归模型的MSE，Bias，AUC，Misclassification均优于Lasso逻辑回归。这正是因为(12)式，我们所提出的SSL逻辑回归模型是自适应于数据的。另外，可以看出，不可分离惩罚的情况下也就是 $\theta ~B\left(1,p\right)$ 时，和可分离惩罚也就是 $\theta =0.5$ 的情况相比，后者的误分类率略优于前者，而前者的MSE，Bias，AUC和步骤数均优于后者。该表说明了在心脏病数据集下spike and slab lasso方法的性能优于Lasso方法。

7. 结论

在本文中，我们提出了Spike and Slab Lasso逻辑回归模型，将两个拉普拉斯密度的混合先验置于单个坐标上，可以自适应地收缩系数，分别讨论了权重θ固定和不固定两种情况。当θ固定时，导致可分离的惩罚，当θ不固定时，令θ服从beta分布，此时导致不可分离的惩罚，并利用预测校正算法分别求解两种惩罚模型。在模拟学习和南非心脏病数据集的研究中，将这两种模型与Lasso逻辑回归模型做对比，说明了本文所提出的Spike and Slab Lasso逻辑回归模型的有效性。

基金项目

本文由国家社会科学基金项目(21BTJ045)资助。

文章引用

齐 琪,张 齐. 基于预测校正算法的Spike and Slab Lasso逻辑回归模型
The Spike and Slab Lasso Logistic Regression Model Based on Prediction Correction Algorithm[J]. 应用数学进展, 2023, 12(01): 292-307. https://doi.org/10.12677/AAM.2023.121032

参考文献