在R中计算加权分组平均值的快速方法？答案

【问题标题】：Fast way to compute weighted group-wise means in R?在R中计算加权分组平均值的快速方法？
【发布时间】：2018-08-30 04:03:28
【问题描述】：

给定纵向数据，我如何计算一个矩阵，其中每列代表给定变量的加权分组平均值？

我开发了一种需要循环的方法，但它太慢了。我认为这可能是矢量化的，但解决方案却让我望而却步。

这是我目前的方法：

library(foreach)

# N is sample size
# g is the number of groups
# p is the number of variables
get_group_mean_matrix <- function(N, g, p){
  X <- matrix(rbinom(N*p, 10, .5), N)
  f <- sort((1:(N)) %% g + 1)
  w <- runif(N)
  dmmat <- foreach(i = unique(f), .combine = rbind) %do% {
    idx <- which(f == i)
    ws <- w[idx]/sum(w[idx])
    t((t(X[idx,]) %*% ws)) %x% rep(1, length(idx))
  }
  dmmat
}

> set.seed(666)
> get_group_mean_matrix(12, 3, 5)
          [,1]     [,2]     [,3]     [,4]     [,5]
 [1,] 5.261103 4.074266 5.828070 4.452703 5.990165
 [2,] 5.261103 4.074266 5.828070 4.452703 5.990165
 [3,] 5.261103 4.074266 5.828070 4.452703 5.990165
 [4,] 5.261103 4.074266 5.828070 4.452703 5.990165
 [5,] 5.560556 4.241942 3.698828 5.572523 4.212532
 [6,] 5.560556 4.241942 3.698828 5.572523 4.212532
 [7,] 5.560556 4.241942 3.698828 5.572523 4.212532
 [8,] 5.560556 4.241942 3.698828 5.572523 4.212532
 [9,] 4.289029 4.771115 5.150607 4.424339 6.346775
[10,] 4.289029 4.771115 5.150607 4.424339 6.346775
[11,] 4.289029 4.771115 5.150607 4.424339 6.346775
[12,] 4.289029 4.771115 5.150607 4.424339 6.346775
> library(microbenchmark)
> microbenchmark(get_group_mean_matrix(1200, 300, 50))
Unit: milliseconds
                                 expr      min       lq     mean   median       uq      max neval
 get_group_mean_matrix(1200, 300, 50) 76.33337 77.39607 80.76586 78.39808 84.46984 93.40047   100

最初，我尝试使用lfe::demeanlist 执行此操作，但它给了我错误的输出！

library(lfe)
get_group_mean_matrix_lfe <- function(N, g, p){
  X <- matrix(rbinom(N*p, 10, .5), N)
  f <- sort((1:(N)) %% g + 1)
  w <- runif(N)
  X - demeanlist(X, list(factor(f)), weights = w)
}
> set.seed(666)
> get_group_mean_matrix_lfe(12, 3, 5)
          [,1]     [,2]     [,3]     [,4]     [,5]
 [1,] 5.138068 4.001781 5.415467 4.722947 5.999827
 [2,] 5.138068 4.001781 5.415467 4.722947 5.999827
 [3,] 5.138068 4.001781 5.415467 4.722947 5.999827
 [4,] 5.138068 4.001781 5.415467 4.722947 5.999827
 [5,] 5.197308 4.067657 3.202478 5.866451 4.066385
 [6,] 5.197308 4.067657 3.202478 5.866451 4.066385
 [7,] 5.197308 4.067657 3.202478 5.866451 4.066385
 [8,] 5.197308 4.067657 3.202478 5.866451 4.066385
 [9,] 4.189951 4.887720 4.953305 4.501874 6.385846
[10,] 4.189951 4.887720 4.953305 4.501874 6.385846
[11,] 4.189951 4.887720 4.953305 4.501874 6.385846
[12,] 4.189951 4.887720 4.953305 4.501874 6.385846
> library(microbenchmark)
> microbenchmark(get_group_mean_matrix_lfe(1200, 300, 50))
Unit: milliseconds
                                     expr      min       lq     mean   median       uq      max neval
 get_group_mean_matrix_lfe(1200, 300, 50) 6.107421 6.202426 6.500411 6.293648 6.582943 8.350876   100

虽然速度要快很多...

我会接受两种答案中的任何一种：

解释lfe::demeanlist 在加权情况下的作用。当我从平均值中减去加权偏差时，我不应该得到加权平均值吗？知道了这一点，我该如何计算加权分组均值矩阵？
不涉及 demeanlist 的方法来计算加权分组均值矩阵。

注意：使用RcppEigen 将%*% 替换为矩阵乘法函数可以加快速度，但还不够。我认为问题在于循环。

这是一些输入示例：

   f X1 X2 X3 X4 X5
1  1  6  5  7  3  6
2  1  6  4  5  5  6
3  1  5  6  3  6  6
4  1  3  5  4  3  5
5  2  5  4  7  7  7
6  2  4  1  4  2  6
7  2  5  6  6  6  5
8  2  6  7  2  5  4
9  3  5  3  4  6  9
10 3  6  6  5  5  6
11 3  5  7  4  6  8
12 3  5  3  7  8  6

其中f 是分组因子。

【问题讨论】：

@Gregor 该函数提供样本输入。只需在全局中设置N、p 和g 并逐行运行。不过，为了清楚起见，我将粘贴一些示例输入。
我现在明白了 - 我只是将输入生成移出函数并给它一个评论，如# x: values, f: groups, w: weights。
是的，我最初就是这样做的，但我想看看速度如何随大小变化。
但是你（大概）不关心数据生成的速度，所以你不应该在你计时的函数中包含它。半相关，我建议f <- rep(1:N, each = g) 而不是你的sort 方法，因为sort 随着向量变长会花费一些时间。
是的，我想是的。但与其他一切相比，它相当微不足道。

标签： r optimization vectorization longitudinal

【解决方案1】：

Hurr durr 我所要做的就是将权重平方根到demeanlist hurr durr

library(foreach)
get_group_mean_matrix <- function(N, g, p){
  X <- matrix(rbinom(N*p, 10, .5), N)
  f <- sort((1:(N)) %% g + 1)
  w <- runif(N)
  dmmat <- foreach(i = unique(f), .combine = rbind) %do% {
    idx <- which(f == i)
    ws <- w[idx]/sum(w[idx])
    t((t(X[idx,]) %*% ws)) %x% rep(1, length(idx))
  }
  dmmat
}

set.seed(666)
A <- get_group_mean_matrix(12, 3, 5)

library(lfe)
get_group_mean_matrix_lfe <- function(N, g, p){
  X <- matrix(rbinom(N*p, 10, .5), N)
  f <- sort((1:(N)) %% g + 1)
  w <- runif(N)
  X - demeanlist(X, list(factor(f)), weights = w^.5)
}

set.seed(666)
B <- get_group_mean_matrix_lfe(12, 3, 5)

> all.equal(A, B)
[1] TRUE

【讨论】：