计算 dplyr 中分组数据的加权中位数的 95% 置信区间

Question

我有一个包含多个组的数据集，我想使用 dplyr 计算每个组的中值。数据经过加权，在计算中位数时需要考虑权重。我从

spatstat

找到了 weighted.median 函数，它似乎工作得很好。考虑以下简化示例：

require(spatstat, dplyr)

tst <- data.frame(group = rep(c(1:5), each = 100))
tst$val = runif(500) * tst$group
tst$wt = runif(500) * tst$val

tst %>%
  group_by(group) %>%
  summarise(weighted.median(val, wt))

# A tibble: 5 × 2
  group `weighted.median(val, wt)`
  <int>                      <dbl>
1     1                      0.752
2     2                      1.36 
3     3                      1.99 
4     4                      2.86 
5     5                      3.45 

但是，我还想为这些值添加 95% 的置信区间，这让我很困惑。我考虑过的事情：

• Spatstat 也有一个

  weighted.var

  函数，但没有文档，我什至不清楚这是围绕中位数还是均值的方差。
• 这篇rcompanion post提出了各种计算中位数附近置信区间的方法，但据我所知，它们都不能处理权重。
• 这篇博文提出了一个计算 CI 和加权数据中位数的函数，这是我能找到的最接近我需要的函数。但是，它不适用于我的 dplyr 分组。我想我可以编写一个循环来一次执行一组操作并构建输出数据帧，但这似乎很麻烦。我也不完全确定我理解帖子中的函数，并对它的结果有点怀疑——例如，测试一下我对

  alpha=0.1

  的估计比对

  alpha=0.05

  的估计更广泛，这对我来说似乎是倒退的。 编辑添加：经过进一步调查，我认为如果我对 95% CI 使用

  alpha=0.95

  ，而不是

  alpha = 0.05

  ，此函数将按预期工作（至少，这会返回直观上感觉正确的值）。我还可以通过编辑返回单个 moe 值而不是一对高/低估计值，使其与 dplyr 一起使用。所以这可能是一个不错的选择 - 但我也在考虑其他选择。

某个库中是否有现有的函数可以执行我想要的操作，或者是否有其他简单的方法来实现此目的？

Answer 1

有几种方法。

您可以使用渐近公式计算样本中位数的标准误差。样本中位数渐近正态，标准误差为 1/sqrt(4 n f(m)^2)，其中 n 是观测值数量，m 是真实中位数，f(x) 是（加权）随机变量的概率密度。您可以使用基本 R 函数

density.default

和

weights

参数来估计概率密度。如果

x

是观测值向量，

w

是对应的权重向量，则

med <- weighted.median(x, w)
f <- density(x, weights=w)
fmed <- approx(f$x, f$y, xout=med)$y
samplesize <- length(x)
se <- 1/(sqrt(4 * samplesize) * fmed)
ci <- med + c(-1,1) * 1.96 * se

这依赖于几个渐近近似，因此可能不准确。样本大小还取决于权重的解释。在某些情况下，样本大小可能等于 sum(w)。

如果每组数据很少，您可以使用更简单的正态参考近似，

med <- weighted.median(x, w)
v <- weighted.var(x, w)
sdm <- sqrt(pi/2) * sqrt(v)
samplesize <- length(x)
se <- sdm/sqrt(samplesize)
ci <- med + c(-1,1) * 1.96 * se

或者，您可以使用 bootstrapping - 生成输入数据的随机重采样（通过选择索引 1, 2, ..., n 的随机重采样），提取相应的加权观测值 (x_i, w_i)，计算每个的加权中位数对数据集进行重采样，并构建 95% 置信区间。 （这种方法隐含地假设样本大小等于 n）