dplyr用于逐行计算分位数

Question

我有一个分层数据框，每个分层都包含1000个样本，这些样本来自该层次估计值的后验分布。

mydf <- as.data.frame(lapply(seq(1, 1000), rnorm, n=100))
colnames(mydf) <- paste('s', seq(1, ncol(mydf)), sep='')

我想为每行添加分布的几个分位数的列。在经典的R语言中，我会这样写。

quants <- t(apply(mydf, 1, quantile, probs=c(.025, .5, .975)))
colnames(quants) <- c('s_lo', 's_med', 's_hi')
mydf <- cbind(mydf, quants)

我怀疑在dplyr中有一种直接的方法可以做到这一点（也许是rowwise？），但我的尝试都失败了。有什么想法吗？

- wylbur

---

plyr::adply may be a better choice: adply(mydf, 1, function(d) quantile(d, c(.025, .5, .975))) - bouncyball

请查看matrixStats包（特别是rowQuantiles函数）。它比使用for/apply要快得多。 - konvas

Answer 1

dplyr并不适用于基于行的计算。虽然你可以使用rowwise()来完成这个任务，但我不建议这样做：性能将非常差。最好的速度可能来自于期望一个matrix并且可以对行进行操作的东西。我建议使用apply。

为了简洁起见，我将使用只有5列的数据框而不是100x1000的数据框。

set.seed(2)
mydf <- as.data.frame(lapply(seq(1, 5), rnorm, n=10))
colnames(mydf) <- paste('s', seq(1, ncol(mydf)), sep='')

将数据框转换为矩阵只有在所有列的class相同时才是合理的。在这种情况下，它们都是numeric类型的，所以我们是安全的。（如果数据框中有非数字列，请在此处仅提取所需列，并稍后绑定回来。）

mymtx <- as.matrix(mydf)
apply(mymtx, 1, quantile, c(0.1, 0.9))
#         [,1]     [,2]     [,3]     [,4]     [,5]       [,6]     [,7]     [,8]     [,9]    [,10]
# 10% 1.028912 1.430939 1.999521 0.305907 1.753824 0.03267599 1.934381 1.270504 2.995816 1.489634
# 90% 4.950067 3.807735 4.881554 6.123989 4.886388 5.55628806 4.207605 4.184460 4.406384 3.782134

使用apply的一个显著特点是，结果以行为基础，可能与预期的不同。只需将其包装在t(...)中，您就会看到您期望的列。

可以使用cbind或类似的函数将其与原始数据框重新组合。

可以像这样在管道中完成此操作:

mydf %>%
  bind_cols(as.data.frame(t(apply(., 1, quantile, c(0.1, 0.9)))))
#            s1         s2        s3       s4       s5        10%      90%
# 1   0.1030855  2.4176508 5.0908192 4.738939 4.616414 1.02891157 4.950067
# 2   1.1848492  2.9817528 1.8000742 4.318960 3.040897 1.43093918 3.807735
# 3   2.5878453  1.6073046 4.5896382 5.076164 4.158295 1.99952092 4.881554
# 4  -0.1303757  0.9603310 4.9546516 3.715842 6.903547 0.30590700 6.123989
# 5   0.9197482  3.7822290 3.0049378 3.223325 5.622494 1.75382406 4.886388
# 6   1.1324203 -0.3110691 0.5482936 3.404340 6.990920 0.03267599 5.556288
# 7   1.7079547  2.8786046 3.4772373 2.274020 4.694516 1.93438093 4.207605
# 8   0.7603020  2.0358067 2.4034418 3.097416 4.909156 1.27050387 4.184460
# 9   2.9844739  3.0128287 3.7922033 3.440938 4.815839 2.99581584 4.406384
# 10  0.8612130  2.4322652 3.2896367 3.753487 3.801232 1.48963385 3.782134

我会把列名的命名交给你处理。

- r2evans

---

我很惊讶这个解决方案比我评论中使用的“adply”解决方案快了多少。 - bouncyball

有时候 plyr，dplyr 等等会更快，通常更易读，也更优雅。我不喜欢在 apply 后进行转置，但是如果没有特定的编译函数，它就不会更快了。矩阵总是比数据框更快（我相信），所以这是有道理的。 - r2evans

如果您愿意，可以将数据框转换为矩阵，将矩阵转换为数据框的强制转换隐式化：cbind(mydf, t(apply(mydf, 1, quantile, c(0.025, .5, .975))))（尽管了解这一点很好）。 - alistaire

管道版本正是我在寻找的，谢谢！我对t()很清楚，它是我最初问题规范的一部分。我怀疑我们这些精通经典R的人总是那样想。 - wylbur

Answer 2

使用类似于data.frame的结构进行逐行操作可能会非常困难，这是由于数据结构的本质所决定的。一种更有效的解决方案可能是重新塑造数据，在列中块状地进行计算，然后再将结果连接回来。通过dplyr+tidyr，可以实现以下操作：

library(dplyr)
library(tidyr)
mydf <- as_data_frame(mydf) %>% 
    mutate(id = row_number())

quants <- mydf %>% 
    gather(sample, value, -id) %>% 
    group_by(id) %>% 
    summarize(q025 = quantile(value, 0.025),
              q500 = quantile(value, 0.5),
              q975 = quantile(value, 0.975)) %>% 
    ungroup()

result <- left_join(quants, mydf)

或者，如果速度特别重要，可以使用 data.table...

library(data.table)
setDT(mydf)
mydf[, id := .I]
mydf_melt <- melt(mydf, id.vars = 'id')
quants <- mydf_melt[, as.list(quantile(value, c(0.025, 0.5, 0.975))), by = id]
setkey(quants, 'id')
setkey(mydf, 'id')
result <- quants[mydf]

- Alexey Shiklomanov

Answer 3

purrr::pmap 可以用于这种情况，它可以在列表中并行迭代，如果使用 data.frame，那么是按照行操作的。但是，如果每个项包含一个参数或者函数接受点，则它更有用；否则您需要使用 c 收集一个向量。

library(tidyverse)
set.seed(47)

mydf <- as.data.frame(lapply(seq(1000), rnorm, n = 100))
names(mydf) <- paste0('s', seq_along(mydf))

# make vector of each row; pass to quantile; convert to list; simplify to data.frame
mydf %>% pmap_df(~as.list(quantile(c(...), c(.025, .5, .975)))) %>% 
    bind_cols(mydf)    # self join to original columns

#> # A tibble: 100 × 1,003
#>      `2.5%`    `50%`  `97.5%`          s1       s2        s3       s4
#>       <dbl>    <dbl>    <dbl>       <dbl>    <dbl>     <dbl>    <dbl>
#> 1  24.52876 501.2313 974.1547  2.99469634 1.857485 4.8062449 5.412425
#> 2  25.96306 501.5381 975.4427  1.71114251 1.534527 5.0045983 4.029735
#> 3  25.36792 499.8048 974.9472  1.18540528 1.575371 2.1515656 4.537178
#> 4  27.15081 500.9932 975.3688  0.71823499 2.747321 0.9841692 3.774623
#> 5  25.77212 498.7223 974.5576  1.10877555 2.659429 4.6865536 5.448446
#> 6  25.43256 501.2437 973.7319 -0.08573747 2.198829 3.7851258 5.769600
#> 7  24.29993 500.8599 975.5050  0.01451784 1.938954 4.1822894 5.205473
#> 8  25.16637 501.8597 974.8636  1.01513086 3.492032 3.2551467 2.570020
#> 9  25.36332 500.3975 973.3588  0.74795410 3.660735 3.3051286 4.270915
#> 10 27.02456 499.8759 974.3890 -0.46575030 2.771156 3.4292355 3.372155
#> # ... with 90 more rows, and 996 more variables: s5 <dbl>, s6 <dbl>,
#> #   s7 <dbl>, s8 <dbl>, s9 <dbl>, s10 <dbl>, s11 <dbl>, s12 <dbl>,
#> #   s13 <dbl>, s14 <dbl>, ...

< p > quantile 生成的名称不是语法上的，但是可以通过在 bind_cols 前插入 set_names(c('s_lo', 's_med', 's_hi')) 来轻松替换。如果您喜欢，还有许多其他重新组装结果的方法。

- alistaire