基于负二项分布的确切推断在微核试验数据分析中的模拟研究*

索瑞鑫 仇玉兰 王 彤△

微核试验是公认的检测染色体或有丝分裂器损伤的一种试验方法，其最大的优点是经济、简单、快速。该实验结果通常以啮齿类动物的骨髓嗜多染红细胞或者人外周淋巴细胞中的微核率来表示〔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．