史海芳, 李 聰, 姬永剛

(1. 中國(guó)民航大學(xué) 理學(xué)院, 天津 300300; 2. 吉林大學(xué) 數(shù)學(xué)學(xué)院, 長(zhǎng)春 130012)

在劑量反應(yīng)實(shí)驗(yàn)中, 時(shí)間可作為解釋變量, 人們可能在收集數(shù)據(jù)前就已獲得了潛在參數(shù)滿(mǎn)足序約束的先驗(yàn)信息. 將這些先驗(yàn)信息應(yīng)用到統(tǒng)計(jì)推斷中可以提供更有效的估計(jì), 特別是在具有低信噪比的小樣本情形下. 文獻(xiàn)[1-2]對(duì)序約束統(tǒng)計(jì)方法進(jìn)行了全面闡述.

目前, 利用Bayes方法解決約束條件下的統(tǒng)計(jì)推斷問(wèn)題已取得許多成果. 文獻(xiàn)[3]借助Savage-Dickey密度比方法[4]討論了序約束下單向分類(lèi)方差分析(ANOVA)模型的假設(shè)檢驗(yàn)問(wèn)題, 利用該方法可近似計(jì)算各潛在的假設(shè)Bayes因子, 并選擇后驗(yàn)概率最大的假設(shè)作為最優(yōu)假設(shè)； 文獻(xiàn)[5]將該方法推廣到其有不等式約束條件的廣義線(xiàn)性模型中, 但Savage-Dickey密度比方法需要在每次Markov鏈Monte Carlo(MCMC)迭代中都近似計(jì)算標(biāo)準(zhǔn)正態(tài)分布的密度和分布函數(shù)值, 計(jì)算時(shí)間較長(zhǎng)； 文獻(xiàn)[6-7]通過(guò)引進(jìn)指標(biāo)變量, 并假設(shè)指標(biāo)變量與參數(shù)相互獨(dú)立, 分別將文獻(xiàn)[8]中的方法應(yīng)用到單向分類(lèi)和雙向分類(lèi)方差分析模型中, 考慮了序約束下的變量選擇問(wèn)題, 其優(yōu)點(diǎn)是沒(méi)有過(guò)多的調(diào)節(jié)參數(shù), 但文獻(xiàn)[9]指出若參數(shù)先驗(yàn)分布的假設(shè)太模糊, 則可能會(huì)導(dǎo)致Markov鏈?zhǔn)諗枯^慢; 文獻(xiàn)[10-11]提出了另一種常用的Bayes變量選擇方法SSVS(stochastic search variable selection)方法, 目前SSVS已被擴(kuò)展到很多模型和應(yīng)用中, 例如遺傳數(shù)據(jù)分析[12]、 向量自回歸模型[13]和分組約束模型[14]等. 該方法的缺點(diǎn)是需要進(jìn)行調(diào)節(jié)的參數(shù)較多. 為解決上述問(wèn)題, 文獻(xiàn)[15]提出了一種改進(jìn)的SSVS方法----NMIG(normal mixture of interse Gamma distributions)方法, 并將其應(yīng)用到基因分析中. 本文將這兩種SSVS方法推廣到單向分類(lèi)方差分析模型中， 并考慮序約束下的變量選擇問(wèn)題. 基于MCMC方法引入指標(biāo)變量, 計(jì)算潛在模型的后驗(yàn)概率.

1 再參數(shù)化ANOVA模型中的變量選擇

令

yij=μi+εij,j=1,2,…,ni,i=1,2,…,k

(1)

表示一個(gè)單向ANOVA混合模型, 其中yij表示第j個(gè)個(gè)體第i次治療的反映變量,μi表示第i次治療的平均效應(yīng),εij表示誤差項(xiàng)且εij～N(0,σ2).實(shí)際應(yīng)用中一般有如下幾種序約束：

(i) 簡(jiǎn)單半序約束μ1≤…≤μk;

(ii) 簡(jiǎn)單樹(shù)序約束μ1≥μi,i=2,…,k-1;

(iii) 傘序約束μ1≤…≤μg≥μg+1≥…≥μk, 其中g(shù)已知.

對(duì)于簡(jiǎn)單半序約束μ1≤…≤μk, 如果令δm-1=μm-μm-1(2≤m≤k), 則第二個(gè)均值μ2可表示為μ1+δ1, 第三個(gè)均值可表示為μ1+δ1+δ2,…, 第m個(gè)均值(2≤m≤k)可表示為μ1+δ1+…+δm-1,μm的簡(jiǎn)單序約束可等價(jià)地表示為δm的非負(fù)約束δm≥0,m=1,2,…,k-1.如果令δ=(δ1,δ2,…,δk-1), 則滿(mǎn)足簡(jiǎn)單半序約束的模型(1)可表示為

y=1nμ1+xδ+ε,δ1≥0, …,δk-1≥0,

(2)

這里1m和0m分別表示元素皆為1和0的m×1維列向量.若均值滿(mǎn)足簡(jiǎn)單樹(shù)序約束或傘序約束, 則本文也可定義相應(yīng)的設(shè)計(jì)矩陣x及再參數(shù)化向量參數(shù)θ.這些不等式約束可參見(jiàn)文獻(xiàn)[3,5].盡管本文給出的方法適用于任何可表示為模型(2)的參數(shù)約束, 但為方便, 本文主要考慮簡(jiǎn)單半序約束.

為能同時(shí)進(jìn)行變量選擇和參數(shù)估計(jì), 本文考慮兩種Bayes變量選擇方法： SSVS方法[10-11]和NMIG方法[15]. 本文將SSVS方法擴(kuò)展到帶有序約束條件的單向ANOVA模型中.引入指示變量γi, 并假設(shè)γi和δi滿(mǎn)足下列先驗(yàn)分布:

δi|γi,φi～TN(0,c(γi)φi,0,+∞),

(3)

(4)

其中:

(5)

π(σ2)∝1,

文獻(xiàn)[16]將該先驗(yàn)分布應(yīng)用于無(wú)約束線(xiàn)性模型的變量選擇問(wèn)題中.

上述先驗(yàn)和超先驗(yàn)分布對(duì)于其各自參數(shù)和超參數(shù)均為共軛的, 因此利用Gibbs抽樣方法易得參數(shù)的后驗(yàn)分布.下面分別給出這兩種方法所有參數(shù)的滿(mǎn)條件分布.

1.1 SSVS方法

1)σ2的條件后驗(yàn)分布.由于

(6)

因此

σ2|y,x,μ1,{γl},{δl}～I(xiàn)Gamma(n/2-1,(y-1nμ1-xθ)′(y-1nμ1-xθ)/2).

2)μ1的滿(mǎn)條件后驗(yàn)分布.由于

3)γi的滿(mǎn)條件后驗(yàn)分布.令γ-i=(γ1,…,γi-1,γi+1,…,γq), 可證明γi的滿(mǎn)條件分布服從如下Bernoulli分布:

(7)

其中

ci=f(y|{δi},σ2,γi=1,γ-i)f(δi|γi=1,γ-i)f(γi=1,γ-i),

di=f(y|{δi},σ2,γi=0,γ-i)f(δi|γi=0,γ-i)f(γi=0,γ-i).

4)πi的滿(mǎn)條件后驗(yàn)分布.由πi和γi的先驗(yàn)可計(jì)算πi的后驗(yàn)分布為

πi|y,x,{βi},σ2,{γi},μ1～Beta(a1+γi,b1-γi+1).

5)δp的滿(mǎn)條件后驗(yàn)分布.由于

1.2 NMIG方法

2 數(shù)值模擬

下面用數(shù)值模擬驗(yàn)證本文方法的有效性, 并與文獻(xiàn)[3]和文獻(xiàn)[6]提出的方法進(jìn)行對(duì)比, 將這兩種對(duì)比方法分別簡(jiǎn)記為Oh方法[3]和Otava方法[6]. 假設(shè)k=4, 且均值滿(mǎn)足簡(jiǎn)單半序約束, 則有下列8個(gè)備選模型：

(8)

根據(jù)文獻(xiàn)[3], 假設(shè)誤差項(xiàng)服從獨(dú)立的標(biāo)準(zhǔn)正態(tài)分布, 考慮以下3種模擬:

模擬1:μ=(0,0,0,0)′; 模擬2:μ=(0,0,0,1)′; 模擬3:μ=(1,2,3,4)′.

設(shè)樣本個(gè)數(shù)n=10,30,100, 每次產(chǎn)生10 000個(gè)Gibbs樣本, 為保證收斂, 丟棄前面的3 000個(gè)抽樣值.每種固定效應(yīng)和樣本個(gè)數(shù)的取值都重復(fù)500次.表1～表3列出了所有模型的平均后驗(yàn)概率.其中標(biāo)*處對(duì)應(yīng)真實(shí)模型的后驗(yàn)概率.

表1 模擬1中不同方法所得模型的平均后驗(yàn)概率

表2 模擬2中不同方法所得模型的平均后驗(yàn)概率

表3 模擬3中不同方法所得模型的平均后驗(yàn)概率

由表1～表3可見(jiàn), 4種方法給出的真實(shí)模型后驗(yàn)概率均為最大. 一般地, NMIG和Otava比其他兩種方法更傾向于提供真實(shí)模型更高的平均后驗(yàn)概率. 為評(píng)價(jià)這些方法在識(shí)別正確模型方面的效果, 類(lèi)似傳統(tǒng)假設(shè)檢驗(yàn)的優(yōu)勢(shì), 本文計(jì)算了每種方法選擇正確模型的頻率P, 結(jié)果列于表4～表6. 由表4～表6可見(jiàn), 在模擬3中Oh方法效果更好, 但在模擬1和模擬2中另外3種方法的正確模型識(shí)別率更高. 數(shù)值模擬結(jié)果表明, 在多數(shù)情形下, NMIG方法比SSVS方法效果更好, 這可能是因?yàn)樵诩僭O(shè)參數(shù)φi服從Gamma分布的條件下調(diào)節(jié)參數(shù)比固定其為常數(shù)更好.

表4 模擬1中不同方法選擇正確模型的頻率

表6 模擬3中不同方法選擇正確模型的頻率

3 實(shí)例分析

下面用本文方法分析一組由文獻(xiàn)[3]分析的實(shí)際數(shù)據(jù). 實(shí)驗(yàn)人員在18個(gè)月內(nèi)測(cè)量了20名男孩在8歲、 8.5歲、 9歲、 9.5 歲時(shí)的支骨高度(ramus bone heights). 令μ1,μ2,μ3，μ4分別表示4次觀(guān)測(cè)的平均支骨高度. 根據(jù)先驗(yàn)信息, Oh方法假設(shè)平均支骨高度滿(mǎn)足簡(jiǎn)單半序約束μ1≤μ2≤μ3≤μ4, 并利用Savage-Dickey密度比方法進(jìn)行了分析. 本文利用SSVS方法和NMIG方法對(duì)該數(shù)據(jù)進(jìn)行分析, 并與Oh方法和Otava方法進(jìn)行比較.考慮式(8)模型H0～HF, 假設(shè)先驗(yàn)的超參數(shù)取值與數(shù)值模擬相同, 并產(chǎn)生40 000個(gè)Gibbs樣本, 將前20 000個(gè)樣本作為初始值. 表7列出了所有模型的后驗(yàn)概率. 由表7可見(jiàn), NMIG方法和Otava方法選擇了相同的最大后驗(yàn)?zāi)Ｐ? 同時(shí)4種方法均認(rèn)為8.5歲時(shí)的平均支骨高度與9歲時(shí)的平均支骨高度有差異.表8和表9分別列出了參數(shù)δi(i=1,2,3)的后驗(yàn)均值和可信區(qū)間.由表8和表9可見(jiàn), Otava方法給出的后驗(yàn)均值更大, 而SSVS方法給出了長(zhǎng)度更短的Bayes可信區(qū)間, 表明SSVS方法在該數(shù)據(jù)集中表現(xiàn)較好. 圖1和圖2分別為SSVS方法和NMIG方法中參數(shù)δi(i=1,2,3)的后驗(yàn)密度估計(jì).由圖1和圖2可見(jiàn), 兩種方法都很好地將參數(shù)δi(i=1,2,3)控制在非負(fù)區(qū)間上.

綜上所述, 本文討論了序約束下單向分類(lèi)方差分析模型中的變量選擇問(wèn)題, 提出了兩種基于SSVS方法序約束下的Bayes變量選擇方法, 并利用數(shù)值模擬和實(shí)例分析驗(yàn)證了本文方法的有效性. 本文提出的后驗(yàn)抽樣方法較簡(jiǎn)單, 且易操作.

表7 支骨高度數(shù)據(jù)中模型的后驗(yàn)概率

表8 支骨高度數(shù)據(jù)中參數(shù)的后驗(yàn)均值

表9 支骨高度數(shù)據(jù)中參數(shù)的95%可信區(qū)間

圖1 SSVS方法參數(shù)δ1,δ2,δ3的后驗(yàn)密度估計(jì)Fig.1 A posteriori density estimation of parameters δ1,δ2,δ3 in SSVS method

圖2 NMIG方法參數(shù)δ1,δ2,δ3的后驗(yàn)密度估計(jì)Fig.2 A posteriori density estimation of parameters δ1,δ2,δ3 in NMIG method