劉瑞銀, 周志慧, 杜 歡

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

0 引 言

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

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

1 保序回歸

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

1.1 常見的半序關系

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

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

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

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

常見的半序關系有以下幾種形式:

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

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

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

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

其中簡單半序是傘型半序當h=k時的一種特殊情況。

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是一個給定的權函數(shù),其中ωi>0,i=1,2,…,k。如果函數(shù)g*∈G且滿足式(1)

(1)

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

2 保序回歸的算法步驟

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

2.1 PAVA算法步驟

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

步驟2 如果g?G,則必定存在一個下標i,有gi-1>gi。令B={i-1,i},則

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

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

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

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

2.2 傘型半序的算法步驟

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

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

例2 設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 傘型半序保序回歸的計算過程(a)

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

以上介紹了簡單半序與傘型半序的算法步驟,其他半序的算法可以通過簡單半序和傘型半序算法的相應變換得到。下面對2種算法進行程序實現(xiàn)。

3 程序實現(xiàn)

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

利用MATLAB軟件,編寫一個名為Iso的函數(shù)解決簡單半序型均值的估計問題[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表示對應的權重。這里繼續(xù)考慮實例1中的問題,在命令窗口中輸入:a=[4 2 10 6 6];w=[10 10 15 20 20]; iso(a,w)

則輸出的結果為

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

這與圖1中的結果是一致的。

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

利用Iso函數(shù)編寫一個名為Umb的函數(shù)計算傘型半序的保序回歸。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表示對應的權重;h表示最小值點。在這里考慮實例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)

則輸出的結果為

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中的結果是一致的。

4 在基因計算中的應用

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

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

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

第1步 選取曲線模式,將表達曲線模式記為C1,C2,…,Ch。對每個基因g=1,2,…,G重復進行下述步驟。

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

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

第5步 進行多重假設檢驗:

原假設為平凡曲線模式,備擇假設為給定曲線模式的并。給定顯著性水平,利用上述步驟得到p值,分別控制不同錯誤度量標準進行多重檢驗,最終挑選出顯著基因。

5 結 語

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