广义线性模型的平方根Lasso选择性推断

Selective Inference of Generalized Linear Models by the Square Root Lasso

摘要: 已经有很多人对线性模型的相关问题做出了大量的选择性推断工作。但是，其适用范围并不全面。实际上，我们会遇到很多非正态非连续的数据，并且随机误差项是异方差等方面的问题，故本文将一些选择性推断的工作推广到广义线性模型上。我们建议使用数据雕琢方法，对数据样本先拆分再进行推断。另外，很多的模型变量选择的工作，都是基于已知噪声水平的。但是，在很多实际情况下，误差方差是未知的，而且在高维数据中，对它的估计存在一定难度。本文中，我们使用平方根lasso，进行变量选择，证明出对于广义线性模型，调优参数的选择不受数据中噪声波动的影响。因此，平方根lasso比lasso调优参数选择更方便，应用范围更广泛。模拟结果表明，使用数据雕琢得到的参数的置信区间更小。

Abstract: In view of the fact that many people have done a lot of work on the selective inference of linear models, however, its scope of application is not comprehensive. In fact, we will encounter many problems such as non normal and non continuous data, and the random error term is heteroscedasticity. Therefore, this paper extends some selective inference work to generalized linear model. We suggest that we use the data carving to split the sample first and then infer. In addition, many of the work of model variable selection are based on the known noise level. However, in practical cases, the error variance is unknown, and in high-dimensional data, it is difficult to estimate. In this paper, we use the square root lasso to select variables, and prove that for the generalized linear model, the selection of tuning parameters is not affected by the noise fluctuations in the data. So the square root lasso is more convenient and widely used than the lasso, and the application is more extensive. Simulation shows that we get the shorter interval length by data carving.

1. 引言

完整的统计调查包括模型选择和推断两个过程。选择过程：我们使用数据来选择模型、进而确定出模型的假设检验、以及关于模型的任何其他问题。推断过程：我们根据所选择出的模型和数据，回答选择过程中提出的问题。简单地说，选择过程决定要问什么问题，推断过程决定回答那些问题。我们以前所做的工作大都是先用所有的数据进行模型选择，然后再进行模型推断。这样两个过程会重复利用数据，造成结果并不那么地精确。以至于，后来统计学家们思考有没有更简单精确的方法进行模型的选择性推断，数据(样本)分割方法 [1] 将数据分成两个独立随机部分，其中一部分数据用于模型选择，另一部分用于推断。Wasserman & Roeder [2] 使用样本分割来解决高维模型中的一致变量选择问题。Meinshausen et al. [3] 通过数据分割处理高维回归模型中的显著性问题，由此产生的p值可以控制错误发现率(FDR)和族错误率(FWER)。样本分割方法解决了控制选择性误差的问题，但这种方法减少了用于模型选择和推断过程中的数据量，成本较高。本文采用数据雕琢的方法，利用其中的一部分数据进行模型选择，另一部分数据和模型选择的剩余数据进行推断，它相对于数据分割方法提高了对于数据的利用率。

针对于线性模型，它仅在响应变量为正态随机变量时适用。Nelder [4] 和Wedderburn [5] 将线性模型推广到广义线性模型。广义线性模型可以适用于以下情况：响应变量为正态和非正态变量，响应变量的分布可以是任意指数族分布。比如，二项分布和伽马分布。可以看出，广义线性模型能很好地解决回归模型中随机误差项的异方差问题。所以，把线性模型推广到广义线性模型上研究问题有重要意义。

近些年来，如生物信息、金融管理等领域产生的高维数据为模型选择带来了更大的挑战。这些领域的实验数据的维数有的甚至超过样本量大小。面对这样的高维数据，我们需要建立一个合适的数学模型，并且还要尽量使用少而有效的数据来进行建模分析，这就产生了一些新的模型选择方法。传统的模型选择方法包括最优子集选择法 [6] [7]，AIC [8] 和BIC [9] 等。当然，还有现在常用的模型选择方法，Robert Tibshirani于1996年提出使用lasso [10]，采用惩罚项压缩一些不重要的变量使其系数变为零。此外，与岭回归相比较，lasso对于重要变量的系数压缩较轻，这样提高了参数估计的精度。自适应lasso [11] 和Relaxed lasso [12] 可以提高参数估计的准确性和一致性。Belloni et al. [13] 提出平方根lasso调优参数γ的选择与误差方差无关。这是相较于lasso方法，用平方根lasso选择模型的一个优点。利用平方根lasso选择变量对正态分布数据进行了大量的研究，但对非正态分布数据的研究还比较少。在本文中，我们尝试使用平方根lasso为广义线性模型选择变量，研究其拥有的特点。

在模型选择后，研究选定模型的推断工作，很多学者已经做了大量的工作。Lockhart et al. [14] 基于lasso拟合值，推导出了原假设下的渐近检验统计量。Lee & Taylor [15]，Lee et al. [16] 提出使用固定调优参数值的lasso进行精确检验。Belloni [17]，Belloni et al. [18]，Zhang & Zhang [19]，重点研究高维稀疏回归模型的推断。Tian et al. [20] 提出在高维线性回归模型中，用平方根lasso进行方差未知的选择性推断。Lehmann & Romano [21] 主要处理了关于指数族的检验问题。Lehmann & Scheffé [22] 提出了一种对于指数族的一致最优势无偏(UMPU)检验的构造方法。Benjamini & Hochberg [23] 提出通过控制FDR来解决多重显著性检验问题。Fithian et al. [24] 提出了通过控制选择性第I型错误来进行模型选择后的推断。Shi et al. [25] 将后选推断推广到广义线性模型，并提出了使用惩罚最小二乘法进行后选推断的新方法。综上所述，大部分工作是在数据服从正态分布时，对线性模型的选择性推断，而对高维异方差数据的研究较少。因此，我们基于平方根lasso方法对广义线性模型的选择性推断。

本篇文章的研究内容主要包括以下几个部分。在第二章中，我们提出广义线性模型的选择性推断包括两个过程：模型选择阶段和模型推断。我们应该注意控制在选定模型下的选择性第I类错误率。第三章，指出了我们使用平方根lasso和lasso选择变量，并且研究了选择事件的特征形式，将选择事件简写成有关于数据的表达式。在第四章中，我们给出了指数族分布的UMPU检验。第五章，给出了指数族分布的参数的充分统计量，以及参数所选模型下的分布形式。第六章，是对于线性回归模型和逻辑回归模型的模拟。

2. 广义线性模型

对于给定的数据 $\left(x_i,y_i\right),i=1,2,\cdots ,n$，是由数据 $x_i=\left(x_{i1},x_{i2},\cdots ,x_{ip}\right)$ 和数据 $y_i,i=1,2,\cdots ,n$ 组成，令 $X=\left\{x_{ij}\right\}$ 是设计矩阵，响应向量 $y\in R^n$。使用下面的广义线型模型来描述设计矩阵和响应向量的关系，

$g\left(y\right)=X\beta +e$ 2-(1)

其中， $\beta$ 是一个p维参数向量，e是一个n维的随机向量，假设它满足以下条件： $E\left(e\right)=0$， $Var\left(e\right)={\sigma }^2{I}_n$。g(·)是连接函数，这里把响应向量和设计矩阵二者进行了连接。

在广义线性模型中，被解释向量y服从的指数族分布可以被表示成如下的形式，

$y~\exp \left\{\frac{yX\beta -B\left(X\beta \right)}{a\left(\phi \right)}\right\}+c\left(y,\phi \right)$ 2-(2)

其中， $\beta$ 是感兴趣参数， $\phi$ 是讨厌参数。

在没有使用选择后推断方法的时候，统计人员在进行分析数据之前，往往根据经验指定模型一个模型M以及要检验的假设。在选定模型M和指定的原假设下，对于原假设的水平 $\alpha$ 检验，必须控制第I型错误率，

$P_{M,H_0}\left(\text{reject}{H}_0\right)\le \alpha$ 2-(3)

从上面的公式可以看到，在计算第I类错误率时，是假定了数据来自模型M，和原假设为真的条件。此时，若我们选定的模型M是不合理的，我们不能保证上面的拒绝率是小于水平 $\alpha$ 的。

因为在大多数的统计工作中，完全排除模型选择这一过程也是不现实的。所以统计专家会利用数据去检查他们的模型，如果诊断出问题，他们就会调整自己的模型。我们认为，如果模型的确定，是通过先用我们手里的数据进行模型选择的，随之，确定原假设。我们就能有效地控制选择性第I型错误率，

$P_{M,H_0}\left(\text{reject}{H}_0|M\text{selected}\right)\le \alpha$ 2-(4)

如果用于模型选择和推断的数据是独立的，则上述条件概率(2-(4))可以写成无条件概率(2-(3))，从而得到如下公式，

$P_{M,H_0}\left(\text{reject}{H}_0|M\text{selected}\right)=P_{M,H_0}\left(\text{reject}{H}_0\right)\le \alpha$ 2-(5)

由此，我们通过数据分割方式将数据y分成两部分，即数据 $y=\left(y_{\text{1}},y_{\text{2}}\right)$，让y₁独立于y₂。那么，将样本数据y₁用于模型选择，样本推断过程仅仅依赖于数据y₂。基于数据y₂的水平 $\alpha$ 的检验满足公式(2-(4))，因此，基于数据y₂的检验可以控制选择性第I类错误率(2-(3))。

数据拆分方法方便了我们的模型选择和推断的工作，也易于让人们理解，这就是当时这种方法被提出后大受欢迎的理由。但是，该方法也同时减少了模型选择和推断过程中的数据量，要付出很大代价。本文采用数据雕刻的方法，利用一部分数据y₁进行模型选择，另一部分数据y₂和选择模型中的数据一起来用于推断。此方法与数据分割方法相比，提高了数据的利用率，进而能够提高参数估计的精度。

3. 通过平方根Lasso选择变量

对于广义线性模型，我们考虑使用平方根lasso选择变量。它是对lasso方法的修正，通过解下面的公式，可以得到 $\beta$ 的最优估计，

$\hat{\beta }=\arg \min {\Vert g\left(y\right)-X\beta \Vert }_2+\gamma {\Vert \beta \Vert }_1$ 3-(1)

其中，参数 $\gamma$ 是调优参数。第一项是广义线性模型的最小二乘目标函数，第二项是平方根lasso的惩罚项，使得许多的变量系数变成零，由此，平方根lasso产生了稀疏解，保留了lasso方法选变量的优点，并且通过定理1可观察到在调优参数的选择上要更加方便。

定理1. 平方根lasso调优参数 $\gamma$ 的选择独立于方差 $\sigma$。

以下集合定义为平方根lasso选择的模型，

$\hat{M}=\left\{j:{\hat{\beta }}_j\ne 0\right\}$

平方根lasso选择的变量符号定义如下，

$\hat{s}=sign\left\{{\hat{\beta }}_j:{\hat{\beta }}_j\ne 0\right\}$

注意到事件 $\left\{\hat{M}=M\right\}$ 是通过平方根lasso进行变量选择产生的，接下来我们将继续注意这个选择事件。并将此选择事件进行具体的公式表达。

对于(3-(1))式的KKT条件的解，其结果等价于下式，

$\begin{array}{l}\frac{X^{\text{T}}\left(g\left(y\right)-X\hat{\beta }\right)}{{\Vert g\left(y\right)-X\hat{\beta }\Vert }_2}=\gamma \hat{s}\\ \text{sign}\left({\hat{\beta }}_{\hat{M}}\right)={\hat{s}}_{\hat{M}}\\ {\Vert {\hat{s}}_{-\hat{M}}\Vert }_{\infty }\le 1\end{array}$ 3-(2)

子集 $\hat{M}=\left\{i\in \left\{1,2,\cdots ,p\right\}:\left|{\hat{s}}_j\right|=1\right\}$ 是使用平方根lasso方法得到的选择模型。 ${\hat{\beta }}_j\ne 0,\left|{\hat{s}}_j\right|=1$ 意味着选定模型的系数是非0的。明显地，从公式(3-(2))可以看出，有些系数是0。但对于几乎所有的参数 $\gamma$，上述集合意味着预测变量具有非零系数。选择事件 $\left\{\hat{M}=M\right\}$ 是通过对对每一个选择事件 $\left\{\hat{M}=M,\hat{s}=s\right\}$ 的中的符号函数s求并集得到的，并将选择事件 $\left\{\hat{M}=M,\hat{s}=s\right\}$ 表达成数据y的表达式。

引理1.X是设计矩阵，M和s表示所选变量及其符号的候选集合。给出以下定义，

$\begin{array}{l}{v}^{*}\left(M,s\right)={\left(X_M^{\text{T}}{X}_M\right)}^{-1}\left(X_M^{\text{T}}g\left(y\right)-\gamma D\left(y\right)s\right)\\ u^{*}\left(M,s\right)=X_{-M}^{\text{T}}{\left(X_M^{\text{T}}\right)}^{+}s+\frac{1}{\gamma D\left(y\right)}{X}_{-M}^{\text{T}}\left(I-P_M\right)g\left(y\right)\end{array}$ 3-(3)

其中， $D\left(y\right)={\Vert g\left(y\right)-X\hat{\beta }\Vert }_2$， $P_M=X_M{\left(X_M^{\text{T}}{X}_M\right)}^{-1}{X}_M^{\text{T}}$ 是 $X_M$ 的投影矩阵。选择事件可以被表示成如下形式，

$\left\{\hat{M}=M,\hat{s}=s\right\}=\left\{\text{sign}\left(v^{*}\left(M,s\right)\right)=s,{\Vert u^{*}\Vert }_{\infty }<1\right\}$ 3-(4)

注意到此时， $\left\{{\Vert u^{*}\Vert }_{\infty }<1\right\}$ 是无效约束，表示当某些变量未被选入模型 $\hat{M}$。然而，当变量被选择进入

模型时，使用条件 $\left\{\text{sign}\left(v^{\text{*}}\right)=s\right\}$ 表示有效约束。

引理2. $v^{*}$ 和 $u^{*}$ 是被定义成上面(3-(3))式。它们可以进一步被重写为如下的形式，

$\left\{\text{sign}\left(v^{*}\left(M,s\right)\right)=s,{\Vert u^{*}\Vert }_{\infty }<1\right\}=\{\left(\begin{array}{c}{C}_0^{*}\left(M,s\right)\left(I-P_M\right)\\ C_1^{*}\left(M,s\right)\end{array}\right)g\left(y\right)<\left(\begin{array}{c}{b}_0^{*}\left(M,s\right)\\ b_1^{*}\left(M,s\right)\end{array}\right)\}$ 3-(5)

其中， $C_0^{\text{*}}\left(I-P_M\right)$， $b_0^{*}$ 对应约束条件 $\left\{{\Vert u^{*}\Vert }_{\infty }<1\right\}$， $C_1^{*}$， $b_1^{*}$ 对应于约束条件 $\left\{\text{sign}\left(v^{*}\left(M,s\right)\right)=s\right\}$，它们被定义成下面的形式，

$C_0^{*}\left(M,s\right)=\frac{1}{\gamma D\left(y\right)}\left(\begin{array}{c}{X}_{-M}\\ -X_{-M}^{\text{T}}\end{array}\right)$

$b_0^{*}\left(M,s\right)=\left(\begin{array}{c}1-X_{-M}^{\text{T}}{\left(X_{-M}^{\text{T}}\right)}^{+}s\\ 1+X_{-M}^{\text{T}}{\left(X_{-M}^{\text{T}}\right)}^{+}s\end{array}\right)$

$C_1^{*}\left(M,s\right)=-\text{diag}\left(s\right){\left(X_M^{\text{T}}{X}_M\right)}^{-1}{X}_M^{\text{T}}$

$b_1^{*}\left(M,s\right)=-\gamma D\left(y\right)\text{diag}\left(s\right){\left(X_M^{\text{T}}{X}_M\right)}^{-1}s$

定理2. 基于引理1和引理2，令 $C=\left(\begin{array}{c}{C}_0^{\text{*}}\left(I-P_M\right)\\ C_1^{\text{*}}\end{array}\right)$， $b=\left(\begin{array}{c}{b}_0^{\text{*}}\left(M,s\right)\\ b_1^{\text{*}}\left(M,s\right)\end{array}\right)$，我们选择事件的形式，

$\left\{\hat{M}=M,\hat{s}=s\right\}=\left\{C^{*}\left(M,s\right)g\left(y\right)<b^{*}\left(M,s\right)\right\}$

所以，选择事件 $\left\{\hat{M}=M\right\}$ 被表示成以下关于数据y的不等式的并集，

$\left\{\hat{M}=M\right\}={\displaystyle \underset{s\in \left\{-1,1\right\}}{\cup }\left\{C^{*}\left(M,s\right)g\left(y\right)<b^{*}\left(M,s\right)\right\}}$

4. 广义线性模型的选择性推断

4.1. 选定模型下的参数推断

在上一章中，已经得到了选定的模型 $\hat{M}$，和它所对应的选择性事件 $\left\{\hat{M}=M\right\}$。在本章中，我们考虑以选定的模型 $\hat{M}$ 作为条件，响应向量y的分布的推断。此分布表示为，

$L\left\{y|\left\{\hat{M}=M\right\}\right\}$ 4-(1)

数据y服从带有干扰参数的p维参数的指数族分布，

$f_{\theta ,\phi }\left(y\right)=\exp \left\{{\theta }^{\prime }T\left(y\right)+{\phi }^{\prime }U\left(y\right)-c\left(\theta ,\phi \right)\right\}{h}_0\left(y\right)$ 4-(2)

其中， $\theta$ 是感兴趣参数， $\phi$ 是讨厌参数。T(y)和U(y)是它们相对应的充分统计量，它们分别是k维和p-k维的向量。

对于任意给定的事件E，在 $y\in E$ 的条件下，y的这个条件分布也是指数族分布，并且与数据y自身的分布拥有相同的参数和充分统计量，即，可以被表达成下面的形式，

$f_{\theta ,\phi }\left(y|y\in E\right)=\exp \left\{{\theta }^{\prime }T\left(y\right)+{\phi }^{\prime }U\left(y\right)-c_E\left(\theta ,\phi \right)\right\}{h}_0\left(y\right)1_E\left(y\right)$ 4-(3)

我们上面提到的模型选择事件 $\left\{\hat{M}=M\right\}$，它作为任意事件E的子集，故样本数据y在选定模型 $\hat{M}$ 下也是服从指数族，且它们的参数和充分统计量都是相同的。例如，在第三章中，我们提到的模型选择事件 $\left\{\hat{M}=M\right\}$，此选择事件被表示成数据y的相关表达式，是事件 $\left\{y\in E\right\}$ 的一中特例，故y的条件密度函数具有类似的形式，也可以写成指数族分布，

$f_{\theta ,\phi }\left(y|\hat{M}=M\right)=\exp \left\{{\theta }^{*}{}^{\prime }{T}^{*}\left(y\right)+{\phi }^{*}{}^{\prime }{U}^{*}\left(y\right)-c^{*}\left({\theta }^{*}\right)\right\}{h}_0^{*}\left(y\right)1_{\hat{M}}\left(y\right)$ 4-(4)

其中 ${\theta }^{*}$， ${\phi }^{*}$ 为选定模型的参数， $T^{*}\left(y\right)$ 和 $U^{*}\left(y\right)$ 为相应的充分统计量。接下来，我们通过消除干扰参数，来对感兴趣的参数进行推断。

如果一维参数 $\theta$ 是感兴趣参数，其它参数 $\phi$ 是干扰参数，T和U是它们分别对应的充分统计量，那么，当讨厌参数的充分统计量 $\left\{U=u\right\}$ 作为已知条件时，以下条件分布仅依赖于感兴趣参数 $\theta$，

$\left(T|U=u\right)~{\psi }_{\theta }\left(t|u\right)=\exp \left\{{\theta }^{\prime }t-c_{\psi }\left(\theta |u\right)\right\}{\psi }_0\left(t|u\right)$ 4-(5)

上述条件分布通过在U上的条件分布消除了讨厌参数φ，当k=1时，我们得到了充分统计量T的单参数指数族。

根据公式(4-(4))，在选择事件 $\left\{\hat{M}=M\right\}$ 下，y的分布仍然是指数族分布。根据公式(4-(5))感兴趣参数 $\theta$ 可由下面的条件分布进行推断，

$L_{\theta }\left(T\left(y\right)|U\left(y\right),\hat{M}=M\right)$

4.2. UMPU选择性检验

给出了原假设 $H_0:\theta \in {\Theta }_0$ 和备择假设 $H_1:\theta \in {\Theta }_1$。如果满足以下不等式，则水平 $\alpha$ 检验 $\varphi \left(y\right)$ 可以成为是选择性无偏的，

$po{w}_{\varphi }\left(\theta |E\right)=E_{\theta }\left[\varphi \left(y\right)|E\right]\ge \alpha ,\theta \in {\Theta }_1$ 4-(6)

在所有满足上述公式的水平 $\alpha$ 的检验中，UMPU水平 $\alpha$ 选择性检验的势函数是一致最优的。

引理3. (UMPU检验)y服从指数族分布(4-(2))，k=1，在水平 $\alpha$ 检验下，考虑以下原假设和备择假设，

$H_0:\theta ={\theta }_0\leftrightarrow H_1:\theta \ne {\theta }_0$，

有UMPU检验 $\varphi \left(y\right)=f\left(T\left(y\right),U\left(y\right)\right)$，其中，

$f\left(t,u\right)=\{\begin{array}{cc}1& t<d_1\left(u\right)或t>d_2\left(u\right)\\ {\pi }_i& t=d_i\left(u\right)\\ 0& d_1\left(u\right)<t<d_2(\; u\; )\end{array}$

其中， $d_i$ 和 ${\pi }_i$ 是下面方程组的解，

$\begin{array}{l}{E}_{{\theta }_0}\left[f\left(T,U\right)|U=u\right]=a\\ E_{{\theta }_0}\left[Tf\left(T,U\right)|U=u\right]=a{E}_{{\theta }_0}\left[f\left(T,U\right)|U=u\right]\end{array}$

由公式(4-(4))可知，当y服从指数族分布时，条件分布 $L\left(y|\left\{\hat{M}=M\right\}\right)$ 也是指数族分布，所以我们得到如下定理。

定理3. (UMPU选择性检验) y遵循指数族分布(4-(2))，k=1，考虑以下在选择事件 $\left\{\hat{M}=M\right\}$ 下的原假设和备择假设，

$H_0:\theta ={\theta }_0\leftrightarrow H_1:\theta \ne {\theta }_0$

有UMPU选择性检验， $\varphi \left(y\right)=f\left(T\left(y\right),U\left(y\right)\right)$，满足，

$f\left(t,u\right)=\{\begin{array}{cc}1& t<d_1\left(u\right)或t>d_2\left(u\right)\\ {\pi }_i& t=d_i\left(u\right)\\ 0& d_1\left(u\right)<t<d_2(\; u\; )\end{array}$

这里的 $d_i$ 和 ${\pi }_i$ 满足下面的方程组，

$\begin{array}{l}{E}_{{\theta }_0}\left[f\left(T,U\right)|U=u,\hat{M}=M\right]=a\\ E_{{\theta }_0}\left[Tf\left(T,U\right)|U=u,\hat{M}=M\right]=a{E}_{{\theta }_0}\left[f\left(T,U\right)|U=u,\hat{M}=M\right]\end{array}$

在下一章中，让我们关注广义线性模型选择性推断的几个特定示例。

5. 例子

例1. 针对于前面讨论的指数族框架下的一些结论，本章中，我们打算举几个具体例子。当数据y来自于多元正态分布时，

$y~N_n\left(\mu ,{\sigma }^2{I}_n\right)$， $\mu =X\beta$

在选定的模型 $\hat{M}$ 的指数族分布形式如下，

$y~\exp \left\{\frac{1}{{\sigma }^2}{\beta }^{\prime }{X}_{M}y-\frac{1}{2{\sigma }^2}{\Vert y\Vert }_2^2-C\left(\beta ,{\sigma }^2\right)\right\}$

若误差方差 ${\sigma }^2$ 已知， ${X^{\prime }}_{j}y$ 是用来推断参数 $\beta$ 的充分统计量，

$L_{{\beta }_j}\left({X^{\prime }}_{j}y|{X^{\prime }}_{M\backslash j}y,\hat{M}=M\right)$.

若误差方差 ${\sigma }^2$ 未知， $\left({X^{\prime }}_{j}y,{\Vert y\Vert }_2^2\right)$ 是参数 $\left(\frac{{\beta }_j}{{\sigma }^2},{\sigma }^2\right)$ 的充分统计量，依据下面的分布进行参数推断，

$L_{\frac{{\beta }_j}{{\sigma }^2}}\left({X^{\prime }}_{j}y|{X^{\prime }}_{M\backslash j}y,{\Vert y\Vert }_2^2,\hat{M}=M\right)$，

同样地，对于参数 ${\sigma }^2$ 的推断，根据下面的分布，

$L_{{\sigma }^2}\left({\Vert y\Vert }_2^2|{X^{\prime }}_{M}y,\hat{M}=M\right)$

例2. 假设数据y由伯努利分布生成，即， $y_{\text{i}}~b\left(1,p_i\right)$。

$y_i~p_i^{y_i}{\left(1-p_i\right)}^{y_i},\ln \frac{p_i}{1-p_i}={\theta }_i,i=1,2,\cdots ,n$

我们把上面这个分布写成指数族的表达形式，

$y~\exp \left\{{\displaystyle \underset{i=1}{\overset{n}{\sum }}{y}_i\ln \frac{p_i}{1-p_i}}+c\left(p_i\right)\right\}$

其中，这个分布有n个参数和相应的充分统计量。如果我们想要对感兴趣参数 ${\theta }_1$ 进行推断，则我们先选定需要的模型 $\left\{\hat{M}=M\right\}=\left\{i=1\right\}$，和讨厌参数 $\left\{{\theta }_i\right\},i=2,3,\cdots ,n$ 及所对应的充分统计量，故对它的推断依赖于以下的条件分布，

$L_{{\theta }_1}\left(y_1|y_2,\cdots ,y_n,\left\{i=1\right\}\right)$，

例3. 若数据 $y=\left(y_1,\cdots ,y_n\right)$ 服从伽马分布 $Gamma\left(\alpha ,\lambda \right)$，即， $y_i~\Gamma \left(\alpha ,\lambda \right)$，

$f\left(y_i\right)=\frac{{\lambda }^{\alpha }}{\Gamma \left(\alpha \right)}{y}_{{}_i}^{\alpha -1}{\exp }^{-\lambda y_i},\alpha ,\lambda >0,i=1,\cdots ,n$

把上述分布进一步写成概率密度的表示形式，

$y~\exp \left\{\left(\alpha -1\right)T\left(y\right)+\lambda N\left(y\right)-c\left(\alpha ,\lambda \right)\right\}$

其中， $T\left(y\right)={\displaystyle {\sum }_{i=1}^n\ln y_i}$， $N\left(y\right)={\displaystyle {\sum }_{i=1}^n{y}_i}$ 分别是参数 $\alpha$， $\lambda$ 的充分统计量。

在任一选择事件E下，我们想要对参数 $\alpha$ 进行推断，则依赖于下面的分布表达式，

$L_{\alpha }\left(T\left(y\right)|N\left(y\right),E\right)$

6. 模拟

6.1. 多元高斯分布

我们用平方根lasso对线性回归模型来选择变量。例如，设置n=100，p=200，X的行是来自变量之间具有成对相关系数为 $\rho =0.\text{4}$ 的等相关多元正态分布，y服从以下的多元正态分布，

$y=X\beta +e,Ee=0,Var\left(e\right)=I_n$

将β设置为200维向量，它的前7个分量为7，调整参数设置成 $E\left(\frac{{\Vert X^{\text{T}}e\Vert }_{\infty }}{{\Vert e\Vert }_2}\right)$。设置置信水平 $\alpha =0.0\text{5}$。

我们考虑将数据分成独立的两部分，以获得Pscreen和参数β的区间长度。其中，Pscreen表示所有的7个变量被成功选择的概率。Carve(n)表示使用样本y的前n个数据进行模型选择，用剩余的数据和选择模型的数据对模型进行推断，以及No carve意味着使用所有的数据进行模型选择，然后再使用所有的数据进行推断。表1展示了我们的模拟运行结果。

Table 1. Pscreen & confidence interval lengths

表1. Pscreen和置信区间长度

此外，Fithian et al. [25] 利用lasso方法对线性模型进行了后选推断。当n = 100和p = 200时，X服从多元正态分布，其变量之间两两相关 $\rho =0.\text{3}$，y服从多元正态分布，列被规范化为长度为1。Pscreen和carve (n)的含义与上面表1中的相同。Power是样本观测值落到拒绝域的概率。表2的模拟运行结果表示如下。

Table 2. Pscreen & power

表2. Pscreen和power

通过对上述结果的分析，用数据雕琢方法推断参数得到的表1的结果与表2的结果是一致的。Pscreen表示通过雕琢选择所有7个变量的概率，概率会随着我们有更多数据用于变量选择而增加，并且它趋于1。上面表1中的Pscreen符合这个理论规律，所得结果与表2的变化趋势一致。在表2中，power随着用于选择模型的数据n的增加而逐渐减小，通过不进行数据雕琢的方式获得的power，比数据雕琢获得的power小。结果表明，通过数据雕琢的得到的power较大。在表1中，通过数据雕琢获得的间隔的长度也随着用于模型选择的数据n的增加而增加。没有用数据雕琢方法得到的置信区间长度介于使用两次不同的数据雕琢后模型选择推断得到的两个区间长度值之间，说明数据雕琢方法可以使得参数区间估计更小。表1中的区间长度和表2中的power变化趋势表明，数据雕琢所得结果与表2的结果一致。

6.2. 伯努利分布

我们用平方根lasso模拟高斯分布，然后将此方法推广到非高斯分布。例如，我们使用数据雕琢方法来选择和推断logistic回归模型，在模型选择阶段，我们使用平方根lasso方法。当n = 100，p = 200，预

测变量y服从伯努利分布，且连接函数为 $g\left(y\right)=\ln \left(\frac{y}{1-y}\right)$ 时，我们有以下广义线性模，

$g\left(y\right)=X\beta +e,Ee=0,Var\left(e\right)=I_n$

将β设置为200维向量，其前7个分量为7。在数据雕琢方法中，我们使用一部分数据进行选择，另一部分数据和选择模型数据一起用于推断。我们用平方根lasso选择模型，并在推断过程中设置置信水平 $\alpha =0.0\text{5}$。我们考虑使用数据雕琢方法，将数据分别拆分成不同的部分，来获得Pscreen，参数β的置信区间和覆盖率。其中，Pscreen和置信区间的含义与第一个模拟中的含义相同，覆盖率表示参数β的覆盖率。表3显示了模拟结果。

Table 3. Pscreen & confidence interval lengths & coverage

表3. Pscreen和置信区间长度和覆盖率

通过分析表3中得到的数据，我们得到以下结论：第一，对于广义线性模型，随着用于模型选择的数据的增加，我们得到的Pscreen会逐渐增大。结果表明，随着模型选择所用数据的增加，模型选择的效果会更好。第二，随着用于模型选择的数据的增加，参数的置信区间将会增大，但是，发现使用数据雕琢会比不使用数据雕琢方法获得的区间长度短。比如，直接用所有数据进行选择性推断获得的置信区间长度值，是在carve (50)和carve (75)之间的。这说明我们可以利用数据雕琢来减小参数的置信区间，从而使我们的选择性推断效果更好。第三，通过数据雕琢获得的参数覆盖率约为0.97~0.98，而不通过数据雕刻获得的覆盖率略大。因此，结合这三个方面，可以利用数据雕琢进行模型选择性推断，这样使得结果更好。

附录

定理1. 考虑广义线性模型，

$g\left(y\right)=X\beta +{\sigma }^2\epsilon ,E\epsilon =0,Var\left(\epsilon \right)=I_n$

$Q\left(\beta \right)={\Vert g\left(y\right)-X\beta \Vert }_2^2$

$\hat{\beta }=\arg \min \frac{1}{2}{\Vert g\left(y\right)-X\beta \Vert }_2+\gamma {\Vert \beta \Vert }_1=\arg \min \frac{1}{2}{Q}^{\frac{1}{2}}\left(\beta \right)+\gamma {\Vert \beta \Vert }_1$

$S:=\frac{\partial Q^{\frac{1}{2}}\left(\beta \right)}{\partial \beta }=\frac{\frac{\partial Q\left(\beta \right)}{\partial \beta }}{2Q^{\frac{1}{2}}\left(\beta \right)}=\frac{Exe}{E{\Vert e\Vert }_2}=\frac{Ex\epsilon }{E{\Vert \epsilon \Vert }_2}$

S是不依赖于误差方差的，故S不受问题中的所有噪声信息的影响。

设置参数 $\gamma$ 来控制噪声，(3-(1))的次梯度是最佳条件， $-S_j+\gamma >0,S_j+\gamma >0$。即， $\left|S_j\right|+\gamma >0,j=1,2,\cdots ,p$。

设置参数 $\gamma$，保证上述公式会成立，当 $\gamma >\max \left|S_j\right|={\Vert S\Vert }_{\infty }$ 时，满足条件。通常，我们将参数设置为 $\gamma =k{\Vert S\Vert }_{\infty }$。

引理1. 重写KKT条件， $X_{-\hat{M}}$ 表示不在模型中的变量。

$\begin{array}{l}{X}_{\hat{M}}^{\text{T}}\left(X_{\hat{M}}{\hat{\beta }}_{\hat{M}}-g\left(y\right)\right)+\gamma D\left(y\right){\hat{s}}_{\hat{M}}=0\\ X_{-\hat{M}}^{\text{T}}\left(X_{\hat{M}}{\hat{\beta }}_{\hat{M}}-g\left(y\right)\right)+\gamma D\left(y\right){\hat{s}}_{-\hat{M}}=0\\ {\Vert {\hat{s}}_{-\hat{M}}\Vert }_{\infty }<1\\ sign\left({\hat{\beta }}_{\hat{M}}\right)={\hat{s}}_{\hat{M}}\end{array}$

由于KKT条件对于我们问题的解是充分必要的，所以当且仅当以下公式成立时，才能得到解，

$\begin{array}{l}{X}_{\hat{M}}^{\text{T}}\left(X_{\hat{M}}{v}^{*}-g\left(y\right)\right)+\gamma D\left(y\right)s=0\\ X_{-\hat{M}}^{\text{T}}\left(X_{\hat{M}}{v}^{*}-g\left(y\right)\right)+\gamma D\left(y\right)u^{*}=0\\ {\Vert u^{*}\Vert }_{\infty }<1\\ sign\left(v^{*}\right)=s\end{array}$

从前两个公式，我们可以得到公式(3-(3))，最后两个公式，我们可以得到公式(3-(4))。

引理2. 我们可以重写两个约束条件，

$\begin{array}{c}\left\{sign\left(v^{*}\left(M,s\right)\right)=s\right\}=\left\{diag\left(s\right)v^{*}>0\right\}\\ =\left\{diag\left(s\right){\left(X_M^{\text{T}}{X}_M\right)}^{-1}\left(X_M^{\text{T}}g\left(y\right)-\gamma D\left(y\right)s\right)>0\right\}\\ =\left\{C_1^{*}\left(M,s\right)g\left(y\right)<b_1^{*}\left(M,s\right)\right\}\end{array}$

$\begin{array}{c}\left\{{\Vert u^{*}\Vert }_{\infty }<1\right\}=\left\{-1<X_{-M}^{\text{T}}{\left(X_M^{\text{T}}\right)}^{+}s+\frac{1}{\gamma D\left(y\right)}{X}_{-M}^{\text{T}}\left(I-P_M\right)g\left(y\right)<1\right\}\\ =\left\{C_0^{*}\left(M,s\right)g\left(y\right)<b_0^{*}\left(M,s\right)\right\}\end{array}$

定理3.

由公式(4-(2))和(4-(3))可知，当数据y服从指数族分布时，条件分布 $L\left(y|y\in E\right)$ 也是指数族分布。选择事件 $\left\{\hat{M}=M\right\}$，作为事件E的特例，自然有 $L\left(y|\hat{M}=M\right)$ 服从指数族分布。我们得到了定理3。

NOTES

^*通讯作者。