R使用dplyr按组计算加权统计

Question

多年来，我一直在使用 Hmisc 包和 base R 来计算加权统计摘要。通常，我使用双重权重，一个是空间影响权重，另一个是数据支持值，例如长度、体积、物理密度等。以“mtcars”数据集为例，其中 mpg 是感兴趣的变量，并且由汽车“wt”和“hp”设计的双重权重，Hmisc + base R 工作流程通常如下所示。

require(Hmisc)

mtcars$Wt2 <- mtcars$wt * mtcars$hp               # double weight
mtcars$Acc <- mtcars$Wt2 * mtcars$mpg             # accumulation

min(mtcars$mpg)                                   # min
sqrt(wtd.var(mtcars$mpg, mtcars$mpg))             # wtd sd
wtd.quantile(mtcars$mpg, mtcars$mpg,0.05)         # wtd 5th 
wtd.quantile(mtcars$mpg, mtcars$mpg,0.50)         # wtd median
wtd.quantile(mtcars$mpg, mtcars$mpg,0.95)         # wtd 95th
max(mtcar$mpg)                                    # max

然后可以使用循环为更大数据帧中的每个感兴趣区域过滤和计算这些加权统计信息。但是，当我尝试学习如何使用 dplyr 时，我想知道如何计算这些加权统计信息。可以使用加权平均选项，但其他选项需要更多工作。下面是我研究出如何从第一原理计算加权平均值的代码（并对照 dplyr 内置函数检查了这一点）。但是，我很困惑如何使用 dplyr 按组计算加权 SD 和分位数，因为我需要以某种方式将平方加权平均差（对于每个组）放入管道链中。

mtcars %>% 
  mutate(Car = row.names(mtcars)) %>%  # variable for car name
  mutate(Wt2 = wt * hp) %>%            # double weight
  mutate(Acc = Wt2 * mpg) %>%          # weighted consumption
  group_by(Car) %>%                    # group by car type
  summarise(n = n(),
            SmWt2 = sum(Wt2),                    # Sum of double weight
            SmAcc = sum(Acc),                    # Sum of accumulations
            WtMn = SmAcc/SmWt2,                  # Weighted mean
            WtMnChk = weighted.mean(mpg, Wt2)    # Check weighted mean
            )

Answer 1

我不确定我是否完全理解您正在使用的方法，但这里有一个使用gear 求加权平均值和加权标准差的示例，使用wt 作为权重：

library(dplyr)
datasets::mtcars %>% 
  group_by(gear) %>%
  summarize(n = n(),
            mpg_weighted_by_weight = sum(mpg*wt) / sum(wt),
            mpg_weighted_by_weight_check = weighted.mean(mpg, wt),
            
            mpg_sd = sqrt(sum(wt * ((mpg - mpg_weighted_by_weight)^2))/(sum(wt)-1)),
            mpg_sd_check = sqrt(Hmisc::wtd.var(mpg, wt)))


# A tibble: 3 x 6
   gear     n mpg_weighted_by_weight mpg_weighted_by_weight_check mpg_sd mpg_sd_check
* <dbl> <int>                  <dbl>                        <dbl>  <dbl>        <dbl>
1     3    15                   15.6                         15.6   3.32         3.32
2     4    12                   23.6                         23.6   4.81         4.81
3     5     5                   19.7                         19.7   5.63         5.63

我不熟悉加权标准差的公式，而是欺骗并依赖 Hmisc::wtd.var 的公式。如果在 RStudio 中按住 Control 键单击公式名称，它会显示函数的底层代码。大部分都是错误处理，直到底部：

#Hmisc::wtd.var
function (x, weights = NULL, normwt = FALSE, na.rm = TRUE, method = c("unbiased", 
  "ML")) 
{
  # ...  skipping error handling
  sw <- sum(weights)
  # ...
  xbar <- sum(weights * x)/sw
  sum(weights * ((x - xbar)^2))/(sw - 1)
}