逐步回归法与PLS-Bootstrap法在肼衍生物与Pu(Ⅳ)反应定量构效关系中的应用

朱晓乐，肖松涛，欧阳应根

中国原子能科学研究院 放射化学研究所，北京 102413

钚是一种重要的战略资源，是制造核武器和核燃料的关键材料。Pu在水溶液中的化学行为比较复杂，以多价态共存的方式存在。其中Pu(Ⅳ)易被TBP萃取到有机相，而Pu(Ⅲ)不易被萃取到有机相，在PUREX流程中可利用这个特点将Pu(Ⅳ)还原为Pu(Ⅲ)实现铀钚的分离和钚的纯化。

王婷[1]、邵开元[2]、石灵娟[3]等分别对有机污染物的混合毒性、催眠类药物的毒性、多环芳烃的荧光强度进行了定量构效关系研究，均得到了理想的结果。其中邵开元[2]为了改善定量构效关系(quantitative structure-activity relationship, QSAR)模型，提出了化学势变化率这一量化参数，最终使QSAR相关性得到提高，降低了预测误差。

Norinder[4]、邓景景[5]和张永红[6]将偏最小二乘法(partial least square, PLS)应用于QSAR研究，在处理多重相关性的问题上表现良好。邓景景还指出，使用PLS方法建模，较多元线性回归(multiple linear regression, MLR)具有更好的拟合能力和稳定性。

PLS-Bootstrap法用于定量构效关系研究的情况较少。金志超[7]提出，使用PLS-Bootstrap法对变量进行筛选，能解决受实验条件限制导致的样本数较少的问题，并可应用于化学药物的活性与亲油性等参数之间相关关系的预测。Bras[8]则将PLS-Bootstrap法进行了拓展，用于筛选波数间隔。虽尚未发现有文献报道将PLS-Bootstrap法用于QSAR研究工作中，但由于该方法可实现变量的筛选，因此本课题组将使用PLS-Bootstrap法进行QSAR研究。本工作将在之前的研究基础上[9]，在剔除错误样本点后重新建立定量构效关系。将分别使用逐步回归法和PLS-Bootstrap法来建立定量构效关系模型并进行交叉检验和外部检验。

1 数据与方法

Pu(Ⅳ)与肼衍生物反应的半反应时间t1/2及它们的量化参数列入表1和表2[9]，肼衍生物结构示于图1[9]。考虑到肼衍生物与Pu(Ⅳ)的半反应时间存在着数量级差异，因此将它们转化成自然对数后再进行定量构效关系研究。肼衍生物与Pu(Ⅳ)反应的半反应时间可用于表示该反应进行的快慢程度。半反应时间越小，该反应进行得越快。量化参数有最高占据轨道能(EHOMO)、最低非占轨道能(ELUMO)、前线轨道能量差(ΔE=ELUMO-EHOMO)、分子总能量(E)、分子偶极矩(μ)、疏水性参数(lgP)、分子折射率(R)、分子摩尔体积(V)、分子表面积(A)、网格化分子表面积(G)、相对分子质量(Mr)、分子极化率(P)以及水合能(EH)等。

表1 肼衍生物与Pu(Ⅳ)反应的半反应时间[9]

为了考察半反应时间与量化参数的相关性、量化参数之间的相关性，求解其Pearson矩阵，列于表3。由于Pearson矩阵为对称矩阵，且主对角元的相关性系数均为1，概率p值均为0，因此仅求出其上半部分的相关性系数及对应的p值。由表3可知，半反应时间(t1/2)与最高占据轨道能(EHOMO)、前线轨道能量差(ΔE)、分子偶极矩(μ)、疏水性参数(lgP)有较强的线性相关性，相关系数均达到0.6以上，且不相关的可能性p均小于0.05，即半反应时间与这些参数均存在相关性，但线性相关性并不强。此外，各参数之间也存在相关关系，如最高占据轨道能(EHOMO)与前线轨道能量差(ΔE)存在强相关性，与偶极矩存在较强的相关性；最低非占轨道能(ELUMO)则与前线轨道能量差(ΔE)、疏水性参数(lgP)、水合能(EH)之外的其他量化参数均存在较强的相关性。即量化参数之间存在着较强的相关性，无法直接进行回归建模，因此将分别使用逐步回归法和PLS-Bootstrap法对量化参数进行筛选，再通过回归分析的方法建立定量构效关系模型，且半反应时间与量化参数之间存在相关性但线性相关性并不强，故还需要考虑非线性定量构效关系模型。

表2 肼衍生物的量化参数[9]

表3 Pearson相关性矩阵

逐步回归法是结合了变量筛选与回归的一种方法，在每一步中，都要分别考查所有的参数与因变量的相关性，将未引入方程中但相关性最显著的参数引入方程，即将p值最小的参数引入方程，而将方程中相关性不够强的参数剔除掉。为防止被剔除的参数能再次被引入而形成死循环，一般设置剔除参数的p临界值(premove)要大于进入方程的p临界值(penter)。一般默认设置为penter=0.05，premove=0.10。逐步回归法的结果是否有意义需要由F检验确定，各参数与因变量的相关性是否显著则应当由t检验确定。

PLS被认为是一种能较好处理自变量共线性和样本数较少的回归方法，在对每一个成分进行回归时都要完成提取主成分、相关性分析和回归三个步骤。首先对自变量与因变量各提取一个主成分，并使自变量的主成分与因变量的主成分相关性最大，其次使用自变量的主成分对因变量进行回归，然后分别求得经过回归之后自变量与因变量的残差，并使用残差分别代替自变量与因变量重复前述步骤进行回归。一般情况下不需要选择所有的成分来建立回归方程，仅需要选用前几个成分即可得到预测能力较好的回归方程。需要注意的是，使用PLS进行回归时所有的参数都会被引入方程，其中可能会包括部分相关性不显著的参数，考虑到Bootstrap作为一种检验方法，能较好地处理样本数较少的情况，通过蒙特卡罗随机抽样的方法建立Bootstrap样本，然后对Bootstrap样本进行回归，将回归系数按大小顺序排列并选定一定的分位数，若该数大于原样本回归方程中的回归系数，则该参数应当从方程中剔除。将PLS与Bootstrap结合起来进行回归分析时，可以将相关性不显著的参数从方程中剔除掉，使方程的意义更明确。由于Bootstrap本身已包含了对各参数的检验，故无需对各参数进行t检验。

本工作使用了Matlab 2014b软件来进行回归处理，在逐步回归法中调用了程序自带的stepwise和stepwisefit函数来进行逐步回归，在PLS-Bootstrap法中调用了程序自带的plsregress函数实现PLS的部分，使用rand函数结合for循环语句，按照文献[7]的步骤来实现Bootstrap的部分。

2 结果与讨论

2.1 逐步回归法对参数的筛选及模型的建立

利用逐步回归法进行回归分析，在完成回归分析任务的同时，也完成了变量的筛选。在默认的设置(penter≤0.05,premove>0.10)下，对其进行逐步回归，得到的结果列入表 4。由表4知，第1、2、6个变量，即最高占据轨道能(EHOMO)、最低非占轨道能(ELUMO)和疏水性参数(lgP)保留在方程中，方程为t1/2=54.939EHOMO+176.841ELUMO+6.363 lgP+27.718。其方程的分子自由度为3，分母自由度为6，F统计量值为37.719，即F(3,6)=37.719>4.76=F0.05(3,6)，即在95%的置信度上该回归方程有意义，对应的调整相关系数即调整r2=0.924。标准化回归方程为y=0.384x1+0.463x2+0.701x6。标准化回归系数的绝对值反映了该自变量对因变量的影响大小，故疏水性参数是该反应的主要特征参数，且疏水性参数越大，半反应时间越大，该反应进行得越慢；最高占据轨道能(EHOMO)和最低非占轨道能(ELUMO)是该反应的次要特征参数，且它们的值越大，半反应时间越大，该反应进行得越慢。

下面考虑非线性关系，由于参数较多，考虑指数形式，假设存在指数关系，则有

y=Ae(Bx1+Cx2+Dx3)=e(Bx1+Cx2+Dx3+A′)

其中A′=lnA。两边同取对数得：

lny=A′+Bx1+Cx2+Dx3

若将半反应时间取自然对数，则其对数值与变量依旧呈线性关系。按此方法处理后进行回归分析，得非线性回归方程t1/2=e(5.543EHOMO+17.201ELUMO+0.560lg P+4.010)，方程的调整r2=0.928，F统计量为39.964>4.76=F0.05(3,6)。标准化回归方程t1/2=e(0.413x1+0.479x2+0.656x3)。疏水性参数(lgP)仍为该反应的主要特征参数，最高占据轨道能(EHOMO)和最低非占轨道能(ELUMO)仍为该反应的次要特征参数，且它们的值越大，半反应时间越大，该反应进行得越慢。

表4 逐步回归结果

注：“In”为变量在方程中出现，“Out”为变量在方程中未出现

2.2 PLS-Bootstrap法对参数的筛选及模型的建立

在使用PLS进行回归分析时，需要确定提取的成分数。一般情况下，当所选取成分包含的信息量达到85%时，便认为信息已经提取完全，因此，把X信息量与Y信息量均达到85%时的最小成分数作为PLS回归分析时提取成分数，对所有样本进行处理，结果列入表5。

表5 PLS回归分析的成分数与信息量关系

注：1) “-”，在第二次回归中，当成分数为4时，X信息量已达到1，继续增加成分数也不会提升信息量，故无需对更多成分数的情况进行计算

在第一次筛选中，成分数为2时，包含的X信息量仅为0.801，而Y信息量达到了0.872；当成分数为3时，包含的X信息量为0.877，Y信息量为0.946，因此，选定PLS提取的成分数为3，进行Bootstrap检验，结果列入表6。在经过第一次筛选后，仅有最高占据轨道能、分子总能量、疏水性参数和相对分子质量被保留。然后进行第二次的PLS-Bootstrap法筛选，首先进行PLS回归分析(表5)，确定提取的成分数为3，然后再进行Bootstrap检验，结果列入表6。由表6可知，最终仅有最高占据轨道能和疏水性参数被保留下来。下面将使用经过两次筛选后仍被保留的参数进行回归分析。使用最高占据轨道能和疏水性参数对半反应时间进行回归分析，得方程t1/2=92.673EHOMO+4.866 lgP+34.956。调整r2=0.766，统计量F=15.730>4.74=F0.05(2,7)，即在95%的置信度上该回归方程有意义，模型示于图2。所有的PLS-Bootstrap法的分位数均选为0.95。

以同样的方法对半反应时间进行处理以考虑其非线性关系，并对其进行回归分析，得非线性回归方程t1/2=e(0.195 8EHOMO+0.138 6lg P+2.371)，调整r2=0.758，模型示于图3。比较其调整r2，显然线性模型较非线性模型略优。

表6 PLS-Bootstrap法对参数进行筛选的结果

注：1) “-”，由于第一次筛选已经将对应的变量排除在方程外，它们将不会参与第二次筛选，故其对应值均不存在

图2 PLS-Bootstrap法线性模型

图3 PLS-Bootstrap法非线性模型

无论是线性模型还是非线性模型，结果均表明，最高占据轨道能(EHOMO)和疏水性参数(lgP)是该反应的特征参数，且它们的值越大，该反应的半反应时间越大，即反应进行得越慢。

2.3 模型的交叉检验与外部检验

为了对模型进行比较，同时也对前人定量构效研究进行补充，将对两种方法建立的模型进行留一法(leave one only, LOO)交叉检验，并对被剔除的样本点叔丁基肼的部分量化参数进行订正后作为外部检验点进行检验[10]。

使用LOO交叉检验后得到的结果列入表 7。由表7可知：两种方法的交叉检验系数(Q2)均大于0.5，其中逐步回归法的Q2达到了0.820，表明逐步回归法建立的线性模型稳定性较优秀，因此就稳定性而言，逐步回归法建立的线性模型更好。由于逐步回归法建立的线性模型与非线性模型相比无明显差别，调整r2分别为0.924和0.928，因此不再单独对其非线性模型进行讨论，即使非线性模型从调整r2上看比线性模型略优。

表7 逐步回归法及PLS-Bootstrap法交叉检验结果

由于两种方法仅用到了最高占据轨道能(EHOMO)、最低非占轨道能(ELUMO)和疏水性参数(lgP)，使用了参考文献[8]的逐步回归法，在对叔丁基肼的结构进行校正后重新计算了这些参数的值，结果如下：最高占据轨道能EHOMO=-0.237，最低非占轨道能ELUMO=-0.006，疏水性参数 lgP=0.370，半反应时间对数值 lnt1/2=16.240。使用两种线性模型分别对lnt1/2进行预测，并与观测值进行对比，得到逐步回归模型预测值为15.991，PLS-Bootstrap模型预测值为14.793，结果列入表8。由表8可知，逐步回归法对外部样本点的预测误差更小，误差不足2%，而PLS-Bootstrap模型的预测误差已接近10%。就叔丁基肼作为外部检验点而言，逐步回归法表现更优秀，但外部检验点样本点数量较少，并不意味着逐步回归法的预测能力一定强于PLS-Bootstrap法。

表8 外部检验结果

总体上来看，逐步回归模型中包含了三个参数，且交叉检验、外部检验结果均优于PLS-Bootstrap模型，而后者虽然只包含了两个参数，但从Pearson相关性矩阵可以看出，这两个参数均与半反应时间有较强的相关性，而差异之处在于相关性相对较弱的最低非占轨道能(ELUMO)，这正好表明了通过提取成分进行相关性分析的PLS更能有效地保证自变量与因变量之间的相关性。

2.4 两种回归方法的比较

为了进一步探讨最低非占轨道能(ELUMO)对模型的影响，以比较两种回归方法的异同，下面将使用相同的变量分别建立回归方程。由于参与回归的变量已经被指定，逐步回归法无需再进行变量筛选，因此可使用最小二乘回归法代替逐步回归法进行比较，而PLS-Bootstrap法无需再进行变量筛选，在回归过程中也无须进行Bootstrap检验。对比实验一指定的变量有疏水性参数、最高占据轨道能和最低非占轨道能，对比实验二指定的变量有疏水性参数和最高占据轨道能。

对比实验一中，逐步回归法回归方程为t1/2=54.939EHOMO+176.841ELUMO+6.363 lgP+27.718，PLS回归方程为t1/2=70.866EHOMO+126.893ELUMO+5.916 lgP+30.992；对比实验二中，逐步回归法回归方程为t1/2=92.673EHOMO+4.866 lgP+34.956，PLS回归方程为t1/2=92.673EHOMO+4.866 lgP+34.956；对比实验一与实验二的各参数列于表9 。由表9可知，当变量数为三个时，PLS在回归时认为仅需提取两个成分数即可建立效果较好的回归方程，当所有的成分数均提取完全即将三个成分全部提取时，其回归结果与逐步回归法的结果相同，即第三个成分对改善模型的拟合优度，即调整r2的贡献不大。对比可知，第三个成分对调整r2的提高仅为0.004，较0.920提高了0.43%，表明了PLS能选取合适的成分数进行回归的同时，舍弃对回归结果作用不明显的部分。但就对比实验一中对逐步回归法与PLS在LOO交叉检验的结果进行对比是不合理的，因为逐步回归法在交叉检验时的变量数为3，而PLS由于仅提取了两个成分，在交叉检验时仅有两个成分进行回归，即变量数为2，故无法简单地将其进行比较。而考虑变量数相同的情况时，在对比实验二中，逐步回归法交叉检验时变量数为2，而PLS在对比实验一中，交叉检验时提取两个成分，变量数也为2，对比交叉检验系数Q2，可得出PLS的模型稳定性更好。综上，对比逐步回归法与PLS，可认为PLS的回归结果更好。

表9 对比实验一与实验二的回归情况

注：1) 回归时提取两个成分，F统计量为F(2,7)=52.776

将对比实验二的逐步回归法与对比实验一中的PLS在交叉检验中的表现进行对比，虽然存在一定的不合理性，但由对比实验二可以看出，PLS在提取所有的成分数的情况下与逐步回归法有相同的结果。将两个对比实验中PLS的表现进行对比，两种情况下PLS均提取了两个成分进行回归，差别在于成分中是否包含ELUMO这一参数。由调整r2的比较可知，该参数的引入对模型有较大的改善，以调整r2的值来看，拟合优度的改善程度为20.10%，但以Q2来看，模型稳定性的提高程度却仅为16.95%。即使这样的比较不准确，却仍能说明ELUMO参数的引入对模型拟合效果的改善作用更大，即ELUMO参数引入的作用存在着过分追求方程拟合效果的可能性。

综上所述，在肼衍生物与Pu(Ⅳ)的反应中，PLS-Bootstrap可使用较少的变量来描述整体的趋势变化，而逐步回归法能获得预测效果更好的模型，两者的结果具有相似性，即两种方法均可得到反映该反应规律的回归方程。就肼衍生物与Pu(Ⅳ)反应而言，ELUMO是否能作为该反应的特征参数，还需进一步研究确定。

3 结 论

通过逐步回归法与PLS-Bootstrap法分别对Pu(Ⅳ)与肼衍生物反应的定量构效关系进行了研究，获得了可描述该反应进行快慢的特征参数：最高占据轨道能(EHOMO)、疏水性参数(lgP)，其中疏水性参数是该反应的主要特征参数，且它们的值越大，半反应时间越大，该反应进行得越慢。而作为次要特征参数的最低非占轨道能(ELUMO)能否作为该反应的特征参数仍需进一步研究确定，因此两种方法得到的结果具有相似性。

对之前的Pu(Ⅳ)与肼衍生物反应的定量构效关系研究的不足之处进行了补充，对错误的样本点进行了校正并作为外部检验点对模型进行了检验，完成了之前模型未进行的交叉检验等检验工作。

本工作使用PLS-Bootstrap法进行定量构效关系研究并进行了交叉检验与外部检验，结果表明，该模型的稳定性较逐步回归模型差，考虑到外部检验点数量少，其预测能力还需进一步研究确定。