我正在尝试模拟每月的数据面板,其中一个变量依赖于 R 中该变量的滞后值。我的解决方案非常慢。我需要大约 1000 个样本,涉及 2545 个人,每个样本在多年内每月观察一次,但第一个样本花了我的电脑 8.5 小时来构建。我怎样才能使它更快?
我首先创建一个不平衡的小组,其中包含不同的出生日期、月龄以及将进行比较以确定Outcome
xbsmall
和error
变量的人。第一个块中的所有代码都只是数据设置。
# 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
(,如果xbsmall > error
,则Outcome == 1
。但是,如果上个月的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 具有非常明显的年龄趋势。 - 仅对较小的样本进行模拟,然后估计样本大小对我需要的值的影响(此处未计算回归系数估计的分布(
一种方法是使用拆分-应用-组合方法。 我取出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 分钟。