我不是调查方法学家或人口统计学家,但我是Thomas Lumley的R Survey软件包的狂热粉丝。我一直在处理一个相对较大的复杂调查数据集,即医疗保健成本和利用项目(HCUP)国家急诊科样本(NEDS)。正如美国医疗保健研究和质量局所描述的那样,这是“947家位于30个州的医院的出院急诊科就诊数据,近似于美国基于医院的急诊科的20%分层样本”。
2006年至2012年的完整数据集共有198,102,435个观测值。我已将数据子集化为40,073,358个与外伤相关的出院记录和66个变量。即使是简单的调查程序也需要极长的时间才能运行这些数据。我曾尝试增加RAM (晚期2013年Mac Pro, 3.7GHz四核, 128GB(!)内存), 在可用时使用multicore, subsetting, 使用MonetDB等out-of-memory DBMS。基于设计的调查程序仍然需要数小时,有时甚至需要多小时。一些适度复杂的分析需要高达15小时。我猜想大部分的计算工作都与一个巨大的协方差矩阵有关?可以预见,使用原始数据的速度要快几个数量级。更有趣的是,根据不同的程序,在这么大的数据集上,未经调整的估计值可能非常接近于调查结果(见下面的例子)。基于设计的结果显然更精确、更可靠,但在时间成本上,数小时对比数秒钟的计算时间并不可忽视,这开始看起来像是一个非常漫长的绕路。
是否有人有这方面的经验?有没有优化R survey程序以适应大数据集的方法?或许可以更好地利用并行处理?使用INLA或Hamiltonian方法(例如Stan)的贝叶斯方法是否是一种可能的解决方案?在调查足够大且代表性足够好的情况下,一些未经调整的估计值是否可以被接受,特别是对于相对度量?
以下是一些未经调整的估计值接近调查结果的例子。
在这个第一个例子中,svymean在内存中花费了不到一小时的时间,而超出内存则需要三个多小时。直接计算只需不到一秒钟。更重要的是,点估计值(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
对于使用svyby(近两个小时)和tapply(约4秒钟)应用于数据子集的结果,svymean与mean之间存在类似的对应关系:
# 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和Anthony借给我你们的神经元。很抱歉迟迟没有跟进,试用了你们两个的建议花费了一些时间。
Severin,你的观察是正确的,Revolution Analytics/MOR版本对于某些操作更快。看起来这与CRAN R附带的BLAS(“基本线性代数子程序”)库有关。它更精确,但速度较慢。因此,我使用了我的机器上的Apple Accelerate vecLib进行了BLAS优化(虽然是专有的,但免费提供给Mac用户),可以进行多线程处理(请参见http://blog.quadrivio.com/2015/06/improved-r-performance-with-openblas.html)。这似乎缩短了一些操作时间,例如从svyby/svymean的3小时到略微超过2小时。
安东尼在使用复制权重设计方法时遇到了一些问题。使用20个复制的type="bootstrap"运行了大约39小时,然后我退出了;当我将选项设置为small="merge",large="merge"时,type="BRR"返回错误信息“无法将奇数数量的PSU分割到层次结构中”,它运行了几个小时后操作系统发出巨大的叹息并耗尽了应用程序内存;type="JKn"返回错误信息“无法分配大小为11964693.8 Gb的向量”。再次感谢您的建议。目前,我将自己限制于分段和长时间运行这些分析。如果我最终想出更好的方法,我会在SO上发布。
tables()
即可。 - Severin Pappadeuxlibrary(survey)
,这就是为什么它的作者尝试使用http://sqlsurvey.r-forge.r-project.org/而不是使用其中之一的原因。 - Anthony Damico