Быстрый способ вычисления взвешенных групповых средств в R?

Question

Учитывая продольные данные, как я могу вычислить матрицу, в которой каждый столбец представляет взвешенное среднее по группе для данной переменной?

Я разработал подход, который требует цикла, и он слишком медленный. Я думаю, что это может быть векторизация, но решение ускользает от меня.

Вот мой текущий подход:

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

Хотя это намного быстрее ...

Я приму любой из двух видов ответов:

1. Объяснение того, что lfe::demeanlist делает в взвешенном случае. Разве я не должен получить взвешенное среднее, когда вычитаю взвешенное отклонение от среднего? И зная это, как я могу вычислить матрицу взвешенных групповых средств?
2. Способы, не включающие список действий для вычисления матрицы взвешенных групповых средств.

Примечание: замена %*% на функцию умножения матриц с использованием 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 - коэффициент группировки.

Answer 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