在 R 中使用 geom_密度_2d() 时出错：“stat_密度2d()”中计算失败：带宽必须严格为正数

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

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      

Also, var(df[,"S1.x"]) > 0 and var(df[,"S1.y"]) > 0.

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

但是，我通过运行以下命令得到了没有错误的密度图：

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

图 2 - 无误差的密度图

如何解决图 1 中的错误？

所以真正的问题是S1.x and S1.y值的列中只有一个非零值。事实证明geom_density_2d仅用一两个值无法真正估计密度。但请继续阅读...

Update:

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

• 看看内部geom_density_2d我们看到它使用了MASS::kde2d计算分布的包函数。
• 看着kde2d我们看到它使用MASS::bandwidth.nrd(df$x)以获得带宽的估计。
• 查看帮助（其中有代码）bandwidth.nrd我们看到它使用了一个经验法则quantile分布，并从第一个分位数减去第二个分位数以获得带宽估计。
• 对原始数据进行分位数，我们发现数据的分位数为零。
• 和跑步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

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

You can let that 0.1 value approach zero, eventually it will not be able to calculate a distribution anymore and you will get your error again.

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

希望有帮助。