比较 PDF 比率的一维核估计：如何设置尾部？

Question

尝试根据数据观察（每周时间点）创建一维内核。我一直在使用

density()

R 包中的

stats

函数。我这样做是为了可以沿 x 轴获取处理 pdf 和对照 pdf 的比率，然后用它执行排列测试。

长话短说，根据我的观察，我知道在第 4 周和第 16 周（函数中的 to 和 from 值），没有任何观察结果。然而，我的尾巴经常超出这些范围。

从统计的角度来看，我不确定如何“强制”函数在这些点处用 0 来估计 pdf。在正确的点上有 0（所有处理之间都相同）对于我能够进行下游分析至关重要。

我已经尝试了

bw

、

adjust

、

n

和

to/from

参数，并且还尝试在已知端手动将密度值设置为 0。我只是不知道从这里开始做什么，但我知道我正在做一些非常天真的事情......

我的代码示例是：

common_x = seq(4,16,by=0.02)

pdf_ctrl = density(control$week, bw=0.7, from=min(common_x),to=max(common_x), n=500)
pdf1 = density(treatment1$week, bw=0.7, from=min(common_x),to=max(common_x), n=500)

然后取 x 的每个点的比率（因为它们具有相同的公共 x 值）

不确定统计数据交换是否有更好的问题...不清楚我是否有统计问题或 R 代码问题。

没有虚拟代码，因为它非常简单，但下面是我正在谈论的内容的绘图示例，如果有帮助，很乐意提供更多代码/数据。

Answer 1

要理解为什么密度不为零，即使您知道计数在某个点为零，您需要了解如何生成一维核估计。

本质上，

density

在 x 轴上有数据的每个点放置一个内核。默认情况下，该内核是正态分布，标准差由

bw

参数控制。

每个点的内核高度按数据点总数缩放，因此，如果有 100 个数据点，钟形曲线将是正态分布高度的 1/100。

density

然后将所有这些钟形曲线相加，得到最终结果。

让我们用一个 5 到 15 之间的整数短向量来演示：

df <- data.frame(x = c(5, 5, 5, 7, 8, 8, 8, 9, 9, 10, 10, 
                       10, 10, 11, 11, 13, 13, 14, 15))

现在让我们绘制该向量的密度曲线：

library(tidyverse)

dens <- density(df$x)

p <- ggplot(as.data.frame(dens[c("x", "y")]), aes(x, y)) +
  geom_line() +
  geom_vline(xintercept = 4, linetype = 2) +
  scale_x_continuous(breaks = 1:20)  +
  theme_bw(base_size = 16) +
  scale_y_continuous(limits = c(0, 0.15), expand = c(0, 0)) +
  guides(fill = "none")

p

请注意，即使我们在 x = 4 处没有值，这里的密度也不为 0。

为了找出原因，让我们重新创建

density

使用的算法。在我们有数据的每个点上，我们都会得到一些正态分布。高度根据每个点的计数进行缩放：

df2 <- df %>% 
  reframe(xval = seq(0, 20, 0.1),
          y = dnorm(seq(0, 20, 0.1), first(x), bw.nrd0(df$x)) * n()/nrow(df), 
          .by = x)

p + geom_area(aes(x = xval, fill = factor(x)), 
              data = df2, position = "identity")

请注意，如果我们将所有这些钟形曲线堆叠在一起，我们就得到了核密度估计：

p + geom_area(aes(x = xval, fill = factor(x)), data = df2, position = "stack")

我们可以立即看到，x = 4 处密度非零的原因是 x = 5 处的内核已“溢出”到 4 上。理论上，我们可以简单地缩小带宽，这样就不会发生这种情况:

plot(density(df$x, bw = 0.3))

但是现在我们当然有一个新问题：我们漂亮的平滑密度内核看起来非常尖锐并且不再平滑。

所以尾巴的大小与密度的平滑程度有着千丝万缕的联系。如果您希望尾部在 4 处达到零，则需要较低的带宽，但这样您就会得到块状密度图。如果你想要一个平滑的密度图，你不能指望它会在计数为 0 的地方神奇地减少。如果你手动压平尾巴 pist-hoc，那么密度曲线下的面积不再是 1，因此它不是一个真正的 p.d.f，除非你手动标准化它。

当然，这里的纯粹主义者会指出，一个简单的计数图在这里是最好的；简单、完整、没有假设、没有浪费墨水：

 plot(table(df$x))

^{创建于 2024-02-06，使用 reprex v2.0.2}