负二项式 GLM 与计数数据的对数转换:增加 I 型错误率

机器算法验证 r 广义线性模型 模拟 负二项分布 类型 i 和 ii 错误
2022-01-30 14:40:42

你们中的一些人可能已经阅读过这篇不错的论文:

O'Hara RB, Kotze DJ (2010) 不要对计数数据进行对数转换。生态与进化方法一:118—122。点击

在我的研究领域(生态毒理学)中,我们正在处理重复性较差的实验,而 GLM 并未被广泛使用。所以我做了一个与 O'Hara & Kotze (2010) 类似的模拟,但模仿了生态毒理学数据。

电源模拟

我用一个对照组模拟了来自因子设计的数据(μC) 和 5 个治疗组 (μ1-5)。处理 1 中的丰度与对照相同(μ1=μC), 处理 2-5 的丰度是对照 (μ2-5=0.5μC)。对于模拟,我改变了样本量(3、6、9、12)和对照组中的丰度(2、4、8、...、1024)。丰度来自具有固定色散参数的负二项式分布(θ=3.91)。使用负二项式 GLM 和高斯 GLM + 对数转换数据生成和分析了 100 个数据集。

结果与预期一致:GLM 具有更大的功效,尤其是在采样动物不多的情况下。 在此处输入图像描述 代码在这里。

第一类错误

接下来我查看了第一类错误。如上所述进行了模拟,但是所有组的丰度都相同(μC=μ1-5)。

然而,结果并不如预期: 在此处输入图像描述 与 LM + 变换相比,负二项式 GLM 显示出更大的 I 型错误。正如预期的那样,随着样本量的增加,差异消失了。 代码在这里。

问题:

与 lm+transformation 相比,为什么 I 型错误会增加?

如果我们的数据很差(样本量小,丰度低(很多零)),我们应该使用 lm+transformation 吗?小样本量(每次处理 2-4 个)是此类实验的典型特征,不能轻易增加。

虽然,否定。斌。GLM 可以证明适用于该数据,lm + 转换可能会阻止我们避免类型 1 错误。

3个回答

这是一个非常有趣的问题。我查看了您的代码,没有立即发现明显的错字。

我希望看到你重做这个模拟,但使用最大似然检验来推断组之间的异质性。这将涉及重新拟合一个空模型,以便您可以获得θs 在组间比率同质性的零假设下。我认为这是必要的,因为负二项式模型不是线性模型(速率是线性参数化的,但是θs 不是)。因此,我不相信drop1论证提供了正确的推论。

大多数线性模型检验不需要您在原假设下重新计算模型。这是因为您可以仅在备择假设下使用参数估计和估计协方差来计算几何斜率(分数检验)和近似宽度(Wald 检验)。

由于负二项式不是线性的,我认为您需要拟合一个空模型。

编辑:

我编辑了代码并得到以下信息: 在此处输入图像描述

此处编辑代码: https ://github.com/aomidpanah/simulations/blob/master/negativeBinomialML.r

O'Hara 和 Kotze 的论文(生态学和进化方法 1:118-122)不是一个很好的讨论起点。我最担心的是摘要第 4 点中的主张:

我们发现转换效果很差,除了 . . .. 准泊松模型和负二项式模型... [显示] 偏差很小。

均值λ对于泊松分布或负二项分布是对于一个分布,对于θ<= 2 和平均值的范围λ被调查的,是高度正偏的。拟合正态分布的均值在 log(y+c) 的范围内(c 是偏移量),并估计 E(log(y+c)]。这种分布比 y 的分布更接近对称.

O'Hara 和 Kotze 的模拟比较了 E(log(y+c)],由 mean(log(y+c)) 估计,与 log(E[y+c])。它们可以是,并且在提到的情况下是,非常不同。他们的图表不会将负二项式与 log(y+c) 拟合进行比较,而是将均值 (log(y+c)] 与 log(E[y+c]) 进行比较。在 log(λ) 在他们的图表中显示的比例,实际上是负二项式拟合更有偏差!

以下 R 代码说明了这一点:

x <- rnbinom(10000, 0.5, mu=2)  
## NB: Above, this 'mu' was our lambda. Confusing, is'nt it?
log(mean(x+1))
[1] 1.09631
log(2+1)  ## Check that this is about right
[1] 1.098612

mean(log(x+1))
[1] 0.7317908

或者试试

log(mean(x+.5))
[1] 0.9135269
mean(log(x+.5))
[1] 0.3270837

估计参数的规模非常重要!

如果一个样本来自泊松,当然,如果根据用于拟合泊松的标准来判断,人们会期望泊松做得更好。负二项式同上。如果比较是公平的,差异可能不会那么大。真实数据(例如,可能在某些遗传背景下)有时可能非常接近泊松。当他们离开泊松时,负二项式可能会也可能不会很好地工作。同样,特别是如果λ对于使用标准正态理论对 log(y+1) 建模,大约为 10 或更多。

请注意,标准诊断在 log(x+c) 范围内效果更好。c 的选择可能无关紧要;通常 0.5 或 1.0 是有意义的。它也是研究 Box-Cox 变换或 Box-Cox 的 Yeo-Johnson 变体的更好起点。[Yeo, I. 和 Johnson, R. (2000)]。请参阅 R 的汽车包中 powerTransform() 的帮助页面。R 的 gamlss 包可以拟合负二项式 I(常见变体)或 II 型,或其他模拟分散和均值的分布,功率变换链接为 0(=log,即对数链接)或更多. 拟合可能并不总是收敛。

示例:死亡与基础损害 数据是针对到达美国大陆的命名大西洋飓风。数据可从最近发布的 R DAAG 包中获得(名称hurricNamed)。数据的帮助页面有详细信息。

稳健的对数线性与负二项式拟合

该图比较了使用稳健线性模型拟合获得的拟合线,以及通过将带有对数链接的负二项式拟合转换为用于图上 y 轴的 log(count+1) 刻度获得的曲线。(请注意,必须使用类似于 log(count+c) 刻度的东西,带有正 c,以在同一张图上显示来自负二项式拟合的点和拟合“线”。)注意大偏差是对数刻度上的负二项式拟合很明显。如果假设计数为负二项式分布,则稳健的线性模型拟合在此尺度上的偏差要小得多。在经典正态理论假设下,线性模型拟合将是无偏的。当我第一次创建基本上是上图的东西时,我发现这种偏见令人惊讶!曲线会更好地拟合数据,但差异在统计变异性的通常标准范围内。稳健的线性模型拟合对于规模低端的计数效果不佳。

注 --- 使用 RNA-Seq 数据进行的研究:比较两种类型的模型对于分析来自基因表达实验的计数数据很有意义。下面的论文比较了稳健线性模型的使用,使用 log(count+1),使用负二项式拟合(如在 Bioconductor 包edgeR 中)。在主要考虑的 RNA-Seq 应用程序中,大多数计数都足够大,以至于适当加权的对数线性模型拟合效果非常好。

Law, CW, Chen, Y, Shi, W, Smyth, GK (2014)。Voom:精确权重解锁了用于 RNA-seq 读取计数的线性模型分析工具。基因组生物学 15,R29。http://genomebiology.com/2014/15/2/R29

注意最近的论文:

Schurch NJ、Schofield P、Gierliński M、Cole C、Sherstnev A、Singh V、Wrobel N、Gharbi K、Simpson GG、Owen-Hughes T、Blaxter M、Barton GJ(2016 年)。RNA-seq 实验需要多少个生物学重复,应该使用哪种差异表达工具?RNA http://www.rnajournal.org/cgi/doi/10.1261/rna.053959.115

有趣的是,使用limma包(如edgeR,来自 WEHI 组)拟合的线性模型相对于具有许多重复的结果而言非常好(在几乎没有偏差的意义上),因为重复的数量是减少。

上图的 R 代码:

library(latticeExtra, quietly=TRUE)
hurricNamed <- DAAG::hurricNamed
ytxt <- c(0, 1, 3, 10, 30, 100, 300, 1000)
xtxt <- c(1,10, 100, 1000, 10000, 100000, 1000000 )
funy <- function(y)log(y+1)
gph <- xyplot(funy(deaths) ~ log(BaseDam2014), groups= mf, data=hurricNamed,
   scales=list(y=list(at=funy(ytxt), labels=paste(ytxt)),
           x=list(at=log(xtxt), labels=paste(xtxt))),
   xlab = "Base Damage (millions of 2014 US$); log transformed scale",
   ylab="Deaths; log transformed; offset=1",
   auto.key=list(columns=2),
   par.settings=simpleTheme(col=c("red","blue"), pch=16))
gph2 <- gph + layer(panel.text(x[c(13,84)], y[c(13,84)],
           labels=hurricNamed[c(13,84), "Name"], pos=3,
           col="gray30", cex=0.8),
        panel.text(x[c(13,84)], y[c(13,84)],
           labels=hurricNamed[c(13,84), "Year"], pos=1, 
           col="gray30", cex=0.8))
ab <- coef(MASS::rlm(funy(deaths) ~ log(BaseDam2014), data=hurricNamed))

gph3 <- gph2+layer(panel.abline(ab[1], b=ab[2], col="gray30", alpha=0.4))
## 100 points that are evenly spread on a log(BaseDam2014) scale
x <- with(hurricNamed, pretty(log(BaseDam2014),100))
df <- data.frame(BaseDam2014=exp(x[x>0])) 
hurr.nb <- MASS::glm.nb(deaths~log(BaseDam2014), data=hurricNamed[-c(13,84),])
df[,'hatnb'] <- funy(predict(hurr.nb, newdata=df, type='response'))
gph3 + latticeExtra::layer(data=df,
       panel.lines(log(BaseDam2014), hatnb, lwd=2, lty=2, 
           alpha=0.5, col="gray30"))    

原始帖子反映了 Tony Ives 的论文:Ives (2015)很明显,显着性检验对参数估计给出了不同的结果。

John Maindonald 解释了为什么估计存在偏差,但他对背景的无知令人讨厌——他批评我们表明我们都同意的方法是有缺陷的。许多生态学家盲目地做对数变换,我们试图指出这样做的问题。

这里有一个更细微的讨论:Warton (2016)

Ives, AR (2015),为了测试回归系数的显着性,继续对数转换计数数据。方法 Ecol Evol,6:828-835。doi:10.1111/2041-210X.12386

Warton, DI, Lyons, M., Stoklosa, J. 和 Ives, AR (2016),为计数数据选择 LM 或 GLM 测试时要考虑的三点。方法 Ecol Evol。doi:10.1111/2041-210X.12552