我不是调查方法学家或人口统计学家，但我是 Thomas Lumley 的 R 调查包的狂热粉丝。我一直在处理一个相对较大的复杂调查数据集，即医疗保健成本和利用项目 (HCUP) 国家急诊室样本 (NEDS https://www.hcup-us.ahrq.gov/nedsoverview.jsp）。正如医疗保健研究和质量局所描述的，这是“来自 30 个州 947 家医院的急诊就诊的出院数据，近似于美国医院急诊的 20% 分层样本”
2006 年至 2012 年的完整数据集包含 198,102,435 个观测值。我已将数据分组为 40,073,358 例与外伤相关的出院，其中包含 66 个变量。即使对这些数据运行简单的调查程序也需要非常长的时间。我尝试过使用 RAM（2013 年末的 Mac Pro，3.7GHz 四核，128GB（！）内存），使用多核 http://r-survey.r-forge.r-project.org/survey/html/surveyoptions.html有空的时候，子集化 http://r-survey.r-forge.r-project.org/survey/html/subset.survey.design.html，与一个内存不足的数据库管理系统 https://faculty.washington.edu/tlumley/tutorials/user-biglm.pdf like MonetDB https://github.com/ajdamico/asdfree/blob/da164016bfdd533b12f40b55045d7a6007a24d12/IPUMS%20International/download%20import%20design%20into%20monetdb.R。基于设计的调查程序仍然需要几个小时。有时要好几个小时。一些不太复杂的分析需要 15 个小时以上。我猜测大部分计算工作都与巨大的协方差矩阵有关？
正如人们所预料的那样，处理原始数据的速度要快几个数量级。更有趣的是，根据程序的不同，对于如此大的数据集，未经调整的估计值可能非常接近调查结果。 （参见下面的示例）基于设计的结果显然更精确且更受欢迎，但几个小时的计算时间与几秒钟的计算时间对于增加的精度来说是一个不可忽视的成本。它开始看起来像是绕着街区走了很长一段路。
有没有人有这方面的经验？有没有办法优化大型数据集的 R 调查程序？也许更好地利用并行处理？贝叶斯方法是否使用INLA https://www.stat.washington.edu/research/reports/2011/tr583.pdf or 哈密顿量 http://www.stat.columbia.edu/~gelman/research/published/Si_et_al-BA14.pdf像斯坦这样的方法可能是解决方案吗？或者，当调查规模足够大且具有足够代表性时，一些未经调整的估计（尤其是相对指标）是否可以接受？
以下是一些未经调整的估计值近似调查结果的示例。
在第一个示例中，内存中的 svymean 花费了不到一个小时，内存不足则需要 3 个多小时。直接计算只需要不到一秒钟的时间。更重要的是，点估计（svymean 为 34.75，未调整为 34.77）以及标准误差（0.0039 和 0.0037）非常接近。
# 1a. svymean in memory
svydes<- svydesign(
id = ~KEY_ED ,
strata = ~interaction(NEDS_STRATUM , YEAR), note YEAR interaction
weights = ~DISCWT ,
nest = TRUE,
data = inj
)
system.time(meanAGE<-svymean(~age, svydes, na.rm=T))
user system elapsed
3082.131 143.628 3208.822
> meanAGE
mean SE
age 34.746 0.0039
# 1b. svymean out of memory
db_design <-
svydesign(
weight = ~discwt , weight variable column
nest = TRUE , whether or not psus are nested within strata
strata = ~interaction(neds_stratum , yr) , stratification variable column
id = ~key_ed ,
data = "nedsinj0612" , table name within the monet database
dbtype = "MonetDBLite" ,
dbname = "~/HCUP/HCUP NEDS/monet" folder location
)
system.time(meanAGE<-svymean(~age, db_design, na.rm=T))
user system elapsed
11749.302 549.609 12224.233
Warning message:
'isIdCurrent' is deprecated.
Use 'dbIsValid' instead.
See help("Deprecated")
mean SE
age 34.746 0.0039
# 1.c unadjusted mean and s.e.
system.time(print(mean(inj$AGE, na.rm=T)))
[1] 34.77108
user system elapsed
0.407 0.249 0.653
sterr <- function(x) sd(x, na.rm=T)/sqrt(length(x)) # write little function for s.e.
system.time(print(sterr(inj$AGE)))
[1] 0.003706483
user system elapsed
0.257 0.139 0.394
svymean 与使用 svyby（近 2 小时）与 tapply（4 秒左右）应用于数据子集的平均值的结果之间存在类似的对应关系：
# 2.a svyby .. svymean
system.time(AGEbyYear<-svyby(~age, ~yr, db_design, svymean, na.rm=T, vartype = c( 'ci' , 'se' )))
user system elapsed
4600.050 376.661 6594.196
yr age se ci_l ci_u
2006 2006 33.83112 0.009939669 33.81163 33.85060
2007 2007 34.07261 0.010055909 34.05290 34.09232
2008 2008 34.57061 0.009968646 34.55107 34.59014
2009 2009 34.87537 0.010577461 34.85464 34.89610
2010 2010 35.31072 0.010465413 35.29021 35.33124
2011 2011 35.33135 0.010312395 35.31114 35.35157
2012 2012 35.30092 0.010313871 35.28071 35.32114
# 2.b tapply ... mean
system.time(print(tapply(inj$AGE, inj$YEAR, mean, na.rm=T)))
2006 2007 2008 2009 2010 2011 2012
33.86900 34.08656 34.60711 34.81538 35.27819 35.36932 35.38931
user system elapsed
3.388 1.166 4.529
system.time(print(tapply(inj$AGE, inj$YEAR, sterr)))
2006 2007 2008 2009 2010 2011 2012
0.009577755 0.009620235 0.009565588 0.009936695 0.009906659 0.010148218 0.009880995
user system elapsed
3.237 0.990 4.186
调查和未调整结果之间的对应关系开始因绝对计数而崩溃，这需要编写一个吸引调查对象的小函数，并使用 Lumley 博士的一些代码来对计数进行加权：
# 3.a svytotal
system.time(print(svytotal(~adj_cost, svydes, na.rm=T)))
total SE
adj_cost 9.975e+10 26685092
user system elapsed
10005.837 610.701 10577.755
# 3.b "direct" calculation
SurvTot<-function(x){
N <- sum(1/svydes$prob)
m <- mean(x, na.rm = T)
total <- m * N
return(total)
}
> system.time(print(SurvTot(inj$adj_cost)))
[1] 1.18511e+11
user system elapsed
0.735 0.311 0.989
结果更难以让人接受。尽管仍在调查程序确定的误差范围内。但同样，为了获得更精确的结果，3 小时与 1 秒相比，成本相当可观。
更新：2016 年 2 月 10 日
感谢塞维林和安东尼允许我借用你们的突触。抱歉延迟跟进，花了很少的时间来尝试您的建议。
Severin，您的观察是正确的，革命分析/MOR 构建对于某些操作来说更快。看起来它与 CRAN R 附带的 BLAS（“基本线性代数子程序”）库有关。它更精确，但速度较慢。因此，我使用允许多线程的专有（但 Mac 上免费）Apple Accelerate vecLib 优化了我的机器上的 BLAS（请参阅http://blog.quadrivio.com/2015/06/improved-r-performance-with-openblas.html http://blog.quadrivio.com/2015/06/improved-r-performance-with-openblas.html）。这似乎减少了一些操作时间，例如从 svyby/svymean 的 3 小时到 2 小时多一点。
安东尼在复制重量设计方法方面运气不佳。 type="bootstrap" withreplicates=20 运行了大约 39 小时，然后我退出了； type =“BRR”返回错误“无法在层中分割奇数个 PSU”，当我将选项设置为small =“merge”，large =“merge”时，它运行了几个小时，然后操作系统抛出了一个错误巨大的叹息，耗尽了应用程序内存； type="JKn" 返回错误“无法分配大小为 11964693.8 Gb 的向量”
再次非常感谢您的建议。现在，我将让自己在很长一段时间内零碎地进行这些分析。如果我最终想出更好的方法，我会发布在SO上
