如何在 R 中创建这个变量？

Question

考虑以下使用 R 的测试数据集：

testdat<-data.frame("id"=c(rep(1,5),rep(2,5),rep(3,5)),
                    "period"=rep(seq(1:5),3),
                    "treat"=c(c(0,1,1,1,0),c(0,0,1,1,1),c(0,0,1,1,1)),
                    "state"=c(rep(0,5),c(0,1,1,1,1),c(0,0,0,1,1)),
                    "int"=c(rep(0,13),1,1))
testdat
   id period treat state int
1   1      1     0     0   0
2   1      2     1     0   0
3   1      3     1     0   0
4   1      4     1     0   0
5   1      5     0     0   0
6   2      1     0     0   0
7   2      2     0     1   0
8   2      3     1     1   0
9   2      4     1     1   0
10  2      5     1     1   0
11  3      1     0     0   0
12  3      2     0     0   0
13  3      3     1     0   0
14  3      4     1     1   1
15  3      5     1     1   1

前 4 个变量是我所拥有的，int 是我想要创建的变量。它类似于treat 和state 之间的交互，但是这将在第 8-10 行中包含 1，这是不希望的。本质上，我只想要在state 在treat 期间发生变化时进行交互，但不是其他情况。关于如何创建这个（尤其是对于具有一百万个观察值的数据集的大规模）有什么想法吗？

编辑：为了澄清我为什么想要这个措施。我想运行类似以下回归的东西：

lm(outcome~treat+state+I(treat*state))

但只有当treat 跨越state 的变化时，我才真正对交互感兴趣。如果我要运行上述回归，I(treat*state) 会汇集我感兴趣的交互的效果，并且当treat 完全为 1 时，state 为 1。理论上，我认为这些会有两种不同的效果，所以我需要分解它们。我希望这是有道理的，我很乐意提供更多详细信息。

Answer 1

另一个基础版本也使用ave。

testdat$treat & c(0, diff(testdat$state))==1 在状态从 0 变为 1 时变为 TRUE，当处理为 1 时。testdat$treat & testdat$state 在两者均为 1 时变为 1。

testdat$int2 <- +ave(testdat$treat & c(0, diff(testdat$state))==1,
  cumsum(c(0, abs(diff(testdat$treat & testdat$state)))),
  FUN=function(x) rep(x[1], length(x)))
testdat
#   id period treat state int int2
#1   1      1     0     0   0    0
#2   1      2     1     0   0    0
#3   1      3     1     0   0    0
#4   1      4     1     0   0    0
#5   1      5     0     0   0    0
#6   2      1     0     0   0    0
#7   2      2     0     1   0    0
#8   2      3     1     1   0    0
#9   2      4     1     1   0    0
#10  2      5     1     1   0    0
#11  3      1     0     0   0    0
#12  3      2     0     0   0    0
#13  3      3     1     0   0    0
#14  3      4     1     1   1    1
#15  3      5     1     1   1    1

或者使用Reduce:

testdat$int2 <- Reduce(function(x,y) {if(y==-1) 0 else if(x==1 || y==1) 1 else 0},
 (testdat$treat & c(0, diff(testdat$state))==1) -c(0, diff(testdat$treat &
  testdat$state) == -1), accumulate = TRUE)

时间安排（继续@Rui-Barradas）：

f3 <- function(testdat) {cbind(testdat, int2=+ave(testdat$treat &
 c(0, diff(testdat$state))==1, cumsum(c(0, abs(diff(testdat$treat &
 testdat$state)))), FUN=function(x) rep(x[1], length(x))))}
f4 <- function(testdat) {cbind(testdat, int2=Reduce(function(x,y) {
 if(y==-1) 0 else if(x==1 || y==1) 1 else 0}, (testdat$treat & c(0,
 diff(testdat$state))==1) -c(0, diff(testdat$treat & testdat$state) == -1),
 accumulate = TRUE))}

microbenchmark::microbenchmark(base = f1(df1), dplyr = f2(df1),
 GKi1 = f3(df1), GKi2 = f4(df1), times = 10)
#Unit: milliseconds
#  expr       min        lq     mean    median        uq       max neval  cld
#  base 1132.7269 1188.7439 1233.106 1226.8532 1293.9901 1364.8358    10   c 
# dplyr 1376.0856 1436.4027 1466.418 1458.7240 1509.8990 1559.7976    10    d
#  GKi1  960.5438 1006.8803 1029.105 1022.6114 1065.7427 1074.6027    10  b  
#  GKi2  588.0484  667.2482  694.415  699.0845  739.5523  786.1819    10 a   

Answer 2

这是使用rle 和ave 的基本R 方式。

r <- rle(testdat$treat)
r$values <- cumsum(r$values) + seq_along(r$values)
int2 <- +(ave(testdat$state, inverse.rle(r), FUN = function(x) x != x[1]) & testdat$treat == 1)
testdat <- cbind(testdat, int2)

testdat
#   id period treat state int int2
#1   1      1     0     0   0    0
#2   1      2     1     0   0    0
#3   1      3     1     0   0    0
#4   1      4     1     0   0    0
#5   1      5     0     0   0    0
#6   2      1     0     0   0    0
#7   2      2     0     1   0    0
#8   2      3     1     1   0    0
#9   2      4     1     1   0    0
#10  2      5     1     1   0    0
#11  3      1     0     0   0    0
#12  3      2     0     0   0    0
#13  3      3     1     0   0    0
#14  3      4     1     1   1    1
#15  3      5     1     1   1    1

时间

由于问题提到性能是一个问题，实际用例数据集有 100 万行，这是我的解决方案的时间安排和r2evans 提供的时间安排。

将两个解决方案都写成函数。

library(dplyr)

f1 <- function(X){
  r <- rle(X$treat)
  r$values <- cumsum(r$values) + seq_along(r$values)
  int2 <- +(ave(X$state, inverse.rle(r), FUN = function(x) x != x[1]) & testdat$treat == 1)
  cbind(X, int2)
}

f2 <- function(X){
  X %>%
    group_by(grp = cumsum(c(FALSE, diff(treat) > 0))) %>%
    mutate(int2 = +(state > 0 & first(state) == 0 & treat > 0)) %>%
    ungroup() %>%
    select(-grp)
}

需要多少份testdat。

log2(1e6/nrow(testdat))
#[1] 16.02468

df1 <- testdat
for(i in 1:15) df1 <- rbind(df1, df1)
nrow(df1)
#[1] 491520

那是半百万，应该足够测试了。

mb <- microbenchmark::microbenchmark(
  base = f1(df1),
  dplyr = f2(df1),
  times = 10
)

rm(df1)    # tidy up
print(mb, unit = "relative", order = "median")
#Unit: relative
#  expr      min       lq     mean   median       uq      max neval
#  base 1.000000 1.000000 1.000000 1.000000 1.000000 1.000000    10
# dplyr 1.283237 1.359772 1.331494 1.369062 1.316815 1.256968    10

基础 R 解决方案的速度提高了大约 36%。

Answer 3

我确信这在基础 R 中是可能的，但这里有一个 tidyversion：

library(dplyr)
testdat %>%
  group_by(grp = cumsum(c(FALSE, diff(treat) > 0))) %>%
  mutate(int2 = +(state > 0 & first(state) == 0 & treat > 0)) %>%
  ungroup() %>%
  select(-grp)
# # A tibble: 15 x 6
#       id period treat state   int  int2
#    <dbl>  <int> <dbl> <dbl> <dbl> <int>
#  1     1      1     0     0     0     0
#  2     1      2     1     0     0     0
#  3     1      3     1     0     0     0
#  4     1      4     1     0     0     0
#  5     1      5     0     0     0     0
#  6     2      1     0     0     0     0
#  7     2      2     0     1     0     0
#  8     2      3     1     1     0     0
#  9     2      4     1     1     0     0
# 10     2      5     1     1     0     0
# 11     3      1     0     0     0     0
# 12     3      2     0     0     0     0
# 13     3      3     1     0     0     0
# 14     3      4     1     1     1     1
# 15     3      5     1     1     1     1

分组的替代逻辑使用游程编码，实际上是相同的（建议您https://stackoverflow.com/a/35313426）：

testdat %>%
  group_by(grp = { yy <- rle(treat); rep(seq_along(yy$lengths), yy$lengths); }) %>%
  # ...

在那个答案中，我希望dplyr 有一个相当于data.table 的rleid。预期的逻辑是能够按列中连续的相同值进行分组，但不是所有行中的相同值。如果你看看这个中间管道（在清理 grp 之前），你会看到

testdat %>%
  group_by(grp = { yy <- rle(treat); rep(seq_along(yy$lengths), yy$lengths); }) %>%
  mutate(int2 = +(state > 0 & first(state) == 0 & treat > 0)) %>%
  ungroup()
# # A tibble: 15 x 7
#       id period treat state   int   grp  int2
#    <dbl>  <int> <dbl> <dbl> <dbl> <int> <int>
#  1     1      1     0     0     0     1     0
#  2     1      2     1     0     0     2     0
#  3     1      3     1     0     0     2     0
#  4     1      4     1     0     0     2     0
#  5     1      5     0     0     0     3     0
#  6     2      1     0     0     0     3     0
#  7     2      2     0     1     0     3     0
#  8     2      3     1     1     0     4     0
#  9     2      4     1     1     0     4     0
# 10     2      5     1     1     0     4     0
# 11     3      1     0     0     0     5     0
# 12     3      2     0     0     0     5     0
# 13     3      3     1     0     0     6     0
# 14     3      4     1     1     1     6     1
# 15     3      5     1     1     1     6     1

但这只是一厢情愿。我想我也可以这样做

my_rleid <- function(x) { yy <- rle(x); rep(seq_along(yy$lengths), yy$lengths); }
testdat %>%
  group_by(grp = my_rleid(treat)) %>%
  # ...