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

Question

我不是调查方法学家或人口统计学家，但我是 Thomas Lumley 的 R 调查包的狂热粉丝。我一直在处理一个相对较大的复杂调查数据集，即医疗保健成本和利用项目 (HCUP) 国家急诊部样本 (NEDS)。正如医疗保健研究和质量机构所描述的那样，它是“来自 30 个州的 947 家医院的急诊就诊的出院数据，接近美国医院急诊的 20% 分层样本”

2006 年至 2012 年的完整数据集包含 198,102,435 个观测值。我已将数据分割为 40,073,358 次与外伤相关的出院，其中包含 66 个变量。对这些数据运行即使是简单的调查程序也需要非常长的时间。我曾尝试使用 multicore 将 RAM(2013 年末 Mac Pro，3.7GHz 四核，128GB(!)内存)投入其中如果可用，subsetting ，使用 out-of-memory DBMS喜欢 MonetDB .基于设计的调查程序仍然需要数小时。有时好几个小时。一些适度复杂的分析需要 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 

使用 svyby(近 2 小时)与 tapply(4 秒左右)应用于数据子集的 svymean 与均值的结果之间存在类似的对应关系:

# 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 小时 vs. 1 秒是一个可观的成本。

更新:2016 年 2 月 10 日

感谢 Severin 和 Anthony 允许我借用你们的突触。很抱歉延迟跟进，几乎没有时间尝试您的两个建议。

Severin ，您的观察是正确的，对于某些操作，Revolution Analytics/MOR 构建速度更快。看起来它与 CRAN R 附带的 BLAS(“基本线性代数子程序”)库有关。它更精确，但速度更慢。因此，我使用允许多线程处理的专有(但在 Mac 上免费)Apple Accelerate vecLib(请参阅 http://blog.quadrivio.com/2015/06/improved-r-performance-with-openblas.html)在我的机器上优化了 BLAS。这似乎减少了一些操作时间，例如从 svyby/svymean 的 3 小时到 2 小时多一点。

Anthony 对复制权重设计方法的运气较差。 type="bootstrap"with replications=20 在我退出之前运行了大约 39 小时； type="BRR"返回错误“不能在层中用奇数个 PSU 拆分”，当我将选项设置为 small="merge", large="merge"时，它运行了几个小时，然后操作系统才开始巨大的叹息并耗尽了应用程序内存； type="JKn"返回错误“无法分配大小为 11964693.8 Gb 的向量”

再次，非常感谢您的建议。现在，我将让自己在很长一段时间内逐步进行这些分析。如果我最终想出更好的方法，我会在 SO 上发布

Answer 1

对于庞大的数据集，线性化设计 ( 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 计算的较大流程在夜间运行