保序回歸算法的MATLAB實(shí)現(xiàn)

劉瑞銀, 周志慧, 杜 歡

(沈陽師范大學(xué) 數(shù)學(xué)與系統(tǒng)科學(xué)學(xué)院, 沈陽 110034)

0 引 言

20世紀(jì)中葉以來,序約束下的統(tǒng)計(jì)推斷在統(tǒng)計(jì)分析中占有十分重要的地位,至今已建立了較為完善的理論[1-2]。在這一領(lǐng)域中,通常研究序約束下求解未知參數(shù)的極大似然估計(jì)。然而,傳統(tǒng)的極大似然估計(jì)所得到的結(jié)果會使總體均方誤差相對較大,缺點(diǎn)明顯。文獻(xiàn)[3-5]中指出,保序回歸和一般的極大似然估計(jì)相比,具有較小的總體均方誤差。因此,研究保序回歸很有意義。

保序回歸在生物醫(yī)學(xué)方面應(yīng)用廣泛。在藥物劑量反應(yīng)中,隨著體內(nèi)藥物劑量或濃度的增加,藥效和藥物的毒性也會隨之增加。由于某些藥物具有毒副作用,因而就會使得藥效呈現(xiàn)一種先升后降的傘型趨勢。可以利用每個(gè)劑量水平下病人毒性反應(yīng)的比率來估計(jì)不同劑量水平下的毒性概率[6],但是這種情況下所估計(jì)出的毒性概率可能不是劑量水平的非減函數(shù),此時(shí)可以采用保序回歸方法來解決這個(gè)問題。

1 保序回歸

假設(shè)某種藥物的使用劑量為X,病人對該藥物的反應(yīng)值為Y,劑量的取值范圍為1,2,…,100,對應(yīng)的真實(shí)反應(yīng)值為Y=y1,y2,…,y100。正常情況下隨著藥物使用劑量的增加,病人對藥物的反應(yīng)值也應(yīng)該增加,但由于病人個(gè)體的原因,Y不一定是一個(gè)非減函數(shù)。如果按照藥物的真實(shí)反應(yīng)排序,對應(yīng)的X將會成為亂序,從而失去了研究意義。本文研究目的就是為了觀察隨著藥物使用劑量的增加,病人的平均反應(yīng)狀況也得到相應(yīng)序約束的變化。在這種情況下,利用保序回歸,可以在不改變X排列順序的條件下,得到滿足序約束關(guān)系的Y。

1.1 常見的半序關(guān)系

定義1[7]設(shè)Θ={θ1,θ2,…,θk}為有限集合,稱集合Θ上的一個(gè)關(guān)系“≤”為半序,如果滿足下列關(guān)系:

1)θi≤θj; 其中θi∈Θ;

2) 如果θi,θj,θk∈Θ,θi≤θj,θj≤θk,那么θi≤θk;

3) 如果θi,θj∈Θ,θi≤θj并且θj≤θi,那么θi=θj。

常見的半序關(guān)系有以下幾種形式:

簡單半序:θ1≤θ2≤…≤θk;

傘型半序:θ1≤…≤θh≥…≥θk;

簡單樹半序:θ1≤θj,j=2,…,k;

簡單環(huán)半序:θ1≤θj≤θk,j=2,…,k-1;

其中簡單半序是傘型半序當(dāng)h=k時(shí)的一種特殊情況。

1.2 保序回歸的定義

定義2 對于定義于Θ上的函數(shù)y=(y1,y2,…,yk)T,其中yi=y(θi),如果θi,θj∈Θ,θi≤θj,有y(θi)≤y(θj),則稱函數(shù)y為對于“≤”的保序函數(shù)[8]。

定義3 假定G是保序函數(shù)的全體,給定函數(shù)g=(g1,g2,…,gk),ω=(ω1,ω2,…,ωk)T是一個(gè)給定的權(quán)函數(shù),其中ωi>0,i=1,2,…,k。如果函數(shù)g*∈G且滿足式(1)

(1)

則稱函數(shù)g*為(g,ω)的保序回歸[9]。

2 保序回歸的算法步驟

本節(jié)主要介紹求簡單半序和傘型半序的保序回歸的算法步驟,其他半序的算法可以通過簡單半序和傘型半序算法的相應(yīng)變換得到。假設(shè)Θ={θ1,θ2,…,θk}是一個(gè)有限集合,G是保序函數(shù)的全體,g和ω為給定函數(shù)和權(quán)函數(shù)。首先考慮簡單半序:θ1≤θ2≤…≤θk,以下求簡單半序保序回歸的算法名為PAVA(pool-adjacent-violators)算法[10-11]。

2.1 PAVA算法步驟

步驟1 如果g∈G,則g*=g。

步驟2 如果g?G,則必定存在一個(gè)下標(biāo)i,有g(shù)i-1>gi。令B={i-1,i},則

圖1 簡單半序保序回歸的計(jì)算過程Fig.1 Calculation process of simple isotonic regression

例1[12]設(shè)i={1,…,5},給定g=(4,2,10,6,6),ω=(20,10,10,15,20),求出g簡單半序的保序回歸g*。計(jì)算過程如圖1所示。

其中3=(10×4+10×2)/(10+10);7.090 9=(35×7.714 3+20×6)/(35+20)。因?yàn)?<7.090 9,所以g*=(3,3,7.090 9,7.090 9,7.090 9)。

下面介紹傘型半序保序回歸的算法[13]。本文中介紹傘型半序的其中一種形式----先降后增,h點(diǎn)為最小值點(diǎn)時(shí)的情況。傘型半序的形式為θ1≥…≥θh≤…≤θk,h對應(yīng)傘型的最小值點(diǎn),當(dāng)h=1或h=k時(shí),傘型半序變?yōu)楹唵伟胄颉?/p>

2.2 傘型半序的算法步驟

步驟1 如果g∈G,則g*=g。

步驟5 重復(fù)步驟4,直到包含h的塊的加權(quán)平均數(shù)比其他所有數(shù)都小為止。

例2 設(shè)k={1,…,9},h=4,給定g=(10,8,9,6,4,7,6,8,9),ω1=…=ωk,求出g的傘型半序的保序回歸g*。

首先對(10,8,9)和(4,7,6,8,9)使用PAVA算法。

表1 傘型半序保序回歸的計(jì)算過程(a)

表2 傘型半序保序回歸的計(jì)算過程(b)

以上介紹了簡單半序與傘型半序的算法步驟,其他半序的算法可以通過簡單半序和傘型半序算法的相應(yīng)變換得到。下面對2種算法進(jìn)行程序?qū)崿F(xiàn)。

3 程序?qū)崿F(xiàn)

3.1 簡單半序的算法實(shí)現(xiàn)

利用MATLAB軟件,編寫一個(gè)名為Iso的函數(shù)解決簡單半序型均值的估計(jì)問題[14-15]。Iso函數(shù)的程序代碼如下:

function is=Iso(a,w)

n=length(a); is=zeros(1,n);uu=0;

while (uu

v=uu+1;b=cumsum(a(v:n).*w(v:n));d=cumsum(w(v:n));b=b./d;u=min(b);

m=find(b==u);m=m(1)+v-1; is(v:m)=u;uu=m;

end

end

其中:a表示待排序的變量;w表示對應(yīng)的權(quán)重。這里繼續(xù)考慮實(shí)例1中的問題,在命令窗口中輸入:a=[4 2 10 6 6];w=[10 10 15 20 20]; iso(a,w)

則輸出的結(jié)果為

ans=3.000 0 3.000 0 7.090 9 7.090 9 7.090 9

這與圖1中的結(jié)果是一致的。

3.2 傘型半序的算法實(shí)現(xiàn)

利用Iso函數(shù)編寫一個(gè)名為Umb的函數(shù)計(jì)算傘型半序的保序回歸。Umb函數(shù)的程序代碼如下:

functiona=Umb(x,w,h)

k=length(x);b=-x(1:(h-1));b=Iso(b,w(1:(h-1)));b=-b;a=x((h+1):k);

a=Iso(a,w((h+1):k));u=[b,a]; [u,INDEX]=sort(u);v=INDEX;ww=[w(1:(h-1)),w((h+1):k)];

fori=1:(k-1)

www(i)=ww(v(i));

end

ww=[w(h),www];a=[x(h),u];ap=cumsum(a.*ww);wp=cumsum(ww);dd=ap./wp;

d=min(dd);m=find(dd==d);

if(m(1)==k)

a=linspace(d,d,m(1));

else

a=[linspace(d,d,m(1)),a((m(1)+1):k)];

fori=1:(k-1)

u(v(i))=a(i+1);

end

a(1:(h-1))=u(1:h-1);a(h)=d;a((h+1):k)=u(h:(k-1));

end

其中:a表示待排序的變量;w表示對應(yīng)的權(quán)重;h表示最小值點(diǎn)。在這里考慮實(shí)例2中的問題,在命令窗口中輸入:

x=[10 8 9 6 4 7 6 8 9];w=[1 1 1 1 1 1 1 1 1];h=4; Umb(x,w,h)

則輸出的結(jié)果為

ans=10.000 0 8.500 0 8.500 0 5.000 0 5.000 0 6.500 0 6.500 0 8.000 0 9.000 0

這與表2中的結(jié)果是一致的。

4 在基因計(jì)算中的應(yīng)用

利用寡核苷酸陣列和c-DNA陣列可以獲得基因的表達(dá)數(shù)據(jù),利用獲取的數(shù)據(jù)可以發(fā)現(xiàn)人們感興趣的信息。本研究小組收集了1 900個(gè)基因在6個(gè)不同時(shí)間點(diǎn)的表達(dá)情況的數(shù)據(jù),研究這1 900個(gè)基因的表達(dá)模式,對于一些表達(dá)模式相似的基因,它們的功能也是大同小異的。反過來思考,如果有2個(gè)基因被同一個(gè)調(diào)控系統(tǒng)控制著,那么可以認(rèn)為它們的表達(dá)模式是相類似的。

為了進(jìn)行基因選擇,在研究過程中,首先要選定幾種感興趣的基因表達(dá)曲線模式,通過對觀測到的數(shù)據(jù)進(jìn)行計(jì)算,把表達(dá)模式屬于選定曲線模式的基因挑選出來。常見的基因曲線模式有平凡曲線模式、簡單序曲線模式、傘狀序曲線模式、循環(huán)曲線模式和分段不等式曲線模式等,然后通常要求出基因在每個(gè)模式下的均值的估計(jì)值。本文所提出的PAVA算法就可以解決這個(gè)問題。

基于本文研究,為了進(jìn)行基因選擇,便于后續(xù)多重假設(shè)檢驗(yàn)技術(shù)的研究,對ORIOGEN算法進(jìn)行了調(diào)整,改進(jìn)后的算法如下:

第1步 選取曲線模式,將表達(dá)曲線模式記為C1,C2,…,Ch。對每個(gè)基因g=1,2,…,G重復(fù)進(jìn)行下述步驟。

第2步 利用PAVA算法求出基因在每個(gè)模式Ci,i=1,2,…,h下的均值的估計(jì)值。

第4步 設(shè)在所有的時(shí)間點(diǎn)上基因真正的均值和方差均相等,每個(gè)基因的bootstrap樣本為N個(gè),抽取方法:將基因的全部觀測值合并成一個(gè)長度為MT的向量,允許有重復(fù)的抽取T個(gè)隨機(jī)樣本,然后通過步驟2和3來計(jì)算抽取樣本的統(tǒng)計(jì)量,最終獲得N個(gè)統(tǒng)計(jì)量求出該基因的p值;對每個(gè)基因重復(fù)上述步驟求出基因的p值。

第5步 進(jìn)行多重假設(shè)檢驗(yàn):

原假設(shè)為平凡曲線模式,備擇假設(shè)為給定曲線模式的并。給定顯著性水平,利用上述步驟得到p值,分別控制不同錯(cuò)誤度量標(biāo)準(zhǔn)進(jìn)行多重檢驗(yàn),最終挑選出顯著基因。

5 結(jié) 語

本文通過對保序回歸算法原理的分析,利用MATLAB軟件編寫出簡單半序和傘型半序的保序回歸函數(shù)并進(jìn)行模擬實(shí)現(xiàn)。在實(shí)際的統(tǒng)計(jì)問題中,許多情況下都要求所估計(jì)的參數(shù)本身滿足某種序關(guān)系,利用Iso和Umb函數(shù)可以簡單快速地估計(jì)出滿足這種序關(guān)系的參數(shù),本文為以后的其他相關(guān)保序回歸內(nèi)容的研究奠定了理論依據(jù)。