Как мне создать эту переменную в 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, но оно будет включать 1 с в строках 8-10, что нежелательно. По сути, я хочу взаимодействия только тогда, когда state изменяется в течение treat, но не иначе. Любые мысли о том, как создать это (особенно в крупном масштабе для набора данных с миллионом наблюдений)?

Редактировать: Для пояснения, почему я хочу эту меру. Я хочу запустить что-то вроде следующей регрессии:

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

Но я действительно заинтересован во взаимодействии, только когда treat перекрывает изменение state. Если бы я запустил вышеупомянутую регрессию, I(treat*state) объединяет эффект взаимодействия, в котором я заинтересован, и когда treat равен 1 полностью, а state равен 1. Теоретически, я думаю, что это будет иметь два разных эффекта, поэтому Мне нужно разбить их на части. Я надеюсь, что это имеет смысл, и я рад предоставить дополнительную информацию.

Answer 1

Я уверен, что это возможно в базе R, но вот тидиверсия:

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

Альтернативный лог c для группировки использует кодирование по длине прогона, фактически то же самое (предложил вам { ссылка }):

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

И, как и в этом ответе, у меня sh dplyr был эквивалент data.table rleid. Ожидаемый лог c должен иметь возможность группировать по последовательным одинаковым значениям в столбце, но не по одному значению во всех строках. Если вы посмотрите на эту среднюю трубу (до очистки 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)) %>%
  # ...

Answer 2

Вот базовый способ R с использованием rle и ave.

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

Времена

Поскольку в вопросе производительность упоминается как проблема, набор данных реального варианта использования имеет 1 миллион строк, вот время моего решения, а число 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

Другая версия base с использованием также ave.

testdat$treat & c(0, diff(testdat$state))==1 переходит на TRUE, когда состояние изменяется с 0 на 1, когда значение равно 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