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

我不是一名调查方法学家或人口统计学家，但却是Thomas Lumley的R调查报告的热心粉丝。我一直在处理一个相对较大的复杂调查数据集，即医疗成本和利用项目（HCUP）国家急诊部样本（NEDS）。如医疗保健研究和质量机构所述，它是“来自位于30个州的947家医院的急诊就诊的出院数据，近似于美国医院急诊科的20％分层样本”R大型复杂调查数据集的方法？

2006年的完整数据集到2012年包括198,102,435个观察。我已经将这些数据分成了40,073,358个与66个变量相关的创伤性损伤相关放电。即使对这些数据运行简单的调查程序也需要非常长的时间。我已经尝试使用multicore（），subsetting，与out-of-memory DBMS（如MonetDB）一起使用RAM（2013年末Mac Pro，3.7GHz四核，128GB（！）内存）。基于设计的调查程序仍需数小时。有时候很多小时。一些适度复杂的分析需要15个小时以上。我猜测大部分的计算工作都与必须是一个巨大的协方差矩阵有关。

正如人们所预料的那样，处理原始数据的速度要快几个数量级。更有意思的是，取决于程序，对于这样大的数据集，未经调整的估计值可能非常接近调查结果。（请参阅下面的示例）基于设计的结果显然更加精确并且更受欢迎，但计算时间与秒数的计算时间对于提高精度而言并不是不小的成本。它开始看起来像在街区周围漫长的漫步。

有没有人有过这方面的经验？是否有方法来优化大数据集的R调查程序？也许更好地使用并行处理？贝叶斯方法是否使用INLA或Hamiltonian方法（如Stan）作为可能的解决方案？或者，是否有一些未经调整的估计值，特别是相对测量值，当调查数据很大且具有代表性时，可以接受？

以下是几个未经调整的估计值的示例，用于近似调查结果。

在这第一个例子中，内存中的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 VS应用于使用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

之间的对应关系调查和未调整的结果开始打破绝对计数，这需要写一个小函数，呼吁调查对象，并使用一些博士拉姆利的代码来衡量计数：

# 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，您的观察结果显示Revolution Analytics/MOR build对某些操作而言速度更快。看起来它与CRAN R附带的BLAS（“基本线性代数子程序”）库有关。它更精确，但更慢。因此，我使用专有（但免费使用Mac）Apple Accelerate vecLib进行多线程优化（请参阅http://blog.quadrivio.com/2015/06/improved-r-performance-with-openblas.html），以优化BLAS。这似乎会缩短操作的时间，例如从3小时的svyby/svymean到2个多小时。

安东尼，与复制重量设计方法运气不多。在退出之前，复制= 20的type =“bootstrap”运行了大约39个小时; type =“BRR”返回错误“不能用分层中的奇数个PSU分割”，当我将选项设置为small =“merge”，large =“merge”时，它在操作系统启动之前运行了几个小时巨大的叹息并且用尽了应用程序内存; type =“JKn”returned he error“can not allocate vector of size of 11964693.8 Gb”

再次，非常感谢您的建议。我现在要辞职，让我自己长时间零碎地进行这些分析。如果我最终想出了一个更好的方法，那么对于庞大的数据集，我将在SO

来源

2016-02-04 charlie

您尝试过http://www.stat.yale.edu/~mjk56/temp/bigmemory- vignette.pdf？ –

感谢Severin的回应。是。已经尝试过bigmemory，ff，data.table。 – charlie

你碰到了交换吗？什么是内存中的对象大小？你提供了时间，但没有记忆信息。对于data.tables只需键入'tables（）' –

上发布，线性化设计（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（尽可能多的可以容忍）。以交互方式运行您的系数并细化您在一天中需要的数字，然后保留需要SE计算的较大过程，并在夜间运行

来源

2016-02-06 12:00:38

感谢，安东尼链接。这是一个很好的建议。没有意识到这种可能的方法。会尝试。 – charlie

不正确的设计，比如'svydesign（ID =〜1，数据= yourdata，权重=〜yourweights）'还是应该给予正确的系数比较迅速，也让你测试你的代码的运行过夜 –

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

回答

相关问题