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大型複雜調查數據集的方法？

回答

相關問題