大数据时代为我们带来了海量的个体信息，这些个体间的联系构成了纵横交错的网络数据。网络数据涵盖了社交网络^[1]、在线营销^[2]、医疗数据^[3]以及气候分析^[4]等一系列当代应用。这些数据通常具有成千上万的维度，这就使得高维数据的作用越来越突出。研究这些高维网络数据间的联系具有非常实际的意义。图模型可以作为研究变量条件关联的一种工具^[5]。在图模型中，图模型的条件独立性可以由它的精度(逆协方差)矩阵来确定^[6]。在高斯图模型中，图模型的条件独立结构完全由精度矩阵的零元素来确定^[7]。因此，挖掘网络数据变量间联系也就是恢复图模型的精度矩阵。

然而由于高维数据的出现，往往会出现变量大于样本容量的情形。在这种情况下，用样本协方差矩阵的逆来估计精度矩阵不再可行。即使可行，求解一个高维矩阵的逆，在内存和时间上的成本也是巨大的。这就使得高维图模型精度矩阵的求解成为各学者研究的热点。当前的研究主要集中在精度矩阵的点估计上^[8-15]。我们知道，通过置信区间来评估估计量的精度是很重要的^[16-17]。然而由于现有估计量抽样分布的复杂性，目前对高维稀疏精度矩阵的置信区间还很少涉及。

近年来，对于高维统计推断的研究主要转向高维线性和广义线性模型回归系数上^[18-21]。特别是van de Geer等^[19]对Lasso估计量进行KKT(Karush-Kuhn-Tucker)反转得到高维线性回归系数的渐近正态统计量。后来他们又将该方法推广到广义线性模型中。遵循文献[19]的思想，Jankov和van de Geer等^[22]将这种消除惩罚偏差的思想运用到精度矩阵的统计推断中。他们通过对Glasso(graphical Lasso)^[11]估计量进行KKT反转得到De-Glasso(De-sparsified graphical Lasso)统计量，并给出渐近正态性的理论保证，进而得到精度矩阵各个元素的置信区间。遵循这种KKT反转的思想，文献[23-24]都通过对精度矩阵点估计的似然函数KKT反转得到不同设定下的渐近正态统计量。还有一类统计推断方法是Ren等^[14]提出的ANT(asymptotically normal estimation and then do thresholding)方法，它是每2个变量对剩余变量进行回归得到一个渐近正态估计量。

本文研究发现，置信区间的计算成本主要来自于精度矩阵的点估计。上述提出的第1类置信区间由于各自点估计具有复杂的似然函数以及需选择调优参数而计算不有效。第2类置信区间由于ANT点估计需要进行O(p²)次scaled Lasso回归而使得计算成本加大。当维度适中时，两类方法都可以有效地计算精度矩阵的置信区间，然而当维度逐渐增加时，上述方法的计算效率开始变低。为解决上述问题，本文受文献[22]渐近正态统计量的启发提出De-ISEE统计量，并给出对应元素的置信区间。

本文提出的De-ISEE统计量是ISEE点估计的简单运算，且ISEE估计量有成熟的算法，这使得De-ISEE统计量易于计算。相比较其他方法，De-ISEE方法由于ISEE点估计具有可伸缩、易调参等优点可处理超高维精度矩阵。从仿真实验中可以看出，De-ISEE方法得出的置信区间不仅覆盖率更接近于理想覆盖率，而且计算高效。并且本文将De-ISEE方法运用到核黄素数据集以及前列腺肿瘤基因表达数据集，发现De-ISEE方法很好地恢复了变量间的联系，这可作为研究基因学的一种辅助工具。

1 精度矩阵的置信区间 1.1 模型建立

选用高斯图模型来模拟网络数据。设X为p维服从多元高斯分布的随机变量，即

$ \boldsymbol{X}=\left(x_{1}, \cdots, x_{p}\right)^{\mathrm{T}} \sim N\left(\boldsymbol{\mu}, \mathit{\pmb{\Sigma}}^{*}\right). $

(1)

定义μ为p维均值向量，Σ^*为协方差矩阵。设G=(V, E)为高斯无向图，V={x₁, …, x_p}为G的顶点集，E={(i, j)}为高斯图模型边的集合。x_i与x_j满足以下的性质

$ x_{i} \perp x_{j} \mid x_{-(i, j)} \Leftrightarrow(i, j) \notin E, $

也就是x_i和x_j之间无边与x_i和x_j条件独立相互等价。

设Θ^*=(Σ^*)^-1为G的精度矩阵。文献[7]提出，在G中，x_i和x_j的边由精度矩阵的元素Θ_ij^*来表示。当Θ_ij^*≠0时，x_i和x_j之间有边。当Θ_ij^*=0时，x_i和x_j之间无边。设υ={1, …, p}，定义S={(i, j)∈υ×υ|Θ_ij^*≠0}为精度矩阵的支撑集。恢复G等价于恢复Θ^*的支撑集。

设X₁，…，X_n∈ $\mathbb{R}^{p}$为独立同分布于(1)的样本。为方便计算，全文假设均值向量μ=0。

1.2 De-ISEE统计量

本节将基于ISEE估计量(innovated scalable efficient estimation)构造De-ISEE统计量。ISEE方法是由Fan和Lyu^[15]提出，是为了高效地估计超高维精度矩阵的点估计。他们是受创新变换的启示，将估计精度矩阵的问题转化为大协方差矩阵估计问题。它的求解如下:

通过创新变换得到Y=Θ^*X，则有Y~N(0，Θ^*)成立。则估计Θ^*的问题可转变为估计Y的协方差。为了估计Y，Fan等将长向量Y分解成小的子向量，即：

$ \boldsymbol{Y}_{A}=\mathit{\pmb{\Theta}}_{A, A}^{*} \boldsymbol{\eta}_{A}, $

这里$\boldsymbol{\eta}_{A}=\boldsymbol{X}_{A}+\mathit{\pmb{\Theta}}_{A, A^{C}}^{*-1} \boldsymbol{X}_{A^{C}}$，A^C为子集$A \subset\{1, \cdots, p\}$的补集。Θ_{A, B}^*为Θ^*的行为A、列为B的子矩阵。

因为X~N(0, Σ*)，则对于任意的A，有XA⊥XAC~N(-ΘA, A*-1ΘA, AC*XACΘA, A*-1)成立。由此得到以下的多元线性回归模型：

$ \boldsymbol{X}_{A}=\boldsymbol{X}_{A^C} \boldsymbol{C}_{A}+\boldsymbol{\eta}_{A}, $

(2)

这里C_A=-Θ_{A^C, A}^*Θ_{A, A}^*-1是回归系数矩阵，η_A~N(0, Θ_{A, A}^*-1为模型的误差向量。

为估计残差向量η_A，用文献[25]提出的scaled Lasso惩罚回归的方法对模型(2)进行拟合。对于(2)中A的每个节点j，有

$ \boldsymbol{X}_{j}=\boldsymbol{X}_{A^C} \boldsymbol{\beta}_{j}+\boldsymbol{\eta}_{j}, $

(3)

这里X_j, β_j, η_j分别是X_A, C_A, η_A的第j个列向量。对(3)进行scaled Lasso惩罚回归：

$ \begin{gathered} \left(\hat{\beta}_{j}, \hat{\sigma}\right)=\arg \min _{\beta \in \mathbb{R}^{p-|A|, }\sigma \geqslant 0} \\ \left\{\frac{\left\|X_{j}-X_{A^C} \beta\right\|_{2}^{2}}{2 n \sigma}+\frac{\sigma}{2}+\lambda \sum\nolimits_{k \in A^C} \frac{\left\|X_{k}\right\|_{2}}{\sqrt{n}}\left\|\beta_{k}\right\|_{1}\right\}, \end{gathered} $

这里‖υ‖_q为任意给定向量υ的L_q范数，σ=var^1/2(η_j)。可得$\hat{\boldsymbol{\eta}}_{j}=\boldsymbol{X}_{j}-\boldsymbol{X}_{A^C} \hat{\boldsymbol{\beta}}_{j}$, 则有$\hat{\boldsymbol{\eta}}_{A}=\left(\hat{\boldsymbol{\eta}}_{j}\right)_{j \in A}$。

由此得到Θ_{A, A}^*的估计为$\hat{\mathit{\pmb{\Theta}}}_{A}=\left(n^{-1} \hat{\boldsymbol{\eta}}_{A}^{\mathrm{T}} \hat{\boldsymbol{\eta}}_{A}\right)^{-1}$。因此得到Y_A的估计$\hat{\boldsymbol{Y}}_{A}=\hat{\mathit{\pmb{\Theta}}}_{A} \hat{\boldsymbol{\eta}}_{A}$，以此得到

$ \hat{\boldsymbol{Y}}=\left(\hat{\boldsymbol{Y}}_{A_{t}}\right)_{1 \leqslant l \leqslant L,} $

这里(A_l)_l=1^L为索引集(1, …, p}的一个分割，即：U_l=1^LA_l={1, …, p}，以及对于任意的1≤l≠m≤L, A_l∩A_m=Ø。

据此得到Y的样本协方差矩阵的初始ISEE估计量为$\hat{\mathit{\pmb{\Theta}}}_{\mathrm{ISEE}, \text { ini }}=n^{-1} \hat{\boldsymbol{Y}}^{\mathrm{T}} \hat{\boldsymbol{Y}}$，给定阈值τ ≥0，定义ISEE估计量为：

$ \hat{\mathit{\pmb{\Theta}}}=T_{\tau}\left(\hat{\mathit{\pmb{\Theta}}}_{\mathrm{ISEE}, \text { ini }}\right), $

这里$T_{\tau}\left(\hat{\mathit{\pmb{\Theta}}}_{\mathrm{ISEE}, \mathrm{ini}}\right)=\left( {{\omega _{jk}}{I_{\{ \mid {\omega _{jk}}\mid≥\tau }}} \right)$表示$\hat{\mathit{\pmb{\Theta}}}_{\text {ISEE, ini }}=\left(\omega_{j k}\right)$在阈值τ下的值。

由于ISEE估计量计算步骤的复杂性，以致很难研究ISEE估计量的分布性质。为了构建精度矩阵Θ^*基于ISEE估计量的置信区间，本文需要一个优良的统计量。为要消除正则化惩罚为ISEE估计量带来的偏差，本文受文献[19, 22]构造去偏统计量的启发，基于文献[22]提出的渐近正态模型

$ \left(2 \hat{\mathit{\pmb{\Theta}}}-\hat{\mathit{\pmb{\Theta}}} \hat{\mathit{\pmb{\Sigma}}} \hat{\mathit{\pmb{\Theta}}}-\mathit{\pmb{\Theta}}^{*}\right)_{i j} \stackrel{d}{\longrightarrow} N\left(0, \sigma_{i j}^{2} / n\right), $

构建De-ISEE统计量$\hat{\boldsymbol{T}}$如下：

$ \hat{\boldsymbol{T}}=2 \hat{\mathit{\pmb{\Theta}}}-\hat{\mathit{\pmb{\Theta}}} \hat{\mathit{\pmb{\Sigma}}} \hat{\mathit{\pmb{\Theta}}}, $

(4)

其中$\hat{\mathit{\pmb{\Sigma}}}=\frac{1}{n} \sum_{i=1}^{n} \boldsymbol{X}_{i} \boldsymbol{X}_{j}^{\mathrm{T}}$为X的样本协方差矩阵，$\stackrel{d}{\longrightarrow}$为依分布收敛，$\sigma_{i j}^{2}=\mathit{\Theta}_{i i}^{*} \mathit{\Theta}_{i j}^{*}+\mathit{\Theta}_{i j}^{* 2}$。以此得出Θ_ij^*的1-α的置信区间为:

$ \begin{gathered} I_{i j, \alpha} \equiv I_{i j}\left(\hat{\mathit{\Theta}}_{i j}, \alpha, n\right)=\left[\hat{T}_{i j}-\Phi^{-1}(1-\alpha / 2) \hat{\sigma}_{i j} / \sqrt{n}, \right. \\ \left.\hat{T}_{i j}+\Phi^{-1}(1-\alpha / 2) \hat{\sigma}_{i j} / \sqrt{n}\right], \end{gathered} $

(5)

此处$\hat{\boldsymbol{\sigma}}_{i j}^{2}=\hat{\mathit{\pmb{\mathit{\Theta}}}}_{i i} \hat{\mathit{\pmb{\mathit{\Theta}}}}_{i j}+\hat{\mathit{\pmb{\mathit{\Theta}}}}_{i j}^{2}$为文献[22]提出的σ_ij²的相合估计。

由式(4)可看出$\boldsymbol{\hat{T}}$的计算成本主要来自于$\hat{\mathit{\pmb{\mathit{\Theta}}}}$。由于ISEE点估计计算高效，这使得我们可以更快速地计算Θ_ij^*的置信区间。而De-Glasso统计量由于Glasso估计量具有复杂的似然函数以及不可表示性条件的约束，在处理高维矩阵时计算效率低下。下一节仿真实验中可以看出，相比De-Glasso方法，由De-ISEE方法计算的置信区间不仅覆盖率更接近于理想覆盖率，计算也更为高效，更适合于处理超高维精度矩阵。

2 仿真实验

通过仿真实验来检验De-ISEE统计量对网络数据变量联系的恢复效果及计算效率。并与De-Glasso统计量^[22]进行对比研究。

采用文献[22]中的模拟设定，用多元高斯分布N_p(0, Σ^*)来模拟网络数据，其中p为维度。设置精度矩阵Θ^*=tridiag(ρ, 1, ρ)为参数ρ＞0的三对角矩阵。样本量都设为n=500，分别使用De-ISEE方法、De-Glasso方法来计算Θ_ij^*的置信区间。为比较方便，本节考虑置信水平为90 %，95 %，99 %的置信区间。

在仿真实验中，当计算De-Glasso时，使用R包(Glasso版本1.10)。它的调节参数使用五折交叉验证从10值的网格中选择。交叉验证的损失函数使用Cai等^[17]提出的L(Σ，Θ)= $\operatorname{trace}\left(\hat{\mathit{\pmb{\mathit{\Theta}}}}^{\mathrm{T}} {\mathit{\pmb{\Sigma}}}\right)-\log \operatorname{det}(\hat{\mathit{\pmb{\mathit{\Theta}}}})$，其中$\operatorname{det}(\hat{\mathit{\pmb{\mathit{\Theta}}}})$为$\hat{\mathit{\pmb{\mathit{\Theta}}}}$的行列式值，$\hat{\mathit{\pmb{\Sigma}}}$为测试集的样本协方差矩阵，$\hat{\mathit{\pmb{\mathit{\Theta}}}}$为用训练集估计的Glasso值。调节参数选择使损失函数最小的λ值。对于De-ISEE统计量中的调节参数，选择Sun等^[25]提出的$\lambda=2 \sqrt{n / \log p}$。

使用文献[22-24]使用的平均覆盖率指标来比较两种方法的覆盖准确性。分别为: AC_S，AL_S，AC_S^c，AL_S^c。其中AC_S为在支撑集S上的平均覆盖率，它的定义为

$ \mathrm{AC}_{S}=\frac{1}{|S|} \sum_{(i, j) \in S} \hat{\alpha}_{i j}, $

这里$\hat{\alpha}_{i j}=\frac{1}{N} \sum{ }_{k}^{N} I_{k}\left\{\mathit{\Theta}_{i j}^{*} \in I_{i j, a}\right\}$表示真实精度矩阵元素Θ_ij^*被置信区间I_{ij, a}覆盖的经验覆盖率，重复迭代次数设为N=100。同理定义，AC_S^c表示在非支撑集S^C上的平均覆盖率。AL_S表示在S上置信区间的平均长度，它的定义为

$ \mathrm{AL}_{S}=\frac{1}{|S|} \sum_{(i, j) \in S} \hat{C}_{i j}, $

这里$\hat{C}_{i j}$定义为$\hat{C}_{i j}=\frac{1}{N} \sum_{k}^{N}\left\{2 \Phi^{-1}(1-\alpha / 2) \hat{\sigma}_{i j} /\right.\sqrt{n}\}_{k}$，重复迭代次数为N。同理，AL_S^c表示在S^C上置信区间的平均长度。设置参数ρ = 0.3时，计算结果展示为表 1。

表 1分别展示维度p=200, 1 000, 2 000, ρ=0.3时，De-ISEE统计量与De-Glasso统计量构成的置信区间平均覆盖率的比较。当p=2 000时，由于De-Glasso统计量计算时间过长，为避免过多的计算损失，只计算De-ISEE统计量的平均覆盖率。

可以看出，当精度矩阵维度p＜n时，De-ISEE统计量与De-Glasso统计量均能很好地恢复图模型。然而当维度p＞n时，De-Glasso方法在S上的平均覆盖率明显变小，甚至随着p的增加无法计算。而De-ISEE方法当p变大时，它的平均覆盖率基本维持不变。

另一方面，De-ISEE方法在S上的平均覆盖率始终高于De-Glasso方法，而De-Glasso方法在S^C上的平均覆盖率始终高于De-ISEE方法。这表明由De-Glasso方法恢复的图模型更容易丢失变量间联系。这不利于在海量信息中挖掘变量间联系。大数据时代促使网络数据的维度激增，例如在经济分析中，影响经济的因素越来越多，然而各个因素之间的联系却是广泛而又稀疏的。这就需要借助工具发现变量间联系，再针对性分析。De-ISEE统计量很好地恢复了图模型，由它得出的平均覆盖率更接近理想覆盖率。这也在一定程度上，弥补了De-ISEE方法平均置信区间长度长于De-Glasso方法平均置信区间长度这一缺陷。另外，平均置信区间长度也在某种程度上反映了样本的随机误差。

为验证De-ISEE统计量的稳定性，设置另一参数ρ= 0.4，计算结果见表 2。

表 2设置了一个更大的三对角矩阵参数值ρ。同理，为了节省时间与空间损失，在表 2中，当p=2 000时，只计算De-ISEE方法的平均覆盖率。可以发现表 2的结论与表 1保持一致。并且相比较于表 1，可发现当参数值ρ变大时，De-Glasso方法在S上的平均覆盖率明显变小，而De-ISEE统计量基本保持不变，可见De-ISEE方法不受参数值变化的影响，具有较高的稳定性。

为了比较计算效率，使用文献[15]使用的平均计算时间指标, 分别计算置信水平为1-α的De-ISEE方法和De-Glasso方法的CPU运行时间(s)对数的平均数。参数设置为n=500, ρ=0.3, α=0.05, N=100。计算结果为图 1，x轴为维度p，y轴为CPU运行时间(s)对数的平均数。

由图 1可见，随着维度p的增加，De-ISEE的平均运行时间增长缓慢，而De-Glasso的平均运行时间增长快速。当维度很高时，De-Glasso方法甚至无法计算。

综上可见，通过比较De-ISEE方法与De-Glasso方法的平均覆盖率，可发现，在高维网络数据中，De-ISEE方法的平均覆盖率不随维度p以及参数值ρ的变化而显著变化，具有较高的稳定性。而De-Glasso方法的平均覆盖率却对维度p和参数值ρ的变化较为敏感。当精度矩阵的维度p变大时，De-ISEE方法依旧计算高效。在处理高维矩阵时，无论从平均覆盖率，还是平均运算时间来看，De-ISEE方法都比De-Glasso方法有效。

3 实例分析

将De-ISEE方法运用于实际网络数据，选用2个实际数据集，分别运用De-ISEE及De-Glasso方法来恢复网络数据间的联系。网络数据间的联系形成图模型，本节的目标是分别运用2种方法计算图模型精度矩阵的置信区间，比较它们恢复图模型的边以及各自运行时间。第1个数据集是由枯草芽孢杆菌产生的核黄素(维生素B2)数据集，可在hdi R包中获得。第2个数据集是前列腺肿瘤基因表达数据集(prostate tumor gene expression)，它来自于spls R包。

对于这2个数据集，为了更方便地分析，选择前500个方差最大的变量进行建模。随后分离样本，用10个随机选择的观察值来估计500个变量的方差。再次利用估计的方差缩放剩余观察值的设计矩阵。在实例分析中，当计算De-ISEE统计量时，调优参数λ选用Fan和L^[15]提出的λ=B/(n-1+B²)^1/2，这里B定义为B=tq(1-n^1/2/(2p log p), n-1)，tq(a, m)为自由度为m、下分位数为α的t分布。在计算De-Glasso统计量时，与仿真实验同理，选用5折交叉验证选择调节参数。显著性水平都设为α=0.05。

第1个数据集包含n=71个枯草芽孢杆菌的转基因突变体，它有p=4 088个基因表达水平的对数观察变量。我们选择前500个方差最大的变量，也就是一个完整的图包含(₂⁵⁰⁰)条边。计算结果显示，De-ISEE统计量得到28条边为显著的，它的CPU运行时间为1.643 s。De-Glasso统计量识别出5条边为显著的^[22]，它的CPU运行时间为3.756 s。对2种方法识别出的边集进行对比发现，由De-Glasso统计量识别的边是De-ISEE识别边的子集。

第2个数据集包含n=102个样本观察值以及p=6 033个变量。本节分别使用De-ISEE方法与De-Glasso方法计算置信区间以此来识别边。通过De-ISEE方法可识别出121条边为显著的，它的CPU运行时间为1.994 s。De-Glasso方法识别出28条边为显著的，它的CPU运行时间为4.013 s。通过对比再次发现，由De-Glasso统计量识别的边是De-ISEE识别边的子集。且De-ISEE方法的计算速度比De-Glasso方法快。

这样的结论与仿真实验类似，De-Glasso方法识别出的边集更稀疏，容易遗漏个体间的联系。而De-ISEE方法估计的边集基本涵盖真实边集。在个体数量逐渐增加的网络数据中，个体与个体之间的联系广泛而又稀疏，通过De-ISEE方法发现个体之间潜在的联系可以帮助研究者从庞大的数据信息中发现有用信息，进而有针对性地进行分析。在第2个数据集中，De-ISEE方法还可以作为分析人类基因的补充工具。

4 结论与展望

本文提出De-ISEE统计量，并且对其进行了仿真实验以及实例分析。与其他方法相比，由于ISEE点估计具有可伸缩、易调优等优点，De-ISEE统计量得出的平均覆盖率更接近于真实覆盖率，且计算更为快速，可以处理超高维矩阵。这对挖掘高维网络数据间的联系具有很重要的作用。例如在实例分析中分析的核黄素数据集以及前列腺肿瘤基因表达数据集，De-ISEE统计量可以精确恢复变量间的联系，从广泛的信息网中帮助研究者提取重要信息，也可作为研究人类基因联系的补充工具。

本文给出了De-ISEE统计量的实用价值，还尚需要一定的理论支撑，这将会是接下来研究的问题。此外，对图模型构建同时置信区间也将是我们进一步研究的问题。