在 R 中使用 geom_density_2d() 时出错：`stat_density2d()` 中的计算失败：带宽必须严格为正

Question

在尝试使用 ggplot2 制作测试 2d 密度图时，我使用了代码片段：

ggplot(df, aes(x = S1.x, y = S1.y)) + geom_point() + geom_density_2d()

我收到错误消息：“计算失败stat_density2d()：带宽必须严格为正”

我的数据框如下所示：

> df

transcriptID S1.x      S1.y      S2.x       S2.y    
DQ459412     0.000000  0.000000  0.000000   0.000000
DQ459413     1.584963  2.358379  4.392317   3.085722    
DQ459415     0.000000  0.000000  0.000000   0.000000    
DQ459418     0.000000  0.000000  0.000000   0.000000    
DQ459419     0.000000  0.000000  4.000000   2.891544    
DQ459420     0.000000  0.000000  0.000000   0.000000      

还有，var(df[,"S1.x"]) > 0和var(df[,"S1.y"]) > 0。

图 1 - 有误差的二维密度图

但是，通过运行，我得到了一个没有错误的密度图：

ggplot(df, aes(x = S2.x, y = S2.y)) + geom_point() + geom_density_2d()

图 2 - 没有错误的密度图

如何解决图 1 中的错误？

Answer 1

所以真正的问题是S1.x和S1.y值在它们的列中只有一个非零值。事实证明，geom_density_2d仅凭一两个值并不能真正估计密度。但是请继续阅读...

更新：

这个问题以前有人问过，答案通常是你的数据列中需要有非零方差。但是您确实有非零方差，那么为什么它不起作用呢？

• 查看内部geom_density_2d我们看到它使用MASS::kde2d包函数来计算分布。
• 查看kde2d我们看到它用于MASS::bandwidth.nrd(df$x)估计带宽。
• 查看帮助（其中包含代码），bandwidth.nrd我们看到它使用经验法则获取quantile分布的 ，并从第 1 个分位数中减去第 2 个分位数以获得带宽估计值。
• 对原始数据进行分位数分析，我们看到数据的分位数为零。
• MASS::kde2d使用bandwidth.nrd该带宽估计值在原始数据上运行会给您带来相同的错误：

library(MASS)
nn <- c("DQ459412","DQ459413","DQ459415","DQ459418","DQ459419","DQ459420")
s1x <- c(0,1.584963,0,0,0,0)
s1y <- c(0,2.358379,0,0,0,0) 
s2x <- c(0,4.392317,0,0,4,0)
s2y <- c(0,3.085722,0,0,2.891544,0) 
df <- data.frame(transcriptID=nn,S1.x=s1x,S1.y=s1y,S2.x=s2x,S2.y=s2y)

> quantile(df$s1x)
      0%      25%      50%      75%     100% 
0.000000 0.000000 0.000000 0.000000 1.584963 
> quantile(df$s1y)
      0%      25%      50%      75%     100% 
0.000000 0.000000 0.000000 0.000000 2.358379 

h <- c(MASS::bandwidth.nrd(df$x), MASS::bandwidth.nrd(df$y))
dens <- MASS::kde2d(df$s1x, df$s1y, h = h, n = n,  lims = c(0,1,0,1))

MASS::kde2d(df$s1x, df$s1y, h = h, n = n, lims = c(0, 1, 0, 1)) 中的错误：带宽必须严格为正

所以使用的真正标准geom_density_2D是 x 和 y 数据都需要在它们的第一个和第二个分位数之间有一个非零的差距。

现在修复它，如果我做一个小的修改 - 用 0.1 替换一个零，如下所示：

nn <- c("DQ459412","DQ459413","DQ459415","DQ459418","DQ459419","DQ459420")
s1x <- c(0,1.584963,0,0,0.1,0)
s1y <- c(0,2.358379,0,0,0.1,0) 
s2x <- c(0,4.392317,0,0,4,0)
s2y <- c(0,3.085722,0,0,2.891544,0) 
df <- data.frame(transcriptID=nn,S1.x=s1x,S1.y=s1y,S2.x=s2x,S2.y=s2y)
print(df)

产生：

  transcriptID     S1.x     S1.y     S2.x     S2.y
1     DQ459412 0.000000 0.000000 0.000000 0.000000
2     DQ459413 1.584963 2.358379 4.392317 3.085722
3     DQ459415 0.000000 0.000000 0.000000 0.000000
4     DQ459418 0.000000 0.000000 0.000000 0.000000
5     DQ459419 0.100000 0.100000 4.000000 2.891544
6     DQ459420 0.000000 0.000000 0.000000 0.000000

然后我得到这个情节而不是你的错误。

您可以让该0.1值接近零，最终它将无法再计算分布，您将再次出现错误。

处理这种情况的一种一般方法是向您的数据添加非常少量的噪声，模拟这样一个事实，即基于连续分布的真实测量的任何有意义的计算都应该不受该少量噪声的影响。

希望有帮助。

Answer 2

@Mike Wise 的答案确实非常可靠，我的答案在某种程度上是对它的补充。实际上，该bandwidth.nrd函数计算第 3 个和第1 个分位数之间的差异，而不是第 2 个和第 1 个（来自函数的代码）：

r <- quantile(distances, c(0.25, 0.75))

我建议不要向数据添加随机噪声，而是建议自己预先计算带宽并将它们传递给函数，测试非零值，如下所示：

kde2d(df$s1x, df$s1y, 
      h = c(ifelse(bandwidth.nrd(df$s1x) == 0, 0.1, bandwidth.nrd(df$s1x)),
            ifelse(bandwidth.nrd(df$s1y) == 0, 0.1, bandwidth.nrd(df$s1y))))

希望这可以帮助。