SAGE的过度分散逻辑回归：多组和协变量建模

摘要

背景

基因表达序列分析（SAGE）数据中的两个主要可识别的变异来源是库内抽样变异性和组内库间异质性。大多数已发表的鉴别差异表达的方法只关注采样变异性。在最近的工作中，通过引入一个明确处理上述两种变异源的β-二项式层次模型，解决了评估两组SAGE库之间差异表达的问题。该模型产生了类似于加权双样本的测试统计t吨-测试。然而，当涉及的小组数量超过两个时，需要一种更通用的方法。

结果

我们描述了过度分散的逻辑回归如何提供这种泛化，并将其他协变量作为副产品合并到模型中。这种方法的优点是，逻辑回归例程在几个常见的统计包中可用。

结论

所述方法为分析SAGE数据提供了一种易于实现的工具，可以正确处理多种类型的变化，并允许更灵活的建模。

背景

SAGE的性质

Velculescu等人介绍的基因表达序列分析（SAGE）方法[1]是一种基于测序的基因表达测量方法。

简言之，将mRNA转录物转化为cDNA，然后进行处理以分离特定的子序列；从poly-A尾部开始，其子序列为10（正常SAGE）或14（长SAGE）bp，紧邻普通限制性内切酶第一次出现裂解位点之前。理想情况下，该子序列或“标签”具有足够的特异性，能够唯一识别其来源的mRNA。对标签进行采样、连接和排序，并组装一个由标签序列及其出现频率组成的表。从给定生物样本中导出的完整表格称为SAGE“库”。由于大多数标记在整个样本中都是稀疏的，因此大多数库都包含数以万计的标记数，以便估计表达式级别。然而，由于当前测序的成本，为给定实验组装的库的总数通常很小：通常为个位数，有时为十位数。

虽然SAGE实验中研究的信息类型，即基因表达，与cDNA或寡核苷酸微阵列实验中的信息类型相同，但在方法上存在一些定性差异。首先，SAGE使用测序而不是竞争杂交。第二，虽然阵列实验报告的表达值是对荧光的测量，并且是松散连续的，但SAGE以计数的形式提供基因表达数据，可能允许进行不同类型的“定量”比较。第三，SAGE是一种“开放”技术，它可以提供样本中所有基因的信息。相比之下，微阵列是“封闭的”，因为我们只能获得关于阵列上打印的基因的信息。

数学上，与样本中特定标签丰度相关的信息总结为两个数字：Y（Y）、库中该标记的计数，以及n个，库中的标记总数。在分析一系列库中的SAGE数据时，兴趣通常集中在评估当我们从一个库转移到另一个库时，基因表达的潜在真实水平是如何变化的。

微分表达式问题的数学公式

当跨一系列库进行调查时，包含单个标记的表达式更改的所有信息的足够统计信息是一组计数{Y（Y） _我 }和库大小集{n个 _我 }，其中下标我表示特定的库。除非另有规定，否则我们将把差分表达式的评估限制在单个标记的情况下。这种方法适用于下面描述的所有程序。在实际分析中，所选测试分别应用于所有标签，并报告显示差异表达的标签列表。不同的测试将提供对单个标签重要性的更改评估，因此提供的列表将取决于所采用的测试。

在大多数感兴趣的问题中，也存在协变量信息X（X） _我 描述库的属性我最常见的情况是比较两组库，如癌症库和对照库。在这种情况下，信息X（X） _我 只需定义哪个组库我属于。如果有两个以上的组，X（X） _我 可以有更多的级别，甚至可以是向量值，但和以前一样，兴趣集中在评估预期比例如何以及是否随X（X）.

关于比较两组之间的表达问题，已经做了很多工作。大多数方法[2——9]处理一个库与另一个库的比较。其中[2,6,7]将他们的考虑扩展到两组库的情况，将库集中在一组中，有效地将足够的统计数据减少到总数

这种方法虽然捕获了数据的计数性质，但由于忽略了组内比例的变化而丢失了信息。如Man等人所述[9]和Ruijter等人[10]在评估显著差异方面，上述大多数测试都给出了相同的结果。相比之下，这两个样本t吨-用于比较两组样本的测试[11]将足够的统计信息减少到一组比例{第页 _我 } = {Y（Y） _我 /n个 _我 }，捕获组成员之间的变化，但忽略了固有的计数采样性质和数据的可变性。两个样本t吨-测试结果可能与合并的测试结果有很大的不同，因为它们关注两种不同类型的变化。这两种方法对一组四个库的影响如表所示1。池将数据减少到右侧的合计计数，而关注比例将数据减少为底部的比例。在这两种情况下，这种减少都会导致信息丢失。当使用池时，我们无法判断其中一个组的比例是大的，另一个是小的，这表明不稳定。当使用比例时，我们无法判断一个库是否比另一个库小得多，因此该比例应该比其他库“不受信任”。

表1数据汇总方法。在四个库中测量的单个标记上的池和比例缩减的影响。池将数据减少为右侧的合计计数，而关注比例将数据减少到底部的比例。在这两种情况下，信息都会丢失。

Baggerly等人[12]提出了SAGE数据的贝塔-二项式层次模型，试图同时对这两种类型的变化进行建模。该模型导致了一个称为加权双样本的测试统计t吨-测试，t吨 _w个 计算此测试统计值需要表主体中的所有8个数字1; 并没有减少足够的统计数字。此测试统计显示了不同的行为，具体取决于给定标记的哪种类型的变化较大。当库内采样变化远大于库间采样变化时，t吨 _w个 给出的结果与池测试提供的结果接近，池测试侧重于库内变体。相反，当库间变量远大于库内变量时，t吨 _w个 给出的结果与双样本t检验的结果非常相似，该检验侧重于库之间的变异。这个t吨 _w个 该模型还允许评估两种变化的相对贡献。Baggerly等人[12]发现对于高计数标签，库间异质性是更大的变异来源，不允许异质性的池方法倾向于发现高计数标签存在显著差异。这可能会导致大量假阳性，这在绘制多个不同标签的结果时变得很明显。

扩展到多个组

文献中描述了两组以上的病例[2,13——15]，目前的分析方法是混合的。已经提出了明确解决多库问题的方法[16,17]，但目前最常见的方法[13,15]似乎涉及将数据的层次聚类与差异表达的成对测试相结合[2]在一组和另一组之间。这种混合方法可以间接捕获这两种类型的变异性，层次聚类侧重于组内比例之间的变异，成对检验侧重于抽样变异。除评估差异表达外，聚类还有其他好处，我们绝对建议使用它来探索数据的结构。然而，聚类往往不提供数值总结，因此将聚类结果与成对比较的结果结合起来可能是一门艺术。另一个缺点是，成对比较可能会由于只关注可用库的子集而错过有关可变性的有用信息。为了评估差异表达，我们认为可以进行更有效的测试。

我们的方法：过度分散的逻辑回归

我们寻求构建一种方法，该方法考虑到数据的计数性质，同时处理多个组，并允许由于单独采样而导致的比例变化。幸运的是，这不是第一次出现这样的问题。

评估多组差异表达的问题对应于方差分析（ANOVA）的经典统计问题。当感兴趣的值是连续的（例如，微阵列对数比）时，测试统计变成F检验，即两个样本的高维推广t吨-测试。当数据是计数（SAGE数据），并且需要处理采样变异性时，可以采用ANOVA检验来给出logistic或Poisson ANOVA。Stekel等人提出的差分表达式的多库测试[17]对应于泊松方差分析，但不考虑过度分散。方差分析处理从两个不同组扩展到更多不同组的问题，但这可以被视为协变量信息连续的特殊情况。对比例对协变量的依赖性建模的一种常见方法是通过逻辑回归或泊松回归，这两种回归都是广义线性模型的特例[18,19]. 这种模型直接包含了抽样可变性的形式。例如，比例的逻辑模型，

根据协变量（比例的logit）和每个测量值的精度定义要建模的数据的功能。该模型参数的最大似然估计可以通过迭代加权最小二乘法（IRLS）得到。

水平内的过度变化或过度分散可以通过多种方式引入逻辑回归框架。最常见和最广泛实施的方法是将上述最大化的二项式似然函数替换为“拟似然”函数，该函数与仅通过引入标度项的初始公式不同，

，到方差方程中，因此

五(Y（Y） _我 ) =n个 _我 第页 _我 （1）-第页 _我 ).

这种方法的优点是，它将每个观测值的方差放大了相同的量，因此

值将是相同的&只是相关的标准错误将被夸大。具有拟似然过分散的Logistic回归在各种统计软件包中实现，包括S-PLUS、R、GLIM和SAS。引入过度分散的另一种方法是假设一个层次模型，在该模型中，协变量在给定水平上的比例是从非退化分布中提取的，并且观察到的计数的分布是以提取的比例值为条件的二项式分布。当比例假设为β分布时，观测计数的最终无条件分布为β二项式。这是Baggerly等人提出的模型[12]用于模拟SAGE数据中的过度分散，也是Crowder使用的模型[20]泛化方差分析以处理过度分散的比例。它可以显示（例如，Collett[18]第201页），β二项式计数的方差的形式

五(Y（Y） _我 ) =n个 _我 第页 _我 (1 -第页 _我 )[1 + (n个 _我 - 1)φ],

这相当于当所有库大小n个 _我 都是一样的。虽然近似相等可能就足够了，但对于SAGE数据来说，即使是这种假设也可能有问题，特别是如果一些库是从不同时间进行的实验中提取的。威廉姆斯[21]显示了如何调整IRLS来处理这种过度分散，并注意到估计涉及φ只需要并且不需要假设beta分布的进一步结构，从而使该过程稍微更通用。这种形式的过度分散是作为dispmod包的一部分在R中实现的。

在逻辑回归框架中，评估差异表达简化为判断一组回归系数是否不同于零。这可能会导致与模型稍有不同的推论，例如t吨-应用于比例的统计数据，假设近似正态性适用于β值而不是比例本身。当我们为第页 _我 我们已经被引导选择了产生相当偏斜分布的参数，这表明logit尺度可能更合适。使用β这些值还有一个额外的优点，即置信区间可以根据折叠变化进行自然解释。

结果

比较两组

我们首先比较最初在Zhang等人[2]. 这8个库包括两个正常结肠（NC1和NC2）、两个原发肿瘤（TU98和TU102）和四个细胞系（CACO2、HCT116、RKO和SW837）。目前，我们重点比较正常结肠和所有肿瘤，原代或细胞系。表中给出了标签计数和相应的库大小2.Aχ²对两组中每一组的合并计数进行测试，得出的测试统计量为444.27；空值的95%截止值分布为3.84，高于此值的值被视为“显著”。两个样本t吨-对两组比例进行测试，得出1.60；空值的95%截止值t吨₆分布为±2.45，因此该测试表明差异不显著，表明专注于不同部分变异性的测试之间可能存在明显分歧。这个t吨 _w个 Baggerly等人提出的统计数据[12]合并了这两种类型的方差，得到了1.60的检验统计量。在这种情况下，此测试统计的空分布近似为t吨₆分布，定性结果与t吨-与合并测试相比，测试反映了患者异质性在驱动该标签总变异方面的相对优势。我们顺便注意到，这两种测试之间的分歧并不是孤立的事件。当我们调查这组库中的所有标签时，我们发现有10个标签|t吨|<2和χ²>200和48个标签|t吨|<2和χ²> 50. In Baggerly等人[12]研究发现，当使用集合测试时，大多数高计数标签出现显著差异，而当使用t吨-测试已经过尝试，在这种情况下t吨-测试更有可能是正确的。

表2 SAGE样本库中的标签计数。Zhang等人[2]的8个结肠文库中标签ATTTGAGAG、TGCTGCCTGT和GCGAAACCT的计数和比例；两个正常结肠（NC）、两个原发性肿瘤（TU）和四个细胞系。

三个logistic回归模型的结果与数据相符，如表所示三在第一个模型中，不允许存在过度分散，在第二个模型中使用了拟似然方法来处理过度分散，而在第三个模型中则使用了分层方法来处理过分散。这里，协变量的值X（X）分别为0或1，因为库位于第一组或第二组中。在型号1和2中，配合比例为=e（电子）^-4.66/（升+e（电子）^-4.66)=0.94%和=e（电子）^{-4.66 - 0.89}/（l）+e（电子）^{-4.66 - 0.89})第一组和第二组分别为0.39%，模型3中的比例仅略有改变。我们注意到，前两个模型的估计系数值完全相同，通常这两种方法也是如此。在不考虑过度分散的情况下对模型进行拟合，可以得到z（z）-的值β₁/s.e.公司。(β₁)=-20.42，这绝对重要。注意，该值的平方与χ得出的值的顺序相同²测试。该模型的皮尔逊残差，然而，这表明了一个问题。如果模型拟合良好，则应将其近似作为标准正态分布，并将一组8个观测值的极值在3或4量级左右。实际值-14.6和19.0太过极端。当模型适合于考虑过度分散的容差时，分散参数的点估计为= 187.57; 如果没有过度分散，这个值应该接近1。有了这个津贴t吨--1.49的值不再显著。这个t吨-值可以从第一个z（z）-值（-20.42）除以= 13.70. 同样，对残差进行缩放，得出的值与标准法线相称得多。我们注意到，由于所用模型的差异，测试统计数据的假定分布发生了变化。如果我们假设不存在过度分散的标准logistic模型成立，则检验统计量具有近似正态分布。这是因为总成功次数将二项分布推向近似正态分布。当我们转向假设存在过度分散的模型时，测试统计现在有一个t吨分配。这是因为我们对方差的估计现在强烈依赖于我们估计过度分散参数的精度，而这种精度取决于库的数量，而不是成功的数量。将该模型与过度分散的分层类型（表中的模型3）进行拟合三，得出的答案略有不同，但大小仍然不重要。不同之处在于之前发现的值是由于在这个模型中，每个比例的过度分散量随库大小略有变化，从而改变了回归模型中使用的权重。分层离散参数的点估计φ是= 3.399e（电子）-03，所以二项式方差的乘数是

表3两组的Logistic回归模型。Logistic回归拟合表2中标签ATTTGAGAG的正常结肠和癌症样本的对比。第一种模型不考虑过度分散，后两种模型以不同的方式引入。引入过分散很重要，因为它会显著影响结果，但选择过分散方法则不那么重要。

• 1 + ({n个 _我 } - 1)=

(169.62, 165.78, 141.62, 190.32,

207.26, 190.12, 175.35, 208.84).

平均这些值得出181.11，接近于准似然色散参数的值。我们注意到，模型2和模型3的系数值的差异在很大程度上是表面的，但模型1和其他模型之间的显著性差异并不是。选择解释过度分散比实现这一点所用的精确模型更重要。

我们注意到，过度分散的逻辑回归方法给出了t吨-值约为-1.49，而两个样本t吨-测试和修改后的版本t吨 _w个 Baggerly等人建议[12]两者都给予t吨-值约为-1.6（如前所述t吨和t吨 _w个 建议对于此标记，between-library变体比within-library变体大得多）。这种差异有两个原因。首先t吨统计适用于比例量表，logistic回归适用于βscale，大致是对数比例标度。其次t吨 _w个 此处使用的统计数据，

不假设所比较的两组中的过度分散系数相同；方差估计值未合并。后一种差异实际上对这种对比更为重要，特别是因为第一组2的方差估计值非常不稳定。这种影响并不总是微妙的；如果我们考虑标签GCGAAACCCCT，其计数在表中给出2，两个样本t吨测试和加权t吨 _w个 检验均为-1.57，logistic回归t吨值为-4.16。

在这两个答案中，我们倾向于选择逻辑回归拟合给出的答案，原因有二。首先，当我们直接将β分布的参数拟合为比例时，我们发现分布非常倾斜。因此，我们发现最好假设β系数标度。其次，当一个组中的库数量非常少（SAGE数据通常是这样的）时，我们更喜欢对方差进行合并估计。这种偏好在很大程度上是由于它通过使用更多自由度而具有更大的稳定性。在逻辑回归中，可以明确地纳入随协变量变化的过度分散水平，但我们在这里并没有追求这一点。

比较三个或更多组

如上所述，我们将结肠文库视为来自两组，但更自然的是将其视为来自三组：正常样本、原发肿瘤和细胞系。当我们有来自多个组的数据时，有两种不同的方式可以改变问题的性质。首先，如果我们只对比较两个组感兴趣，那么将其他组的数据合并到模型中通常是值得的。原因是，当过度分散导致方差时，我们的结果的重要性在很大程度上取决于估计过度分散参数的精度。未直接参与兴趣比较的组中的库仍然可以提供关于过度分散参数的信息，并增加相关t检验的自由度。其次，通过检查所有组的拟合比例，可以评估过渡的相对大小。

我们首先查看Zhang等人在论文中标记为感兴趣的单个标签的结果[2]即TGCTGCCTGT，我们假设最感兴趣的对比是正常结肠和原发肿瘤之间的对比。表中给出了此标签的计数和相应的库大小2.

我们首先尝试比较正常结肠和原发性肿瘤中的水平，而忽略细胞系（即仅使用四个库），然后使用包含所有三组的模型。使用具有分层过度分散的逻辑回归的结果如表所示4.

表4从两组扩大到三组的对比。Logistic回归模型检验正常结肠和原发肿瘤之间差异的意义(β₁)用于表2中的标签TGCTGCCTGT。在第一个模型中，只使用直接涉及的四个库中的数据。在第二个模型中，还包括来自四个细胞系库的数据，从而对过度分散参数提供了更稳定的估计φ.

在只有两组的模型中，我们有一个协变量向量x个₁=（0,0,1,1）表示该库属于两组中的哪一组。该模型产生的过度分散估计为= 8.938e（电子）-05，通货膨胀系数

• 1 + ({n个 _我 }-1）=（5.43、5.33、4.70、5.98）。

事实上，这些因素明显大于1，这表明组内异质性是模型未解释的方差的主要组成部分。在三组模型中，我们不能使用单个协变量向量x个₁，因为这不适合以无序的方式指示3个或更多组（分别对这三个组使用0、1和2将通过说原发性肿瘤是正常样本和细胞系之间的中间来强制排序）。一般来说，如果我们有k个组，我们需要使用k个-1个协变量向量。在这里，我们使用x个₁=（0,0,1,1,0,0,0）和x个₂= (0, 0, 0, 0, 1, 1, 1, 1). 所有0s的集合(x个₁＝0时，x个₂=0）对应于第一组，这里是正常冒号，其他组由其他协变量中的一个非零定义：组2（原色）(x个₁= 1,x个₂=0），组3（细胞系）(x个₁＝0时，x个₂=1）由于我们仍然关注正常结肠和原发肿瘤之间的差异，因此logit值为β₀和β₀+β₁主要关注的是β₁与零和细胞系组的预测logit显著不同，β₀+β₂，不直接输入问题。拟合此模型会产生以下过度分散估计= 1.160e（电子）-04，对于通货膨胀系数

• 1 + ({n个 _我 } - 1)=

(6.76, 6.62, 5.80, 7.46,

8.04, 7.45, 6.95, 8.09).

在这两种情况下（考虑两个组或三个组），正常结肠和原发肿瘤之间的对比，表示为

一旦考虑到过度分散，就显得很重要，但有一点值得注意。即使包括细胞系时过度分散的点估计值增加t吨-统计(/瑞典())与差异下降相关的p值表明显著性增加。在不使用细胞系的情况下，我们只有4个库，在估计各组的平均比例后，估计的自由度只有2度φ。当我们使用细胞系时，我们有8个库和5个自由度用于估计φ因此t吨-测试从2变为5。这个t吨₂分布有非常广泛的截止点，并且t吨₅更接近正常。一般来说，包含相关组可以通过提高我们对过度分散估计的精度来改进估计。

当然，在将模型与三组进行拟合时，我们也获得了查看其他对比的能力。例如，我们可以比较正常结肠和细胞系，其中logit是β₀和β₀+β₂分别通过检查同样，我们可以观察原发性肿瘤和细胞系之间的差异，logits是β₀+β₁和β₀+β₂通过测试差异的显著性虽然表中没有直接列出这种显著性，但我们可以计算出这种对比的标准误差，s.e.公司。 ，将估计除以其标准误差，得到具有所列自由度的t统计量（此处为5），并相应地计算p值。

还可以对各组之间是否存在任何显著差异进行综合测试，即比例的逻辑方差分析。常规ANOVA测试查看模型中感兴趣的术语解释的方差量，并将其与剩余方差量进行比较。根据每组的自由度进行调整，可以得到F类-测试。在处理广义线性模型时，数量-2*log（似然比），即偏差，在ANOVA中起着类似于方差的作用，因此我们可以说是偏差分析。由于模型中包含了过度分散，偏差分析变得复杂，需要一种多步骤的方法，其中几个不同的模型可以连续拟合。表中列出了这些型号5首先，使用所有可用的协变量对模型进行拟合，并估计过分散参数。这里，可用的协变量是x个₁和x个₂，并将完整模型与现有模型进行拟合，β₀+β₁+β₂，给出= 1.160e（电子）-04，如上所述。其次，子模型与输入的过分散值相匹配。在这种情况下，子模型为β₀+β₂，使用x个₂作为唯一的协变量，β₀+β₁，使用x个₁作为唯一的协变量，以及β₀，不使用协变量，只对所有数据进行单一比例拟合。结果如表所示5通过将偏差的比例减少与比例残差与适当的F类分配。例如，在这里测试整个模型是否包括β₁和β₂解释事情要比完全符合相同比例要好得多(β₀)减少到

表5偏差分析。各种子模型的偏差表与表2中给出的标签TGCTGCCTGT的数据相符。所有这些模型都使用最广泛模型的过度分散值，= 1.160e（电子）- 04.

表明各组之间的总体差异在5%水平上不显著。可以注意到，子模型包括β₁除了常数似乎解释得很少；这是由于我们选择参赛作品的方式X（X），从而包括β₁分离原发肿瘤组的影响，但排除β₂仍将正常结肠组与细胞系组进行对比。后一种分类模糊了正常结肠与原发肿瘤的区别，发现前者稍大一些。

合并其他协变量

可以使用逻辑回归方法在多个利益影响之间划分方差。例如，在上面的部分中，我们考虑了一个从原发肿瘤和细胞系中提取结肠文库的案例。这些数据也适用于其他器官，如胰腺。如果我们有兴趣确定原发性肿瘤和细胞系之间的一致性差异，那么使用这两种器官类型的库是很自然的。然而，如果将其作为两组进行比较，即原代细胞系与细胞系，则由于原代细胞组和细胞系组内的组织类型之间存在较大差异，因此很难分离出差异。解决方法是用两个协变量拟合一个模型x个₁样本分别为结肠或胰腺时为0或1，以及x个₂样本为0或1分别是原发肿瘤或细胞系。推理简化为测试β₂，只有在组织类型变化的影响之后才能评估自然变化的规模，β₁，已被排除在外。

在上面的例子中，我们考虑了另一种效果，即组织类型。原则上，可以通过包含其他协变量来考虑多种因素。同样，尽管上述示例中的两个协变量都是具有有限数量无序水平的“因子”，但也可以在建模过程中包括连续协变量。

为了说明这一点，我们使用表中GCGAAACCCT标记的计数给出了两个假设示例2在第一个例子中，我们假设我们正在尝试评估正常组织和原发性肿瘤之间的差异，如图所示，前4个库来自正常结肠和原发肿瘤，其余4个库不是来自细胞系，而是来自正常组织（库5和库6）和原发癌（库7和库8）来自其他器官。如上所述，这导致了一种情况，即我们希望用两个协变量拟合模型：x个₁=（0,0,1,1,0,0,1,1），表示库是正常的（0）还是原发肿瘤（1），以及x个₂=（0,0,00,1,1,1），表示库的来源器官。在第二个例子中，我们假设除了上述信息外，我们还可以获得潜在预测生存率的生物标记物水平。这些水平作为第三个协变量向量的值提供，x个_三= (0.89, 0.35, 0.66, 0.23, 0.30, 0.54, 0.90, 0.90). 的值x个_三从均匀分布中随机抽取。在拟合模型方面，力学与前面介绍的相似。模型拟合如表所示6.

表6将协变量纳入模型。处理表2中标签GCGAAACCCT计数拟合的模型，假设细胞系分配给正常组织B（库5和6）和癌组织B（库7和8）。本部分旨在说明两种不同的影响，正常与癌症以及组织A与组织B(β₁和β₂可根据其重要性进行划分。在模型2中，我们进一步引入了一个连续的协变量(β_三)与生物标志物的水平相对应，以显示如何将其计算在内。

当逻辑回归中断时

逻辑回归拟合过程可能会失败，或表现出缺乏收敛性。通常，这意味着其中一组中的所有比例为零或一；只有前者在SAGE数据的背景下才是现实的。这是很自然的，因为组比例的最大似然点估计为0β涉及到第二组中比例的折叠变化，导致被零除。当比例如此之小时，二项式变异性主导了异质性，并且这些值对于过度分散的估计是完全无信息的。

我们提出了一种迭代性质的修复方法，因为它需要运行三次逻辑拟合例程。为了说明这个过程，我们将使用表中标签ATTTGAGAG的数据2，前两个标签计数（来自组1）设置为零。

第一次运行拟合程序是为了估计过分散参数。这种拟合只使用计数不为零的组，忽略了有问题的组。在这里，这涉及到将一个比例与第2组中的六个库相匹配。拟合比例为0.40%，超差估计为

= 3.71e（电子）- 03.

第二次拟合过程采用给定的过分散参数，并在用相同的非零小比例替换组中的零比例后对数据进行拟合，为我们提供了一个保守的折迭变化估计值。这种类型的替换是常用的，并且最常见的是通过对比例的模糊先验分布的假设来证明，其中点估计被导出为后验平均值或模式。在处理比例时，先验的一个常见假设是均匀分布。成功0次后的后验平均值n个 _我 试验为1/(n个 _我 + 1); 如果进行多次试验，则为1/（∑n个 _我 ) + 1). 这就是我们使用的价值。由于两个原因，这个值实际上非常保守。首先，均匀分布使比例超过百分之几的可能性大大增加，这在SAGE数据中是无法观察到的。将分布限制在[0,0.02]范围内是均匀的，这应该是足够的。其次，过度分散的存在意味着汇集样本低估了观测比例的零方差所提供的小比例证据。虽然我们可以追求一个更优化的比例，但在这种情况下，我们选择简单地使用上面提到的简化界限。在这里，由于第一组中的库大小为49610和48479，比例为1/（49610+48479），伪造计数分别为0.506=49610/（49610/48479）和0.494。表中显示了此拟合的一些重新格式化结果7（模型1）。

表7拟合嵌套偏差模型。将嵌套模型拟合到数据中，以获得偏差分数。模型之间偏差的差异更好地表明了相关影响的重要性(β₁)当logistic回归拟合接近空间边界时，比例接近于零。

这种拟合的结果是荒谬的“微不足道”。问题在于使用t吨-值（Wald检验）依赖于最大值附近的似然函数的近似正态性，如果一组中的计数数量较少，则此形状假设将严重破坏。基于标度偏差变化的测试（与似然比测试相对应）更好。

第三次拟合过程适用于一个更简单的子模型，在这种情况下，所有八个库的比例都是单一的，使用相同的过度分散估计来测量偏差的变化。该拟合结果如表所示7（模型2）。显著性偏差检验分析

在这里，我们不能得出结论（考虑到过度分散的程度），这种差异是真实的。注意，分母中使用的自由度为5；这是因为只有6个库用于估计过分散参数，而估计比例需要6个自由度中的一个。

一般来说，当任何组的计数非常少时，检查偏差的变化是一个好主意。

讨论

过度分散的Logistic回归解决了SAGE数据的三个问题：同时建模多种类型的方差，同时处理多个组，并允许纳入协变量。此过程在可用的软件中得到了广泛的实现。此外，最重要的是，在逻辑回归设置中查看SAGE数据为思考描述此类数据的模型提供了框架。

处理多种类型的方差会产生显著性估计，我们认为其优于从集合计数或t吨-测试。回归设置还具有其他优点，例如关于模型检查、残差分析和异常值检测的完善工作。例如，可以评估任何给定库标记计数对整体分析的影响，并且可以通过限定这些函数使方法更加稳健，这样就不会有单个库驱动结果。

我们可以在一些领域发现困难并看到改进的余地。

首先，我们用于误差的模型可能会得到改进。对于SAGE数据，与特定标签关联的比例很少是百分数，因此

罗吉特(第页 _我 )≈对数(第页 _我 )

如果愿意的话，我们可以说使用日志而不是logit转换。假设方差在对数尺度上稳定，则会导致处理微阵列数据时通常假设的对数正态分布。假设对数正态分布等同于以另一种方式引入过度分散，即作为作用于β比例尺。这里是图书馆的真实比例我形式为

罗吉特(第页 _我 ) =β₀+β₁x个 _我 +ε _我 ,

哪里ε _我 是均值为0且方差为0的正态随机变量这里描述的模型是广义线性混合模型（GLMM）的一个特例，其中“混合”指的是我们既有固定的利益效应，也有协变量的变化，还有需要估计和考虑其方差的随机冲击。威廉姆斯[21]建议如何使用Taylor系列类型扩展来适应此模型，再次调用IRLS。然而，如Collett所述[18]第272页，“这种方法对于将此类模型拟合到二进制数据中并不完全令人满意，因为在某些情况下估计值可能会有偏差。此外，此类模型计算的偏差通常非常近似，不能推荐用于比较替代模型。”在SAS和S-PLUS中存在用于拟合GLMM的基于最大似然的方法，但在将混合效应模型拟合到具有少量聚类或库的二进制数据方面存在已知的问题。更准确地解决这个问题的一种方法是通过模拟（例如通过BUGS[22]). 我们现在正在研究这些不同的错误模型。

其次，上面开发的方法一次只能处理一个标记。这样做并没有充分利用SAGE数据的独特功能。此类开发的示例包括通过查看序列相似性用于定义邻域网络的邻域来纠正测序错误，以及通过使用诸如φ过度相似的群体。有关这些问题的工作正在进行中（例如，科林奇和费格[23]，N.Blades（2002），未发表的论文，Johns Hopkins），我们认为这些功能可以与此处介绍的方法有效结合。

方法

数据

这里使用的数据最初在Zhang等人中进行了描述[2]. 实际使用的数字库是从SAGE精灵Boon等人介绍的网络资源[24,25]. 这些库已删除链接器标记。

过度分散逻辑回归

这里只对该方法进行了粗略的描述；Collett中给出了更详细的处理方法[18]还有McCullagh和Nelder[19]等等。

我们希望符合观察到的比例，第页 _我 =Y（Y） _我 /n个 _我 ，作为协变量的函数X（X） _我 。这个过程的第一步是指定关系将采取何种形式。如果关系是线性的，那么第页 _我 =β₀+β₁X（X） _我 +ε，然后我们可以潜在地得到区间[0,1]以外的拟合比例，因此我们通常选择拟合转换版本的第页 _我 s在协变量中呈线性。比例的常见选择是逻辑变换，罗吉特(第页 _我 ) =日志(第页 _我 /1 -第页 _我 ) =β₀+β₁X（X） _我 +ε这种特殊的选择是由二项式数据的似然函数形式提出的（参见McCullagh和Nelder[19]，第28–32页），我们将按照这里的假设，但要注意的是，虽然logit可以覆盖所有实际值，但相应的比例都在0到1之间。此时，我们正在将一条直线拟合到数据的转换版本；这类似于标准线性回归，它通过最小化观测值与其拟合值之间的平方偏差之和来拟合：最小二乘法。现在，最小二乘法中的默认假设是，所有观测值都具有相同的精度，因此获得相同的权重。这里不是这样，因为比例的方差是五(第页 _我 ) =第页 _我 (1 -第页 _我 )/n个 _我 因此，已知观测值的精度取决于该观测值和总量的大小n个 _我 从中导出比例。如果观测值具有不同的精度，则标准平差是拟合加权最小二乘，将观测值与其拟合值之间的平方差的加权和最小化，其中权重与观测值的方差成反比。因此，在第一步中，我们使用加权最小二乘法拟合逻辑曲线，其中权重与我们对比例的初始估计相关的方差成反比(Y（Y） _我 + 0.5)/(n个 _我 + 1). 在第一次拟合之后，我们现在有了每个观测值的预测值，而这些预测值又为方差和权重提供了新的值。因此，第二步是使用新的权重重新调整数据。该过程将被迭代（迭代重加权最小二乘法，IRLS），直到预测值从一次拟合到下一次拟合的变化足够小，以至于该过程被称为收敛。

即使在该过程收敛后，如果方差的形式与上面给出的形式完全相同，则平方偏差的大小通常也会大大大于预期。在这种情况下，数据相对于假设模型表现出过度分散，我们试图估计过度分散的规模。我们在这里处理过分散的拟似然情况，其中方差实际上是形式上的五(第页 _我 ) =n个 _我 第页 _我 （l）-第页 _我 )，用于> 1. 计算层次结构形式的新增机制有些复杂，我们建议读者参考Williams[21]了解详细信息。使用过分散的拟似然模型，最佳拟合模型的实际参数不会改变，因为加权最小二乘例程中使用的权重都与方差的倒数成正比，并且用相同的因子缩放所有方差，使权重的相对大小保持不变。改变的是与这些参数相关的假定精度；参数的方差也将乘以，显著性检验需要相应调整。为了估计，我们回到上面提到的观测值和预测值之间的加权平方偏差。理想情况下，加权残差的平方和将具有具有k个-第页自由度，其中k个是库的总数第页是的数字β正在估计的条款。由于齐方分布的平均值等于其自由度，我们得到了以下的初始估计通过将加权残差的平方和除以假定的自由度：

鉴于

，测试统计数据按并重新计算了重要性。在下面的例子中，我们概述了过程，并将描述与免费软件包的脚本结合起来R（右）在每种情况下，该方法首先加载与标签计数对应的数据Y（Y） _我 和库大小n个 _我 ，用于提供观察到的比例。案例之间的主要区别在于协变量如何X（X）值已定义。所有模型都假设存在一个常数向量X（X）₀所有这些中的一个；这就产生了相应的估计β₀我们的讨论同样会将此协变量视为所有建模步骤中存在的协变量。

带注释的R代码

#本文中使用的模型的源代码

#“过度分散的Logistic回归

#SAGE：建模多个组和

#协变量”，Baggerly等人。

##########################################

#首先，我们处理两人的案件

#组，并介绍以下方法

#拟合logistic回归模型。

##########################################

如果（0）{

#加载ATTTGAGAG的标签计数（y）

#来自Zhang等人的8个图书馆。

# [2]，相关的库大小（n）

#协变量向量表示

#两组中哪一组是图书馆

#属于、正常或癌症（x）。

y<-c（320、600、312、549、，

246, 65, 41, 52);

n<-c（49610、48479、41371、55700、，

606825564125129461148）；

x<-c（0，0，1，1，1,1，1）；

#现在拟合一个标准的logistic回归

#根据数据建模，不留任何余量

#过度分散。完成了

#通过调用广义

#线性模型（glm）例程。帮助（glm）

#提供了有关

#这里论点的性质。

fit1<-glm（cbind（y，n-y）~x，

家族=二项式）；

#检查结果

总结（fit1）；

#接下来，我们在

#允许过度分散

#拟似然型；所有差异

#由于一个共同的因素而膨胀。这个

#调用与第一个调用的区别仅在于

#glm“家族”的定义

#使用。

fit2<-glm（cbind（y，n-y）~x，

族=拟二项式）；

#检查结果

总结（fit2）；

#理想情况下，平方和

#皮尔逊残差应该有一个chi-

#均数相等的平方分布

#自由度。分割

#自由度总和

#给出了我们对

#过分散参数。

varQL<-总和（残差（fit2，

“皮尔逊”）^2）/fit2$df.residual；

#最后，我们使用

#建议的过度分散方法

#威廉姆斯[21]，其中方差为

#因一些轻微因素而膨胀

#取决于基础

#库大小。这个例行程序是

#在R包“dispmod”中实现

#可在以下网址购买

#网址：http://cran.r-project.org

图书馆（“dispmod”）；

fit3<-glm.二项分布（fit1）；

#检查结果

总结（fit3）；

φ<-fit3$离散；

#注意，报告的p值来自

#这种配合不正确。这是由于

#假设test-stats

#具有正态分布，即使

#我们不得不估计

#过分散参数。当我们有

#要执行此估计

#正确的测试是t测试，带有

#自由度

#对应于

#图书馆数量减少估计

#参数。作为库的数量

#通常不是很大，这可以

#创造巨大的差异。

sumfit3<-摘要（fit3）；

t值<-摘要(

fit3）$系数[，“z值”]；

p.值<-2*pt（-abs（t.值），

fit3$df.剩余）；

}

##########################################

#接下来，我们处理三组

##########################################

如果（0）{

#我们首先关注收益

#当有多个组时可用

#即使其他组

#不是直接对比的一部分

#利息，由于额外

#添加的组可以

#提供关于

#过度分散。

#这里，我们使用标签中的数据

#TGCTGCCTGT，这次我们注意到

#有三组库：

#法线（库1–2），主

#肿瘤（库3-4）和细胞系

#（库5-8）。如果我们感兴趣

#法线和

#原发性肿瘤，我们可以用

#只有这两组的数据，

#或者使用这三个数据。

#首先，将模型拟合为

#只有两组人出席。

y<-c（0，1，1，15）；

n<-c（49610、48479、41371、55700）；

x<-c（0，0，1，1）；

fit1<-glm（cbind（y，n-y）~x，

家族=二项式）；

fit2<-glm.二项式.disp（fit1）；

#得到正确的p值

fit2.t.值<-摘要(

fit2）$系数[，“z值”]；

fit2.p.值<-2*pt（-abs(

fit2.t.values），fit2$df.residual）；

#接下来，假设

#有三组。在这种情况下，

#我们不能使用单一的协变

#向量x，因为这不适合

#表示一个

#无序时尚（使用0、1和2

#分别针对三组

#通过这样说来强迫命令

#原发性肿瘤是中间性的

#正常样本和细胞系之间）

#一般来说，如果我们有k个组，我们

#需要使用k-1协变量向量。

#在这里，我们使用

#x1<-c（0、0、1、1、0、0，0）；

#x2<-c（0、0、0，1、1、1，1）；

#所有0s的集合（x1=0，x2=0）

#对应于第一组，此处

#法线和其他组是

#由其中一个定义

#协变量非零：

#第2组（初级），（x1=1，x2=0），

#第3组（细胞系），（x1=0，x2=1）

y<-c（0、1、1、15、9、1、12、27）；

n<-c（49610、48479、41371、55700、，

60682, 55641, 51294, 61148);

x1<-c（0、0、1、1、0、0，0）；

x2<-c（0、0、0，1、1、1，1）；

fit3<-glm（cbind（y，n-y）~x1+x2，

家族=二项式）；

fit4<-glm.二项分布（fit3）；

#获得正确的p值

fit4.t.values<-摘要(

fit4）$系数[，“z值”]；

fit4.p.值<-2*pt（-abs(

fit4.t.values），fit4$df.residual）；

#上述方法符合模型

#所有的协变量都可用，

#但为了进行分析

#我们想要适应各种偏差

#使用相同估计值的子模型

#这里发现的过度分散。在这个

#在这种情况下，有3个子模型：

fit5<-glm（cbind（y，n-y）~x1，

族=二项式，

重量=fit4$disp.weights）；

fit6<-glm（cbind（y，n-y）~x2，

族=二项式，

重量=fit4$disp.weights）；

fit7<-glm（cbind（y，n-y）~1，

族=二项式，

重量=fit4$disp.weights）；

#或者，anova函数可以

#使用，但这仅考虑

#通过添加项获得的子模型

#顺序。因此，我们得到

#beta0（空模型）的偏差，

#beta0+beta1（添加x1

#仅协变量），和beta0+beta1+

#β2（将x2协变量添加到

#我们已经拥有的。

fit4.anodev<-anova（fit4）；

}

##########################################

#接下来，我们处理另一个

#协变量，可能是连续的。

##########################################

如果（0）{

#这里，我们使用的是来自

#GCGAAACCCT标签，但我们正在治疗

#来自组织的8个图书馆

#类型1（库1-4）和组织类型

#2个（文库5-8），具有正常组织

#两种类型（库1–2、5–6）和

#两种类型的原发肿瘤（库

# 3–4, 7–8). 在这个假设中

#例如，我们可以对

#更改为与相关的效果

#正常/主要差异（x1）或

#组织1/组织2差异（x2）。

y<-c（167、566、64、98、33、47、40、27）；

n<-c（49610、48479、41371、55700、，

60682, 55641, 51294, 61148);

x1<-c（0、0、1、1、0、0，1、1）；

x2<-c（0、0、0，1、1、1，1）；

fit1<-glm（cbind（y，n-y）~x1+x2，

家族=二项式）；

fit2<-glm.二项分布（fit1）；

#获得正确的p值

fit2.t.values<-摘要(

fit2）$系数[，“z值”]；

fit2.p.值<-2*pt（-abs(

fit2.t.values），fit2$df.residual）；

#接下来，再次使用上述标记，我们

#假设我们也可以访问

#潜在生物标志物水平

#生存预测，作为

#另一个协变量x3的水平。

#此处提供的值为

#作为随机抽取从

#均匀（0,1）分布

x3<-c（0.89、0.35、0.66、0.23、，

0.30, 0.54, 0.90, 0.90);

fit3<-glm（cbind（y，n-y）~x1+x2+x3，

家族=二项式）；

fit4<-glm.二项分布（fit3）；

#获得正确的p值

fit4.t.values<-摘要(

fit4）$系数[，“z值”]；

fit4.p.值<-2*pt（-abs(

fit4.t.values），fit2$df.residual）；

}