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"))