data.table：每组过去 24 小时内的观察子集答案

【问题标题】：data.table: subset by observations in last 24 hours, per groupdata.table：每组过去 24 小时内的观察子集
【发布时间】：2016-12-04 15:02:42
【问题描述】：

假设我有一个如下所示的 data.table：

dt = data.table(group = c(1,1,1,2,2,2,3,3,3),time = c("2016-03-09T08:31:00-05:00","2016-03-08T11:31:00-05:00","2016-03-06T08:31:00-05:00",
                                                "2016-04-04T23:28:00-04:00","2016-04-10T23:28:00-04:00","2016-04-09T23:28:00-04:00",
                                                "2016-05-11T19:52:00-04:00","2016-05-10T20:52:00-04:00","2016-04-11T19:52:00-04:00"))

dt
   group                      time
1:     1 2016-03-09T08:31:00-05:00
2:     1 2016-03-08T11:31:00-05:00
3:     1 2016-03-06T08:31:00-05:00
4:     2 2016-04-04T23:28:00-04:00
5:     2 2016-04-10T23:28:00-04:00
6:     2 2016-04-09T23:28:00-04:00
7:     3 2016-05-11T19:52:00-04:00
8:     3 2016-05-10T20:52:00-04:00
9:     3 2016-04-11T19:52:00-04:00

对于此 data.table 中的每个组，我只想保留最近日期 24 小时内的观察结果。我为此制定了一个令人讨厌的解决方案，但它并没有我在大型数据集上所需的速度那么快。

library(lubridate)
set(dt,j = "time",value = ymd_hms(dt[["time"]]))
dt[,.(mostRecent = max(time),time),by = group][
  time > (mostRecent - days(1)),.(group,time)]

   group                time
1:     1 2016-03-09 13:31:00
2:     1 2016-03-08 16:31:00
3:     2 2016-04-11 03:28:00
4:     3 2016-05-11 23:52:00
5:     3 2016-05-11 00:52:00

有没有人有关于如何更优雅/更快地完成的提示？

【问题讨论】：

标签： r data.table data-manipulation

【解决方案1】：

为每个组创建截止时间的简单解决方案（假设时间已经转换）：

dt[, cutoff_time := max(time) - 24*60*60, by = group]
dt[time > cutoff_time]

编辑：

关于“GForce 优化最大值”的评论让我很好奇，所以我创建了一些更大的假数据来比较速度。请注意，integer 与 max 和 >= 配合得很好：

require(data.table)

require(microbenchmark)

N = 100000
N_g = 100

all_times = seq(from = as.POSIXct('2016-01-01 10:00:00'),
                to = as.POSIXct('2016-06-30 10:00:00'),
                by = 60)

all_times_int = as.integer(all_times)

idx = sample(seq.int(length(all_times)), N, replace = TRUE)

dt = data.table(group = sample(seq.int(N_g), N, replace = TRUE),
                time = all_times[idx],
                time_int = all_times_int[idx])

f1a = function (x) {
  x[, cutoff_time := max(time) - 24*60*60, by = group]
  x[time >= cutoff_time, list(group, time)]
}

f1b = function (x) {
  x[, cutoff_time := max(time_int) - 24*60*60, by = group]
  x[time_int >= cutoff_time, list(group, time)]
}

f2 = function (x) {
  thresh_dt = x[, .(time = max(time)), by=group][, time := time - 24*60*60]
  thresh_dt[x, on=c("group", "time"), roll=TRUE, nomatch=0][, list(group, time)]
}

microbenchmark(f1a(dt),
               f1b(dt),
               f2(dt))

Unit: milliseconds
   expr       min        lq      mean   median        uq       max neval
 f1a(dt)  9.842106 10.593243 11.593148 11.62311 12.478853 14.335338   100
 f1b(dt)  3.391178  3.763598  4.403264  4.00142  5.018182  8.335717   100
  f2(dt) 14.422669 15.701397 17.090674 16.56990 17.695653 52.926897   100

identical(f1a(dt), f1b(dt)) # TRUE
identical(f1a(dt), f2(dt)) # TRUE

编辑 2： 还有一个 N = 1,000,000 和 N_g = 10,000 组：

> microbenchmark(f1a(dt),
+                f1b(dt),
+                f2(dt),
+                times = 10)
Unit: milliseconds
    expr       min       lq      mean    median        uq      max neval
 f1a(dt) 634.91473 647.5662 670.74597 663.28238 694.29595 728.2481    10
 f1b(dt)  64.61488  67.3692  76.68925  68.42335  72.36862 113.1407    10
  f2(dt) 205.67688 208.6491 229.65610 213.59476 249.16703 278.7713    10

> microbenchmark(f1a(dt),
+                f1b(dt),
+                f2(dt),
+                times = 10)
Unit: milliseconds
    expr       min        lq     mean    median        uq       max neval
 f1a(dt) 620.11090 624.33587 645.0220 642.13648 657.74347 697.27674    10
 f1b(dt)  64.80214  67.43851  67.9140  67.99647  68.63552  69.74466    10
  f2(dt) 198.39200 199.56088 209.6908 204.60183 216.23255 241.76792    10

> microbenchmark(f1a(dt),
+                f1b(dt),
+                f2(dt),
+                times = 10)
Unit: milliseconds
    expr      min        lq      mean    median        uq      max neval
 f1a(dt) 619.2903 645.22617 656.58883 660.99508 664.82678 682.7618    10
 f1b(dt)  63.2454  67.31781  72.10255  68.19679  71.91441 106.7493    10
  f2(dt) 195.9335 210.06171 222.19868 215.75979 241.74100 245.9022    10

【讨论】：

我对 all_times_int 感到困惑。从某种意义上说，那些仍然是时代吗？
AFAIK 您可以将时间（以秒为单位）转换为integer，只需从某个起点开始数秒即可。如果您关心时差，您可以简单地使用integers，因为您不关心实际的YYYY-mm-dd HH:MM:SS 表示...使用integer 可能比使用POSIX 更快。也许在这方面有更多经验的人可以谈谈data.table::IDateTime...
好的，谢谢，听起来你只是失去了几毫秒和一个很好的表现。
我研究了一些，发现组数对您在基准测试中看到的内容很重要，而不仅仅是行数。我在下面添加了它。这在 data.table 中很常见。
哦，刚刚注意到你在 N_g 中添加了，但是，它仍然没有在你的代码中使用（仍然有 1:100 硬编码）。

【解决方案2】：

首先，将阈值放在一个表格中：

thresh_dt = dt[, .(time = max(time)), by=group][, time := time - 24*60*60][]

max 与减去一天的秒数以利用 the "GForce" optimized max 分开计算。另见?datatable.optimize。

接下来，进行滚动或非等值连接：

thresh_dt[dt, on=c("group", "time"), roll=TRUE, nomatch=0]

# or, on data.table 1.9.7+
thresh_dt[dt, on=.(group, time <= time), nomatch=0]

   group                time
1:     1 2016-03-09 13:31:00
2:     1 2016-03-08 16:31:00
3:     2 2016-04-11 03:28:00
4:     2 2016-04-10 03:28:00
5:     3 2016-05-11 23:52:00
6:     3 2016-05-11 00:52:00

基准测试。 GForce 最大和滚动的好处只有在您有足够的组时才会显现出来。我的示例数据扩展了@sbstn，因此组数是一个参数：

N = 5e6
ng = 1e5

all_times = seq(from = as.POSIXct('2016-01-01 10:00:00'),
                to = as.POSIXct('2016-06-30 10:00:00'),
                by = 60)
all_times_int = as.integer(all_times)    
idx = sample(seq.int(length(all_times)), N, replace = TRUE)
dt = data.table(group = sample(ng, N, replace = TRUE),
                time = all_times[idx],
                time_int = all_times_int[idx])

# sbstn, no gmax
system.time({
  dt[, cutoff_time := max(time) - 24*60*60, by = group]
  dt[time >= cutoff_time]
})
#    user  system elapsed 
#    8.50    0.01    8.47 

# sbstn, with gmax
system.time({
  dt[, maxtime := max(time), by = group][, cutoff_time := maxtime - 24*60*60]
  dt[time >= maxtime]
})
#    user  system elapsed 
#    4.98    0.01    4.99

# gmax and roll
system.time({
  thresh_dt = dt[, .(time = max(time)), by=group][, time := time - 24*60*60]
  thresh_dt[dt, on=c("group", "time"), roll=TRUE, nomatch=0][, list(group, time)]
})
#    user  system elapsed 
#    1.29    0.06    1.36 
# (Caveat: I didn't verify that these results match.)

我的答案涉及对行进行两次分组（一次计算最大值，再次与原始数据连接）。通过将其减少到一个分组操作，Clayton Stanley 的答案也设法很快（至少我认为这是正在发生的事情）：

system.time(dt[order(group, -time)
    ][, groupP := shift(group, type='lag')
    ][, head := is.na(groupP) | group != groupP
    ][, copy(.SD)[.SD[head == T], rTime := i.time, on=c(group='group')]
    ][time > (rTime - 24*60*60)
    ][, .(group, time)
    ][order(group, -time)
    ])
#    user  system elapsed 
#    1.32    0.25    1.14

【讨论】：

【解决方案3】：

通过计算，瓶颈可能在 max(*) 中。如果是这样：

dt[order(group, -time)
   ][, groupP := shift(group, type='lag')
   ][, head := is.na(groupP) | group != groupP
   ][, copy(.SD)[.SD[head == T], rTime := i.time, on=c(group='group')]
   ][time > (rTime - 24*60*60)
   ][, .(group, time)
   ][order(group, -time)
   ]

   group                time
1:     1 2016-03-09 13:31:00
2:     1 2016-03-08 16:31:00
3:     2 2016-04-11 03:28:00
4:     3 2016-05-11 23:52:00
5:     3 2016-05-11 00:52:00
>

【讨论】：

仅供参考，当max 是j 中唯一的东西时，会进行一些额外的优化，如dt[, max(time), by=group, verbose=TRUE] 所示
@Frank 我使用了你的测试代码，发现使用 days(*) 的减法是这段代码的瓶颈。好奇这个在你的机器上与你测试的其他人相比如何公平，因为它可能很接近。
1.14 秒，比滚动快一点，是的。（非 equi 连接，在上面的基准测试中没有显示，类似于滚动。）