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

索瑞鑫 仇玉蘭 王 彤△

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

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

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

基于負二項分布的確切推斷

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

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

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

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

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

模擬研究

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

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

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

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

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

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

當離散參數φ→0時，檢驗功效結果見表3。

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

當離散參數φ→1時，檢驗功效結果見表4。

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

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

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

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

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

討 論

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