R: быстро симулировать несбалансированную панель с переменной, которая зависит от самих запаздывающих значений

Question

Я пытаюсь смоделировать месячные панели данных, где одна переменная зависит от запаздывающих значений этой переменной в R. Мое решение очень медленное.Мне нужно около 1000 образцов из 2545 человек, каждый из которых наблюдается ежемесячно в течение многих лет, но первый образец потребовался моему компьютеру за 8,5 часов.Как я могу сделать это быстрее?

Я начинаю с создания несбалансированной группы людей с разными датами рождения, месячным возрастом и переменными xbsmall и error, которые будут сравниваться для определения Outcome.Весь код в первом блоке является просто настройкой данных.

# Setup:
library(plyr)

# Would like to have 2545 people (nPerson). 
#Instead use 4 for testing.
nPerson = 4
# Minimum and maximum possible ages and birth dates
AgeMin = 10
AgeMax = 50
BornMin = 1950
BornMax = 1963

# Person-specific characteristics
ind = 
  data.frame(
    id = 1:nPerson,
    BornYear = floor(runif(length(1:nPerson), min=BornMin, max=BornMax+1)),
    BornMonth = ceiling(runif(length(1:nPerson), min=0, max=12))
  )

# Make an unbalanced panel of people over age 10 up to year 1986
# panel = ddply(ind, ~id, transform, AgeMonths = BornMonth)
panel = ddply(ind, ~id, transform, AgeMonths = (AgeMin*12):((1986-BornYear)*12 + 12-BornMonth))

# Set up some random variables to approximate the data generating process
panel$xbsmall = rnorm(dim(panel)[1], mean=-.3, sd=.45)
# Standard normal error for probit
panel$error = rnorm(dim(panel)[1])

# Placeholders
panel$xb = rep(0, dim(panel)[1])
panel$Outcome = rep(0, dim(panel)[1])

Теперь, когда у нас есть данные, вот часть, которая работает медленно (около секунды на моем компьютере всего за 4 наблюдения, а часы за тысячинаблюдений).Каждый месяц человек получает два розыгрыша (xbsmall и error) из двух разных нормальных распределений (это было сделано выше) и Outcome == 1, если xbsmall > error.Однако, если Outcome равно 1 в предыдущем месяце, то Outcome в текущем месяце равно 1, если xbsmall + 4.47 > error.Я использую xb = xbsmall+4.47 в приведенном ниже коде (xb - это «линейный предиктор» в пробитной модели).Я игнорирую первый месяц для каждого человека для простоты.К вашему сведению, это симуляция пробитного DGP (но это не обязательно знать, чтобы решить проблему скорости вычислений).

# Outcome == 1 if and only if xb > -error
# The hard part: xb includes information about the previous month's outcome
start_time = Sys.time()
for(i in 1:nPerson){
  # Determine the range of monthly ages to loop over for this person
  AgeMonthMin = min(panel$AgeMonths[panel$id==i], na.rm=T)
  AgeMonthMax = max(panel$AgeMonths[panel$id==i], na.rm=T)
  # Loop over the monthly ages for this person and determine the outcome
  for(t in (AgeMonthMin+1):AgeMonthMax){
    # Indicator for whether Outcome was 1 last month
    panel$Outcome1LastMonth[panel$id==i & panel$AgeMonths==t] = panel$Outcome[panel$id==i & panel$AgeMonths==t-1] 
    # xb = xbsmall + 4.47 if Outcome was 1 last month
    # Otherwise, xb = xbsmall
    panel$xb[panel$id==i & panel$AgeMonths==t] = with(panel[panel$id==i & panel$AgeMonths==t,], xbsmall + 4.47*Outcome1LastMonth)
    # Outcome == 1 if xb > 0
    panel$Outcome[panel$id==i & panel$AgeMonths==t] =
      ifelse(panel$xb[panel$id==i & panel$AgeMonths==t] > - panel$error[panel$id==i & panel$AgeMonths==t], 1, 0)
    }
  }
end_time = Sys.time()
end_time - start_time

Мои мысли по сокращению компьютерного времени:

• Что-то с cumsum()
• Какая-то замечательная функция данных панели, о которой я не знаю
• Найдите способ заставить цикл t проходить одинаковые начальную и конечную точки для каждого отдельного человека изатем каким-то образом используйте plyr::ddpl() или dplyr::gather_by()
• Итеративное решение: сделайте обоснованное предположение о значении Outcome в каждом месячном возрасте (скажем, в режиме) и каким-то образом скорректируйте значения, которые не соответствуют предыдущемумесяц.Это будет работать лучше в моем реальном приложении, потому что xbsmall имеет очень четкую тенденцию в возрасте.
• Выполните моделирование только для небольших выборок, а затем оцените влияние размера выборки на нужные мне значения (распределения коэффициента регрессии).оценки здесь не рассчитываются)

Answer 1

Один из подходов заключается в использовании метода разделения-применения-объединения.Я вынимаю цикл for(t in (AgeMonthMin+1):AgeMonthMax) и помещаю содержимое в функцию:

generate_outcome <- function(x) {
  AgeMonthMin <- min(x$AgeMonths, na.rm = TRUE)
  AgeMonthMax <- max(x$AgeMonths, na.rm = TRUE)
  for (i in 2:(AgeMonthMax - AgeMonthMin + 1)){
    x$xb[i] <- x$xbsmall[i] + 4.47 * x$Outcome[i - 1]
    x$Outcome[i] <- ifelse(x$xb[i] > - x$error[i], 1, 0)
  }
  x
} 

, где x - это кадр данных для одного человека.Это позволяет нам упростить конструкцию panel$id==i & panel$AgeMonths==t.Теперь мы можем просто сделать

out <- lapply(split(panel, panel$id), generate_outcome)
out <- do.call(rbind, out)

и all.equal(panel$Outcome, out$Outcome) возвращает TRUE.Для вычисления 100 человек потребовалось 1,8 секунды, используя этот метод, по сравнению с 1,5 минутами в исходном коде.