数据统计分析方法

无论2D还是3D方法，多是通过各自的方式计算出不同的参数，而用这些参数建立模型还需要相应的统计方法。如Co MFA法采用PLS、Compass采用神经网络、GERM采用遗传算法等。

6.5.1 主成分分析与偏最小二乘回归

由于分子参数和生物活性的真实关系未知，因此许多时候采用最简单的线性回归形式，其中最常用的回归统计方法是多元线性回归（multip lelinear regres-sion，MLR）。但是通常参数与生物活性数据之间不会恰好呈线性关系，其相关性系数常常较低。为了提高相关性系数，经常需要增加用于建模的参数个数，参数用得越多，模型的拟合能力越好。但是单纯增加参数个数提高模型拟合往往造成其推广能力的下降，即拟合能力很好而预测能力很差的情况，称之为过拟合（overfit-ting）。因此不能不加限制地增加建模所用的参数个数，有人认为建模的样本数目至少应是建模所用参数个数的3至5倍，也有人提出该比例应大于2n，n为选取的参数个数。

在QSAR中评价一个模型的好坏有一些统计学标准。一个QSAR模型的拟合能力有以下几个评价标准：相关性系数r，标准偏差s，Fisher检验值F。r越接近于1，s越小，F越大，则模型的拟合能力越好。

上面三式中：Yexp为生物活性实验值；Ycal为生物活性的模型拟合值；n为样本数；k为建模所用的参数个数。

其次是模型的稳定性，即模型是否会因为少量样本的变化而剧烈变化。评价一个模型的稳定性的标准是交叉验证相关性系数（rcv）。交叉验证就是每次从所有样本中抽出k个，用剩余的样本建模预测被抽取的k个样本的活性，重复这一步骤直至所有样本都被抽取和预测过一遍，最后以类似于求取相关系数r的方法求交叉相关系数rcv。特别地，当k=1时，这种方法称为留一法交叉验证（leave-one-out，LOO），也叫抽一法交叉验证。

式中：Ypred表示交叉验证的预测值。

最后一个方面是模型的预测能力，也称泛化能力或推广能力，这也是定量构效关系模型最根本的指标。所有的样本分成训练集或预测集，如果可能的话预测集应当相比训练集更大，但是由于样本的缺乏，在实际操作中常常相反，用2/3左右的样本作为训练集。预测能力的评价指标是预测集的相关性系数。

采用多元线性回归的方法快速简便，不过需要各个参数彼此间线性无关，但是随着各种参数计算软件能够得到的分子参数越来越多，它们彼此间的相关性也越来越高，为了解决这一问题，需要采用主成分分析的方法。

主成分分析方法（principle component analysis）中的主成分是原来变量的某种线性组合，它们能够最大限度地解释原变量的方差，同时彼此之间又保持线性无关。主成分分析的第一步就是将原始参数化为主成分，使原始参数降维（主成分数常常远远少于原参数个数）；然后再用主成分代替原始参数建立线性回归模型。例如自变量是X，因变量是Y，X中有些变量严重线性相关，那么首先计算得分因子矩阵T=XW，并用回归因子矩阵作回归模型：Y=TQ+E（E是噪声）。但是将原始参数转化为主成分时不考虑生物活性数据的因素，有可能保留噪声而去除有用的主成分，因此更常用的方法是偏最小二乘回归。

相比于主成分分析，偏最小二乘法（partialleastsquare，PLS）在考虑自变量的同时也折中了因变量的作用，是现在进行线性回归的主流方法。主成分回归产生的权重矩阵W反映的是预测变量X之间的协方差，偏最小二乘回归产生的权重矩阵W反映的是自变量X与因变量Y之间的协方差。

6.5.2 人工神经网络

人工神经网络（artificial neural network，ANN）是模仿人类神经元及人脑神经网络结构而建立的信息处理系统。神经网络主要有以下优点。

（1）非线性。一个神经元本身可以是线性的，也可以是非线性的。有非线性的神经元参与组成的神经网络也是非线性的，并且可以通过神经元之间的相互制约形成输入和输出之间复杂的非线性映射。当我们进行QSAR分析时，很难指望在所有体系中参数和生物活性数据之间的内部物理机制能够恰好符合线性关系，如果线性较差或者根本不存在线性关系，那么采用传统的线性回归方法无论如何都无法得到理想的模型。因此对于QSAR分析而言，神经网络的这种非线性性质是很重要的。

（2）并行性。人工神经网络通常包含若干层，某些层中的神经元并行地排列，虽然每个单个神经元的功能非常简单，但是由于这种并行结构的存在，整个系统处理数据的能力大大提高。

（3）适应性。每个神经元均有权值与阈值，当神经元的输入发生变化时，可以通过调整权值阈值进行适应，有的网络还能够调整网络的结构。这种特性使得神经网络具有强大的学习能力。

（4）容错性。人工神经网络通过自身网络结构将记忆的信息分散保存在所有神经元之间的权值中，而不是某几个特定的权值中，这种分布式的储存方式使得人工神经网络具有良好的容错性。当某个神经元损坏时，整个人工神经网络可能逐渐恶化，但不会立即崩溃。当我们要将一个神经网络以硬件形式实现时，这一优点便得以充分体现。

鉴于以上优点，以及神经网络善于处理因果关系不明确问题的能力，神经网络经常用于线性回归无法处理的复杂的QSAR问题的分析。

但是神经网络也存在一些问题。首先，神经网络的非线性拟合能力很强，但是相比于传统线性方法更容易出现过拟合，当训练集较小时这种现象尤其明显；其次，神经网络最终得到的输入与输出的关系相当复杂，不如传统线性模型那样简单明了，易于解释，所以神经网络也就难于指导活性更好的药物的设计；最后，神经网络的层数和节点个数如何选取，神经元选取何种非线性函数，也都没有普遍适用的方法来确定。

人工神经网络的种类多种多样，在2D-QSAR中采用的通常是三层误差反向传播网络（三层BP网，back-propagation neuralnet work），因为根据Kolmogorov定理，三层BP神经网络已经能够以任意精度逼近任意函数。BP网络分输入层、隐藏层和输出层三层。隐藏层的层数和节点数均可调，在三层BP网中，隐藏层仅一层。此外，偶尔也有将自组织人工神经网络、连续型连接网络等应用于QSAR的尝试。

神经网络的基本单元是神经元，每个神经元可以接受多个输入信号，经过某种非线性映射（对于BP网络而言，通常是某种可微的单调递增函数，如对数函数或正切函数），再加上阈值的影响，得到（通常是单个）输出信号。

虽然每个神经元的功能十分简单，然而联结成为网络之后便可以产生十分复杂的模型。三层BP网的结构如图6-18所示。

图6-18 三层BP神经网络

输入层神经元的个数等于参数的个数，每个输入的维数等于样本的个数；隐含层神经元的个数须自行确定；由于通常只须研究药物的某一种生物活性数据，因此输出层一般只设一个神经元。

Aoyama率先将人工神经网络引入到QSAR研究中。1989年和1990年， Aoyama分别在降冰片烯衍生物分类和抗肿瘤药物卡波醌衍生物QSAR中采用了三层BP网进行分析，均取得了成功。

6.5.3 支持向量机

支持向量机有时也被归为人工神经网络中的一种，不过与之前QSAR中最常采用的三层BP网有明显的差别，因此在这里单独介绍。

支持向量机出现的时间较晚，由Cortes和Vapnik于1995年提出，在2000年之后才运用于QSAR。SVM方法的主要优点是：①需要的样本较少。虽然无论任何方法，样本总是多多益善，但是处理相同复杂程度的问题，SVM方法需要的样本比其他方法要少。②非线性。SVM方法可以通过松弛变量和核函数处理线性不可分的分类问题。③高维分类及回归。由于SVM方法只用到少部分称为支持向量的样本，因此就能节省相当多的计算和存储资源，从而可以处理相对更高维的样本。SVM方法的主要功能是分类，不过也能够推广应用到函数拟合中。

SVM的原理首先要从SVM分类说起。

SVM算法的目标是最小化结构风险。结构风险包括经验风险和置信风险（见图6-19）。经验风险是模型在训练集样本上的误差，我们在分类或回归时会使这个值尽量地小，而置信风险则是模型的可靠程度，由两个因素决定：样本数量越多，分类函数的VC维越小，则置信风险越低，模型越可靠。SVM的基本思想就是使VC维最小化。

图6-19 各类风险

所谓VC维是函数容量（复杂程度）的一种度量。假设有p个点，显然这p个点若分别属于正类和负类，就有2p种可能的情况。假设有一个最大的p，使每一种情况的正类和负类都能被某一函数集分开，则称这个可被散列的点的最大数量为该函数集的最大VC维。图6-20给出了一个简单的例子，以便于理解。

图6-20中，函数集是线性函数，那么函数集可分开平面上的3点，但无法分开4点，因此在这个例子中VC维就等于3。

图6-20 VC维示意（引自Steve R.Gunn）

首先考虑样本线性可分的情况。图6-21中的任意一个分类超平面都可以将正类与负类样本正确地分开，但是SVM方法则需要找出其中的最佳分类超平面（深色线）。

图6-21 最佳分类面（引自Steve R.Gunn）

最佳分类超平面的标准是，距离分类超平面最近的样本点（也就是所谓的“支持向量”），距离分类超平面的距离最远。显然，对于最佳分类超平面，正类和负类的支持向量距离该分类超平面的距离相等。

之所以要这样设立分类超平面，是因为样本的错误分类次数的上限与支持向量到超平面的距离的平方成反比（见图6-22）。

图6-22 支持向量与边界距离（引自Burges）

但是，有的时候样本并非线性可分的，SVM方法通过核函数和松弛变量来解决这一问题。

对于低维空间的线性不可分的情况，如果我们将其映射到高维空间，便可能变得线性可分，而SVM需要的则是高维空间中的内积值。所谓核（kernel）函数就是这样一些函数，它们能够接受低维空间输入，并得出高维空间的内积。比较常用的核函数有以下几个：

线性核函数：K（x，xi）=（x·xi）

多项式核函数：K（x，xi）=[（x·xi）+1]q

S型核函数：K（x，xi）=tanh[v（x·xi）+c]

其中高斯核函数（径向基核函数）是相对较为通用的核函数，不过总的来说，对何种样本该采用何种核函数，并没有绝对的标准。

但是有的时候，即使采用核函数也无法完全将正负类分开（见图6-23），因此SVM还同时采用了另一种处理方法：松弛变量。

图6-23 无法完全分类的情况（引自Steve R.Gunn）

改为：

式中：ξi就是松弛变量；而C称为惩罚因子，是一个需要提前设定的量。不同样本的C可以设为不同，如果C很大，则说明这个样本很重要，不希望它被错分。如果正负类样本数量严重失衡，有时候可以通过给予样本少的类较大的惩罚因子来弥补。

SVM方法同样可以应用到回归问题方面，因此才可以用于QSAR研究。用SVM方法处理回归问题必须要引入损失函数。损失函数包括：二次函数、Laplace函数、Huber函数、ε-不敏感函数。通常SVM方法采用ε-不敏感函数，此时SVM回归的解的展开形式采用的支持向量较少。

图6-24 SVR原理与ε-不敏感损失函数（引自Pao-Shan Yu等）

SVM回归（SVR）的方程形式上和SVM分类（SVC）十分接近，不过除了C之外，SVM回归还要预先设置ε参数：

假设我们要构造一个函数g（x），使它以ε为精度逼近函数f（x），那么实际上我们就要使f（x）处于以g（x）为中心的ε管道内，由于管壁趋向于平缓，因此管壁将和f（x）上的某些点重合，这些点就是支持向量。ε管道越宽，支持向量的个数便越少。我们通过控制C和ε参数来调控SVM的推广能力。

6.5.4 遗传算法

遗传算法（genetic algorithm，GA）是一种模拟自然选择和遗传变异等生物进化过程求解极值问题的全局性概率搜索算法，在形式上是一种迭代的方法。所以遗传算法通常不单独使用，而要结合PLS，ANN，SVM等方法运用。在遗传算法当中，子代的个体能够保留父带的大部分模式，同时还有一定的概率得到更好的模式。

遗传算法最早起源于20世纪60年代初期。1975年，美国密歇根大学的Hol-land首先系统阐述了遗传算法的基本理论和方法，并且提出了模式（schema）理论。同年De Jong提出了衡量遗传算法性能的在线（online）和离线（off-line）指标，并使用了5种试验函数，沿用至今。

遗传算法在2D-QSAR中的主要作用是进行用于建模的分子参数的筛选。当可供选取的分子参数个数较少时，可以采用完全枚举的方法。但是随着目前各类相关软件功能的完善，可以计算的分子参数的数量日益增长（例如DRAGON2007版可计算的分子参数为3224个，那就将是23224-1种参数组合可用于建模，显然不能采用穷举的方法）。使用适当的方法从数量庞大的候选参数中选取合适的参数建模是非常重要的。而遗传算法正是一种相对容错性强、简便、并且相对快速的算法，1994年，Roger和Kubinyl等率先分别将遗传算法用于2D-QSAR领域。

6.5.4.1 基本流程

遗传算法的基本流程如图6-25所示。

图6-25 遗传算法的一般流程

首先，将随机产生的初始参数组合编码成某种字符串的格式。随后，根据一定格式计算出其适应度，随后通过选择、交叉、逆转（很少采用）、突变等步骤产生下一代个体，重复以上过程，直至达到某种终止条件（如：种群适应度达到某个标准，种群适应度不再改变，或单纯地设定代数上限）。

6.5.4.2 遗传编码

遗传算法的遗传编码方式多种多样，如二进制编码、多进制数编码、序列编码、实数编码、树编码、自适应编码、乱续编码、二倍体编码等。在2D-QSAR中，通常采用二进制编码或整数编码。

在用二进制数列表示建模所用参量时，有两种常用方法。第一种是使每个个体的长度等于所有候选参量的个数，“1”表示该位置对应的参数参与建模，“0”表示不参与建模，这种编码方式最终建模的参量个数不固定，由于遗传算法的适应性函数常直接采取偏最小二乘的相关性系数，采用这种方式编码容易造成过拟合，需要根据建模的参数个数调整个体的适应值；第二种方法是固定建模采用的长度，而个体序列的每一位数字所代表的是建模所选参数的序号。例如候选参数个数为10个，编码方法2的参数个数固定为4个，则第1，3，4，9号参数建模用两种方式编码的区别如图6-26所示。

图6-26 建模所用参量的两种编码方法

同时，以上两种编码方式也可以稍作改动，将每个个体序列长度加倍。奇数位保留原序列信息，偶数位则用一系列整数表示相应的数学变化（例如：1表示不变， 2表示平方，3表示平方根，4表示对数），这样可以将模型扩展成某些非线性模式，如图6-27所示。

图6-27 编码的非线性扩展

遗传算法对于初始值不敏感。因此只要编码方式确定后，初始种群可以随机生成，只要遗传算法各个参数设置得当，初始值的变化不会对最后的结果产生影响。

6.5.4.3 适应函数

遗传算法对适应函数的选择比较灵活，只要满足以下两点便可作为适应函数：①越好的个体适应性函数值越大；②所有个体的适应性函数值非负。

在2D-QSAR中，适应性函数可以直接采用QSAR模型（一般是MLR或PLS方法，其他如ANN或SVM也可，不过速度较慢）的回归系数r、留一法交叉验证的回归系数q，或者这两者之和。有些时候模型很差，其回归系数r可以是负数，这时需对其作归一化或标准化处理。

6.5.4.4 遗传算子

遗传算法当中需要设置的遗传算子有选择因子、交叉因子、逆转因子、变异因子等。

（1）遗传算法选择方式多种多样，主要的有以下几种：①适应值比例选择，即简单地根据适应值的大小依比例选择，也叫轮盘赌选择。②Boltzmann选择，在开始时选择压力较小，个体间选择概率差别较小，保证初期的个体多样性，随后逐步升高选择压力。③联赛选择（或称锦标赛选择），从所有个体中随机选择一定数量的个体参加“联赛”，保留参赛个体中的最好个体，这种方法容易使搜索过程早熟或停顿。④其他的选择方式还有排序选择、稳态选择、精英选择等。

（2）交叉因子（或称杂交算子）：用于控制基因重组的步骤。最简单的是单点交叉，即在每个个体中随机位置交叉，但是单点交叉有较强的模式破坏能力，所以通常采用两点交叉。

（3）逆转因子：逆转是遗传算法中很少采用的一个步骤，由Holland提出，指的是个体的某一段区域中的基因左右倒转。在2D-QSAR中通常不需要设置这一步骤。

（4）变异因子：每一个位点突变的概率称为个体变异概率，整个个体突变的概率称为基因变异概率。通常采用个体变异即可，其变异概率一般取0.005～0.01。

6.5.4.5 部分替换

为了防止遗传算法陷入局部最小值，每隔数百代，可以采用一次部分替换的方法。可以将种群中适应度最低的个体，或者适应度低于某一域值的个体全部去除，重新补充随机序列。

6.5.4.6 精英种群

当遗传算法收敛后，最后一代个体便是我们所需要的模型。不过也可以额外建立一个精英种群，每一代结束后，将最优的部分个体与精英种群中的最差个体比较，适应度超出便进行替换，这样使精英种群中的个体适应度不断提高，最后的模型也从精英种群中选取，这种方法可以提高遗传算法的效率。

6.5.4.7 种群大小

遗传算法的群体大小建议设置在20～200之间。群体规模越大，包含的模式越多，求得全局最优解的概率也就越高，但同时计算速度变慢，效率降低。同时，如果每个个体的字符串长度很长，群体也要相应扩大。

有时种群的大小并不固定，如meta-GA。这是一种二元遗传算法，有一个种群专门用于控制另一个种群的规模。