微阵列研究中的功率和样本量估计

摘要

背景

DNA微阵列技术为同时研究数百或数千个不同基因的表达谱提供了工具。微阵列研究的一个基本目标是确定在感兴趣的实验条件下差异表达的基因子集。在进行微阵列实验之前，需要确定的一个重要问题是为了有足够的能力识别差异表达基因所需的阵列（重复）数量。

针对不同的I类错误规范，如家庭错误率（FWE），已经开发了许多样本量估计方法[1-三]，错误发现率（FDR）[4-8]，以及假阳性的数量[7，9]。微阵列研究的样本量通常计算为达到指定功率所需的阵列数量平均值（例如[三-6，9，10]). 这种能力，即预期检测到的真正差异表达基因的比例，被称为敏感性。通过计算样本量估计值以达到规定的灵敏度平均值，检测到的真正差异表达基因的比例往往低于平均值。因此，计算的样本量往往会给出过于乐观的结果。或者，Wang和Chen[2]，Tsai等人[7]邵和曾[8]提出了一种替代公式：计算样本量以确保以指定的概率检测到至少指定的灵敏度水平。这将被称为（置信）概率公式。

当制定样本量问题以达到规定的灵敏度时平均值，我们将表明，可以使用单变量样本量公式简单地计算所需的样本量，而不考虑基因之间的依赖性。另一方面，如果问题的形式是以特定的概率达到特定的灵敏度，则需要估计灵敏度分布的百分位。在这种情况下，需要考虑基因之间的依赖性。Tsai等人[7]提出了一种在所有基因具有恒定功率的独立或等相关正态分布模型下控制比较误差率（CWER）的方法。邵和曾[8]提出了一种在正态分布下估计一般相关矩阵的无模型方法。他们使用了72个样本的数据集来说明相关性矩阵的估计。然而，先导数据的大小通常很小，每组10个或更少，并且真阳性的估计方差通常为负（零），导致我们的模拟研究中对样本大小的估计较差。提比什拉尼[10]提出了一种用于估计特定样本量的FDR和平均灵敏度的置换方法。Tibshirani的方法只需要少量的试验数据集，并且完全无模型，即不需要对测试统计的分布、效应大小和相关性进行假设。然而，检验统计量的标准差估计值（标准误差）取决于样本量。小样本量的测试统计值的变化比大样本量的变化大。由于试验数据集的样本量通常很小，因此基于小试验数据集得出的截止水平通常会超过所需样本的实际截止水平，从而导致对所需样本量的过高估计。

本文概述了幂规范和参数规范，并提出了概率公式下样本容量确定的置换方法([2，7，8]). Tibshirani的方法[10]通过加入一个调整因子来改进以获得更正确的排列分布。该方法使用每组4到6个样本的小型试验数据集；该方法需要的样本比Tibshirani少[10]方法。当试验数据集的样本量较大时，所提出的方法和Tibshirani[10]方法是等效的。

方法

让米表示在数组中研究的基因数量，其中米₀和米₁分别是非差异表达基因和差异表达基因的数量。给定显著性水平α（按比较误差率）米测试可以总结为一个2×2的表（表​（表11).

V（V）/米₀是宣布显著的未差异表达基因的比例，其期望值是每次比较的错误率E(V（V）)/米₀=α.V（V）/R（右）是已声明的重要基因在已声明的、实际上没有差异表达的重要基因总数中的比例。它的期望值是错误发现率E(V（V）/R（右）) =q个，给定R（右）> 0.U型/米₁是正确声明的真正差异表达基因的比例。在诊断问题中，这个比例通常被称为真阳性率或灵敏度。通过取期望值，我们得到了“平均灵敏度”E(U型)/米₁，表示为λ.

样本量估计

在样本量估计中，米，米₁和（标准化）效果大小δ= (δ₁,...,δ_米1)因为差异表达基因是由研究者预先指定的。估计达到规定灵敏度所需的样本量λ₀一般来说，这很简单。自米₁和λ₀是预先指定的，给定FDR级别q个*每次比较误差率的相应显著性水平α很容易计算。设置α= [米₁λ₀q个*]/[米₀(1-q个*)]，FDR将控制在q个*足够大的米₁和米₀.

如果δ_我=δ₀对所有人来说都是不变的我，然后是比较功率（1-β)单变量检验的结果与λ₀.给定α，δ₀、和（1-β) =λ₀，样本量可以基于单变量样本量计算，如下所示

(1)

哪里t吨_α和t吨_β是a的百分位数t吨-分配。如果δ_我那么，是不同的β_我=t吨^-1()根据方程式（1）。样本量n个*可以通过以下公式计算

（2）

所需的样本量为n个= ⌈n个*É，其中n个*⌉是大于或等于的最小整数n个*. 给定样本大小n个经计算，一个真正差异表达基因的单变量测试结果可以用贝努利随机变量建模，成功概率至少为（1-β_我)自n个≥n个*. 实际检测的预期数量至少为米₁λ₀无论基因之间的相关结构如何，因此平均可以达到期望的灵敏度。大多数样本量估计方法要么基于此方法，要么基于扩展[三-6，9-11]。然而，根据该公式计算的样本量不足；下面显示了一个独立模型下的简单演示。

鉴于米，π₁(=米₁/米)，效果大小恒定δ_我=δ₀，q个*,λ₀，以及计算的样本量n个（基于方程1），在独立模型下，检测到的真正差异表达基因的总数U型是具有成功概率（1）的二项式随机变量-β) (≥λ₀自从n个≥n个*). 概率ϕ_λ0至少确定λ₀的分数米₁差异表达基因可以计算为二项式概率之和[2，7]:

(3)

使用方程（1）估计样本量的方法称为单变量方法。表第3-5列​表22显示估计的样本大小n个，平均灵敏度λ和概率ϕ_λ0在λ₀= 0.6, 0.7, 0.8, 0.9. 计算中使用的参数为：米= 2,000,π₁= 5%, 10%, 20%,δ₀=2和q个* = 0.05. 可以看出ϕ_λ0可以小于60%。也就是说，使用此公式计算所需阵列可能会导致实验的灵敏度低于指定的λ₀概率超过40%。

表2

独立模型下的平均公式与95%概率公式。^一

		平均配方： 单变量法			95%概率公式： 二项式方法
π1	λ0	n个b条	λ	ϕλ0	n个c（c）	λ	ϕλ0
5%	60%	9	0.70	0.985	9	0.70	0.985
	70%	9	0.70	0.576	10	0.81	0.997
	80%	10	0.81	0.681	11	0.88	0.993
	90%	12	0.92	0.866	13	0.95	0.992
10%	60%	8	0.70	0.999	8	0.70	0.999
	70%	8	0.71	0.687	9	0.82	1
	80%	9	0.82	0.841	10	0.89	1
	90%	11	0.93	0.977	11	0.93	0.977
20%	60%	7	0.72	1	7	0.72	1
	70%	7	0.74	0.975	7	0.74	0.975
	80%	8	0.85	0.996	8	0.85	0.996
	90%	9	0.91	0.792	10	0.95	1

在单独的窗口中打开

a.估计样本量n个，平均灵敏度λ和概率ϕ_λ0对于规定的灵敏度λ₀=独立模型下的60%、70%、80%、90%。计算中使用的参数为：米= 2,000,π₁= 5%, 10%, 20%,δ₀=2和q个* = 0.05.

b.样本量n个由方程式（1）中的单变量方法计算得出，以获得灵敏度λ₀平均而言。

c.样本量n个使用Tsai等人[7]确保概率的方法ϕ_λ0至少检测到λ₀差异表达基因的比例至少为95%。

或者，王和陈[2]将问题公式化为：达到指定灵敏度所需的阵列数量λ₀有可能ϕ_λ0.在此配方中λ₀和ϕ_λ0需要指定，但不一定相等。这个ϕ_λ0设置为95%，因为它符合使用95%置信概率的常见统计实践。在此配方下，对于指定的λ₀计算所需的阵列数量，使平均灵敏度大于λ₀和5^第个百分位数，λ₅灵敏度分布U型/米₁大于λ₀:

在独立和恒定效应规模模型中，Tsai等人[7]使用方程（1）和（3）估计所需的样本量，即二项式方法。表第6-8列​表22显示估计的样本大小n个，平均灵敏度λ，以及概率ϕ_λ0对于λ₀= 0.6, 0.7, 0.8, 0.9. 第8列中的概率均高于95%，原因是n个≥n个*. 该程序将确保以至少95%的概率检测到特定比例的差异表达基因。

在表中​表2，2，理论结果表明，这两种方法给出了相当接近的样本量估计。估计值的差异反映了两种公式的差异；什么时候δ₀=2，差值最大为1。对于给定的灵敏度，所需的样本大小随着效应大小的增加而增加δ₀减少，估计值中两个公式的差异较大。我们使用与表中相同的参数计算样本大小​表22对于δ₀= 1. 样本量差异大约是表的四倍​表22（未显示数据）。

结果

进行了两次模拟分析，以评估上述两种样本量估计公式。第一次分析在独立和恒定效应大小模型下评估了这两种配方。两种配方的理论结果如表所示​表2。2模拟分析提供了经验验证。第二个分析在相关模型下评估了四种方法：1）单变量方法（例如Jung[4]); 2） 邵和曾[8]无模型方法，3）Tibshirani[10]置换法；以及4）所提出的置换方法。单变量方法是为平均公式设计的，而其他三种方法是使用先导数据集以95%的概率考虑的。表中的相同模型参数​表22用于评估。I类错误率基于将FDR设置为q个* = 0.05. 请注意，有许多具有不同策略的多重测试FDR程序。例如，Storey的FDR程序[13]包括对非差异表达基因数量的估计米₀然而，为了最大限度地减少估计中的变化带来的混淆影响米₀，我们只是使用了米₀在我们的仿真分析中。计算给定参数值的样本大小。FDR的经验估计，平均灵敏度λ和概率ϕ_λ0然后进行计算和评估。使用true米₀在FDR的控制下，对提议的程序进行直接验证。

第一次模拟研究的目的是验证表中所示两种方法的样本量、灵敏度和95%概率的理论结果​表22在独立模型下。我们从表的第3列或第6列中生成了1000个模拟样本，每组样本大小​表2。2。对于空模型，米₀=米× (1 -π₁)基因来自独立标准正常N（0,1）；对于替代模型，米₁=米×π₁基因是基于独立的正常N生成的(δ₀, 1). 对于每个模拟样本集t吨-计算统计学和相应的p值，以及FDR水平的假阳性和真阳性的数量q个*=0.05。FDR的经验估计，平均灵敏度λ和概率ϕ_λ0然后进行计算。估计ϕ_λ0是1000个模拟中真阳性数不小于的时间比例米₁×λ₀.

表​表3三显示了两种方法的实证结果。在这两种方法中，经验FDR似乎接近名义水平。对于单变量方法，经验平均灵敏度λ除了π₁=0.05和λ₀= 70%. 概率ϕ_λ0小于50%，对于π₁=0.05和λ₀= 70%. 对于二项式方法，经验平均灵敏度λ的都大于指定的级别。大多数可能性ϕ_λ0的值超过95%，除了π₁= 0.05,λ₀= 60%,π₁= 0.10,λ₀=90%和π₁= 0.20,λ₀= 70%. 表的实证结果​表3三通常与表中所示的理论值一致​表2。2也就是说，使用单变量方法计算的样本量通常会达到规定的平均灵敏度；然而，达到规定灵敏度的概率可以低于50%。

为了进行比较，表的最后一列还提供了使用组大小为4的试验数据集的拟议排列方法的样本大小估计值的平均值和标准偏差​表3。三。引导数据是从每次模拟的正态分布中随机生成的。提出的方法往往会高估多达五个阵列所需的样本量。

第二个分析是评估四种方法，即单变量方法（Jung[4])、Shao和Tseng[8]，提比什拉尼[10]，并在使用已知结肠癌数据集的相关模型下提出置换方法[14]。结肠癌数据集[14]由22个正常结肠肿瘤组织和40个结肠肿瘤组织样本组成，共有2000个基因。分析包括两个步骤。第一步根据每组样本量为4和6的试验数据集，评估了三种95%概率公式方法获得的样本量估计值。第二步将第一步提出的方法估计的样本量与单变量方法估计的样本量进行比较。

在第一步中，从结肠数据集中随机选择4个样本，每组不进行替换，以形成试点数据集。上述算法用于估计所提方法和Tibshirani的样本大小[10]方法。例如，对于π₁= 5%,q个*=0.05和λ₀=90%，初始样本量为n个=13（表第6列​表2）2)和α= 0.00249. 恒定的效果大小δ_我=δ₀已考虑=2。对于建议的排列方法（f）是（f）₁=0.6777和（f）₂==1.155，而Tibshirani未进行调整[10]方法。对于邵和曾[8]无模型方法，相关矩阵t吨-统计数据是通过使用来自先导数据集的所有可能的置换数据集进行估计的。然而，邵氏和曾氏[8]大多数情况下，无模型方法存在计算困难。详情稍后给出。

重复该过程1000次，从每组中选择大小为4的不同先导数据集，以说明先导数据集的变化。Tibshirani样本量估计的平均值和标准偏差[10]和建议的方法进行了计算，如表的第4列和第5列所示​表4。4一元方法被视为标准方法，估计值列于第3列。Tibshirani估计的所需样本量[10]或者在每种情况下，所提出的方法都大于单变量方法。在每种情况下，单变量方法和提出的方法之间的差异小于每组5个数组。来自Tibshirani的平均值和标准差估计[10]该方法比所提出方法的估计值大得多。差异随着λ₀增加或π₁减少。注意，在独立模型下，拟议方法的样本量和标准偏差估计值较小（表​（表3三).

对初始试验数据集的6个样本重复该过程。估计值如第6列和第7列所示。建议的程序从两个试点样本大小得出了一致的结果；然而，Tibshirani的结果[10]方法有很大不同。Tibshirani方法没有充分考虑试点样本大小。当试点样本量远小于所需样本量时，Tibshirani对样本量的高估[10]方法变得严重。随着试点研究规模越来越接近所需的样本规模，Tibshirani[10]所提出的方法将给出类似的结果。

在我们的仿真中，Shao和Tseng中的算法B[8]无法在所有1000次复制中成功为组大小为4的试点数据生成解决方案。当组大小增加到6时，该算法仅在以下情况下有效π₁= 20%,λ₀=60%和70%；样本量估计值的平均值（标准偏差）分别为6.4（0.012）和6.8（0.012”）。估计值似乎太小而不正确。该方法似乎不适用于小样本。使用整个结肠癌数据集[14]在62个样本中，样本量估计值如第8列所示。与单变量方法相比，估计通常需要多出一个或两个数组，但比提出的方法少。自从Tibshirani[10]方法给出了更大的估计值，Shao和Tseng[8]给出的估计值小于建议的方法。在第二步分析中，对单变量方法和建议的方法进行了评估。

两种方法的性能比较与表中所示类似​表3。三这些数据是从结肠癌数据集而不是从独立模型下的正常随机变量中取样的，没有进行替换。样本大小基于表的第3列或第4列​表4。4。然后对数据进行随机排列，以消除两组之间的差异，以及一个共同的效应大小δ₀=2被添加到一组随机选择的米₁肿瘤组中的基因。对于每个重新采样的数据集，使用置换测试生成p值，并使用以下公式计算假阳性和真阳性的数量q个* = 0.05. 计算置换测试的重复次数为10000次。FDR的经验估计，λ和ϕ_λ0已计算。整个过程重复1000次。

表​表55显示了对q个*,λ、和ϕ_λ0对于这两种方法。这两种方法都能很好地控制FDR并达到所需的灵敏度。因此，可以预期这两种方法在实践中具有令人满意的性能。然而，对于单变量方法，经验ϕ_λ0估计值在55%到75%之间，只有一个是80%。人们必须承担敏感度低于规定水平的风险。

的效果大小δ₀=2（表​（表4）4)使用结肠癌数据集在相关模型下验证所提出的置换方法[14]。实际上，效果大小可以小得多。我们使用与表中相同的参数计算样本大小​表44具有效果大小δ₀=1，对于两个试验样品尺寸4和6。样本量估计如表所示​表6。6.建议的程序对两个试点样本大小给出了类似的结果，这与以下结果一致δ₀=表中的2​表4。4。单变量方法和建议方法之间的差异约为每组15个数组。Tibshirani家族[10]对于4个和6个导频样本，每个组需要多达67个和35个额外阵列。邵和曾的估算[8]只有当试点研究规模大于或等于所需样本规模时，才能对方法进行评估。

讨论和结论

在进行微阵列实验之前，确定所需的样本量是一个重要问题。样本量问题通常表示为达到规定灵敏度所需的阵列数量λ平均而言。本文证明了在该公式下计算的样本量可能具有灵敏度λ平均在规定的水平，但概率ϕ_λ由于灵敏度分布的差异，达到规定的灵敏度可能很低（小于50%）。此外，在这个公式下，本文表明，无论基因表达水平之间的相关结构如何，样本量都可以用单变量方法计算；解释相关性的程序，如Li等人[6]不需要（表​（表5）。5). 这些发现与Jung报告的结果一致[4]多宾和西蒙[11]。然而，本文为这种方法提供了理论解释。

在置信概率公式下，估计样本大小时需要考虑基因表达之间的依赖性，因为敏感性分布的百分位数不仅取决于单个基因的影响大小，还取决于它们的相关性。我们基于Tibshirani提出的方法提出了一种置换方法[10]，但包括一个调整系数和一个要求，以95%的概率达到特定的灵敏度。调整系数可以更准确地估计功率和样本大小。邵和曾[8]还根据置信概率制定了所需的样本量。在正态假设下，邵和曾[8]使用初步数据集提出了基因间轻度相关性的算法。他们表明，他们的方法对72个样本的示例数据集很有效。然而，在我们的结肠数据集模拟中使用他们的算法B（结肠数据集的平均相关性约为0.4），当初始样本量为4或6时，真阳性的估计方差可能为负。他们的程序对于样本量较小的小样本试验数据集没有很好的表现。实际上，试点数据的样本量通常很小。我们的模拟研究表明，我们的程序可以很好地处理每组4到6个样本。然而，当相关性很小时，我们的程序似乎过高估计了所需的样本量，尤其是在小效应量的情况下。在这种情况下，我们的模拟结果表明（f）₂可能不需要（未显示数据）。

用于数据分析的特定多重测试程序的选择会影响样本量估计中的错误率和功率。在数据分析中使用保守的程序可能会降低研究的“力量”；有时，计算出的样本量可能具有低于指定水平的灵敏度。例如，在本文中，计算基于非差异基因的真实数量米₀然而，如果数据分析使用高估值米₀比如Benjamini和Hochberg程序[15]，则功率可能低于所需水平。另一种方法是使用基因总数米而不是非差异基因的数量米₀估计样本量。无论使用哪种多重测试程序进行数据分析，该程序均应产生适当的样本量，以达到规定概率的期望灵敏度。

作者的贡献

JJC构思了这项研究并撰写了手稿。JJC和WJL开发了该方法并证明了理论结果。WJL实现了这些算法。HMH改进了平均和95%置信概率公式的概念。JJC、HMH和WJL进行了分析。所有作者阅读并批准了最终手稿。

补充材料

工具书类