当向量化不可行时，tidyverse如何迭代数据帧中的行？

Question

我希望知道迭代数据框行的最佳方法，当第n行的变量的值依赖于第n-1行和/或n-2行的变量时。理想情况下，我想以“tidyverse”方式完成此操作，可能使用purrr::pmap()函数。

例如，假设我有以下数据框:

library(dplyr)

x <- tibble(t = c(1:10),
            a = c(seq(100, 140, 10), rep(NA_real_, 5)),
            b = c(runif(5), rep(NA_real_, 5)),
            c = c(runif(5), rep(NA_real_, 5)))

x
#> # A tibble: 10 x 4
#>        t     a      b         c
#>    <int> <dbl>  <dbl>     <dbl>
#>  1     1   100  0.750  0.900   
#>  2     2   110  0.898  0.657   
#>  3     3   120  0.731  0.000137
#>  4     4   130  0.208  0.696   
#>  5     5   140  0.670  0.882   
#>  6     6    NA NA     NA       
#>  7     7    NA NA     NA       
#>  8     8    NA NA     NA       
#>  9     9    NA NA     NA       
#> 10    10    NA NA     NA

我已经知道时间(t)=5时的数值，现在想要预测未来的数值，使用下面的公式：

a = lag(a) * 1.1
b = a * lag(b)
c = b * lag(a, 2)

这段代码实现了预期的输出，但是它采用的是笨重、可怕的循环方式，适应大型数据集的能力不佳。

for(i in 1:nrow(x)) {
  x <- x %>%
    mutate(a = if_else(!is.na(a), a, lag(a, 1) * 1.1),
           b = if_else(!is.na(b), b, a * lag(b, 1)),
           c = if_else(!is.na(c), c, b * lag(a, 2)))
}

x
#> # A tibble: 10 x 4
#>        t     a        b        c
#>    <int> <dbl>    <dbl>    <dbl>
#>  1     1  100  7.50e- 1 9.00e- 1
#>  2     2  110  8.98e- 1 6.57e- 1
#>  3     3  120  7.31e- 1 1.37e- 4
#>  4     4  130  2.08e- 1 6.96e- 1
#>  5     5  140  6.70e- 1 8.82e- 1
#>  6     6  154  1.03e+ 2 1.34e+ 4
#>  7     7  169. 1.75e+ 4 2.45e+ 6
#>  8     8  186. 3.26e+ 6 5.02e+ 8
#>  9     9  205. 6.68e+ 8 1.13e+11
#> 10    10  225. 1.51e+11 2.80e+13

- Matt Cowgill

---

解决方案是否需要在NAs散布在数据框中时工作，或者它们总是在列的底部？如果可以依赖某一列一旦成为NA就“保持”NA，那么有些事情会变得更简单。 - Peter Ellis

Answer 1

我认为对于这种本质上迭代的过程，使用 for 循环确实是最好的选择。@Shree 提出的方法依赖于缺失值连续且从已知位置开始。

这是我对你的循环的小改进，我认为它更易读且速度约为你的方法的2.5倍，而且可能比结合矢量化操作和循环的方法更加可扩展。通过完全摆脱 tidyverse 并采用逐行循环，我们可以在这两个方面都获得一些效率提升：

method_peter <- function(x){
  for(i in 2:nrow(x)){
    x[i, "a"] <- ifelse(is.na(x[i, "a"]), x[i - 1, "a"] * 1.1,       x[i, "a"])
    x[i, "b"] <- ifelse(is.na(x[i, "b"]), x[i, "a"] * x[i - 1, "b"], x[i, "b"])
    x[i, "c"] <- ifelse(is.na(x[i, "c"]), x[i, "b"] * x[i - 2, "a"], x[i, "c"])
  }
  return(x)
}

毫无疑问还存在更多的效率提升空间，当然这是一个可以用C++重写的理想候选人:)

从以下数据可以看出，这种方法的速度大约是您方法的两倍：

method_matt <- function(x){
  for(i in 1:nrow(x)) {
    x <- x %>%
      mutate(a = if_else(!is.na(a), a, lag(a, 1) * 1.1),
             b = if_else(!is.na(b), b, a * lag(b, 1)),
             c = if_else(!is.na(c), c, b * lag(a, 2)))
  }
  return(x)
}

set.seed(123)
x <- tibble(t = c(1:10),
            a = c(seq(100, 140, 10), rep(NA_real_, 5)),
            b = c(runif(5), rep(NA_real_, 5)),
            c = c(runif(5), rep(NA_real_, 5)))

stopifnot(identical(method_matt(x), method_peter(x)))

library(microbenchmark)
microbenchmark(
  method_matt(x),
  method_peter(x)
)

这将返回:

Unit: milliseconds
            expr     min       lq     mean   median      uq     max neval
  method_matt(x) 24.1975 25.50925 30.64438 26.33310 31.8681 74.5093   100
 method_peter(x) 10.0005 10.56050 13.33751 11.06495 13.5913 42.0568   100

@Shree的方法再次证明要快得多，并且非常适合示例数据，但我不确定它是否足够灵活，以适用于您所有的用例。

如果有的话，我想看到tidyverse解决方案。

- Peter Ellis

Answer 2

编辑：添加了tidyverse方法

这是一种易读且灵活的tidyverse方法。缺点是它非常慢。

accumutate <- function(df, ...){
  df %>% group_by(row_number()) %>%
    nest() %>%
    pull(data) %>%
    accumulate(function(x,y){bind_rows(x,y) %>% mutate(!!!enquos(...)) }) %>%
    .[[length(.)]]
}

x %>%
  accumutate(a = ifelse(is.na(a), 1.1 * lag(a,1), a)) %>%
  accumutate(b = ifelse(is.na(b), a * lag(b), b)) %>%
  accumutate(c = ifelse(is.na(c),b * lag(a, 2), c))

#> # A tibble: 10 x 4
#>        t     a        b        c
#>    <int> <dbl>    <dbl>    <dbl>
#>  1     1  100  2.88e- 1 4.56e- 2
#>  2     2  110  7.88e- 1 5.28e- 1
#>  3     3  120  4.09e- 1 8.92e- 1
#>  4     4  130  8.83e- 1 5.51e- 1
#>  5     5  140  9.40e- 1 4.57e- 1
#>  6     6  154  1.45e+ 2 1.88e+ 4
#>  7     7  169. 2.45e+ 4 3.43e+ 6
#>  8     8  186. 4.57e+ 6 7.04e+ 8
#>  9     9  205. 9.37e+ 8 1.59e+11
#> 10    10  225. 2.11e+11 3.94e+13

^{由 reprex package (v0.3.0) 于2020年10月07日创建}

---

这里有另一种方法，你可能会觉得有趣。它不够简洁易读，但它受到了tidyverse（或至少是功能上）的启发。而且它表现得相当不错。它使用了半群模式，将mutate表达式转换为二元函数，创建相应的列表，然后使用accumulate。

library(tidyverse)
library(dplyr)
library(microbenchmark)
options(width =100)

set.seed(123)

# Create the data frame
x <- tibble(t = c(1:100),
            a = c(seq(100, 140, 10), rep(NA_real_,100- 5)),
            b = c(runif(5), rep(NA_real_, 100-5)),
            c = c(runif(5), rep(NA_real_, 100-5)))

a_mappend <- function(a1, a2) {
  ifelse(is.na(a2), a1 * 1.1, a2)
}

b_mappend <- function(ab1, ab2) {
    list(a = ab2$a, b =  ifelse(is.na(ab2$b), ab2$a * ab1$b,ab2$b))
}

c_mappend <- function(abc12, abc23) {
  list(abc1 = list(a = abc12$abc2$a, b = abc12$abc2$b, c = abc12$abc2$c),
       abc2 = list(a = abc23$abc2$a, b = abc23$abc2$b, c = ifelse(is.na(abc23$abc2$c),abc12$abc1$a * abc23$abc2$b,abc23$abc2$c)))
}

method_ian <- function(x) {
  x %>%
    mutate(a = accumulate(a, a_mappend)) %>%
    mutate(b = list(a, b) %>% 
                  pmap(~ list(a = .x, b = .y)) %>% 
                  accumulate(b_mappend) %>% map_dbl(~ .x$b)) %>%
    mutate(c = list(a, b, c, c(a[-1], NA), c(b[-1], NA), c(c[-1], NA)) %>%
                  pmap(~ list(abc1 = list(a = ..1, b = ..2, c = ..3),
                              abc2 = list(a = ..4, b = ..5, c = ..6))) %>% 
                  accumulate(c_mappend) %>% map_dbl(~ .x$abc1$c))
}


method_matt <- function(x){
  for(i in 1:nrow(x)) {
    x <- x %>%
      mutate(a = if_else(!is.na(a), a, lag(a, 1) * 1.1),
             b = if_else(!is.na(b), b, a * lag(b, 1)),
             c = if_else(!is.na(c), c, b * lag(a, 2)))
  }
  return(x)
}

method_peter <- function(x){
  for(i in 2:nrow(x)){
    x[i, "a"] <- ifelse(is.na(x[i, "a"]), x[i - 1, "a"] * 1.1,       x[i, "a"])
    x[i, "b"] <- ifelse(is.na(x[i, "b"]), x[i, "a"] * x[i - 1, "b"], x[i, "b"])
    x[i, "c"] <- ifelse(is.na(x[i, "c"]), x[i, "b"] * x[i - 2, "a"], x[i, "c"])
  }
  return(x)
}

stopifnot(identical(method_matt(x), method_ian(x)))

microbenchmark( method_matt(x), method_peter(x), method_ian(x))
#> Unit: milliseconds
#>             expr       min        lq      mean    median        uq       max neval
#>   method_matt(x) 324.90086 330.93192 337.46518 334.55447 338.38461 426.30457   100
#>  method_peter(x) 208.27498 211.60526 213.59438 212.66088 214.36421 242.59854   100
#>    method_ian(x)  13.06774  13.43105  14.30003  13.86428  14.32263  19.54843   100

^{此内容创建于2020年10月6日，使用reprex软件包(v0.3.0)}

- Ian

---

这非常巧妙。它在性能方面得分很高，但可读性不太好！如果我需要提取性能，并且有选择的话，我宁愿用Rcpp重新编写我的代码。 - Peter Ellis

@PeterEllis 可能在可读性方面存在一些小问题... ;) 同意 - Rcpp 将是性能和清晰度的最佳组合。顺便说一下，如果性能根本不是问题，我已经添加了另一个备选项。 - Ian

Answer 3

我认为在tidyverse中没有任何简单的方法可以进行基于行的计算。使用Reduce或gather + spread可能是可行的，但我不希望它们在可读性方面得分。

无论如何，好消息是，您可以使用dplyr和zoo软件包对您的计算进行向量化处理 -

"最初的回答"

x %>% 
  mutate(
    a = ifelse(is.na(a), na.locf(a) * 1.1^(t-5), a),
    b = ifelse(is.na(b), na.locf(b) * c(rep(1, 5), cumprod(a[6:n()])), b),
    c = ifelse(is.na(c), b * lag(a, 2), c)
  )

 # A tibble: 10 x 4
 t     a        b        c
 <int> <dbl>    <dbl>    <dbl>
 1     1  100  1.85e- 1 9.43e- 1
 2     2  110  7.02e- 1 1.29e- 1
 3     3  120  5.73e- 1 8.33e- 1
 4     4  130  1.68e- 1 4.68e- 1
 5     5  140  9.44e- 1 5.50e- 1
 6     6  154  1.45e+ 2 1.89e+ 4
 7     7  169. 2.46e+ 4 3.45e+ 6
 8     8  186. 4.59e+ 6 7.07e+ 8
 9     9  205. 9.40e+ 8 1.59e+11
10    10  225. 2.12e+11 3.95e+13

Data -

set.seed(2)
x <- tibble(t = c(1:10),
            a = c(seq(100, 140, 10), rep(NA_real_, 5)),
            b = c(runif(5), rep(NA_real_, 5)),
            c = c(runif(5), rep(NA_real_, 5)))

- Shree

---

该方法依赖于知道所有NAs都在第5个观察值之后开始，并且彼此连续（即不停止和开始）。 - Peter Ellis