基于負(fù)二項分布的確切推斷在微核試驗數(shù)據(jù)分析中的模擬研究*

索瑞鑫 仇玉蘭 王 彤△

微核試驗是公認(rèn)的檢測染色體或有絲分裂器損傷的一種試驗方法，其最大的優(yōu)點是經(jīng)濟(jì)、簡單、快速。該實驗結(jié)果通常以嚙齒類動物的骨髓嗜多染紅細(xì)胞或者人外周淋巴細(xì)胞中的微核率來表示〔1〕。

國際上有多位生物統(tǒng)計學(xué)家對微核實驗的統(tǒng)計分析方法進(jìn)行過討論，但并未形成統(tǒng)一的推薦方法〔2〕。同時我們還注意到，EPA、OECD、ECC、UKEMS等國際上關(guān)于微核試驗的指南里也沒有明確統(tǒng)一的數(shù)據(jù)分析方法，大多都是類似于“應(yīng)使用合適的統(tǒng)計方法分析數(shù)據(jù)”這類陳述〔2〕。為了解目前對這類數(shù)據(jù)所采用的統(tǒng)計分析方法現(xiàn)況，我們利用PubMed-NLM文獻(xiàn)數(shù)據(jù)庫進(jìn)行了簡單的檢索。限定題目中包含“micronucleus”，年限從1995年到2010年，語言為英文，得到文獻(xiàn)總數(shù)959篇，可得到的全文文章152篇，涉及52種雜志。除去其中關(guān)于微核檢測技術(shù)、關(guān)聯(lián)性研究或僅作統(tǒng)計描述而非組間比較問題的文獻(xiàn)，最終103篇文章中使用t檢驗35篇，方差分析23篇，非參數(shù)秩和檢驗36篇，卡方檢驗13篇，F(xiàn)isher確切概率法7篇，Kastenbaum and Bowman檢驗4篇。這一匯總結(jié)果雖然不夠全面，但確實反映出在統(tǒng)計分析方法選擇上的不一致。

通常而言參數(shù)方法的統(tǒng)計功效會高于非參數(shù)方法，但需正確指定其分布假設(shè)。而關(guān)于微核數(shù)據(jù)的分布假設(shè)也有一些不同的討論，如二項分布、泊松分布或負(fù)二項分布等〔3〕。由于微核數(shù)的發(fā)生較少(尤其是包含對照組時)，所以多數(shù)學(xué)者傾向于可以用描述稀有事件發(fā)生數(shù)的泊松分布來假定。然而實際數(shù)據(jù)中泊松分布均數(shù)等于方差的條件有時不一定會被滿足，而負(fù)二項分布可以用離散參數(shù)來描述這種過離散現(xiàn)象，并且當(dāng)離散參數(shù)取極端值的時候，負(fù)二項分布就退化為泊松分布，所以將微核數(shù)據(jù)假定為負(fù)二項分布可以合理地考慮到動物個體間變異的客觀存在〔4〕。雖然利用隨機(jī)效應(yīng)模型也可以處理這種動物個體間變異，然而，較少的微核數(shù)并不滿足大樣本的漸近性條件(包括卡方檢驗等方法和基于廣義線性混合模型的方法)，故而考慮一些確切推斷方法是必要的?；诓此煞植肌⒍椃植嫉拇_切推斷方法存在已久，且在不同方面也有所進(jìn)展，但基于負(fù)二項分布的確切推斷方法最近才被實現(xiàn)〔5〕。故以下我們將介紹基于負(fù)二項分布的條件確切概率方法并通過對文獻(xiàn)中常見方法進(jìn)行統(tǒng)計模擬來比較其Ⅰ類錯誤和檢驗功效。

基于負(fù)二項分布的確切推斷

若Y服從均數(shù)為μ，離散參數(shù)為φ的負(fù)二項分布，記為Y～NB(μ，φ)，則其概率函數(shù)為:

Y的期望E(Y)=μ，方差 var(Y)=μ+φμ2。當(dāng) φ→0時，均數(shù)等于方差，負(fù)二項分布則退化為泊松分布。在微核數(shù)據(jù)中，這里的離散參數(shù)φ可描述動物個體間的變異。

在微核實驗中假設(shè)每個動物的微核數(shù)Y1，…，Yn相互獨立，服從負(fù)二項分布NB(μ=miλ，φ)，其中mi表示每個動物的細(xì)胞數(shù)，λ表示微核率，如果每個動物細(xì)胞數(shù)相等(如mi都等于2000)，則Y獨立同分布，此時λ的最大似然估計為微核總數(shù)除以細(xì)胞總數(shù)，與φ無關(guān)。

如果所有動物的細(xì)胞數(shù)都是一樣的(即mi≡m)，那么它們的計數(shù)之和也服從負(fù)二項分布，即Z=Y1+…+Yn～NB(nmλ，φn－1)。這里 λ 是充分統(tǒng)計量，λ與φ彼此獨立，可以通過讓(2)式的條件似然值取最大來確定φ值，從而通過(1)式得到推斷的確切概率〔6〕:

若ZtA和ZtB分別是A、B兩組的計數(shù)總和，動物數(shù)分別為 nA和 nB。則在無效假設(shè) ZtA～NB(nkmλ1，)，k∈A，B下，可以建立Fisher精確概率的列聯(lián)表，計數(shù)為ZtA+ZtB也是服從負(fù)二項分布的隨機(jī)變量，我們可以運用公式(1)負(fù)二項分布的精確概率計算公式得到比觀測值更極端的概率，即雙側(cè)的P值。

模擬研究

我們按照毒理學(xué)微核實驗有關(guān)的指南文獻(xiàn)提供的依據(jù)對涉及到的模擬參數(shù)進(jìn)行如下設(shè)定:動物數(shù)n的模擬參數(shù)設(shè)定為5到10，細(xì)胞數(shù)的模擬參數(shù)設(shè)定為2000，微核陽性數(shù)的總體差值δ為0到18。所比較的方法包括Poisson精確概率法、Fisher精確概率法、負(fù)二項精確概率法、t檢驗及其變換、確切概率秩和檢驗。

由于我們研究的是mi≡m，即每組細(xì)胞數(shù)相等的情況，離散參數(shù)φ是獨立的變量，在模擬研究中我們采用了固定離散參數(shù)φ的方法，其值分別設(shè)置為1和0.001，當(dāng)φ→1時，模擬數(shù)據(jù)服從負(fù)二項分布;φ=0.001時，模擬數(shù)據(jù)近似服從Poisson分布。

對Ⅰ類錯誤的模擬估計結(jié)果見表1，2。

表1 Ⅰ類錯誤及排序(φ→0)

表2 Ⅰ類錯誤及排序(φ→1)

從表1與表2可以看到無論φ為1還是接近0，Poisson精確概率法的Ⅰ類錯誤是最大的，且遠(yuǎn)大于規(guī)定的檢驗水準(zhǔn)α=0.05。φ→0，即數(shù)據(jù)來自于Poisson分布時，平方根變換t檢驗與平方根反正弦變換t檢驗以及t檢驗的Ⅰ類錯誤大致相同，都僅次于Poisson精確概率法;經(jīng)log變換t檢驗Ⅰ類錯誤是最小的;Fisher精確概率法與負(fù)二項精確概率法(NB)，及秩檢驗的Ⅰ類錯誤大致相同，均略高于log變換的t檢驗，但遠(yuǎn)遠(yuǎn)小于Poisson精確概率法的Ⅰ類錯誤。φ→1時，F(xiàn)isher精確概率法的Ⅰ類錯誤較高，僅次于Poisson精確概率法;經(jīng)平方根變換的t檢驗與平方根反正弦變換得到t檢驗的Ⅰ類錯誤大致相同，僅次于Fisher精確概率法;經(jīng)log變換的t檢驗Ⅰ類錯誤僅次于平方根變換的t檢驗;而秩檢驗、t檢驗、負(fù)二項精確概率法的Ⅰ類錯誤均小于規(guī)定的檢驗水準(zhǔn)α=0.05。由此結(jié)果可見，在φ→0時，Poisson精確概率法的Ⅰ類錯誤遠(yuǎn)遠(yuǎn)超過檢驗水準(zhǔn)，故不可接受;在φ→0時，Poisson精確概率法、Fisher精確概率法的Ⅰ類錯誤也遠(yuǎn)遠(yuǎn)超過檢驗水準(zhǔn)α=0.05，為不可接受的方法;其他方法尚可。

當(dāng)離散參數(shù)φ→0時，檢驗功效結(jié)果見表3。

當(dāng)φ→0時，分布近似于Poisson分布，此時Poisson精確概率法的檢驗功效最高;Fisher精確概率法、負(fù)二項精確概率法、t檢驗、平方根反正弦變換的t檢驗(arcsin變換t)、平方根變換的t檢驗、秩檢驗、log變換的t檢驗都隨著n和δ的增高呈現(xiàn)一定的波動，但總得來講，在δ=3到18時各個統(tǒng)計方法的檢驗功效排序趨于穩(wěn)定(圖1)，我們從表3中可以得到Poisson，F(xiàn)isher，NB方法的檢驗功效較大，其值也較相近;t檢驗、平方根反正弦變換t檢驗、平方根變換t檢驗、以及l(fā)og變換t檢驗的檢驗功效居中，其大小近似;而秩檢驗的檢驗功效較小。

當(dāng)離散參數(shù)φ→1時，檢驗功效結(jié)果見表4。

當(dāng)φ→1時，分布近似于負(fù)二項分布，此時Poisson精確概率法的檢驗功效最高;Fisher精確概率法的檢驗功效僅次于Poisson精確概率法;負(fù)二項精確概率法、t檢驗、平方根反正弦變換的t檢驗(arcsin變換t)、平方根變換的t檢驗、秩檢驗、log變換的t檢驗都隨著n和δ的增高有一定的波動性。我們以總體差別δ為橫坐標(biāo)，功效為縱坐標(biāo)作圖，發(fā)現(xiàn)除了Poisson和Fisher在δ為5時趨于穩(wěn)定外，其他方法都隨δ的變化較大，直到δ為12時才相對穩(wěn)定(圖1)，故我們以δ=12為例，從表4，圖2中可以發(fā)現(xiàn):負(fù)二項精確概率法(NB)的檢驗功效較大;平方根反正弦變換t檢驗、平方根變換t檢驗、以及秩檢驗的檢驗功效居中;而log變換t檢驗、t檢驗的檢驗功效較小。

表3 檢驗功效及排序(φ→0，δ=3)

表4 檢驗功效及排序(φ→1，δ=12)

圖1 δ與功效的關(guān)系圖(n=10)

圖2 不同δ值時n與功效示意

討 論

雖然 David J．Kirkland 等〔2〕已經(jīng)明確指出過未經(jīng)變換的t檢驗方法很多時候不適用于微核數(shù)據(jù)的研究分析，但用PubMed-NLM檢索近3年的關(guān)于微核實驗可獲全文的16篇文章中仍能查到5篇文獻(xiàn)使用此方法。我們的模擬研究結(jié)果表明功效隨著動物數(shù)n和容許誤差δ的變化而變化，其中δ的影響略大，但當(dāng)δ達(dá)到一定值時，功效值便趨于穩(wěn)定。為了說明功效的增高是由于δ增大引起還是由于不同統(tǒng)計方法導(dǎo)致，我們選取了功效穩(wěn)定時的結(jié)果作為討論依據(jù)。在φ→0時，模擬數(shù)據(jù)服從Poisson分布，此時Poisson分布的精確概率法功效最高，但由于其Ⅰ類錯誤非常大，遠(yuǎn)遠(yuǎn)超出檢驗水準(zhǔn)故不推薦使用;負(fù)二項精確概率法與Fisher精確概率法的功效很相近，這是由于二者的算法類似，其值也高且Ⅰ類錯誤均小于0.05;t檢驗在樣本量n小于7時功效很小，在大于7時，功效尚能達(dá)到0.8左右，但其Ⅰ類錯誤也增大了，故不可接受，而log變換的t檢驗功效很低故不推薦使用。arcsin變換的t檢驗、平方根變換的t檢驗及秩檢驗在樣本量n為10時，功效較大，他們的Ⅰ類錯誤也小于檢驗水準(zhǔn)0.05;故在φ→0時，若動物數(shù)很少，建議使用負(fù)二項精確概率法和Fisher精確概率法，若動物數(shù)較大(大于10時)使用Arcsin變換的t檢驗、平方根變換的t檢驗及秩檢驗尚可。在φ→1時，Poisson分布的精確概率法與Fisher精確概率法的效能均很高，但他們的Ⅰ型錯誤很大，遠(yuǎn)大于檢驗水準(zhǔn)α=0.05，故不推薦使用;t檢驗、log變換的t檢驗功效較低，故不推薦使用;arcsin變換的t檢驗、平方根變換的t檢驗、秩檢驗在動物數(shù)大于9時功效尚可接受，但其Ⅰ型錯誤此時也大于檢驗水準(zhǔn)0.05，故不建議使用;負(fù)二項分布精確概率法的效能較高，僅次于Fisher精確概率法，其Ⅰ型錯誤也小于檢驗水準(zhǔn)0.05。故在φ→1時，建議使用負(fù)二項分布精確概率法。綜上表明，在φ→0時，動物數(shù)很少時，建議使用負(fù)二項精確概率法和Fisher精確概率法;若動物數(shù)較大(大于10時)使用arcsin變換的t檢驗、平方根變換的t檢驗及秩檢驗尚可。在φ→1時，只推薦使用負(fù)二項分布精確概率法。這里需要強(qiáng)調(diào)的是，很多時候由于動物間的個體變異存在，負(fù)二項分布假定比Poisson分布更合理，應(yīng)用者應(yīng)謹(jǐn)慎檢查數(shù)據(jù)分布假設(shè)再進(jìn)行恰當(dāng)?shù)慕y(tǒng)計分析〔3〕。

1．International Conference on Harmonisation Guidance on Genotoxicity Testing and Data Interprerpretation for Pharmaceuticals Intended for Human Use．Draft January 28th Version 5．3，2008．

2．Kirkland DJ．Statistical evaluation of mutagenicity test data:recommendations of the U．K．Environmental Mutagen Society．Environmental Health Perspectives Supplements vol．102 Suppl，1994，1:43-47．

3．Ludwig Hothorn．Biostatistical analysis of the micronucleus mutagenicity assay based on the assumption of a mixing distribution．Environmental Health Perspectives Supplements vol．102 Suppl，1994，1:121-125．

4．Hothorn LA，Gerhard D．Statistical evaluation of the in vivo micronucleus assay．Arch Toxicol，2009 Jun，83(6):625-634．

5．Robinson MD，McCarthy DJ，Smyth GK．edgeR:a Bioconductor package for differential expression analysis of digital gene expression data，Bioinformatics，2010 Jan 1，26(1):139-140．

6．Robinson MD，Smyth GK．Small Sample Estimation of Negative Binomial Dispersion，with applications to SAGE data．Biostatistics，2008，9(2):321-332．