《1. 引言》

1. 引言

生物学和医学早已进入大数据时代,测序、蛋白质分析、代谢分析等分析方法和诊断技术的发展则加速了这一时代的到来。例如,相较于国际上对人类基因组计划[1]多年的努力和工作,近10年来下一代测序技术可以对个体基因组和宏基因组进行测序,甚至在单个实验室也是如此。微生物组的研究也得益于测序技术加速发展[2],它揭示了微生物群落在人类健康和疾病[3]等领域的重要性。这些发展产生了前所未有的海量高维数据,从而促进了人们对多元回归分析的研究兴趣[45]。多元回归旨在对一系列响应和一系列特征之间的关系进行建模,而普通回归通常描述的是一对一的关系[67]。响应变量(或因变量)是研究者希望能够解释的实验结果,预测变量(或自变量)是可能引起响应变化的受控输入。例如,在基因组学研究中,此类回归中的响应变量可能是人的性状,其特征可能是遗传基因或环境因素。因此,多元回归可以应用于我们日常生活的各个方面。例如,多元回归在经济学中被广泛应用于研究影响股票收益的因素[8]。其在生物学领域也很常见[9]。例如,将其应用于临床试验,以帮助研究者解释药物成分与农药效应之间的关系[10]。近年来,在基因组学研究中,包括宏基因组学研究,以及与代谢组学、蛋白质组学等相结合的研究中,多元回归在理解重要性状的关联和潜在因果关系方面发挥了重要作用[11]。

人们在运用多元方法应对具体挑战方面进行了各种尝试。主成分分析法(PCA)是最古老、最著名的基于特征向量的多变量分析技术之一[12]。当变量个数较多时,它被广泛用于利用正交变换找到描述方差最大的变量的线性组合。通过将数据投影到低维空间可以显示其主导梯度,PCA可以揭示数据的内部结构[1314]。在实际应用中,主成分回归法(PCR)是一种利用PCA估计回归系数矩阵的线性回归模型。典范对应分析(CCA)是另一种常用的方法,它通过降维来解释两组变量之间的关系[15]。它旨在找到一个能够描述预测变量和响应变量之间最大相关性的线性组合。另一种常用于寻找两个矩阵之间关系的方法是偏最小二乘(PLS)回归法[16]。通过将响应变量和预测变量投影到新的空间中构建线性回归模型,从而总结出协方差结构。虽然上述方法在研究中得到广泛的应用,但其存在三个主要的统计学问题。第一个问题是传统方法往往忽略了观测数据和响应变量之间可能存在的相互关系。第二个问题是,有些方法不允许变量选择,但这在预测变量数量较大的探索性实验中是必不可少的。第三,一些真实数据库往往变量总数大并且样本量小,导致出现不可靠的解决方案[1718]。

基于这些考虑,我们分析了一类新的方法[降秩回归(RRR)及其扩展],这类方法通过考虑响应变量之间的关联来提高回归模型的可解释性[1920]。RRR是一种类似于PCA的数据缩减方法,它创建新的变量来总结原始数据中的大量信息[21]。尤其是,它通过定义一组预测变量的线性组合,以最好地解释响应变量中的总方差,并且具有简单、计算效率高、预测性能好等许多可取的优点。当预测变量数量较多时,如何选择重要变量是研究者感兴趣的另一个问题。尽管一些传统的方法仍然可用——如随机森林,它通过构建大量的决策树来进行预测——但它们更适用于响应变量单一的情况。因此,一些基于RRR的方法采用了群体选择方法的思想,包括群最小绝对值收敛和选择算子(群lasso)方法等,其用于在确定变量间的群体结构时用于选择变量。在营养流行病学和遗传学方面,许多研究使用了RRR方法;但对不同基于RRR方法的性质以及其对真实数据(尤其是以宏基因组为中心数据)的适用性的全面分析还有待进行[2225]。在此,我们使用不同维度的仿真数据,比较了基于RRR的各种方法与其他性质相似的多元回归方法的性能,考察了它们在不同场景下的优势和局限性,并将其应用于大规模公共宏基因组数据集。

《2. 方法》

2. 方法

《2.1. 本研究中的测试方法》

2.1. 本研究中的测试方法

在本研究中,我们以基本的多元线性回归方法为出发点,然后与一些基于RRR和其他多元回归的方法进行比较。每一种方法的定义和原理解释如下。

《2.1.1. 多元线性回归》

2.1.1. 多元线性回归

多元线性回归模型由多个预测变量X1, X2, …, X和多个响应变量Y1, Y2, …, Y组成。每个响应变量由预测因子的线性回归来表示,即

Yj=k=1pXkckj+ϵj, j=1, 2,, q(1)

式中,c是与XY相关的回归系数;是均值为零的误差项。

另外,这个公式可以用n个观测值改写如下:

Y=XC+E(2)

式中,Xn×p的预测矩阵;Yn×q的响应矩阵;Cp×q的回归系数矩阵;En×q的误差矩阵。

我们基于最小二乘准则估计出系数矩阵C,即

minC Y-XC2(3)

式中,·表示弗罗贝尼乌斯范数。

利用普通最小二乘法(OLS)得到C, Ĉ的估计值,计算公式如下:

Ĉ =XTX-1XTY(4)

《2.1.2. 降秩回归》

2.1.2. 降秩回归

然而,OLS方法忽略了响应变量之间可能存在的相互关系,而只对每个响应变量单独进行估计,因此它只能提供一个粗略的估计。因此,我们引入了RRR方法,它限制了系数矩阵C的秩。我们假设C的秩较低,r = rank(C) ≤ min(p,q),C可以表示为两个秩为r的矩阵乘积,C = BAT,其中,BA的维数分别是p×rq×r。多元回归模型(2)可以改写为:

Y=XBAT+E(5)

此外,XB有n×r维数,表示X的一组r个线性组合,可以解释为驱动Y变化的潜在因素,因此,RRR有助于降低预测变量的维度并且提高计算效率。

我们可以重写优化函数(3),即

minA, B  Y-XBAT2(6)

ÂB̂解的集合如下形式:

Â=V(7)

B̂=ΣXX-1ΣXYV(8)

式中,ΣXX= (1/n)XTX, ΣXY= (1/n)XTY, ΣYX= (1/n)YTXV表示ΣYXΣXX-1ΣXY所对应特征值的特征向量[26]。

《2.1.3. 稀疏降秩回归》

2.1.3. 稀疏降秩回归

稀疏降秩回归(SRRR)是一种广泛的RRR方法,该方法不仅关注降维,而且关注变量选择[2627]。它通过在最小二乘法估计中加入惩罚项来增加系数矩阵的稀疏性,从而具有独特的性质。与RRR利用所有的预测变量来构建潜在因子相比,SRRR可以从大量变量中选择有用变量,并通过引入群lasso惩罚来排除冗余变量[28]。因此,优化公式(6)可以改写如下:

minA, B  Y-XBAT2+j=1pλiBi s.t. ATA=I(9)

式中,λ是惩罚参数。该约束ATA = I用于满足可辨识条件,其中,I表示单位矩阵。此外,如果将Bi设为零,则整个第i行的矩阵B将为零,第i个预测变量将不再活跃。

利用次梯度法或变分法可以求解该优化问题,但通过交叉验证(CV)定义p惩罚参数(λ)可能会消耗时间。为了减少调优参数的数量,通常使用如下两种策略[26]。

(1)群lasso惩罚:设置所有λ值为λ

(2)自适应加权lasso惩罚:根据原始数据结构λi= 1/C˜iγ·λ计算每一个λ,这里C˜C的根-n一致估计,γ是正整数[29]。

《2.1.4. 具有行稀疏性的子空间辅助回归》

2.1.4. 具有行稀疏性的子空间辅助回归

具有行稀疏性的子空间辅助回归(SARRS)也侧重于解决系数矩阵中低秩性和稀疏性的问题[30]。这种新的估计方法可以用于自适应稀疏降秩多元回归,并达到降维和变量选择的目的。此外,与前面讨论的SRRR相比,当预测变量个数超过样本容量时,SARRS能够改善其性能。

在利用群体稀疏性优化回归的过程中,可以使用两个惩罚函数。

(1)群lasso惩罚:ρB;λ= λ B,其中,λ为惩罚参数,B为待优化的参数矩阵。

(2)群极大极小凹惩罚(MCP):ρB;λ=λ·0B1-t/γλ+dt,其中,γ是大于1的正整数,(1-t/γ)+表示其正部分,即1-t/γλ+= 1-t/γλ·11-t/γλ0 [31]。

《2.1.5. 稀疏偏最小二乘回归》

2.1.5. 稀疏偏最小二乘回归

稀疏偏最小二乘回归(SPLS)方法是在PLS的基础上,进一步增强预测空间多维方向上的稀疏性,从而实现变量选择[17]。它首先选取与响应具有较强相关性的预测变量,然后再加入具有较强偏相关性的预测变量。与SRRR相比,SPLS采用了不同的降秩结构,并且没有直接聚焦于响应变量的预测,从而形成了其预测方面的弱点。

《2.1.6. 识别主预测变量的正则化多元回归》

2.1.6. 识别主预测变量的正则化多元回归

识别主预测变量的正则化多元回归(REmMap)方法与上述方法不同,它不仅假设只有部分预测变量与响应相关,而且这些预测变量可能只影响部分响应[32]。这一做法是合理的,因为在实际情况中,研究者往往比其他人更关注具体的响应。REmMap能够拟合高维和小样本规模的多元回归模型,并且能够在系数矩阵中同时引入总体稀疏性和群体稀疏性来检测主预测因子。

《2.1.7. 总结》

2.1.7. 总结

表1总结了本文所讨论方法的特点。

《表1》

表1 方法的比较

MethodsData reduction method (low rankness)Variable selection (sparsity)Explains interrelation between responsesFeatures
RRRRestricts the rank of the regression coefficient matrix
SRRRUses a group lasso penalty to allow row-wise sparsity
SARRSSuitable when the number of predictors exceeds the sample size
SPLSUses PLS to impose sparsity
REmMapEach response has different relevant predictors
PCRProjects the predictors into a lower-dimensional space
Group lassoEnables variable selection considering group structure
Random forestUsed to rank the importance of predictor variables

《2.2. 基于仿真数据的方法性能测试》

2.2. 基于仿真数据的方法性能测试

《2.2.1. 仿真设置》

2.2.1. 仿真设置

为了解释和比较SRRR、SARRS、SPLS、REmMap以及一些传统方法(PCR、群lasso和随机森林)的性能,我们首先引入一个仿真研究,利用上述方法生成数据并进行分析。我们采用与Chen和Huang [26]类似的仿真设置。仿真研究的核心思想是分析一些与响应变量相关的预测变量和一些不相关的预测变量。然后,我们使用这些方法来检验其中哪种方法能够最准确地判定其关系以及取得良好的预测性能。

我们用多元线性方程Y=XBAT+E来生成数据。在该模型中,n×p的设计矩阵X服从多元正态分布N0, ΣX,其中,ΣX的对角元素为1,非对角元素为ρx。矩阵B和矩阵A构成了模型的系数矩阵。在p×r的矩阵B中,第一行p0服从N(0,1),其余p-p0行为零。q×r的矩阵AN(0,1)生成。矩阵E是由N0, σ2ΣE定义的随机噪声矩阵,其中,σ2是噪声的大小,ΣE表示对角线元素为1和非对角线元素为ρ的矩阵。然后,由上述模型计算出n×q的矩阵Y

我们生成三组数据,包括训练集、验证集和测试集。训练集用于拟合基于各种方法的模型。验证集用于调整模型内部的参数和估计噪声方差。最后,我们使用测试集的数据来评估所构建模型的性能。

为探索在不同情况下这些方法的适用性,我们采用几种不同的情况来进行仿真研究。首先,由于研究人员有时很难获得足够的样本进行试验,因此我们希望同时在小样本容量和大样本容量的情况下测试这些方法的性能;其次,我们还对变量数量对实验造成的影响感兴趣。真实数据如微生物学数据和遗传学数据往往包括高维预测变量或响应变量。基于以上考虑,我们仿真了以下6种情况,其中,n为数据的样本容量,pq分别为XY中变量的个数。

案例1:小样本容量,n < p

案例1a:n = 20;p = 100;q = 25

案例1b:n = 20;p = 25;q = 25

案例1c:n = 20;p = 25;q = 100

案例2:大样本容量,n > p

案例2a:n = 200;p = 100;q = 25

案例2b:n = 200;p = 25;q = 25

案例2c:n = 200;p = 25;q = 100

在R中使用spls、rrpack、remMap、pls、glmnet和randomForest软件包对所讨论的仿真过程和方法进行了编码。SARRS方法的代码由各软件包的作者提供。附录A中指定和列出了计算过程。

《2.2.2. 各种方法的评价准则》

2.2.2. 各种方法的评价准则

在每种情况下,我们重复了仿真过程500次,使用以下指标来衡量和比较上述多元回归方法的性能。

(1)均方误差(MSE):MSE用来表示这些方法的预测精度,具体定义如下:

MSE= 1ni=1n(Ŷi-Yi)2(10)

式中,ŶiY的预测值。

(2)该指数是指可由预测变量解释的响应变量中方差所占的比例。它通常用来描述模型与数据的拟合程度。R2数值越大表明模型的拟合程度越好。

(3)泰尔不平等系数(TIC):TIC是另一个用来反映拟合值与真实值差异的指标。它的定义如下:

TIC=1ni=1n(Ŷi-Yi)21ni=1nŶi2+1ni=1nYi2(11)

TIC取值范围为0 ~ 1,TIC越小,模型的预测精度越高。

(4)敏感性(TPR)和特异性(SPC):TPR和SPC常用来评价变量选择的准确性。TPR是选择真正相关变量的能力,由正确选择次数与相关输入变量总数的比值计算得出。SPC是选择真正无关变量的能力,它由正确选择次数与不相关输入变量总数的比值计算得出。一种方法如果既具有高TPR又具有高SPC就意味着它可以准确地选取出相关变量。

(5)曲线下的面积(AUC):AUC也用于衡量正确选择真实相关变量的比率[33]。

(6)总体评分:采用上述评价指标计算总体评分指数。性能较好的方法整体得分较高。

《2.3. 应用于真实的种群级微生物数据》

2.3. 应用于真实的种群级微生物数据

为了说明上述方法对于现实世界问题的实用性,我们将其应用于Falony等[34]工作的比利时佛兰芒肠道菌群项目(FGFP;发现群体:n = 1106)的数据。本研究旨在发现微生物菌群变异与宿主特征、地理和药物摄入等环境因素之间的关系。讨论了将66个基于临床和问卷的变量作为预测变量,74个微生物群作为选择后的响应变量。

我们通过CV将数据随机分成训练集和测试集。我们还建立了一个模型来拟合数据,并通过上述指标来评价其性能。我们将这个过程重复50次来检验变量选择的稳定性。

《3. 结果》

3. 结果

《3.1. 仿真研究揭示了每种方法的不同性质》

3.1. 仿真研究揭示了每种方法的不同性质

在仿真中,我们对不同的案例应用SRRR(带有群lasso惩罚和自适应加权群lasso惩罚)、SARRS(带有群lasso惩罚和群MCP惩罚)、SPLS、REmMap、PCR、群lasso和随机森林,并使用CV来调整每种方法的低秩参数。它们的总体表现如图1所示。

《图1》

图1 所有方法的总体评价,以热图的形式显示。x轴表示不同的案例,y轴表示不同的方法。每个单元格的颜色代表相应的总体评分。总体评分越高,表现越好。glasso:群lasso惩罚;adglasso:自适应加权群lasso惩罚;gMCP:群MCP惩罚。

热图显示,所有的方法在案例1的性能都比案例2差,这与我们的预测一致。此外,很明显,在所有案例中都最为适用的是SARRS(带有群MCP惩罚),当样本容量增多时,SRRR(带有自适应加权群lasso惩罚)和SPLS同样适用,且性能都很好。

每种方法的效果都是通过上述标准来衡量的,案例1的结果如图2所示。

《图2》

图2 对案例1a、1b和1c中所有方法的效果评估,以雷达图形式显示。圆心为0。对于R2、TIC、TPR、SPC和AUC,圆周为1。因此,如果一种方法在响应矩阵中的R2、TPR、SPC和AUC较高,TIC和MSE较低,我们认为该方法效果较好。NO.Var指标给出了每种方法选择的变量数目,其中圆心表示0,边表示预测变量的数量。

在案例1a中,样本容量非常小,预测变量数目大于样本容量;因此,大多数方法都没有很好的预测和变量选择效果。除PCR无法选择相关变量外,各方法的SPC约为0.75,TPR约为0.55,表明选择不足。然而,与传统方法(PCR、群lasso和随机森林)相比,本文讨论的新方法都具有更好的效果,特别是带有群MCP惩罚的SARRS方法。该方法的MSE和TIC最低,R2、SPC和AUC最高。这一结果与方法论部分的论述一致,其中特别强调了当预测变量数目远远大于样本容量时,SARRS是最合适且最准确的方法。

在案例1b和1c中,预测变量数目更接近样本容量。我们发现,由于R2较高和TIC较低,所有模型对仿真数据的拟合效果都优于案例1a。从图中还可以看出SRRR在预测精度方面的优势,因为SPLS和REmMap相比SRRR和SARRS具有更大的MSE。此外,在变量选择方面,我们可以看到SARRS(带有群MCP惩罚)、SRRR(带有自适应加权群lasso惩罚)和SPLS的效果较好,其SPC值要高得多,表明其在选择真正的相关变量和避免过度选择之间取得了平衡。

在案例2中,我们研究了大样本容量的情况;如图3所示,很明显所有的方法都比案例1的效果更好。

《图3》

图3 对案例2a、2b和2c中所有方法的效果评估,以雷达图形式显示。圆心为0。对于R2、TIC、TPR、SPC和AUC,圆周为1。因此,如果一种方法在响应矩阵中的R2、TPR、SPC和AUC较高,TIC和MSE较低,我们认为该方法效果较好。NO.Var指标给出了每种方法选择的变量数目,其中圆心表示0,边表示预测变量的数量。

我们首先讨论案例2a,在这种情况下,预测变量远多于响应变量。在预测效果方面,我们感兴趣的新方法在系数矩阵中的MSE都非常低,甚至接近0,平均TIC约为0.27,平均R2约为0.72,表明模型拟合良好。此外,传统方法在这方面的表现仍然不佳。在变量选择方面,所有方法的TPR和AUC均高于案例1。SARRS(带有群MCP惩罚)、SRRR(带有自适应加权群lasso惩罚)和SPLS也具有极高的SPC,说明它们能够准确地选择所有相关变量,剔除所有无关变量。

在案例2b中,我们将预测变量的数目减少到与响应变量相同,从而增大了样本容量与变量数目的差值。在此情况下,所有新方法的预测效果都很好。然而,在变量选择方面,这些方法的结果呈现两极化。表现最好的仍然是SRRR(带有自适应加权群lasso惩罚),其次是SARRS(带有群MCP惩罚)和SPLS,其选择精度接近1。而对于其他方法,SPC普遍较低,说明存在过度选择问题。由于案例1c类似于案例1b,因此案例2c也类似于案例2b。但是,与案例1b和1c相比,样本容量的增加使得案例2b和2c在预测精度和变量选择方面表现得更好。

《3.2. 应用于实际种群水平的微生物组数据》

3.2. 应用于实际种群水平的微生物组数据

案例研究的数据特征与上述仿真研究中的案例2b一致,其中样本容量(n = 1106)远大于变量数目(p = 66,q = 74)。基于之前的讨论,我们知道最适合应用于该数据集的方法是SRRR(带有自适应加权群lasso惩罚)。因此,我们构建了一个SRRR模型来分析环境指标与细菌组成之间的关系,并讨论其预测精度和变量选择结果。结果如图4所示。

《图4》

图4 对SRRR(带有自适应加权群lasso惩罚)应用于实际种群水平微生物组数据的效果评估。(a)使用TIC指数来展示SRRR模型对每个响应变量的预测效果的条形图。图中显示了20个TIC值最低的变量。(b)显示每个预测变量在50次交叉验证中选择的百分比图。

模型的平均TIC为0.56,高于在仿真研究中案例2b的0.26。但由于我们知道真实的数据比仿真的数据噪声更大,因此我们认为,虽然这一TIC是可以接受的,但就充分的模型预测而言其并不令人信服。然而,在对每个响应变量进行更仔细的检查后,我们发现TIC较低的变量在以往许多研究报告中就出现过,包括FGFP [34],如Faecalibacterium(TIC为0.24)、Blautia(0.32)、Bacteroides(0.33)、Roseburia(0.35)和Ruminococcus(0.40)。这些是生产丁酸的关键细菌,它们的低丰度与许多疾病有关,因为微生物组产生的低丁酸会引发较高水平的炎症和代谢紊乱[3536]。因此,SRRR选择的预测因子可以很好地解释这些变量,并且也可以通过SRRR计算出的系数矩阵来预测。

最后,我们检验了变量选择的鲁棒性。由于我们将CV流程重复了50次,因此在80%以上的案例中选择的预测变量是最有意义的。图4(b)显示了最常选择的34个变量,包括性别、吸烟者、红细胞计数(RBC)、肌酸酐、粪便评分、平均红细胞血红蛋白浓度(MCHC)和多种药物。这一结果与先前一项研究中所讨论的药物作用的重要性一致[37]。

《4. 结论》

4. 结论

随着各个研究领域的数据量和维度的增长,生物医学研究仍然是最重要和发展最快的领域之一。当从数据中提取最大值时,在不同度量或组学水平之间获取正确、有用和有意义的关联是一项重大挑战[38]。在此,我们研究了一些有代表性的方法,包括RRR和其他多元回归方法的扩展形式,并使用仿真数据和以微生物组为中心的真实数据来说明这些方法的优势和局限性,这可能对微生物组和其他相关组学数据的未来应用具有指导意义。

我们总共研究了9种方法-参数组合,其中包括7种方法;此外,对其中两种组合使用了两种不同惩罚。我们仿真了不同样本容量/维度的数据,并比较了在维度上存在/不存在较大差异的预测变量和响应变量。从案例1和案例2的对比结果可以看出大样本容量的重要性,这将大大提升所有方法的效果。特别地,与案例1相比,SPLS在案例2中具有更好的预测精度,表明其在样本量较大时更适用。在类似于案例1中小样本容量的情况下,最好的方法是带有群MCP惩罚的SARRS,该方法在预测和变量选择方面都有优异的性能。当样本量较大时,如案例2,SARRS(带有组MCP惩罚)、SRRR(带有自适应加权群lasso惩罚)和SPLS均表现良好;通过进一步观察发现,SRRR(带有自适应加权群lasso惩罚)的效果略好于其他两种方法。

我们利用公布的FGFP数据来为一个真实场景选择最好的方法,从中可获得微生物组数据和研究中确定的环境因素。由于两类变量的维度大致相似,我们决定使用带有自适应加权群lasso惩罚的SRRR。首先,我们确定了环境变化最能解释的菌群。所确定的菌群证实了先前的论断,即产生丁酸盐的细菌对人体健康非常重要,且有可能与这些环境因素有关。另外,由于环境因素被认为是预测变量(即选定的与菌群有关的特征),我们还设法重现了已发表的研究中最重要的特征,这再次证明了所选方法的可靠性和鲁棒性。

研究人员在拟合模型时应谨慎选择合适的惩罚。例如,在SRRR方法中,当n > p时,自适应加权群lasso惩罚提高了预测精度和变量选择。当n < p时,与未加权组相比,自适应加权组的TPR较低,但SPC较高。这可以解释为,当我们在SRRR计算过程中引入权重时,之前被过滤掉的变量在其惩罚项中具有较大的权重,并且不再包含在模型中。因此,带有未加权惩罚的SRRR会选择更多的变量,导致SPC较高。

综上所述,我们检验了几种多元回归方法的适用性,并检测了它们在不同的组学情景下的效果,而在现实中,这些情景可能在样本量和维度上存在巨大差异。基于此,我们能够推荐最佳方法。诚然,我们的初步分析在现阶段无法进一步扩展,无法将不同指标(如物种)之间的系统信息纳入多组学数据中,因为这需要关于这些指标之间连通性和相似性的先验信息。我们还使用了一个著名的微生物组数据集,证明我们的选择方法可以在很大程度上概括由单变量分析获得的结果,并促进了在变量和组合特征选择方面的思考。这些发现将有助于在未来更大规模的组学研究中方法的选择,包括以微生物组为中心的研究。