基因组学研究的一个关键问题是如何确定与疾病相关的基因及其生物途径, 常见的做法是通过将高维基因组数据(如微阵列基因表达数据)与各种临床结果联系起来构建疾病诊断预测模型. 然而, 迄今为止, 虽然许多基因生物标志物研究^[1,2]已经完成, 但目前提出的众多相关方法在临床应用中均难以得到令人满意的结果. 其原因主要在于研究样本量太小^[3,4], 从而导致统计效能降低, 进而得到可信度较低甚至错误的结论. 因此, 充足的样本是产生有效统计分析和结论的必要条件. 另一方面, 数据收集技术的进步促使现行可用生物数据日益增多, 于是有学者提出了数据融合的思想, 即综合多个数据集或有关结果. 然而, 尽管一些基因表达研究有着相同的目标, 但所用数据集通常是来自不同的处理设备、不同的数据平台, 甚至彼此之间具有不同的数值尺度, 从而导致批量效应的存在. 因此, 直接整合不同的基因表达数据将会给统计分析带来巨大挑战.

为解决上述问题, 研究者们做了大量的工作, 主要分为以下两类: 元分析和融合分析^[5] 方法. 元分析即利用统计的概念与方法去收集、整理以及分析之前学者针对某个主题所做的众多实证研究. 然而, 元分析对一些必要条件较为敏感, 稍加违反就可能造成错误性结论^[6]. 融合分析是对不同的数据集进行整合并以此作为研究数据集. 相比元分析, 融合分析具有更多的样本从而更具统计效用. 近年来, 基于融合分析的方法层出不穷, 如Benito等^[7]提出的距离加权判别法(DWD), Johnson等^[8]提出的经验贝叶斯方法(EB), Shabalin等^[9]提出的跨平台标准化方法(XPN), Deshwar等^[10]提出的PLIDA方法以及Deng等^[11]提出的WaveICA方法. 然而, 由于批量效应的存在, 且其来源复杂无法消除, 导致以上方法均可能给融合数据集带来新的系统误差, 使其变得更加复杂. 因此, 直接分析融合后的数据可能会引起一些问题^[12,13], 需要提出一种新的方法来解决数据融合问题.

Kumar等^[14]提出的自主学习(Self-Paced Learning, SPL)方法可以根据模型已经学习的内容自适应地识别简单和困难样本, 并且随着模型训练的不断进行, 越来越多的困难样本进入模型. SPL方法可以在很大程度上克服批量效应, 并且其应用较为广泛, 目前已成功应用于各种机器学习问题^[15]. 此外, Ma等^[16]还对SPL方法的收敛性质进行了补充和讨论, 使其在理论上更加丰富.

除样本规模问题之外, 样本维度是另一研究热点. 许多研究中的样本维数通常远远大于样本数量, 即常见的高维度低样本问题. 这在生物统计中尤为常见, 如基因表达数据. 为解决该问题, 研究者们提出了许多正则化方法, 用于在回归框架中识别与临床表型相关的基因, 如Lasso^[17]、SCAD^[18]、Elastic-Net ^[19]、Fused Lasso^[20]、Lars^[21]、adaptive Lasso^[22]、Group Lasso^[23]以及 $ {L_{1/2 + 2}}$ 混合正则化方法^[24-26]. 然而, 以上正则化方法都存在共同的局限性, 即这些方法仅是从计算或算法的角度出发, 没有利用任何先验知识或信息. 但对于许多复杂的疾病尤其是癌症, 许多生物学途径信息对于了解治疗疾病具有较大的效用, 并且该信息可以从多年的生物医学研究中获得, 故将此种先验信息纳入模型考虑应该会有更好的预测效果.

本文将基于SPL方法构建一个更精确的基因表达预测模型. 首先我们将不同的基因表达数据集融合到一个统一的数据集中, 紧接着在线性回归的背景下将SPL方法与SCAD网络惩罚相结合得到最终的回归预测模型. 具体来说, 该模型由3部分组成: (1) SCAD罚函数. 利用SCAD惩罚来增强模型的稀疏性, 该惩罚不仅为大系数提供了无偏估计, 并且具有较高的理论价值, 例如Oracle性质^[18]; (2)基于网络的惩罚, 利用网络惩罚来实现基因调控网络上相邻节点系数之间的平滑; (3) SPL方法, 促使模型自适应地从简单样本(高置信度样本)向复杂样本(低置信度样本)上过渡. SPL方法对于分析融合数据是至关重要的, 因为融合数据往往存在较大的噪声以及异常值点.

本文接下来内容安排如下: 第1节提出了一个基于SCAD网络惩罚的线性回归模型, 紧接着介绍了自主学习(SPL)方法并将其与SCAD网络惩罚相结合从而得到最终的预测模型; 第2节首先对SCAD网络惩罚函数的理论性质进行简单分析, 包括群组效应以及渐近性质; 然后给出一种求解本文所提出模型的有效算法; 在第3节中, 通过不同情形下的模拟数据以及在乳腺癌细胞系数据集上的分析结果来评估本文所提出模型的预测效果. 第4节是结论与展望.

1 SCAD网络正则化与自主学习方法 1.1 SCAD网络正则化

假设数据集 $D = \{ ({{{x}}_{{1}}},{y_1}),({{{x}}_{{2}}},{y_2}), \cdots ,({x_n},{y_n})\}$ , 其中 ${{{x}}_{{i}}} = {({x_{i1}},{x_{i2}}, \cdots ,{x_{ip}})^{\rm{T}}}$ 表示第 $i$ 个样本, ${y_i}$ 为对应的响应变量, 记 ${{X}} = ({{{x}}_{\rm{1}}},{{{x}}_2}, \cdots ,{{{x}}_p}),{{Y}} = {({y_1},{y_2}, \cdots ,{y_n})^{\rm{T}}}$ . 进一步, 假设各个预测因子 ${{{x}}_{{i}}}\;(i = 1,2, \cdots ,p)$ 经过标准化处理, 响应变量 ${{y}}$ 经过去中心化处理, 从而有:

$\sum\limits_{i = 1}^n {{y_i}} = 0,\sum\limits_{i = 1}^n {{x_{ij}}} = 0,\sum\limits_{i = 1}^n {{x_{ij}}^2} = 1,\;j = 1,2, \cdots ,p$

本文考虑最简单的线性回归模型:

${y_i} = {{{x}}_{{i}}}^{\rm{T}}{{\beta }} + {\varepsilon _i}$

式中, $\;{\bf{\beta }} = {({\beta _1},{\beta _2}, \cdots ,{\beta _p})^{\rm{T}}}$ 为要估计的参数, ${\varepsilon _i}$ 表示均值为0, 方差为 ${\sigma ^{\rm{2}}}$ 的误差项. 上述模型的平方损失函数可以表示为:

$l({{\beta }}) = \sum\limits_{i = 1}^n {{{({y_i} - {{{x}}_i}^{\rm{T}}{{\beta }})}^2}} $

在许多研究当中, 样本维度通常远远大于样本数量, 即高维度低样本问题. 在这种情形下, 线性回归并不能够直接用来估计回归参数. 由此, 引入了正则化方法, 即:

$\ell ({\bf{\beta }},\lambda ) = l({\bf{\beta }}) + P({\bf{\beta }})$

其中, $P({\bf{\beta }})$ 表示正则化项. 高维变量选择中常用的正则化方法为 ${L_1}$ 约束, 即Lasso方法, 具体可以表示为 ${P_{\lambda ,{\rm{Lasso}}}}({\bf{\beta }}) = \displaystyle\sum\limits_{i = 1}^p {{{\left| {{\beta _j}} \right|}_1}}$ , 这里 $\lambda $ 表示任意非负数, 一般可使用 $k$ 折交叉验证方法确定. 由于 ${L_1}$ 罚函数具有奇异性, 故基于 ${L_1}$ 惩罚的线性回归模型可以将一些系数较小的参数压缩为0从而达到变量选择的效果. 但当 $\lambda $ 过大时, $\;{\bf{\beta }}$ 估计量中系数较大的参数会存在较大偏差, 而当 $\lambda $ 过小时, $\;{\bf{\beta }}$ 估计量则不够稀疏. 为克服这一问题, Fan等^[18]提出了SCAD惩罚函数, 其具体形式为:

$ {P}_{\lambda }(\beta )=\left\{\begin{array}{ll} \lambda \left|\beta \right|, & \left|\beta \right|\le \lambda\\ \dfrac{({a}^{2}-1){\lambda }^{2}-{(\left|\beta \right|-a\lambda )}^{2}}{2(a-1)}, & \lambda <\left|\beta \right|\le a\lambda \\ \dfrac{(a+1){\lambda }^{2}}{2}, & \left|\beta \right|>a\lambda \end{array}\right.$

(1)

其中, $a$ 为一个大于2的超参数, 根据文献[18], 可将其设为3.7, 当然也可以通过交叉验证的方法加以确定. 从表达式(1)可以发现, 当 $\left| \beta \right|$ 较小时, 惩罚函数为线性函数; 当 $\left| \beta \right|$ 较大时, 惩罚函数为二次惩罚; 当 $\left| \beta \right|$ 很大时, 惩罚项为常数. SCAD惩罚函数关于 $\;\beta $ 的一阶导函数为:

$\begin{split} {P}_{\lambda }{}^{'}(\beta )=&\left\{\begin{array}{ll} \lambda ,& \left|\beta \right|\le \lambda \\ \dfrac{a\lambda -\beta }{a-1}, & \lambda <\left|\beta \right|\le a\lambda \\ 0, & \left|\beta \right|> a\lambda \end{array}\right.\\ =&\lambda \{I(\beta \le \lambda )+\dfrac{{(a\lambda -\beta )}_+}{(a-1)\lambda }I(\beta > \lambda )\}\end{split}$

(2)

在非0处, 对任意 $\;{\beta _j} \approx {{\textit{z}}_j}$ , 由二阶泰勒展开可得:

${{P}_{\lambda }}({{\beta }_{j}})\approx {{P}_{\lambda }}(| {{{\textit{z}}}_{j}} |)\text+\frac{\text{1}}{\text{2}}\text{ }\left\{\text{ }\frac{{{P}_{\lambda }}^{'}(| {{{\textit{z}}}_{j}} |)}{| {{{\textit{z}}}_{j}} |}\text{ }\right\}\text{ (}\beta _{j}^{2}-\text{{\textit{z}}}_{j}^{2}\text{)}$

(3)

此外, 对于协变量之间存在高相关性的问题, Zou等^[19]提出了Elastic-Net惩罚函数, 其具体表达为 ${P_{{\lambda _1},{\lambda _2},{\rm{enet}}}}({\bf{\beta }}) = \displaystyle\sum\limits_{}^{} {{\lambda _1}{{\left| {{\beta _j}} \right|}_1} + {\lambda _2}} {\left| {{\beta _j}} \right|_2}$ . Zeng等^[27]提出了SCAD- ${L_2}$ 惩罚, 该惩罚同时结合SCAD和 ${L_2}$ 惩罚. 以上几种方法都可以实现群组效应, 即具有强相关性的预测因子要么同时被选中, 要么同时被剔除. 然而, 上述几种方法都是从计算的角度出发, 都没有考虑先验信息, 如网络拓扑信息.

我们知道, 生物统计中的基因交互信息对于识别基因组模式具有重要价值. 该先验信息可以用一个加权图 $G = (N,E,W)$ 来表示, 其中 $N$ 是网络节点集, 表示 $p$ 个预测因子, $E = \{ u\sim v\} $ 是图中边的集合, 表示节点 $u$ 和 $v$ 之间有边相连, $W = \{ w(u,v)\} $ 则表示边上的权重. 近年来, 网络惩罚出现在大量现实应用中, 例如, Li等^[28], Chen等^[29]以及Wang等^[30]利用基于网络的 ${L_1}$ 惩罚对基因组数据进行回归分析并进行变量选择. 在这些研究当中, 网络惩罚函数被定义为拉普拉斯矩阵的二次型. 然而, 在某些情况下, ${L_1}$ 惩罚存在偏差并且可能导致结果不够稀疏, 而SCAD惩罚可以避免过度惩罚并且具有良好的统计性质. 因此, 本文给出基于网络的SCAD-Net惩罚函数(SCAD Network-based penalized fuction, SCAD-Net):

${P_{{\lambda _1},{\lambda _2},\text{scad-net}}}({\bf{\beta }}) = {P_{{\lambda _1},{\rm{scad}}}}({\bf{\beta }}) + {\lambda _2}{{\bf{\beta }}^{\rm{T}}}{{L\beta }}$

(4)

其中, L表示拉普拉斯矩阵, 根据文献[31], 可将其定义为:

$ L(u,v)=\left\{\begin{array}{ll} 1-w(u,v)/{d}_{u}, & u=v,{d}_{u}\ne 0\\ -w(u,v)/\sqrt{{d}_{u}{d}_{v}}, & u,v相邻接\\ \rm{0,} &其他 \end{array}\right.$

其中, ${d_u} = \displaystyle\sum\limits_{u\sim v} {w(u,v)}$ 表示节点 $u$ 的度, ${d_u} = 0$ 说明节点u是孤立的点. L融合了网络的结构信息, 且 $\;{{{\beta }}^{\rm{T}}}{{L\beta }}$ 项使得 $\;{{\beta }}$ 的估计更加光滑. 进一步, 根据L的性质, 表达式(4)可以写成:

$ {P_{{\lambda _1},{\lambda _2},{\rm{scad} \text{-} \rm{net}}}}({\bf{\beta }}) = {P_{{\lambda _1},{\rm{scad}}}}({\bf{\beta }}) + {\lambda _2}\sum\limits_{i\sim j} {{{\left(\frac{{{\beta _i}}}{{\sqrt {{d_i}} }} - \frac{{{\beta _j}}}{{\sqrt {{d_j}} }}\right)}^2}w(i,j)} $

(5)

其中, 参数 ${\lambda _{\rm{1}}}$ 与 ${\lambda _{\rm{2}}}$ 分别控制参数估计的稀疏度与光滑度. 进一步, 在线性回归背景下, 可以得到基于SCAD-Net惩罚的线性回归模型(SCAD-Net penalized Linear regression, SNL):

$ {{\hat \beta }} = \mathop {\min }\limits_{\bf{\beta }} \left\{ l({\bf{\beta }}) + {P_{{\lambda _1},{\rm{scad}}}}({\bf{\beta }}) + {\lambda _2}\sum\limits_{i\sim j} {{{\left(\frac{{{\beta _i}}}{{\sqrt {{d_i}} }} - \frac{{{\beta _j}}}{{\sqrt {{d_j}} }}\right)}^2}w(i,j)} \right\} $

(6)

其中, 第1项表示线性回归的损失函数; 第2项表示SCAD惩罚函数, 保证参数估计的稀疏性, 并且强化结果的可解释性; 最后一项表示基于网络的惩罚函数, 保证参数估计的光滑性, 并且将网络结构信息与基因表达数据相融合.

1.2 自主学习策略

尽管正则化方法在基因数据分析, 变量选择等方面扮演非常重要的角色, 但最终得到的结论却鲜少在临床中得到应用. 这是因为上述结论都是基于小样本数据, 导致结果的可信度较低. 为解决这一问题, 有人提出通过整合不同的数据集来生成人工大样本数据. 然而, 这些数据整合的方法并不能消除内部偏差, 甚至可能会增加新的误差.

受人类学习机制的启发, Kumar等^[14]提出了自主学习(SPL)方法, 该方法首先通过低噪声样本学习一个基础模型, 然后通过高噪声样本学习使模型变得更加稳健, 该方法可以显著提高融合数据集的统计分析效用. 并且Kumar表示, 通过引入一个惩罚项, 可以将自主学习方法视为优化模型, 具体可以表示为:

$\ell ({{\beta }},{{v}}) = \sum\limits_{i = 1}^n {\bigg\{ {v_i}{l_i}({y_i},d({{{x}}_i},{{\beta }}))} + f({v_i},\tau )\bigg\} $

(7)

其中, $\ell ({{\beta }},{{v}}) = \displaystyle\sum\limits_{i = 1}^n {\bigg\{ {v_i}{l_i}({y_i},d({{{x}}_i},{{\beta }}))} + f({v_i},\tau )\bigg\}$ 为每个样本的加权损失, ${l_i}({y_i},d({{{x}}_i},{{\beta }})$ 表示一般损失函数, $d({{{x}}_i},{{\beta }})$ 为决策函数, ${{v}} = {({v_1},{v_2}, \cdots ,{v_n})^{\rm{T}}}$ 为样本集上的加权向量, $\tau $ 为超参数, 用来调整学习步长. $ f\left( {{v}_{i}},\tau \right)$ 表示施加在样本权重上的惩罚, 一般情况下, 我们可以取 $f(v,\tau ) = {\rm{ - }}v\tau $ . 显然, 当 $\tau $ 较小时, 为使式(7)最小, 训练集中只能包括少量的样本, 而当 $\tau $ 增大时, 训练集中又必须加入更多的样本, 样本由少变多, 模型从简单到复杂.

为了加强对融合数据分析的准确性与鲁棒性, 本文将SPL方法与SCAD-Net正则化在线性回归的背景下相结合, 从而得到最终的回归模型(Self-paced learning and SCAD-Net penalized Linear regression, SSNL):

$\begin{split} \ell ({{\beta }},{{v}}) =& \frac{{\rm{1}}}{{\rm{2}}}\sum\limits_{i = 1}^n {({v_i}{l_i}({{\beta }})} - {v_i}\tau ) + \sum\limits_{j = 1}^p {{P_{{\lambda _1}}}(\left| {{\beta _j}} \right|)} + \frac{{{\lambda _2}}}{2}{{{\beta }}^{\rm{T}}}{{L\beta }} \\ = &\frac{{\rm{1}}}{{\rm{2}}}\sum\limits_{i = 1}^n {({v_i}{{({y_i} - {{{x}}_i}^{\rm{T}}{{\beta }})}^2}} - {v_i}\tau ) + \sum\limits_{j = 1}^p {{P_{{\lambda _1}}}(\left| {{\beta _j}} \right|)} + \frac{{{\lambda _2}}}{2}{{{\beta }}^{\rm{T}}}{{L\beta }} \\ \end{split} $

(8)

其中, 第1项表示加权的线性回归模型, 最后两项表示SCAD-Net惩罚函数.

2 理论性质及求解方法 2.1 理论性质

本小节我们给出与SCAD-Net正则化方法相关的性质, 包括群组效应以及在 $p$ 固定且 $n \to \infty $ 情形下的渐近性质.

2.1.1 群组效应

Huang等^[32]证明了SCAD-Net惩罚函数具有群组效应, 如引理1与引理2所示, 其具体证明过程见文献[32].

引理1. 若 $\; {\bf{\hat \beta }}$ 是由表达式(6)计算得到, 对任意 ${{\lambda }_{\text{2}}} > \dfrac{1}{2(a-1)}$ , 若满足 $ {{x}_{i}}={{x}_{j}}$ , 则有:

${\hat \beta _i} = {\hat \beta _j}$

引理1可看作Zou等^[19]中引理2的进一步结果, 其保证在两个预测变量相等时, 估计参数具有群组效应.

引理2. 若 $\; {\hat \beta _i}{\hat \beta _j} > 0$ 且 ${{\lambda }_{\text{2}}} > \dfrac{1}{2(a-1)}$ , 定义:

$ D(i,j) = \frac{{\left| {{{\hat \beta }_i} - {{\hat \beta }_j}} \right|}}{{{{\left| {{y}} \right|}_1}}} $

则有:

$ D(i,j) \le \dfrac{1}{{2{\lambda _2} - \dfrac{1}{{a - 1}}}}\sqrt {2(1 - \rho )} $

其中, $\;\rho = {{x}}_i^{\rm{T}}{{{x}}_j}$ 表示样本相关系数.

引理2给出SCAD-Net惩罚函数群组效应的量化描述, 即在满足以上条件的前提下, 两个参数的差异具有上界约束. 进一步, 若样本相关系数 $\;\rho $ 趋于1, 则两个估计参数几乎相同.

2.1.2 渐近性

SCAD-Net惩罚线性回归的目标函数为:

$ \begin{split} & \sum\limits_{i = 1}^n {{{({y_i} - {{{x}}_i}^{\rm{T}}{{\beta }})}^2}} + \sum\limits_{j = 1}^p {{P_{\lambda _n^{(1)}}}(\left| {{\beta _j}} \right|)} + \frac{{_{\lambda _n^{(2)}}}}{2}{{{\beta }}^{\rm{T}}}{{L\beta }} \\ & =\sum\limits_{i = 1}^n {{{({y_i} - {{{x}}_i}^{\rm{T}}{{\beta }})}^2}} + \sum\limits_{j = 1}^p {{P_{\lambda _n^{(1)}}}(\left| {{\beta _j}} \right|)} \\ & \;\;\;\;+ \lambda _n^{(2)}\sum\limits_{u\sim v} {{{\Bigg(\frac{{{\beta _u}}}{{\sqrt {{d_u}} }} - \frac{{{\beta _v}}}{{\sqrt {{d_v}} }}\Bigg)}^2}w(u,v)} \end{split} $

其中, $\lambda _n^{({\rm{1}})}$ 与 $\lambda _n^{(2)}$ 为样本大小 $n$ 的函数. 定理1给出了估计量的渐近性质.

定理1. 若 $\lambda _n^{(m)}/\sqrt n \to \lambda _0^{(m)} \ge 0(m = 1,2)$ , 并且 ${{C}} = $ $ {\lim _{n \to \infty }}\Bigg(\dfrac{1}{n}\displaystyle\sum\limits_{i = 1}^n {{{{x}}_i}^{\rm{T}}} {{{x}}_i}\Bigg)$ 非奇异, 则: $\sqrt n ({{\mathbf{\hat \beta }}_n} - {\bf{\beta }})\xrightarrow{d}\arg \min (V)$

其中,

$\begin{split} V({{u}}) = & - 2{{{u}}^{\rm{T}}}{{W}} + {{{u}}^{\rm{T}}}{{Cu}} + \lambda _0^{(1)}\sum\limits_{j = 1}^p {\operatorname{sgn} ({\beta _j})} {u_j}\left\{ I(\left| {{\beta _j}} \right|) \le \lambda _0^{(1)} + \frac{{{{(a\lambda _0^{(1)} - \left| {{\beta _j}} \right|)}_ + }}}{{(a - 1)\lambda _0^{(1)}}}I(\left| {{\beta _j}} \right|) > \lambda _0^{(1)}\right\} \\ & + 2\lambda _0^{(2)}\sum\limits_{i\sim j} {\Bigg(\frac{{{\beta _i}}}{{\sqrt {{d_i}} }} - \frac{{{\beta _j}}}{{\sqrt {{d_j}} }}\Bigg)\Bigg(\frac{{{u_i}}}{{\sqrt {{d_i}} }} - \frac{{{u_j}}}{{\sqrt {{d_j}} }}\Bigg)w(i,j)},\;且{{W}}\sim N(0,{\sigma ^2}{{C}}) \end{split}$

证明: 定义

$\begin{split} {V_n}({{u}}) =&\sum\limits_{i = 1}^n {\Bigg\{ {{\Bigg({\varepsilon _i} - \frac{{{{{u}}^{\rm{T}}}{{{x}}_i}}}{{\sqrt n }}\Bigg)}^2} - {\varepsilon _i}^2\Bigg\} } + \sum\limits_{j = 1}^p {\left\{ {P_{\lambda _n^{(1)}}}({\beta _j} + {u_j}/\sqrt n )\right\} } \\ & + \lambda _n^{(2)}\sum\limits_{i\sim j} {\Bigg\{ {{\Bigg(\Bigg(\frac{{{\beta _i}}}{{\sqrt {{d_i}} }} - \frac{{{\beta _j}}}{{\sqrt {{d_j}} }}\Bigg) + \Bigg(\frac{{{u_i}}}{{\sqrt {{d_i}} }} - \frac{{{u_j}}}{{\sqrt {{d_j}} }}\Bigg)\Bigg/\sqrt n \Bigg)}^2}}w(i,j)- {\Bigg(\frac{{{\beta _i}}}{{\sqrt {{d_i}} }} - \frac{{{\beta _j}}}{{\sqrt {{d_j}} }}\Bigg)^2}w(i,j)\Bigg\} \end{split}$

其中, ${{u}} = {({u_1},{u_2}, \cdots ,{u_p})^{\rm{T}}}$ , 注意到 ${V_n}({{u}})$ 在 $\sqrt n ({{\mathbf{\hat \beta }}_n} - {\bf{\beta }})$ 处取得最小值, 且在有限维收敛的情况下有:

$\sum\limits_{i = 1}^n {\Bigg\{ {{\Bigg({\varepsilon _i} - \frac{{{{{u}}^{\rm{T}}}{{{x}}_i}}}{{\sqrt n }}\Bigg)}^2} - {\varepsilon _i}^2\Bigg\} } \to - 2{{{u}}^{\rm{T}}}{{W}} + {{{u}}^{\rm{T}}}{{Cu}}$

进一步, 根据式(2)和式(3)可知:

$ \begin{split}&\sum\limits_{j=1}^{p}{\left\{{{P}_{\lambda _{n}^{(1)}}}({{\beta }_{j}}+{{u}_{j}}/\sqrt{n})\text-{{P}_{\lambda _{n}^{(1)}}}({{\beta }_{j}})\right\}}\\ &=\frac{1}{2}\sum\limits_{j=1}^{p}{\frac{{{P}_{\lambda _{n}^{(1)}}}^{'}(\left| {{{\textit{z}}}_{j}} \right|)}{\left| {{{\textit{z}}}_{j}} \right|}\left\{{{({{\beta }_{j}}+{{u}_{j}}/\sqrt{n})}^{2}}-{{\beta }_{j}}^{2}\right\}}\to\lambda _{0}^{(1)}\sum\limits_{j=1}^{p}{\frac{{{\beta }_{j}}{{u}_{j}}}{\left| {{{\textit{z}}}_{j}} \right|}}\left\{I(\left| {{{\textit{z}}}_{j}} \right|)\le \lambda _{0}^{(1)}+\frac{{{(a\lambda _{0}^{(1)}-\left| {{{\textit{z}}}_{j}} \right|)}_+}}{(a-1)\lambda _{0}^{(1)}}I(\left| {{{\textit{z}}}_{j}} \right|)>\lambda _{0}^{(1)}\right\} \end{split} $

对任意 $\;{\beta _j} \approx {{\textit{z}}_j}$ , 有:

$ \sum\limits_{j = 1}^p {\left\{ {P_{\lambda _n^{(1)}}}({\beta _j} + {u_j}/\sqrt n ) - {P_{\lambda _n^{(1)}}}({\beta _j})\right\} } \to \lambda _0^{(1)}\sum\limits_{j = 1}^p {\operatorname{sgn} ({\beta _j})} {u_j} \Bigg\{ I(\left| {{\beta _j}} \right|) \le \lambda _0^{(1)} + \frac{{{{(a\lambda _0^{(1)} - \left| {{\beta _j}} \right|)}_ + }}}{{(a - 1)\lambda _0^{(1)}}}I(\left| {{\beta _j}} \right|) > \lambda _0^{(1)}\Bigg\} $

同样地, 关于第三项有:

$\begin{split} &\lambda_{n}^{(2)} \sum\limits_{i\sim j} \Bigg\{{ {\left(\left(\frac{\beta_i}{\sqrt {d_i}}-\frac{\beta_j}{\sqrt {d_j}}\right)+\left(\frac{u_i}{\sqrt {d_i}}-\frac{u_j}{\sqrt {d_j}}\right)/\sqrt n\right)}^2w(i,j)-\left(\frac{\beta_i}{\sqrt {d_i}}-\frac{\beta_j}{\sqrt {d_j}}\right)^2 w(i,j)}\Bigg\}\\ &\to 2\lambda_0^{(2)}\sum\limits_{i\sim j}\left(\frac{\beta_i}{\sqrt {d_i}}-\frac{\beta_j}{\sqrt {d_j}}\right)\left(\frac{u_i}{\sqrt {d_i}}-\frac{u_j}{\sqrt {d_j}}\right)w(i,j) \end{split}$

因此, 在有限维收敛的情况下有:

$ {V_n}({{u}})\xrightarrow{d}V({{u}}) $

又 ${V_n}$ 为凸函数且 $V$ 有最小值, 可得:

$\arg \min ({V_n}) = \sqrt n ({{\mathbf{\hat \beta }}_n} - {\bf{\beta }})\xrightarrow{d}\arg \min (V).$

证毕.

2.2 求解方法

本节给出模型SSNL的求解算法, 具体如下:

(1) 固定v更新 $\;{{\beta }}$ 时, 相当于解决SCAD-Net正则化回归问题, 即SNL. 本文利用坐标下降法进行求解, 具体来说, 式(8)关于 $\;{\beta _j}(j = 1,2, \cdots ,p)$ 求导, 可得:

$\begin{split} \frac{\partial \ell }{\partial {{\beta }_{j}}}&=-\sum\limits_{i=1}^{n}{{{v}_{i}}{{x}_{ij}}}({{y}_{i}}-{{{x}}_{i}}^{\rm T}{\beta })+{{P}_{{{\lambda }_{1}}}}^{'}(\left| {{\beta }_{j}} \right|)+{{\lambda }_{2}}{{{\beta }}^{\rm T}}{{{L}}_{j}} \\ & = -\sum\limits_{i=1}^{n}{{{v}_{i}}{{x}_{ij}}}({{y}_{i}}-{{{x}}_{i,-j}}^{\rm T}{{{\beta }}_{-j}}-{{x}_{ij}}{{\beta }_{j}})+\frac{{{P}_{{{\lambda }_{1}}}}^{'}(\left| {{{\textit{z}}}_{j}} \right|)}{\left| {{{\textit{z}}}_{j}} \right|}{{\beta }_{j}} \\ &\;\;\;\;+{{\lambda }_{2}}{{{\beta }}_{-j}}^{\rm T}{{{L}}_{-j,j}}+{{\lambda }_{2}}{{L}_{jj}}{{\beta }_{j}} \end{split} $

令其等于0, 有:

$ {{\hat{\beta }}_{j}}=\dfrac{\displaystyle\sum\limits_{i=1}^{n}{{{v}_{i}}{{x}_{ij}}}({{y}_{i}}-{{{x}}_{i,-j}}^{\rm T}{{{\beta }}_{-j}})-{{\lambda }_{2}}{{{\beta }}_{-j}}^{\rm T}{{{L}}_{-j,j}}}{\displaystyle\sum\limits_{i=1}^{n}{{{v}_{i}}{{x}_{ij}}^{2}}+\dfrac{{{P}_{{{\lambda }_{1}}}}^{'}(\left| {{{\textit{z}}}_{j}} \right|)}{\left| {{{\textit{z}}}_{j}} \right|}+{{\lambda }_{2}}{{L}_{jj}}} $

(9)

具体更新算法如算法1所示.

算法1. SNL

1. 令 $\scriptstyle t = 0,\;{\beta _j}(t) = {\beta _j}^0$ , 其中 $\;\scriptstyle{{\bf{\beta }}^0}$ 为OLS估计量, 各调优参数可利用交叉验证方法得到;

2. 通过式(9)依次更新 $\; \scriptstyle{\beta _j}(t),j = 1,2, \cdots ,p$ ;

3. 令 $\scriptstyle t \leftarrow t + 1$ , 重复步骤2直至:

　　　　 $\scriptstyle\sum\limits_{j = 1}^p {\left| {{\beta _j}(t) - {\beta _j}(t - 1)} \right|} < {10^{ - {\rm{6}}}} $

(2) 固定 $\;{{\beta }}$ 更新v时, 式(8)关于 ${v_i}$ 求导, 可得:

$\begin{split} \frac{{\partial \ell }}{{\partial {v_i}}} = {l_i}({{\bf{\beta }}^t}) - \tau \end{split}$

进而有:

$ {\hat v_i} = \left\{ {\begin{array}{*{20}{c}} 1, &{{l_i}({{\bf{\beta }}^t}) \le \tau } \\ 0, &{{l_i}({{\bf{\beta }}^t}) > \tau } \end{array}} \right. $

(10)

对于样本i, 若其损失小于超参数 $\tau $ , 则可将其视为高质量样本, 相对应的 ${v_i}$ 设为1, 否则设为0. 显然, 对于样本损失小于 $\tau $ 的样本会被纳入模型中. 一旦得到v, 我们进一步增大 $\tau $ 的值, 这样具有更大损失的样本将会进入模型当中, 重复上述步骤直至收敛, 完整算法如算法2所示.

算法2. SSNL

输入: 训练集 $\scriptstyle \{ {{{X}}_{n \times p}},{y_{n \times 1}}\}$ , 超参数 $\scriptstyle\tau$ , $\scriptstyle \mu > 1$

输出: 模型参数 $\;\scriptstyle {{\beta }}$

1. 初始化 $\scriptstyle v_i^0 = 1\;(i = 1,2, \cdots ,n)$ 及 $\scriptstyle\tau$

2. 基于算法1更新β^m;

3. 基于式(10)更新v^m

4.  $\scriptstyle\tau \leftarrow \mu \tau $ ;

5. 令 $\scriptstyle m \leftarrow m+1$ , 重复步骤2–4, 直至:

　　　　 $\;\scriptstyle \sum\limits_{j=1}^p |\beta_j(m)-\beta_j(m-1)| \lt 10^{-6}$

3 数值结果 3.1 模拟数据分析

为检验本文所提出SSNL模型的预测表现, 我们首先按照以下方式模拟出一个简单的基因调控网络: 假设有200个转录因子(TFs), 每个转录因子调控10个基因,由此产生由2200个基因(节点)组成的生物基因调控网络, 转录因子之间以及与其调控的基因之间形成网络的边. 为了简单起见, 我们进一步假设模型中只有4个转录因子以及其调控的基因与响应变量y有关. 对于第一个模型, 我们按照以下方式来生成相关数据:

(1) $y = {{X\beta }} + {{\varepsilon }}$

(2) $\;{\bf{\beta }} = \left( {{5,\underbrace {\frac{{5}}{{\sqrt {10} }}\;, \ldots ,\frac{5}{{\sqrt {10} }}}_{10}}, - 5,\underbrace {\frac{{{\rm{ - }}5}}{{\sqrt {10} }}\;, \ldots ,\frac{{ - 5}}{{\sqrt {10} }}}_{10},\;3,} \right. $ $ {\left. \underbrace {\frac{{3}}{{\sqrt {10} }}\;, \ldots ,\frac{3}{{\sqrt {10} }}}_{10}{, - 3,\underbrace {\frac{{ - 3}}{{\sqrt {10} }}\;, \ldots ,\frac{{ - 3}}{{\sqrt {10} }}}_{10}} \right)^{\rm{T}}} \in {\mathbb{R}^{2200}}$

(3) ${\mathbf{\varepsilon }}$ 为误差项, 且 ${\varepsilon _i}\sim 5\times N(0,1).$

(4) 200个转录因子服从标准正态分布, 即 ${{{x}}_{T{F_j}}}\sim $ $ N(0,1),\;j = 1,2, \cdots ,200.$

(5) 每个TF与其调控的单个基因均服从二元正态分布, 且相关系数为 $\;\rho $ .

对于模型2, 我们假设

$\begin{split} {\bf{\beta }} =& (5,\underbrace {\frac{{{\rm{ - }}5}}{{\sqrt {10} }}\;, \ldots ,\frac{{ - 5}}{{\sqrt {10} }}}_3,\underbrace {\frac{5}{{\sqrt {10} }}\;, \ldots ,\frac{5}{{\sqrt {10} }}}_7, \\& - 5,\underbrace {\frac{5}{{\sqrt {10} }}\;, \ldots ,\frac{5}{{\sqrt {10} }}}_3,\underbrace {\frac{{ - 5}}{{\sqrt {10} }}\;, \ldots ,\frac{{ - 5}}{{\sqrt {10} }}}_7, \\& 3,\underbrace {\frac{{{\rm{ - }}3}}{{\sqrt {10} }}\;, \ldots ,\frac{{ - 3}}{{\sqrt {10} }}}_3,\underbrace {\frac{3}{{\sqrt {10} }}\;, \ldots ,\frac{3}{{\sqrt {10} }}}_7, \\& - 3,\underbrace {\frac{3}{{\sqrt {10} }}\;, \ldots ,\frac{3}{{\sqrt {10} }}}_3, \underbrace {\frac{{ - 3}}{{\sqrt {10} }}\;, \ldots ,\frac{{ - 3}}{{\sqrt {10} }}}_7{)^{\rm{T}}} \in {\mathbb{R}^{2200}} \end{split} $

其他设置与模型1完全一致. 该模型假设同一转录因子所调控的基因既可以对响应变量y产生正影响, 也可以对其产生负影响.

我们将模拟数据分为训练集和测试集, 其中训练集所占比例为70%. 在实际应用中, 人们通常使用k (k=3, 5或10)折交叉验证的方法来选择调优参数, 然而, 不同的k折交叉验证的方法很可能产生非常相似的预测结果^[33,34]. 此外, 研究结果表明将交叉验证的折数从10减少到3, 可以使算法的计算时间减少一半以上. 鉴于本文所提出的模型存在多个调优参数, 因此, 最终我们采用3折交叉验证的方法来选择最佳调优参数.

将基因相关系数 $\rho$ 分别设为0.2、0.5以及0.7. 每种情况都独立重复模拟50次, 并计算得到相应的均方误差(PMSE). 我们还进一步给出表征模型特征选择效果的两个指标, 分别是P和TP. 其中P表示模型中非零系数的个数, TP表示实际模型中非零系数的个数. 表1给出了各模型在不同情形下的模拟结果, 其中Lasso-Net表示Lasso和网络正则化; SCAD-Net表示不使用自主学习方法的SCAD和网络正则化.

如表1所示, 在不同模型不同相关系数下, 本文提出的SSNL方法相比于Lasso-Net和SCAD-Net方法均给出最小的PMSE. 此外, 在识别真正相关基因时, SSNL相比于其他两种方法具有最高的准确性(TP). 例如, 在 $\; \rho = 0.7$ 的情形下, 利用SSNL模型计算得到的TP值均超过40, 几乎达到模型的真实值44. 上述结果表明SSNL方法在处理高维度低样本、高噪声、高相关性的复杂数据集时具有良好的表现.

3.2 实际数据分析

为进一步论证SSNL模型的预测效果, 我们收集得到了乳腺癌细胞系数据集. 该数据集共有56个样本, 其中每个样本都隶属于一个确定的细胞亚型, 通过对其进行编码可以得到一个数值型响应变量. 此外, 每个乳腺细胞样本包含39653个基因, 并且这些基因之间存在交互关系. 通过加权基因共表达网络分析, 我们可以得到相应的加权网络. 然后将基因表达数据与该调控网络相结合, 得到最终的研究数据集. 我们旨在探索基因网络与关注的表型之间的关联关系以及网络中的核心基因.

我们将数据集随机打乱, 使约70%的样本成为训练样本, 剩余30%的样本作为测试样本. 类似于上文模拟中的情形, 我们采用3折交叉验证来估计得到最佳的调优参数. ${\lambda _1}$ 与 ${\lambda _{\rm{2}}}$ 的候选值均来自于{0.01: 0.1: 5}(起始值: 步长: 终值), $\; \mu $ 来自于{1.1: 0.1: 3}以及 $\tau $ 来自于{0.1: 0.05: 0.5}. 独立重复10次, 计算得到相应的均方误差(PMSE)以及模型中非零系数的个数P, 具体结果如表2所示.

从表2可以看出, 本文提出的SSNL方法给出了最小的PMSE, 其表现显著优于Lasso-Net方法, 且优于不使用自主学习的SCAD-Net方法. 此外, 在特征选择方面, 尽管3种方法的数值表现效果相当, 但SSNL方法仍优于其他两种对比方法. 上述结果再次说明本文所提出的SSNL模型在处理高维复杂网络数据集时具有良好的表现.

4 结论与展望

融合分析为基因组研究提供了一种有效的分析角度. 传统的融合分析方法是将多个数据集组合成一个集成的数据集, 然后直接对数据进行分析. 然而, 这种集成方法非但不能消除内部偏差, 甚至可能给融合数据集增加新的随机噪声和估计误差, 从而降低融合分析的统计功效. 本文提出了一种新的融合分析模型SSNL, 该模型融合了自主学习(SPL)和SCAD-Net正则化方法. 一方面, SPL方法能够先从低噪声样本中学习出一个基本模型, 然后通过高噪声样本学习使得模型更加稳健. 另一方面, 特征选择是SSNL模型的重要组成部分. SCAD罚函数是一种常见的特征选择方法, 但SCAD罚函数仅是从计算的角度出发, 没有利用任何先验信息. 故在已有研究的基础上, 本文给出了结合网络结构信息的SCAD-Net惩罚, 并对这一问题进行了一些理论探究, 包括群组效应和渐近性质. 不同情形下的模拟分析结果以及在乳腺癌细胞系数据集上的分析结果均表明, SSNL方法在处理高维复杂网络数据集时具有良好的预测表现.

本文使用3折交叉验证(CV)方法来选择SSNL模型中出现的惩罚参数. 然而, 当遇到多个超参数时, 使用CV方法进行网格搜索需要消耗大量的时间与内存. 最近, 一种进化计算(EC)方法被用来调整惩罚参数, 并且表现良好^[35]. 针对本文情形, EC方法可能是一个更好的选择. 此外, 我们还考虑将SPL+SCAD-Net方法拓展到其他回归模型中, 如广义线性回归等.