代码之家 › 专栏 › 技术社区 › randy

R:用依赖于自身滞后值的变量快速模拟不平衡面板

randy · 技术社区 · 8 年前

我尝试模拟每月的数据面板,其中一个变量依赖于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个观测值,而数千个观测值只有几个小时)。每个月,一个人会得到两张抽奖单( XB小 和 错误 Outcome == 1 如果 xbsmall > error . 但是,如果 结果 上个月等于1,然后 结果 当月等于1,如果 xbsmall + 4.47 > error . 我用 xb = xbsmall+4.47 在下面的代码中( xb 是probit模型中的“线性预测因子”)。为了简单起见,我忽略了每个人的第一个月。对于您的信息,这是模拟probit 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()
迭代解:对 结果 在每个月的年龄(比如说模式),并以某种方式调整与前一个月不匹配的值。这在我的实际应用程序中会更好,因为XBSill在年龄上有一个非常明显的趋势。
只对较小的样本进行模拟,然后估计样本大小对我需要的值的影响(此处不计算回归系数估计的分布)

1 回复 | 直到 8 年前

Weihuang Wong 8 年前

一种方法是使用拆分应用组合方法。我拿出了 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分钟。