r - 有效地创建变量，指示日期变量是否在事件之前（按组）答案

【问题标题】：r - Efficiently create variable indicating if date variable precedes event (by group)r - 有效地创建变量，指示日期变量是否在事件之前（按组）
【发布时间】：2018-08-10 07:23:50
【问题描述】：

我在 data.frame 中有两个日期（date1 和 date2）和一个 id 变量：

dat <- data.frame(c('2014-02-11', '2014-05-04', '2014-05-22'), c('2014-04-12', '2014-09-22', '2014-07-04'), c('a', 'a', 'b'))
names(dat) <- c('date1', 'date2', 'id')
dat$date1 <- as.character.Date(dat$date1, format = '%Y-%m-%d')
dat$date2 <- as.character.Date(dat$date2, format = '%Y-%m-%d')
> dat
       date1      date2 id
1 2014-02-11 2014-04-12  a
2 2014-05-04 2014-09-22  a
3 2014-05-22 2014-07-04  b

我想创建一个新变量 var 来指示 any date2 日期值是否在该行的 date1 日期值之前（不仅仅是紧接在前面的 date2 值它）：

> dat
       date1      date2 id var
1 2014-02-11 2014-04-12  a   0
2 2014-05-04 2014-09-22  a   1
3 2014-05-22 2014-07-04  b   0

我已经能够通过以下循环实现这一点：

ids <- as.vector(unique(unlist(dat$id)))
dat$var <- as.numeric(0)
for (i in ids) {
  date2s <- as.vector(unlist(filter(dat, id == i)$date2))
  for (j in date2s) {
    dat <- dat %>% mutate(var = replace(var, (j < date1) & (id == i), 1)) # if any cdate precedes rdate
  }
}

但是，我的数据集非常大，如果可能的话，我想使用data.table 来实现这一点，但如果有有效的方法，我很乐意使用dplyr 来解决这个问题。

【问题讨论】：

你的意思是date1 > min(date2)?
我的意思是如果给定行中存在任何小于date1 值的date2 值。

标签： r date group-by dplyr data.table

【解决方案1】：

data.table 和 dplyr 都不是，但首先要编写一个函数，假设列未分组

function(x, y)
    as.Date(x) > min(as.Date(y))

然后使用split()将数据分组，Map()将函数应用于每个组，split<-()分配新值

answer <- logical(nrow(dat))
split(answer, dat$id) <-
    Map(fun, split(dat$date1, dat$id), split(dat$date2, dat$id))

这将是相对有效的，即使有大量数据，只要没有太多的组。不确定为什么日期在示例数据中被转换为字符； fun() 可以用其他方式概括。

对于使用@chinsoon12 中的数据进行计时（实际上只有几组），我有

df <- as.data.frame(dat)
mtm1 <- function(df) {
    answer <- logical(nrow(dat))
    split(answer, df$id) <-
        Map(fun, split(df$date1, df$id), split(df$date2, df$id))
    answer
}

与

> identical(mtm1(df), frankMtd()$v)
[1] TRUE
> microbenchmark::microbenchmark(frankMtd(), mtm(df), times=5L)
Unit: milliseconds
       expr        min        lq       mean     median         uq        max
 frankMtd() 1917.95697 1927.2548 1928.65821 1928.45893 1933.34159 1936.27878
   mtm1(df)   47.00293   47.0198   48.02849   47.10012   47.18432   51.83523
 neval cld
     5   b
     5  a

如果有 1000 个组 (id = sample(1000, N, replace = TRUE))，那么时间会更均匀

Unit: milliseconds
       expr       min        lq      mean    median        uq      max neval
 frankMtd() 140.87859 140.88647 141.97093 141.86977 142.28619 143.9336     5
   mtm1(df)  61.82032  64.55505  64.61313  65.53642  65.53768  65.6162     5
 cld
   b
  a

通过将日期值向量化强制转换为数字可以显着提高速度

mtm2 <- function(df) {
    answer <- logical(nrow(df))
    split(answer, df$id) <- Map(
        function(x, y) x > min(y),
        split(as.numeric(df$date1), df$id),
        split(as.numeric(df$date2), df$id)
    )
    answer
}

在 1e4 组中有 1e5 个值，使用 id 一个因子（），与最快的 frank_*() 相比，结果是

> identical(frank_any()$v, mtm1(df))
[1] TRUE
> identical(frank_any()$v, mtm2(df))
[1] TRUE

和

Unit: milliseconds
        expr       min        lq      mean    median        uq       max neval
 frank_any()  79.90262  80.43112  81.79228  81.18565  83.18963  84.25236     5
    mtm1(df) 237.00027 241.40299 244.83638 246.26495 249.47713 250.03658     5
    mtm2(df)  44.11074  46.17133  51.26976  47.03285  52.77204  66.26184     5
 cld
  b 
   c
 a

【讨论】：

如果您使用的是chinsoon12的数据，则不需要fun中的as.Date转换。顺便说一句，看起来如果我计算 min，它会缩短时间以支持非 equi 连接：chat.stackoverflow.com/transcript/message/41465718#41465718
@Frank 它是 data.table 应该快速完成的那种操作，所以我并不惊讶有比基本 R 更快的实现。我实际上无法读取数据。表语法，以及基本 R 中的总体策略——编写一个函数来完成这项工作，然后弄清楚如何在组上执行它（通常是 split() <- lapply(split(), fun) 或在这种情况下使用 Map() 的扭曲）似乎很一般并且相当快，因此是一个很好的工具。
是的，它是一个非常好的工具并且易于阅读。我总是忘记split<-。我只想到我是否可以用ave 干净地做到这一点，当我做不到时，放弃并转到 data.table。

【解决方案2】：

在目前其他三个答案的基础上...

library(data.table)

frank_first = function() dat[, v0 := as.logical(copy(.SD)[copy(.SD), on=.(id, date2 < date1), mult="first", .N, by=.EACHI]$N)]

frank_which = function() dat[, vw := !is.na(copy(.SD)[copy(.SD), on=.(id, date2 < date1), mult="first", which=TRUE])]

frank_any = function() dat[, v1 := .SD[copy(.SD), on=.(id, date2 < date1), .N, by=.EACHI]$N > 0L]

frank_min = function() dat[, v := as.logical(.SD[, min(date2), by=id][copy(.SD), on=.(id, V1 < date1), .N, by=.EACHI]$N)]

fun = function(x, y) x > min(y)
mtm <- function(df) {
    df$var <- NA  # new column, to be updated
    split(df$var, df$id) <-
        Map(fun, split(df$date1, df$id), split(df$date2, df$id))
    df
}

由于an open issue/bug，需要copy 的东西。

以 chinsoon + Martin Morgan 的数据为基准：

set.seed(2L)
N <- 1e5
ng = 1e4
dat <- data.table(date1=sample(seq(as.Date("1970-01-01"), Sys.Date(), by="1 day"), N, replace=TRUE), 
    date2=sample(seq(as.Date("1970-01-01"), Sys.Date(), by="1 day"), N, replace=TRUE),
    id=sample(ng, N, replace=TRUE))

df = data.frame(dat)

microbenchmark::microbenchmark(frank_first(), frank_which(), frank_any(), frank_min(), mtm(df), times=5L)

Unit: milliseconds
          expr       min        lq      mean    median        uq       max neval cld
 frank_first()  70.38654  70.72610  80.37284  73.33607  86.87363 100.54186     5  a 
 frank_which()  55.90631  57.16385  62.89525  61.82535  64.63895  74.94178     5  a 
   frank_any()  38.56254  39.42893  40.53816  39.85976  41.47074  43.36885     5  a 
   frank_min()  36.73850  36.90551  62.55768  45.44839  55.41056 138.28545     5  a 
       mtm(df) 186.44924 190.26654 209.38918 219.73829 224.06300 226.42884     5   b

因此，最小方式（受 Martin Morgan 回答的启发）在此示例数据中获胜。

【讨论】：

这只是 thelatemail 答案的一个变体，但我是应 OP 的要求发布的。
谢谢！是的，但这保留了 data.table 的其余部分（假设还有其他变量），因此更实用。
是的，那个错误很不幸......我认为它与我链接的问题/错误有关。根据我的经验，包装 copy(.SD) 而不是 .SD 通常可以修复它。是的，另一个 var 应该与 on= 中的定义一起使用。您可能还需要/想要使用cbind，如下所示：stackoverflow.com/a/45663715 如果我们无法在 cmets 中弄清楚，您当然可以发布一个新问题。
@Frank，明白了。非常感谢！
FWIW 将日期强制转换为数字有助于基本 R 实现，我用这个更新了我的回复。

【解决方案3】：

按照@thelatemail 的建议，在自加入之后使用.EACHI 的建议如下

dat[dat, .(date1=i.date1, date2=i.date2, var=any(date2 < i.date1)), by=.EACHI, on=.(id)]

#   id      date1      date2   var
#1:  a 2014-02-11 2014-04-12 FALSE
#2:  a 2014-05-04 2014-09-22  TRUE
#3:  b 2014-05-22 2014-07-04 FALSE

编辑：一些时间供参考

set.seed(2L)
N <- 1e5
dat <- data.table(date1=sample(seq(as.Date("1970-01-01"), Sys.Date(), by="1 day"), N, replace=TRUE), 
    date2=sample(seq(as.Date("1970-01-01"), Sys.Date(), by="1 day"), N, replace=TRUE),
    id=sample(letters, N, replace=TRUE))

dt1 <- copy(dat)
tlmMtd <- function() {
    dt1[, rownum := .I]
    dt1[dt1[dt1, on="id", rownum[i.date2 < date1], allow.cartesian=TRUE], hit := 1]
}

dt2 <- copy(dat)
csMtd <- function() dt2[dt2, .(date1=i.date1, date2=i.date2, var=any(date2 < i.date1)), by=.EACHI, on=.(id)]


dt3 <- copy(dat)
frankMtd <- function() dt3[, v := .SD[copy(.SD), on=.(id, date2 < date1), .N, by=.EACHI]$N > 0L]

microbenchmark::microbenchmark(
    tlmMtd(),
    csMtd(),
    frankMtd(),
    times=5L)

# Unit: milliseconds
#       expr        min         lq       mean     median         uq       max neval
# tlmMtd()   18528.9799 18652.2217 23486.4213 19116.8014 21140.5923 39993.511     5
# csMtd()     3801.2146  3943.6201  4984.6274  5341.4322  5673.6878  6163.182     5
# frankMtd()   176.4477   177.5576   191.9636   178.9564   182.0311   244.825     5

【讨论】：

我认为这会比我的努力更快，因为.EACHI 不应该在进行计算之前返回连接的结果。我只是为了让它工作而陷入纠结。
@thelatemail 我在尝试使用.EACHI 时陷入了困境。以前从未使用过它并返回很多东西并弄清楚我是否应该这样做。或者我不。甚至：=。需要了解更多
感谢@chinsoon12，但它会删除 data.table 中的所有其他变量，因为您正在选择特定的列（如果将其应用于非玩具数据会出现问题）。
点赞数具有误导性...请参阅stackoverflow.com/questions/49063024/…

【解决方案4】：

我很确定这可以通过data.table 中的自加入来实现。例如：

library(data.table)

setDT(dat)
dat[, rownum := .I]
dat[dat[dat, on="id", rownum[i.date2 < date1]], hit := 1]
dat

#        date1      date2 id rownum hit
#1: 2014-02-11 2014-04-12  a      1  NA
#2: 2014-05-04 2014-09-22  a      2   1
#3: 2014-05-22 2014-07-04  b      3  NA

我基本上创建了一个行参考号，然后将表格加入到自身on"id"，找到日期比较符合预期的行，然后使用这些行号分配最终的hit变量。

【讨论】：

一个不同的自我加入：dat[, v := .SD[copy(.SD), on=.(id, date2 < date1), .N, by=.EACHI]$N > 0L] 由于一个未解决的问题，需要copy：github.com/Rdatatable/data.table/issues/1926
当我尝试@thelatemail 的建议时，我收到以下错误：Error in vecseq(f__, len__, if (allow.cartesian || notjoin || !anyDuplicated(f__, : Join results in 546303250 rows; more than 2304900 = nrow(x)+nrow(i). Check for duplicate key values in i each of which join to the same group in x over and over again. If that's ok, try by=.EACHI to run j for each group to avoid the large allocation. If you are sure you wish to proceed, rerun with allow.cartesian=TRUE.
鉴于该消息，我尝试了@Frank 提出的建议，该建议有效。但是，结果变量是合乎逻辑的。我如何将其转化为指标变量？
@kathystehl 好吧，逻辑 is 是一个指示变量 :) 但是，如果您想要 0/1 代替 false/true，as.integer 应该这样做，例如dat[, v := as.integer(.SD[copy(.SD), on=.(id, date2 < date1), .N, by=.EACHI]$N > 0L)]
哈！是的，是的，我想是的。谢谢@弗兰克！我没有意识到您可以使用data.table 生成类似的新变量（我正在慢慢尝试从dplyr 转换）。您能否写下您的评论作为答案，以便我接受？