logo标签列表

R中用于大型复杂调查数据集的方法？

rsurveylarge-data
9

我不是调查方法学家或人口统计学家，但我是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上发布。
1
你试过访问http://www.stat.yale.edu/~mjk56/temp/bigmemory-vignette.pdf吗？ - Severin Pappadeux
感谢您的回复，Severin。是的，我已经尝试过bigmemory、ff和data.table。 - charlie
你是否执行了交换操作？内存中的对象大小是多少？你提供了时间信息，但没有内存信息。对于data.tables，只需键入 tables() 即可。 - Severin Pappadeux
1
嗨@SeverinPappadeux，这些都不兼容library(survey)，这就是为什么它的作者尝试使用http://sqlsurvey.r-forge.r-project.org/而不是使用其中之一的原因。 - Anthony Damico
2
我们正在考虑在数据库中进行调查分析，只有对于非常大的数据集才值得实施（https://github.com/chrk623/svydb），因此了解这样的示例是有用的。 - Thomas Lumley
显示剩余2条评论
2个回答
6
对于大数据集，线性化设计（svydesign）比复制设计（svrepdesign）慢得多。请查看survey::as.svrepdesign中的加权函数，并使用其中一个来直接生成复制设计。您不能使用线性化来完成此任务。最好不要使用as.svrepdesign，而是使用其中的函数。
例如，直接将cluster=、strata=和fpc=应用到复制加权设计中，请参见https://github.com/ajdamico/asdfree/blob/master/Censo%20Demografico/download%20and%20import.R#L405-L429 请注意，您还可以在此处查看每个事件的时间戳的分分秒秒速度测试：http://monetdb.cwi.nl/testweb/web/eanthony/ 另外，请注意replicates=参数几乎占设计运行速度的100%。因此，可能需要创建两种设计，一种用于系数（只需几个副本），另一种用于标准误差（尽可能多）。在白天运行系数并优化您需要的数字，然后将需要计算SE的更大进程留在夜间运行。
谢谢，安东尼。这是一个很好的建议。我之前不知道这种可能的方法。我会试一下。 - charlie
1
一个不正确的设计，比如 svydesign( id = ~ 1 , data = yourdata , weights = ~ yourweights )，仍然可以相对快速地给出正确的系数，并且让您测试代码是否能够在过夜运行中正常工作。 - Anthony Damico
2
网页内容由stack overflow 提供, 点击上面的
可以查看英文原文，
原文链接