Создать последовательный счетчик, начиная с события и нулей перед событием для групп на панели

Для набора панельных данных (GSOEP) мне нужно создать счетчик времени, который дает мне дельту t после события, которое имеет фиктивный код 1 для этого конкретного года для каждого человека. Например. есть наблюдения для отдельного человека в течение случайного диапазона лет, например 1990-2006, с отдельной переменной, указывающей 1 для определенного события в году, например. 1996. Счетчик должен начинаться в следующем году, заканчиваться следующим человеком (id) и должен быть равен нулю до того, как произойдет событие для этого человека.

На данный момент данные выглядят так:

df <- data.frame(id= rep(c("1","2","3"), each=6), year=rep(1998:2003, times=3), event=c(0,0,1,0,0,0,0,0,0,0,1,0,0,1,0,0,0,0), stringsAsFactors=FALSE)

   id year event
1   1 1998     0
2   1 1999     0
3   1 2000     1
4   1 2001     0
5   1 2002     0
6   1 2003     0
7   2 1998     0
8   2 1999     0
9   2 2000     0
10  2 2001     0
11  2 2002     1
12  2 2003     0
13  3 1998     0
14  3 1999     1
15  3 2000     0
16  3 2001     0
17  3 2002     0
18  3 2003     0

Что нужно, так это:

df <- data.frame(id= rep(c("1","2","3"), each=6), year=rep(1998:2003, times=3), event=c(0,0,1,0,0,0,0,0,0,0,1,0,0,1,0,0,0,0),delta=c(0,0,0,1,2,3,0,0,0,0,0,1,0,0,1,2,3,4), stringsAsFactors=FALSE)

   id year event delta
1   1 1998     0     0
2   1 1999     0     0
3   1 2000     1     0
4   1 2001     0     1
5   1 2002     0     2
6   1 2003     0     3
7   2 1998     0     0
8   2 1999     0     0
9   2 2000     0     0
10  2 2001     0     0
11  2 2002     1     0
12  2 2003     0     1
13  3 1998     0     0
14  3 1999     1     0
15  3 2000     0     1
16  3 2001     0     2
17  3 2002     0     3
18  3 2003     0     4

Как я могу этого добиться? Ближайшее, что я получил, было здесь: Создать последовательный счетчик, который перезапускается по условию в группах данных панели

Но я не знаю, как изменить его, чтобы он запускался только после того, как событие произошло один раз, и ставил нули перед событием. Также есть некоторые индивидуумы, для которых нет событий, где счетчик должен выдавать нули. Количество лет (наблюдений) для каждого человека разное, поэтому некоторые идентификаторы варьируются от 1984-1999 годов, а другие - от 1995-2015.

Вы очень помогли бы мне, и я хочу заранее поблагодарить вас за ваше время и усилия.

С уважением,

Юлий


r counter dplyr panel timedelta
person Julius    schedule 07.03.2018    source источник

Ответы (2)


arrow_upward
2
arrow_downward

Вы можете использовать group_by(id) и cumsum(cummax(event)), чтобы приблизиться - производит 1...N, начиная с event==1. Я оборачиваю его в ifelse(...), чтобы вычесть 1 из тех значений, которые равны > 0.

library(tidyverse)
df %>%
  group_by(id) %>%
  mutate(delta = ifelse(cumsum(cummax(event)) > 0, cumsum(cummax(event)) - 1, 0)) %>%
  ungroup()

# A tibble: 18 x 4
   # id     year event delta
   # <chr> <int> <dbl> <dbl>
 # 1 1      1998    0.    0.
 # 2 1      1999    0.    0.
 # 3 1      2000    1.    0.
 # 4 1      2001    0.    1.
 # 5 1      2002    0.    2.
 # 6 1      2003    0.    3.
 # 7 2      1998    0.    0.
 # 8 2      1999    0.    0.
 # 9 2      2000    0.    0.
# 10 2      2001    0.    0.
# 11 2      2002    1.    0.
# 12 2      2003    0.    1.
# 13 3      1998    0.    0.
# 14 3      1999    1.    0.
# 15 3      2000    0.    1.
# 16 3      2001    0.    2.
# 17 3      2002    0.    3.
# 18 3      2003    0.    4.
person CPak    schedule 07.03.2018
comment
Спасибо за ваш ответ, это то направление, в котором я хотел идти. Однако, когда я пропускаю его через свои данные, векторная дельта равна NULL и, похоже, не создается, но ошибки нет, хотя я удалил все NA ранее. У вас есть идеи, в чем может быть проблема? - person Julius; 07.03.2018
comment
Не назначил фрейм данных, моя беда. Еще раз спасибо за полезный ответ! - person Julius; 07.03.2018
comment
Я думаю, что этот ответ очень хорош, поскольку он очень краток, поэтому никакой критики CPak, хороший подход (+1)! Просто хотел отметить в отношении @Julius, что мой ответ дает те же результаты и в 4 раза быстрее, см. Сравнительный анализ, который я добавил в свой ответ ниже. Убедительно прошу вас указывать ваши требования в будущих вопросах. Вы не указали, что будете принимать только dplyr / tidyverse решения, что, опять же, совершенно нормально, просто укажите такие требования более четко, чтобы сэкономить время других, предлагающих другие подходы. Спасибо. - person Manuel Bickel; 07.03.2018
comment
Уважаемый Мануэль, дело в том, что я не мог заставить ваше решение работать, поскольку оно приводило к тому, что счетчик не перезапускался для новых идентификаторов, и мне было не так легко справиться, как очевидному новичку в R (моя ошибка). В связи с этим я полностью уверен, что ваш ответ, по крайней мере, одинаково хорош, но я взял CPak, поскольку он был наиболее удобным с моей точки зрения, извините, что я не могу выбрать оба ответа одновременно. Спасибо еще раз! - person Julius; 08.03.2018
comment
@Julius Хорошо, спасибо за отзыв и извините за то, что вы так разборчивы. Я согласен с тем, что для начинающих пакетов, таких как dplyr, часто есть очень хорошие готовые, краткие и понятные решения, только для повышенных требований к скорости и функциональности могут быть интересны самодельные решения, подобные моему. Желаю успехов в изучении R! - person Manuel Bickel; 11.03.2018

arrow_upward
1
arrow_downward

Возможно, не самая элегантная версия, но если ваш набор данных не слишком велик, следующие строки могут быть началом.

library(data.table)
df <- data.frame(id= rep(c("1","2","3"), each=6), year=rep(1998:2003, times=3), event=c(0,0,1,0,0,0,0,0,0,0,1,0,0,1,0,0,0,0), stringsAsFactors=FALSE)
DT <- as.data.table(df)

get_delta <- function(x) {
  if (all(x == 0)) {
    return(x)
  } else {
    event_position <- which(x == 1)
    x[event_position] <- 0
    if (event_position == length(x)) {
     return(x) 
    } else {
     x[(event_position+1):length(x)] <- seq(length(x)-event_position)
     return(x)
    }
  }
}


DT[, delta:= get_delta(event), by = c("id")]
DT
# id year event delta
# 1:  1 1998     0     0
# 2:  1 1999     0     0
# 3:  1 2000     1     0
# 4:  1 2001     0     1
# 5:  1 2002     0     2
# 6:  1 2003     0     3
# 7:  2 1998     0     0
# 8:  2 1999     0     0
# 9:  2 2000     0     0
# 10:  2 2001     0     0
# 11:  2 2002     1     0
# 12:  2 2003     0     1
# 13:  3 1998     0     0
# 14:  3 1999     1     0
# 15:  3 2000     0     1
# 16:  3 2001     0     2
# 17:  3 2002     0     3
# 18:  3 2003     0     4

n_rows <- 1e6
DT_large <- data.table(id= as.character(rep(c(1:n_rows), each=6))
                       ,year=rep(1998:2003, n_rows), 
                       event = as.vector(sapply(1:n_rows, function(x) {
                         x <- rep(0, 6)
                         x[sample(6, 1)] <- 1  
                         x
                       }))
                       ,stringsAsFactors=FALSE)

system.time(DT_large[, delta:= get_delta(event), by = c("id")])
# User      System     elapsed 
# 9.30        0.02        9.35

#some benchmarking...
library(tidyverse)
library(data.table)
library(microbenchmark)

df <- data.frame(id= rep(c("1","2","3"), each=6), year=rep(1998:2003, times=3), event=c(0,0,1,0,0,0,0,0,0,0,1,0,0,1,0,0,0,0), stringsAsFactors=FALSE)

CPak_approach <- function() {
  df %>%
    group_by(id) %>%
    mutate(delta = ifelse(cumsum(cummax(event)) > 0, cumsum(cummax(event)) - 1, 0)) %>%
    ungroup()  
}

manuelbickel_approach <- function(x) {
  DT <- as.data.table(df)
  get_delta <- function(x) {
    if (all(x == 0)) {
      return(x)
    } else {
      event_position <- which(x == 1)
      x[event_position] <- 0
      if (event_position == length(x)) {
        return(x) 
      } else {
        x[(event_position+1):length(x)] <- seq(length(x)-event_position)
        return(x)
      }
    }
  }
  DT[, delta:= get_delta(event), by = c("id")]
}


microbenchmark(
  (dplyr_approach()),
  (manuelbickel_approach())
)

# Unit: microseconds
#       expr                      min        lq     mean   median       uq       max neval
# (dplyr_approach())         3731.146 3872.6625 4098.923 3985.363 4194.183  6441.475   100
# (manuelbickel_approach())   803.705  829.5605 1148.891 1014.105 1049.829 13993.372   100
person Manuel Bickel    schedule 07.03.2018
comment
Спасибо за быстрое предложение! К сожалению, это взорвало мой барана, поскольку вектор стал слишком большим, я попытаюсь найти решение с помощью dplyr. - person Julius; 07.03.2018
comment
Насколько велики ваши данные. Я обновил пример и использовал DT с 10e6 строками. У меня это работает. (Я немного обновил функцию). Другое дело, что ваша проблема программирования не связана в первую очередь с конкуренцией между dplyr и data.table. - person Manuel Bickel; 07.03.2018