如何使用统计数据来加速 data.frame 上的逐行计算？

Question

我有一个包含 10,000 行和 40 列的数据框。我正在尝试对这些行中的每一行应用一个函数。对于每一行，我期望返回一个标量，它是我在这个函数中计算的统计数据的值。以下是我到目前为止所做的；

library(sandwich)

# Creating example data #

nrows=10000
ncols=40
n1=20
n2=20
df=data.frame(t(replicate(nrows, rnorm(ncols, 100, 3))))
cov=data.frame(group=as.factor(rep(c(1,2),c(n1,n2))))

# Function to evaluate on each row of df #

get_est= function(x){
mod = lm(x~cov$group)
vcov = vcovHC(mod)
coef = as.numeric(mod$coefficients[2])
se = sqrt(as.numeric(diag(vcov)[2]))
stats = coef/se
return(stats)
}

# Applying above function to full data #

t1=Sys.time()
estimates=apply(df, 1, function(x) get_est(x))
t2=Sys.time()-t1

# Time taken by apply function

Time difference of 32.10623 secs

有没有办法显着减少对完整数据实施 get_est() 的时间？我需要在单个 df 上加快计算速度的主要原因是因为我还有 1000 个具有相同维度的数据帧，并且我必须将此函数同时应用于每个数据帧的每一行。为了说明，下面是我正在处理的更广泛的情况；

# Creating example data

set.seed(1234)
nrows = 10000
ncols = 40
n1 = 20
n2 = 20
df.list = list()
for(i in 1:1000){
  df.list[[i]] = data.frame(t(replicate(nrows, rnorm(ncols, 100, 3))))
}

# Applying get_est() to each row and to each of data frame in df.list #

all.est = foreach(j = 1:length(df.list), .combine = cbind, .packages = 'sandwich') %dopar% {
  cov=data.frame(group=as.factor(rep(c(1,2),c(n1,n2))))
  est = apply(df.list[[j]], 1, function(x) get_est(x))
  return(est)
}

即使在并行化之后也需要数小时才能完成。我的最终目标是显着缩短获取“all.est”的时间，该“all.est”将包含 10000 行和 1000 列，其中每列都有相应数据集的统计估计值。任何帮助深表感谢！！提前致谢！

Answer 1

您的函数get_est 使用了一些“昂贵”的函数，例如lm、vcovHC 等。如果你想到 OLS 方程，

$$ \hat{\beta} = (X^TX)^{-1}X^Ty, $$

然后你可以看到第一部分 $(X^TX)^{-1}X^T$ 在你的模拟中没有改变，所以设计矩阵是常数。为了利用这一点，我在开始模拟之前计算了 $(X^TX)^{-1}X^T$。这种方法还需要使用公式手动计算 HC3 标准误差（参见例如here）

$$ \widehat{\text{Cov}}_{\text{HC3}}(\hat{\beta}) = (X^TX)^{-1}X^T \text{diag} \left[ \frac{ e_i^2}{(1-h_{ii})^2} \right] X(X^TX)^{-1}。 $$

在您的模拟迭代中，除了残差之外的所有内容都是恒定的，因此可以预先计算。一旦我实施了这些技巧，我的速度大约提高了 50 倍。

（注意：lm 使用 QR 分解，在这里也可以类似地实现。也许您可以通过并行化代码来提高速度。）

nrows = 10000
ncols = 40
n1 = 20
n2 = 20
df = data.frame(t(replicate(nrows, rnorm(ncols, 100, 3))))
cov = data.frame(group=as.factor(rep(c(1,2),c(n1,n2))))

# old function
get_est_old = function(x){
  mod = lm(x~cov$group)
  vcov = sandwich::vcovHC(mod)
  coef = as.numeric(mod$coefficients[2])
  se = sqrt(as.numeric(diag(vcov)[2]))
  stats = coef/se
  return(stats)
}

# new function
# first construct design matrix
X = matrix(c(rep(1, ncols), rep(0, ncols / 2), rep(1, ncols / 2)), ncol = 2)
# these quantities will be used below
inv = solve(crossprod(X)) %*% t(X)
h = diag(X %*% inv)

get_est_new= function(x){
  coef = (inv %*% x)
  resid = x - (X %*% coef)
  bread = (resid^2 / (1 - h)^2)[,1]
  hc3 = inv %*% diag(bread) %*% t(inv)
  se = sqrt(hc3[2,2])
  stats = coef[2,1]/sqrt(hc3[2,2])
}

# Applying above function to full data #

system.time({
  estimates_old = apply(df, 1, function(x) get_est_old(x))  
})
#>    user  system elapsed 
#>   7.876   0.042   7.929

system.time({
estimates_new = apply(df, 1, function(x) get_est_new(x))
})
#>    user  system elapsed 
#>   0.141   0.016   0.158

# check
all.equal(estimates_old, estimates_new)
#> [1] TRUE

^{由reprex package (v2.0.1) 于 2021-09-04 创建}

这些帖子可能很有趣：

• https://grantmcdermott.com/efficient-simulations-in-r/
• https://grantmcdermott.com/simulations-remix-turn-up-the-base/