化學(xué)建模與模型集群分析

云永歡　鄧百川　梁逸曾

摘 要 本文簡單介紹了化學(xué)建模與模型集群分析的思想，并列舉了基于模型集群分析的思路與框架。近年來，應(yīng)用于化學(xué)建模各個(gè)方面的許多新算法包括奇異樣本診斷、變量選擇、模型參數(shù)與評(píng)價(jià)、穩(wěn)健與模型應(yīng)用域。本文通過應(yīng)用于不同的數(shù)據(jù)類型，包括近紅外光譜、定量構(gòu)效關(guān)系及代謝組學(xué)數(shù)據(jù)，舉例闡述模型集群分析方法的可行性與應(yīng)用性，為未來開發(fā)化學(xué)建模新算法提供一個(gè)好的思路和框架。

關(guān)鍵詞 化學(xué)建模； 模型集群分析； 采樣； 統(tǒng)計(jì)分析； 綜述

1 引 言

隨著化學(xué)量測(cè)數(shù)據(jù)的不斷累積和大數(shù)據(jù)信息處理技術(shù)，包括數(shù)據(jù)發(fā)掘和機(jī)器學(xué)習(xí)各種新方法的不斷涌現(xiàn)，采用化學(xué)建模（Chemical modeling）方法進(jìn)行化學(xué)知識(shí)規(guī)律發(fā)現(xiàn)及建立定量模型等研究得到了飛速發(fā)展。此外，在分析化學(xué)的發(fā)展過程中，由于儀器分析的飛速發(fā)展，復(fù)雜體系的快速儀器分析，包括近紅外和拉曼光譜無損分析及各類波譜如質(zhì)譜、激光誘導(dǎo)擊穿光譜（LIBS）等的分析、代謝組學(xué)中核磁共振譜及各種色譜分析、中藥色譜指紋圖譜分析等，現(xiàn)都已成為了分析化學(xué)的重要研究方向[1]。值得提出的是，這樣的化學(xué)建模的共同特點(diǎn)是它們的模型都可由下述簡單算式給出： y=f（X）。式中， y為含n個(gè)元素的列矢量，每個(gè)元素都表征一個(gè)樣本的定性特征或定量指標(biāo)，而矩陣X則為含n行的矩陣，每行為一系列表征樣本屬性特征（含p個(gè)元素）或一個(gè)樣本的測(cè)量譜（波譜或色譜）； f（.）為不定的函數(shù)關(guān)系，它可以是線性的，如主成分回歸（PCR）或偏最小二乘（PLS）；也可以是非線性的，如支撐向量機(jī)（SVM）或人工神經(jīng)網(wǎng)絡(luò)（ANN）等。其關(guān)系見圖1。

由圖1可見，此類數(shù)據(jù)體系（包括紫外、近紅外、拉曼光譜分析、定量構(gòu)效關(guān)系和代謝組學(xué)數(shù)據(jù)）是一類極具復(fù)雜性的體系，由于其函數(shù)關(guān)系f（.）是未知的，線性或非線性無法確定，變量與y的關(guān)系不明確，沒有任何物理或化學(xué)定理可作為基礎(chǔ)，解空間類似美國著名統(tǒng)計(jì)學(xué)家George E. P. Box所說的那樣，即“所有模型都是錯(cuò)誤的，但其中有些是有用的（All models are wrong， and some are useful.）”。所以，對(duì)于這樣的復(fù)雜體系，找到盡量逼近的基空間，并通過有效模型評(píng)價(jià)方法及其可靠應(yīng)用域的定義方法十分重要。

近年來，化學(xué)與生物領(lǐng)域引入大量高通量分析技術(shù)，使得上述模型中的 xi 這個(gè)行矢量變得很長，而且其中還有很多變量與yi無關(guān)，甚至還有干擾作用[2～5]。此外，由于目前樣本數(shù)（n）相對(duì)較少，出現(xiàn)了在統(tǒng)計(jì)學(xué)稱為維數(shù)災(zāi)禍的“大p，小n”問題，這是目前統(tǒng)計(jì)學(xué)及其應(yīng)用領(lǐng)域研究的重大挑戰(zhàn)[6～8]。對(duì)于這樣的體系，很容易出現(xiàn)模型過擬合，建模須謹(jǐn)慎[9，10]。

2 化學(xué)建模與模型集群分析

化學(xué)計(jì)量學(xué)和化學(xué)信息學(xué)研究的一個(gè)主要目標(biāo)就在于建立一個(gè)有效并可靠的化學(xué)模型，以對(duì)未知的化學(xué)樣本的濃度/性質(zhì)等進(jìn)行預(yù)測(cè)。從上述的分析可知，這個(gè)任務(wù)不簡單，由于模型完全未知，建模有點(diǎn)類似“瞎子摸象”的任務(wù)。而模型集群分析（Model population analysis，MPA）[9，11] 打破傳統(tǒng)一次性建模思路，力求最大限度地利用已有樣本集的信息，通過隨機(jī)采樣，從不同角度考察數(shù)據(jù)集的內(nèi)在性質(zhì)，通過對(duì)所得結(jié)果進(jìn)一步統(tǒng)計(jì)分析，獲得數(shù)據(jù)集的內(nèi)在結(jié)構(gòu)。從這個(gè)角度看來，模型集群分析與貝葉斯統(tǒng)計(jì)分析的追求后驗(yàn)分布有些類似。而且，模型集群分析中主要是強(qiáng)調(diào)集群分析，強(qiáng)調(diào)所得的各種不同結(jié)果的分布，與一次性建模分析形成了強(qiáng)烈對(duì)比。

基于模型集群分析的化學(xué)建模算法之構(gòu)建框架示于圖2。它的構(gòu)建框架主要包括3個(gè)基本要素：（1）通過隨機(jī)采樣獲取子數(shù)據(jù)集；（2）針對(duì)每個(gè)子數(shù)據(jù)集，建立一個(gè)子模型；（3）從樣本空間、變量空間、參數(shù)空間或模型空間對(duì)所有建立的集群子模型的感興趣的參數(shù)進(jìn)行統(tǒng)計(jì)分析，獲取有用的信息。

3 模型集群分析的3個(gè)基本要素

3.1 隨機(jī)采樣獲取子數(shù)據(jù)集

采樣是數(shù)據(jù)進(jìn)行統(tǒng)計(jì)分析中常用的有力工具[12]。通過隨機(jī)采樣，可以從給定一個(gè)數(shù)據(jù)集的樣本或變量空間中獲取N個(gè)子數(shù)據(jù)集。如果從樣本空間獲取，子數(shù)據(jù)集由部分樣本組成，從變量空間獲取，則子數(shù)據(jù)集由部分變量組成； 如果同時(shí)從樣本和變量空間獲取，子數(shù)據(jù)集則由選取部分樣本和部分變量組成。目前，應(yīng)用于MPA中常用的隨機(jī)采樣方法有4種：（1）蒙特卡洛采樣（Monte Carlo sampling）[13]、（2）自助法采樣（Bootstrap sampling）[14]、（3）二進(jìn)制采樣（Binary matrix sampling）[15～19]、（4）重排技術(shù)（Permutation）[20]。

蒙特卡洛采樣，無放回采樣，也稱刀切法（Jack knife），隨機(jī)選取一定比例的樣本作為訓(xùn)練集，剩下的樣本作為獨(dú)立測(cè)試集。

自助法采樣，有放回采樣，每個(gè)樣本被選中的概率相同，因此會(huì)出現(xiàn)有些樣本被多次選中，這些樣本可作為重復(fù)樣本，也可以一次使用。隨機(jī)性選取一定比例的樣本作為訓(xùn)練集，其余的樣本為獨(dú)立測(cè)試集。

二進(jìn)制采樣是新近提出的一種采樣方法，首先產(chǎn)生二進(jìn)制矩陣，行代表采樣次數(shù)，列代表對(duì)應(yīng)數(shù)據(jù)的變量，矩陣只是含0和1，0代表變量被選中，1代表變量沒被選中，每一列0和1的比例統(tǒng)一設(shè)定，接著每一列自主打亂，根據(jù)每行有1的位置選取變量，由于每列1的數(shù)目是固定的，這種方法能夠保證按行選取變量時(shí)，變量被選擇的概率相同。

重排技術(shù)是對(duì)樣本矩陣或響應(yīng)值矢量進(jìn)行重排打亂，然后再建模，普遍應(yīng)用于檢查過擬合風(fēng)險(xiǎn)。

3.2 建立子模型

對(duì)所有產(chǎn)生的N個(gè)數(shù)據(jù)集，采用選定的建模方法建立模型，將得到N個(gè)子模型。由于每個(gè)子模型是建立在相對(duì)應(yīng)的子數(shù)據(jù)集上，僅反映了原數(shù)據(jù)集的局部信息，建立N個(gè)子模型可較全面地反映原始數(shù)據(jù)集的信息。目前比較常用的建模方法有線性方法和非線性方法。（1）線性方法包括多元線性回歸（MLR）、偏最小二乘（PLS）、主成分回歸（PCR）、嶺回歸（RR）、Lasso回歸、線性判別分析（LDA）；（2）非線性方法包括支撐向量機(jī)（SVM）、神經(jīng)網(wǎng)絡(luò)（ANN）、分類回歸樹（CART）、隨機(jī)森林（Random forest）。endprint

3.3 統(tǒng)計(jì)分析

MPA的核心思想是對(duì)獲得的由N個(gè)子模型構(gòu)成的集群模型的某個(gè)感興趣的參數(shù)進(jìn)行統(tǒng)計(jì)分析，通過統(tǒng)計(jì)分布獲取對(duì)解決實(shí)際問題有價(jià)值的信息。實(shí)際上，由于應(yīng)用的復(fù)雜性與多樣性，對(duì)感興趣的參數(shù)進(jìn)行統(tǒng)計(jì)分析的策略需要根據(jù)具體情況進(jìn)行具體分析與設(shè)計(jì)，不同的策略與設(shè)計(jì)將會(huì)產(chǎn)生不同的算法。而這些對(duì)所有建立的集群子模型的感興趣的參數(shù)是從樣本空間＼，變量空間＼，參數(shù)空間或者模型空間4個(gè)空間獲取的。樣本空間：與樣本相關(guān)的模型輸出，如回歸模型里樣本的擬合誤差，預(yù)測(cè)誤差；分類模型里，樣本類別的預(yù)測(cè)準(zhǔn)確率。變量空間：與變量相關(guān)的模型輸出，如變量的回歸系數(shù)。參數(shù)空間：與模型自身相關(guān)的參數(shù)，如PLS模型主成份的個(gè)數(shù)，隨機(jī)森林與支撐向量機(jī)模型的間隔。模型空間：模型相對(duì)于其它模型的參數(shù)，是由集群模型共同決定，如模型在空間的位置＼，模型的分布密度。

繼通過統(tǒng)計(jì)分布分析所有建立的集群子模型的感興趣的參數(shù)，如：（1）對(duì)正常樣本與奇異樣本的兩類預(yù)測(cè)誤差分布診斷奇異樣本；（2）不同的變量組合的交互檢驗(yàn)預(yù)測(cè)誤差分布來找出最優(yōu)的變量子集；（3）比較某個(gè)變量組合中某個(gè)變量存在和不存在模型里時(shí)的兩個(gè)交互檢驗(yàn)預(yù)測(cè)誤差分布來獲得變量重要性；（4）比較每個(gè)變量被重排前后的預(yù)測(cè)誤差分布來獲得每個(gè)變量的重要性。

針對(duì)這些分布，利用其最大最小值、均值、標(biāo)準(zhǔn)差、均值標(biāo)準(zhǔn)差比、95%置信區(qū)間、t檢驗(yàn)（有參數(shù)檢驗(yàn)）和Mann-Whitney U檢驗(yàn)[21]（無參數(shù)檢驗(yàn)）獲得有價(jià)值的信息。4 基于MPA的新算法在化學(xué)建模中的應(yīng)用

MPA是基于建立集群模型的一種數(shù)據(jù)分析思路。它是數(shù)據(jù)分析的一般性框架，為系統(tǒng)研究數(shù)據(jù)結(jié)構(gòu)、建立模型及算法設(shè)計(jì)等提供了新的思維方式?；贛PA的3個(gè)基本要素和4個(gè)空間，近來已經(jīng)開發(fā)了應(yīng)用于化學(xué)建模的許多新算法，包括奇異樣本診斷、變量選擇、模型參數(shù)與評(píng)價(jià)、穩(wěn)健與模型應(yīng)用域。下面對(duì)化學(xué)建模的這幾個(gè)方面的應(yīng)用進(jìn)行簡單舉例介紹。

4.1 奇異樣本診斷

構(gòu)建穩(wěn)健的化學(xué)模型主要取決于訓(xùn)練集數(shù)據(jù)樣本。如果訓(xùn)練集數(shù)據(jù)中包括一些遠(yuǎn)離數(shù)據(jù)主體的奇異樣本，它們將會(huì)破壞整個(gè)數(shù)據(jù)結(jié)構(gòu)，從而影響模型的建立以及預(yù)測(cè)。因此，奇異樣本診斷是穩(wěn)健化學(xué)建模的一個(gè)關(guān)鍵步驟[22]。Cao等[23]提出的基于模型集群分析的奇異樣本診斷方法（Monte Carlo sampling， MCS）主要研究了基于模型特征分布診斷奇異樣本，其步驟如下：（1）采用蒙特卡洛采樣從原始數(shù)據(jù)總樣本中選取一定比列的樣本作為訓(xùn)練樣本，如r=80%，剩下的20%樣本作為獨(dú)立測(cè)試集樣本。這個(gè)過程重復(fù)N次，即可得到N個(gè)子訓(xùn)練集和與之對(duì)應(yīng)的N個(gè)子測(cè)試集；（2）每個(gè)子訓(xùn)練集建立模型并對(duì)相應(yīng)的測(cè)試集樣本進(jìn)行預(yù)測(cè)；（3）設(shè)每個(gè)樣本被等概率采樣，則其被選進(jìn)測(cè)試集的次數(shù)接近N（1-r）。因此，每個(gè)樣本將約有N（1-r）個(gè)預(yù)測(cè)誤差，其預(yù)測(cè)誤差的統(tǒng)計(jì)分布特征可用于診斷奇異樣本。

舉一個(gè)使用MCS方法診斷奇異樣本的例子：選取了常用的一組近紅外量測(cè)玉米的光譜數(shù)據(jù)，光譜測(cè)量采用mp5儀器，該量測(cè)數(shù)據(jù)x包含80個(gè)玉米樣本，玉米的淀粉含量作為響應(yīng)變量y，光譜波長區(qū)間為1100～2498 nm，間隔為2 nm，總獲得700個(gè)量測(cè)波長（該數(shù)據(jù)可從http：//www.eigenvector.com/data/Corn/index.html免費(fèi)下載）。PLS作為校正模型的方法，蒙特卡洛采樣次數(shù)N=10000，每次采樣，80%樣本作為訓(xùn)練集，剩下的20%樣本作為獨(dú)立測(cè)試集，PLS潛變量個(gè)數(shù)由10折交互檢驗(yàn)選取為9，采用MCS方法對(duì)該數(shù)據(jù)的奇異樣本診斷結(jié)果如圖3所示，MCS判斷出了三類樣本分別為正常樣本（圖3a），X方向奇異樣本（圖3b）和Y方向奇異樣本（圖3c和d）。圖3A顯示了對(duì)應(yīng)的a， b， c 三類樣本的預(yù)測(cè)誤差分布圖，可以看出正常樣本A，預(yù)測(cè)誤差分布在原點(diǎn)附近，均值接近0，分布高而窄，說明其有很小的不確定性。對(duì)于X方向的奇異樣本，由于其遠(yuǎn)離數(shù)據(jù)主體，用不同的樣本得到的不同的模型將會(huì)產(chǎn)生一個(gè)很寬的預(yù)測(cè)誤差分布，標(biāo)準(zhǔn)差較大。而對(duì)于Y方向的奇異樣本，預(yù)測(cè)誤差分布的均值遠(yuǎn)離原點(diǎn)和標(biāo)準(zhǔn)差也較大。因此，這些結(jié)果都表明，僅利用一次模型得到的一次預(yù)測(cè)誤差診斷奇異樣本是不充分，采用預(yù)測(cè)誤差分布來診斷奇異樣本才更加可靠和穩(wěn)健，這是采用模型集群分析獲得多個(gè)模型的重要原因。

圖3 A： 根據(jù)預(yù)測(cè)誤差的均值與標(biāo)準(zhǔn)差畫出的奇異值診斷圖， 三類樣本分別為正常樣本（a），X方向奇異樣本（b）和Y方向奇異樣本（c和d）； B： 三類樣本（a，b和c）的預(yù)測(cè)誤差分布圖

Fig.3 （A） The diagnostic plot for outlier detection based on the mean and STD value of prediction errors. Three representative samples are a normal sample （a）， an X-outlier （b）， and a Y-outlier （c and d ）； （B） The distributions of prediction errors of these a， b and c samples

4.2 變量選擇

現(xiàn)代高通量分析儀器的成千上萬個(gè)分析通道可提供豐富的測(cè)量數(shù)據(jù)，但常遇到到“樣本少，變量多”問題[24，25]。而變量選擇無疑是解決此類問題的有效方法[26]。Yun等[27]也證明了復(fù)雜分析體系中變量選擇的重要性與必要性。選擇變量有3個(gè)目的：（1）提高預(yù)測(cè)能力；（2）降低數(shù)據(jù)維數(shù)并選擇更有效的變量；（3）增強(qiáng)模型的可解釋性[28]。然而，變量選擇是一個(gè)NP問題，隨著變量個(gè)數(shù)的增加，變量空間成指數(shù)增大，找到一個(gè)最佳變量組合非常具有挑戰(zhàn)性?；谀Ｐ图悍治龅目蚣芩悸?，新近提出了眾多變量選擇方法，這里簡單介紹一個(gè)代表性方法并舉出相關(guān)應(yīng)用例子。基于變量組合的變量重要性分析（VIAVC）[29]是基于模型集群分析思路對(duì)每個(gè)變量進(jìn)行重要性分析的方法。具體步驟如下：（1）采用二進(jìn)制采樣從原始數(shù)據(jù)總樣本中產(chǎn)生N個(gè)變量組合，每個(gè)變量組合含有一組隨機(jī)變量；（2）每個(gè)變量組合建立一個(gè)子模型并計(jì)算其交互檢驗(yàn)預(yù)測(cè)誤差或準(zhǔn)確率，即可獲得交互檢驗(yàn)預(yù)測(cè)誤差或準(zhǔn)確率的分布；（3）對(duì)每個(gè)變量，觀察其存在或不存在某個(gè)固定變量組合時(shí)前后的差別，因有N個(gè)變量組合，每個(gè)變量都有存在與不存在某個(gè)固定變量組合的分布，采用統(tǒng)計(jì)檢驗(yàn)對(duì)對(duì)其進(jìn)行評(píng)價(jià)，得出的p值即可作為評(píng)價(jià)變量重要性的標(biāo)準(zhǔn)；（4）只保留p<0.05的變量，重復(fù)上述步驟1～3，直至無p>0.05的變量。endprint

選取一組代謝組學(xué)數(shù)據(jù)[30]作為此方法的應(yīng)用例子，該數(shù)據(jù)兩類樣本采自中南大學(xué)湘雅醫(yī)學(xué)院的16例正常兒童血漿樣本和13例超重兒童的血漿樣本。通過島津GCMS-QP2010氣相色譜與質(zhì)譜聯(lián)用儀分析并采用NIST質(zhì)譜庫檢索定性定量分析了30個(gè)代謝產(chǎn)物。VIAVC目的是找出重要的代謝物，這些代謝物用于建模時(shí)能夠達(dá)到變量選擇的3個(gè)目的，即（1）提高兩類樣本的預(yù)測(cè)準(zhǔn)確率；（2）選擇少并有效的變量來建模；（3）變量的可解釋性。根據(jù)VIAVC原理，以受試者工作特征曲線 （Receiver operating haracteristic curve， ROC曲線）下面的面積（Area under roc curve，AUC）作為模型評(píng)價(jià)指標(biāo)[31， 32] ，結(jié)合統(tǒng)計(jì)t檢驗(yàn)，挖掘出4類變量（圖4），分別為強(qiáng)有信息變量、弱有信息變量、無信息變量和干擾變量。經(jīng)過3次迭代后，VIAVC最終保留了13個(gè)有信息變量，根據(jù)統(tǒng)計(jì)檢驗(yàn)對(duì)兩個(gè)分布計(jì)算的p值來排序這13個(gè)變量，最后再利用10折雙層交互檢驗(yàn)按照排序向前選擇，找出最佳的變量組合，前3個(gè)變量。這3個(gè)代謝物分別為β-羥基丁酸、甘油酸和棕櫚酸，他們的交互檢驗(yàn)預(yù)測(cè)準(zhǔn)確率為86.21%，與全部變量的交互檢驗(yàn)預(yù)測(cè)準(zhǔn)確率65.52%相比，選擇變量大大提高了預(yù)測(cè)能力，所選擇的3個(gè)代謝物也被驗(yàn)證與肥胖疾病有關(guān)[33～35]。

圖4 4類變量分布圖，a為強(qiáng)有信息變量，存在該變量時(shí)，交互檢驗(yàn)AUC值（AUCCV）顯著性提高，t檢驗(yàn)p0.05；b為弱有信息變量，存在該變量時(shí)， AUCCV稍微有提高但不顯著，t檢驗(yàn)p>0.05；c為無有信息變量，存在該變量時(shí)，AUCCV稍有降低但不顯著，t檢驗(yàn)p>0.05；d為干擾變量，存在該變量時(shí)， AUCCV顯著性降低，t檢驗(yàn)p<0.05。

Fig.4 Four kinds of variable， a is the strongly informative variable， when inclusion of this variable， area under roc curve （AUC） value of cross validation （AUCCV） has improved significantly， and p value of t test is much less than 0.05； b is the weakly informative variable， when inclusion of this variable， AUCCV has improved but not significantly， and p value of t test is more than 0.05； c is the uninformative variable， when inclusion of this variable， AUCCV has decreased a little， and p value of t test is more than 0.05； d is the interfering variable， when the inclusion of this variable， AUCCV has decreased significantly， and p value of t test is much less than 0.05

除了以上方法，近年來有很多基于模型集群分析思路開發(fā)的變量選擇新方法，Monte Carlo based uninformative variable elimination （MC-UVE）[36]， Competitive adaptive reweighted sampling （CARS）[37，38]，Margin influence analysis （MIA）[39]，Iteratively retaining informative variables （IRIV）[40]，Random frog[41，42]，Variable combination population analysis（VCPA）[17]，Variable iterative space shrinkage approach （VISSA）[15，16]，Modified mutual information-based feature selection algorithm （MMIFS）[43]，Randomization test（RT）[44]，Variable complementary network（VCN）[45]， Subwindow permutation analysis （SPA）[4，46]。在圖5中，每種方法的采樣技術(shù)、采樣空間、參數(shù)輸出、統(tǒng)計(jì)分析均通過連接線畫出。如MC-UVE方法，首先采用“蒙特卡洛采樣”從“樣本空間”里產(chǎn)生子數(shù)據(jù)集，對(duì)每個(gè)變量的“回歸系數(shù)”進(jìn)行“均值方差比”統(tǒng)計(jì)分析來評(píng)價(jià)變量重要性。用于光譜波段選擇及QSAR描述符選擇的方法有：MC-UVE，CARS，IRIV，VISSA，VCPA，RT，Random Frog。用于代謝組學(xué)的生物標(biāo)記物選擇的方法有VIAVC，VISSA，SPA，MMIFS，CARS，Random frog，VCN。用于基因組學(xué)及蛋白組學(xué)特征選擇的方法有VIAVC，Random frog，MIA。

4.3 模型參數(shù)與評(píng)價(jià)

模型參數(shù)與評(píng)價(jià)是化學(xué)建模研究的一個(gè)基礎(chǔ)問題[47]，任何模型的研究都離不開模型評(píng)價(jià)。目前，有很多模型評(píng)價(jià)信息標(biāo)準(zhǔn)，如AIC，BIC，DIC和Cp統(tǒng)計(jì)量等[48]。交互檢驗(yàn)是比較常用的模型評(píng)價(jià)方法，只得到一個(gè)值用來評(píng)價(jià)，如交互檢驗(yàn)預(yù)測(cè)誤差。然而，僅用預(yù)測(cè)誤差評(píng)價(jià)模型是不足的[49～51]。Deng等[52]提出了一種基于模型集群分析并結(jié)合模型預(yù)測(cè)能力與模型穩(wěn)定性評(píng)價(jià)模型的方法。該方法應(yīng)用于PLS回歸模型，以確定PLS潛變量個(gè)數(shù)這個(gè)參數(shù)。對(duì)于PLS回歸模型，模型穩(wěn)定性可以根據(jù)回歸系數(shù)的方差判斷。具體步驟如下：（1）采用蒙特卡洛采樣從原始數(shù)據(jù)總樣本中選取80%的樣本作為訓(xùn)練樣本，剩下20%的樣本作為獨(dú)立測(cè)試集樣本。這個(gè)過程重復(fù)N次，將去獲取N個(gè)子訓(xùn)練集和與之對(duì)應(yīng)的N個(gè)子測(cè)試集；（2）每個(gè)子訓(xùn)練集建立模型并對(duì)相應(yīng)的測(cè)試集樣本進(jìn)行預(yù)測(cè)。同時(shí)，記錄每個(gè)子模型的PLS回歸系數(shù)；（3）根據(jù)模型預(yù)測(cè)誤差來獲得模型預(yù)測(cè)能力，同時(shí)根據(jù)PLS回歸系數(shù)來獲得模型穩(wěn)定性；（4）根據(jù)不同的PLS潛變量個(gè)數(shù)，重復(fù)步驟1～3，選擇同時(shí)具有好的預(yù)測(cè)能力及穩(wěn)定性的PLS模型。endprint

在此，選取常用的一組近紅外光譜測(cè)量角叉膠的數(shù)據(jù)[53]，該量測(cè)數(shù)據(jù)包含128個(gè)樣本，每條近紅外光譜包含701個(gè)數(shù)據(jù)點(diǎn)。PLS作為校正模型的方法。留一交互檢驗(yàn)（LOOCV），五折交互檢驗(yàn)（5-fold CV）和蒙特卡洛交互檢驗(yàn)（MCCV）的結(jié)果顯示，最優(yōu)的潛變量數(shù)很難確定，因?yàn)椴煌瑵撟兞繑?shù)的模型有非常接近的交互檢驗(yàn)均方根誤差（圖6A）。然而，從模型穩(wěn)定性的角度我們可以發(fā)現(xiàn)潛變量為6的模型穩(wěn)定性明顯高于其它潛變量數(shù)的模型，如圖6B所示。模型的穩(wěn)定性在這里用回歸系數(shù)之間的歐式距離衡量。潛變量數(shù)為6的模型回歸系數(shù)之間的歐式距離的值明顯小于潛變量數(shù)為1和20的模型，歐式距離的分布也更集中。值得注意的是，潛變量數(shù)為20的模型比潛變量數(shù)為6的模型預(yù)測(cè)誤差稍小。但是，結(jié)合模型穩(wěn)定性，選擇的最優(yōu)潛變量數(shù)是6。

圖6 （A）不同方法確定潛變量數(shù)（nLV）的結(jié)果，藍(lán)色代表留一法交互檢驗(yàn)（LOOCV），紫紅色代表5折交互檢驗(yàn)（5-fold CV），淺藍(lán)色代表蒙特卡洛交互檢驗(yàn)（MCCV），紅色代表基于模型集群分析的方法（MPA-based）； （B）不同模型復(fù)雜度的回歸系數(shù)之間的歐氏距離

Fig.6 （A）Determination of the number of latent variables by different methods， blue represents leave one out cross validation， purplish red represents five-fold cross validation， light blue represents Monte Carlo cross validation and red represents MPA-based. （B） Euclidean distance between PLS regression coefficient on different model complexities （The number of latent variable， nLV=1， 6 and 20）

4.4 穩(wěn)健與模型應(yīng)用域

奇異樣本診斷往往應(yīng)用于建立模型前去除奇異樣本，而模型應(yīng)用域則是在模型建立后在應(yīng)用上需要定義的，是化學(xué)建模中至關(guān)重要的一步，決定著建立好的模型的應(yīng)用范圍。給定一個(gè)建立好的模型，對(duì)于需要預(yù)測(cè)的外來新樣本，其與模型應(yīng)用域的關(guān)系存在3種情況：（1）新樣本在模型應(yīng)用域內(nèi)，即所建模型考慮到了該樣本的信息，可被很好的預(yù)測(cè)，即預(yù)測(cè)誤差??；（2）新樣本處在應(yīng)用域邊緣，即模型只考慮了該樣本的部分信息，其可被預(yù)測(cè)但精度不高，預(yù)測(cè)誤差較大；（3）新樣本完全處在模型應(yīng)用域外，即所建模型完全沒有考慮了該樣本的任何信息，因此該樣本無法被準(zhǔn)確預(yù)測(cè)，預(yù)測(cè)誤差極大。而目前有很多模型應(yīng)用域的方法[54]，基于范圍和幾何原理的方法[55]，基于主成份分析的方法[56]，基于凸包原理（Convex Hulls）的方法[57]，基于化學(xué)相似性的方法[58]，基于概率密度的方法[59]和基于模型集群分析的方法[60，61]。

基于模型集群分析的方法，通過在樣本空間或者變量空間隨機(jī)得到多個(gè)子模型，并且對(duì)子模型的預(yù)測(cè)方差進(jìn)行統(tǒng)計(jì)分析，從而確定模型的應(yīng)用域[60]，其步驟如下：（1）采用蒙特卡洛采樣從原始數(shù)據(jù)總樣本中選取一定比列的樣本作為訓(xùn)練樣本，如r=80%，剩下20%的樣本作為獨(dú)立測(cè)試集樣本。這個(gè)過程重復(fù)N次，將去獲取N個(gè)子訓(xùn)練集和與之對(duì)應(yīng)的N個(gè)子測(cè)試集；（2）每個(gè)子訓(xùn)練集建立模型并對(duì)相應(yīng)的測(cè)試集樣本進(jìn)行預(yù)測(cè)，即可得到N個(gè)測(cè)試結(jié)果；（3）統(tǒng)計(jì)每個(gè)測(cè)試樣本的預(yù)測(cè)誤差的分布，計(jì)算每個(gè)樣本預(yù)測(cè)誤差的標(biāo)準(zhǔn)差。在此選取Hou等[62]報(bào)道的一組QSAR數(shù)據(jù)進(jìn)行分析。該數(shù)據(jù)包含1290個(gè)化合物 （樣本），324個(gè)分子描述符 （變量）。首先，將樣本劃分為訓(xùn)練集（411個(gè)），測(cè)試集1（410個(gè)）和測(cè)試集2 （466個(gè)）。訓(xùn)練集和測(cè)試集1的化合物只包含C， H， O和N元素，而測(cè)試集2的化合物除了含有C， H， O和N元素外，還含有其它元素。 用訓(xùn)練集進(jìn)行建模，可以預(yù)計(jì)測(cè)試集1的樣本在模型應(yīng)用域里面，而測(cè)試集2的樣本則在模型應(yīng)用域之外。結(jié)果如圖7a所示，測(cè)試集2的樣本（紫紅色）預(yù)測(cè)誤差的標(biāo)準(zhǔn)差明顯大于訓(xùn)練集（藍(lán)色）和

圖7 應(yīng)用域描述圖：a，預(yù)測(cè)誤差值與其標(biāo)準(zhǔn)差值之間的關(guān)系，藍(lán)色圓圈代表訓(xùn)練集，綠色菱形代表測(cè)試集1，紫紅色星號(hào)代表測(cè)試集2； b，兩個(gè)典型的樣本（a圖的A和B）的預(yù)測(cè)誤差分布圖，A樣本在模型應(yīng)用域里，B樣本在模型應(yīng)用域外

Fig.7 Description of applicability domain： （a） The relationship between the standard deviation values and the values of prediction error， green rhombus denote training data； blue circles denote test data 1； and purplish red star denote test data 2； （b） The distrubution of prediciton errors for two selected sampels marked in the left panel. Sample A is inside the applicability domain and sample B is outside the applicability domain測(cè)試集1的樣本（綠色），說明該方法能夠很好地劃分模型的應(yīng)用域。圖7b呈現(xiàn)的分別是測(cè)試集1和測(cè)試集2中的兩個(gè)典型樣本A和B的預(yù)測(cè)誤差的分布。樣本A為模型應(yīng)用域內(nèi)的樣本，而樣本B為模型應(yīng)用域外的樣本。樣本B的預(yù)測(cè)誤差明顯大于樣本A，同時(shí)樣本B的預(yù)測(cè)的變化范圍也遠(yuǎn)大于樣本A。endprint

5 結(jié)論與展望

值得指出的是，模型集群分析的主要思路是從多個(gè)視點(diǎn)看待化學(xué)建模，并通過多次建模以盡量逼近建模基空間；同時(shí)，通過多個(gè)子模型比較，以避免模型的過擬合或其他建模陷阱，為化學(xué)建模提供了新思路。另一方面，模型集群分析實(shí)際是提供了一種一般性研究框架，可以從3個(gè)基本要素4個(gè)空間里選擇改變某一點(diǎn)或幾點(diǎn)作為切入口，開發(fā)一種應(yīng)用于化學(xué)建模的新算法。也就是，在樣本空間、變量空間、參數(shù)空間或模型空間的任何一個(gè)空間或多個(gè)空間，改變隨機(jī)采樣方法，改變建立子模型的方法，以及用不同統(tǒng)計(jì)方法分析不同的感興趣的輸出參數(shù)。同時(shí)，模型集群分析為化學(xué)建模在奇異值診斷，變量選擇，模型參數(shù)與評(píng)價(jià)以及模型應(yīng)用域的算法設(shè)計(jì)上提供了一種新的思維方式，為解決復(fù)雜多組分體系的高通量分析提供了新手段。在這里，我們討論的基于模型集群分析的方法及應(yīng)用是有限的，實(shí)際上，模型集群分析的泛化性很強(qiáng)，可以延伸到其它領(lǐng)域，如基因組學(xué)、蛋白組學(xué)等組學(xué)分析，以及生物信息學(xué)領(lǐng)域。未來可基于模型集群分析從基因和蛋白組學(xué)等大數(shù)據(jù)挖掘方向開發(fā)化學(xué)建模新算法，解決當(dāng)今大數(shù)據(jù)時(shí)代急需的計(jì)算應(yīng)用，但開發(fā)化學(xué)建模新算法同時(shí)應(yīng)注重算法的應(yīng)用性，應(yīng)以能解決實(shí)際問題為出發(fā)點(diǎn)。

References

1 LIANG Yi-Zeng， XU Qing-Song. Instrumental Analysis of Complex Systems-White， Grey and Black Analytical Systems and Their Multivariate Methods. Beijing： Chemical Industry Press， 2013： 1-18

梁逸曾， 許青松. 復(fù)雜體系儀器分析─白、灰、黑分析體系及其化學(xué)計(jì)量學(xué)算法. 北京： 化學(xué)工業(yè)出版社， 2013： 1-18

2 Cawley G C， Talbot N L C. Bioinformatics， 2006， 22（19）： 2348-2355

3 Chen T， Martin E. Anal. Chim. Acta.， 2009， 631（1）： 13-21

4 Wang Q， Li H D， Xu Q S， Liang Y Z. Analyst， 2011， 136（7）： 1456-1463

5 Yeung K Y， Bumgarner R E， Raftery A E. Bioinformatics， 2005， 21（10）： 2394-2402

6 Candes E， Tao T. Ann. Stat.， 2007： 2313-2351

7 Johnstone I M， Titterington D M. Philos. Trans. A. Math. Phys. Eng. Sci.， 2009， 367（1906）： 4237-4253

8 Zou H， Hastie T. J. Roy. Stat. Soc. B.， 2005， 67（2）： 301-320

9 Li H D， Liang Y Z， Xu Q S， Cao D S. J. Chemometr， 2010， 24： 418-423

10 Li H D， Liang Y Z， Long X X， Yun Y H， Xu Q S. Chemometr. Intell. Lab.， 2013， 122： 23-30

11 Li H D， Liang Y Z， Cao D S， Xu Q S. TRAC-Trend. Anal. Chem.， 2012， 38： 154-162

12 Efron B， Efron B. The jackknife， the bootstrap and other resampling plans. SIAM， 1982： 1-92

13 Miller R G. Biometrika， 1974， 61（1）： 1-15

14 Efron B， Tibshirani R J. An Introduction to the Bootstrap.， Boca Raton： CRC Press， 1994： 1-404

15 Deng B C， Yun Y H， Liang Y Z， Yi L Z. Analyst， 2014， 139（19）： 4836-4845

16 Deng B C， Yun Y H， Ma P， Lin C C， Ren D B， Liang Y Z. Analyst， 2015， 140（6）： 1876-1885

17 Yun Y H， Wang W T， Deng B C， Lai G B， Liu X B， Ren D B， Liang Y Z， Fan W， Xu Q S. Anal. Chim. Acta.， 2015， 862： 14-23

18 Yun Y H， Wang W T， Tan M L， Liang Y Z， Li H D， Cao D S， Lu H M， Xu Q S. Anal. Chim. Acta.， 2014， 807： 36-43

19 Zhang H Y， Wang H Y， Dai Z J， Chen M S， Yuan Z M. BMC Bioinformatics.， 2012， 13（1）： 298-317

20 Edgington E， Onghena P. Randomization tests. Boca Raton： CRC Press， 2007： 1-998endprint

21 Mann H B， Whitney D R. Ann. Math. Statist.， 1947， 18（1）： 50-60

22 Egan W J， Morgan S L. Anal. Chem.， 1998， 70（11）： 2372-2379

23 Cao D S， Liang Y Z， Xu Q S， Li H D， Chen X. J. Comput. Chem.， 2010， 31（3）： 592-602

24 Fan J， Li R： Statistical Challenges with High Dimensionality： Feature Selection in Knowledge Discovery. In： Proceedings of the Madrid International Congress of Mathematicians： 2006； Madrid

25 Cai T T， Shen X. High-Dimensional Data Analysis. Beijing： Higher Education Press， 2010： 119-145

26 Spiegelman C H， McShane M J， Goetz M J， Motamedi M， Yue Q L， Coté G L. Anal. Chem.， 1998， 70（1）： 35-44

27 Yun Y H， Liang Y Z， Xie G X， Li H D， Cao D S， Xu Q S. Analyst， 2013， 138（21）： 6412-6421

28 Guyon I， Elisseeff A. J. Mach. Learn. Res.， 2003， 3： 1157-1182

29 Yun Y H， Liang F， Deng B C， Lai G B， Vicente Gonalves C， Lu H M， Yan J， Huang X， Yi L Z， Liang Y Z. Metabolomics.， 2015， doi：10.1007/s11306-015-0803-x

30 Zeng M M， Liang Y Z， Li H D， Wang M， Wang B， Chen X， Zhou N， Cao D S， Wu J. J. Pharmaceut. Biomed.， 2010， 52（2）： 265-272

31 Marrocco C， Duin R P W， Tortorella F. Pattern. Recogn.， 2008， 41（6）： 1961-1974

32 Zweig M H， Campbell G. Clin. Chem.， 1993， 39（4）： 561-577

33 Hulver M W， Berggren J R， Cortright R N， Dudek R W， Thompson R P， Pories W J， MacDonald K G， Cline G W， Shulman G I， Dohm G L et al. Am. J. Physiol. Endocrinol. Metab.， 2003， 284（4）： 741-747

34 Kien C L， Bunn J Y， Ugrasbul F. Am. J. Clin. Nutr.， 2005， 82（2）： 320-326

35 Proenza A M， Roca P， CrespíC， Llad I， Palou A. J. Nutr. Biochem.， 1998， 9（12）： 697-704

36 Cai W S， Li Y K， Shao X G. Chemometr. Intell. Lab.， 2008， 90（2）： 188-194

37 Li H D， Liang Y Z， Xu Q S， Cao D S. Anal. Chim. Acta.， 2009， 648（1）： 77-84

38 Zheng K Y， Li Q Q， Wang J J， Geng J P， Cao P， Sui T， Wang X， Du Y P. Chemometr. Intell. Lab.， 2012， 112： 48-54

39 Li H D， Liang Y Z， Xu Q S， Cao D S， Tan B B， Deng B C， Lin C C. Ieee. Acm. T. Comput. Bi.， 2011， 8（6）： 1633-1641

40 Deng B C， Yun Y H， Liang Y Z， Yi L Z. Analyst.， 2014， 139（19）： 4836-4845

41 Li H D， Xu Q S， Liang Y Z. Anal. Chim. Acta.， 2012， 740： 20-26

42 Yun Y H， Li H D， E. Wood L R， Fan W， Wang J J， Cao D S， Xu Q S， Liang Y Z. Spectrochim. Acta. A.， 2013， 111： 31-36

43 Long X X， Li H D， Fan W， Xu Q S， Liang Y Z. Chemometr. Intell. Lab.， 2013， 121： 75-81

44 Xu H， Liu Z C， Cai W S， Shao X G. Chemometr. Intell. Lab.， 2009， 97（2）： 189-193endprint

45 Li H D， Xu Q S， Zhang W， Liang Y Z. Metabolomics.， 2012， 8（6）： 1218-1226

46 Li H D， Zeng M M， Tan B B， Liang Y Z， Xu Q S， Cao D S. Metabolomics.， 2010， 6（3）： 353-361

47 Gramatica P. Qsar. Comb. Sci.， 2007， 26（5）： 694-701

48 Akaike H. IEEE. T. Automat. Contr.， 1974， 19（6）： 716-723

49 Breiman L. Mach. Learn.， 2001， 45（1）： 5-32

50 Varma S， Simon R. BMC Bioinformatics.， 2006， 7（1）： 91-98

51 Hawkins D M， Basak S C， Mills D. J. Chem. Inf. Comp. Sci.， 2003， 43（2）： 579-586

52 Deng B C， Yun Y H， Liang Y Z， Cao D S， Xu Q S， Yi L Z， Huang X. Anal. Chim. Acta.， 2015， 880： 32-41

53 Dyrby M， Petersen R V， Larsen J， Rudolf B， Nrgaard L， Engelsen S B. Carbohyd. Polym.， 2004， 57（3）： 337-348

54 Dimitrov S， Dimitrova G， Pavlov T， Dimitrova N， Patlewicz G， Niemela J， Mekenyan O. J. Chem. Inf. Model.， 2005， 45（4）： 839-849

55 Sahigara F， Mansouri K， Ballabio D， Mauri A， Consonni V， Todeschini R. Molecules.， 2012， 17（5）： 4791-4810

56 Wold S， Esbensen K， Geladi P. Chemometr. Intell. Lab.， 1987， 2（1-3）： 37-52

57 Preparata F， Shamos M： Convex Hulls： Basic Algorithms. In： Computational Geometry. Springer New York； 1985： 95-149

58 Netzeva T I， Worth A P， Aldenberg T， Benigni R， Cronin M T， Gramatica P， Jaworska J S， Kahn S， Klopman G， Marchant C A. ATLA. 2005， 33： 155-173

59 Jaworska J， Nikolova-Jeliazkova N， Aldenberg T. ATLA-NOTTINGHAM-. 2005， 33（5）： 445-459

60 Kaneko H， Funatsu K. J. Chem. Inf. Model.， 2014， 54（9）： 2469-2482

61 Yan J， Zhu W W， Kong B， Lu H B， Yun Y H， Huang J H， Liang Y Z. Mol. Inform.， 2014， 33（8）： 503-513

62 Hou T， Xia K， Zhang W， Xu X. J. Chem. Inf. Comp. Sci.， 2004， 44（1）： 266-275endprint