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

Posted 2023-03-11

【中文标题】R中用于大型复杂调查数据集的方法？

【英文标题】：Methods in R for large complex survey data sets?

【发布时间】：2016-05-14 15:17:19

【问题描述】：

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

从 2006 年到 2012 年的完整数据集包含 198,102,435 个观测值。我已将数据子集化为 40,073,358 次与外伤相关的出院，其中包含 66 个变量。即使对这些数据进行简单的调查程序也需要非常长的时间。我试过在它上面扔 RAM（2013 年末 Mac Pro，3.7GHz 四核，128GB（！）内存），使用multicore（如果可用），subsetting，使用out-of-memory DBMS，如MonetDB。基于设计的调查程序仍然需要数小时。有时好几个小时。一些适度复杂的分析需要 15 个小时以上。我猜大部分计算工作都与必须是巨大的协方差矩阵有关？

正如人们所预料的那样，处理原始数据的速度要快几个数量级。更有趣的是，根据程序的不同，对于这么大的数据集，未经调整的估计值可能与调查结果非常接近。 （见下面的例子）基于设计的结果显然更精确和更受欢迎，但是几个小时的计算时间与几秒钟的时间相比，对于增加的精度来说是一个不可忽视的成本。它开始看起来像是绕着街区走了很长一段路。

有没有人有这方面的经验？有没有办法优化大型数据集的 R 调查程序？也许更好地利用并行处理？使用像 Stan 这样的 INLA 或 Hamiltonian 方法的贝叶斯方法是一种可能的解决方案吗？或者，当调查规模大且具有足够代表性时，是否可以接受一些未经调整的估计值，尤其是相对测量值？

以下是一些近似调查结果的未调整估计值示例。

在第一个示例中，内存中的 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 和 Anthony 允许我借用你们的突触。很抱歉延迟跟进，没有时间尝试您的两个建议。

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

Anthony，对复制权重设计方法的运气不太好。 type="bootstrap" with replicates=20 在我退出之前运行了大约 39 小时； type="BRR" 返回错误 "Can't split with odd number of PSUs in a stratum", 当我将选项设置为 small="merge", large="merge" 时，它运行了几个小时，然后操作系统启动了叹了口气，应用程序内存不足； type="JKn" 返回错误“无法分配大小为 11964693.8 Gb 的向量”

再次感谢您的建议。现在，我将放弃自己，在很长一段时间内零碎地运行这些分析。如果我最终想出更好的方法，我会在 SO 上发帖

【问题讨论】：

你试过stat.yale.edu/~mjk56/temp/bigmemory-vignette.pdf吗？

感谢您的回复，Severin。是的。试过bigmemory、ff、data.table。

你打过交换吗？内存中对象的大小是多少？您提供了时间，但没有提供内存信息。对于 data.tables 只需输入 tables()

嗨@SeverinPappadeux 没有一个与library(survey) 兼容，这就是为什么它的作者尝试使用sqlsurvey.r-forge.r-project.org 而不是使用其中之一

我们正在考虑在数据库中进行调查分析，仅对非常大的数据集 (github.com/chrk623/svydb) 进行实施工作才值得，因此了解此类示例很有用

【参考方案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 计算的更大流程在一夜之间运行

【讨论】：

谢谢，安东尼。这是一个很好的建议。不知道这种可能的方法。会试试的。

像svydesign( id = ~ 1 , data = yourdata , weights = ~ yourweights ) 这样不正确的设计仍然可以相对较快地给出正确的系数，并且还可以让您测试您的代码以备不时之需

【参考方案2】：

已经有一段时间了，但关闭了这个循环。正如 Lumley 博士在上面最近的评论中提到的那样，Charco Hui 将实验性 sqlsurvey 包重新命名为“svydb”，我发现它是在 R 中处理非常大的调查数据集的好工具。请参阅此处的相关帖子：How to get svydb R package for large survey data sets to return standard errors