Методы в R для больших комплексных наборов данных опроса?

Question

Я не методолог опроса или демограф, но ярый поклонник опроса Р. Томаса Ламли. Я работаю с относительно большим набором комплексных данных обследования, Образцом национального плана неотложной помощи в здравоохранении (HCUP) (NEDS). Как описано Агентством по исследованиям и качеству в здравоохранении, «данные о выбросах для посещений ED из 947 больниц, расположенных в 30 штатах, приближаются к 20-процентной стратифицированной выборке американских ЭД в больницах»

Полный набор данных за 2006 год до 2012 года состоит из 198 102 425 наблюдений. Я подвел данные к 40 073 358 травматическим повреждениям, связанным с травмой, с 66 переменными. Выполнение даже простых процедур обследования по этим данным занимает чрезмерно длительное время. Я пробовал бросать в него RAM (в конце 2013 года Mac Pro, 3,7 ГГц Quad Core, память 128 ГБ (!)), Используя multicore, когда есть, subsetting, работающий с out-of-memory DBMS, как MonetDB. Процедуры обследования на основе дизайна все еще занимают часы. Иногда много часов. Некоторые скромно сложные анализы занимают более 15 часов. Я предполагаю, что большая часть вычислительных усилий связана с тем, что должно быть сложной ковариационной матрицей?

Как и следовало ожидать, работа с необработанными данными на несколько порядков выше. Более интересно, в зависимости от процедуры, с набором данных это большое, нескорректированные оценки могут быть весьма близки к результатам опроса. (См. Примеры ниже). Результаты, основанные на дизайне, явно более точны и предпочтительны, но несколько часов вычислительного времени в секундах - это не несущественная стоимость для этой добавленной точности. Он начинает напоминать очень длинную прогулку вокруг квартала.

Есть ли у кого есть опыт работы с этим? Существуют ли способы оптимизации процедур опроса R для больших наборов данных? Возможно, лучше использовать параллельную обработку? Используются ли байесовские подходы с использованием методов INLA или Hamiltonian, таких как Stan, возможное решение? Или, какие-то нескорректированные оценки, особенно для относительных мер, приемлемые, когда опрос является большим и достаточно представительным?

Вот несколько примеров нескорректированных оценок, приближающих результаты опроса.

В этом первом примере svymean в памяти занял чуть меньше часа, из памяти потребовалось более 3 часов. Прямой расчет занял второе место. Что еще более важно, точные оценки (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 

Соответствия между обзор и нескорректированные результаты начинают разрушаться с абсолютным отсчетом, что требует написания небольшой функции, которая обращается к объекту съемки, и использует небольшой код кода доктора Ламли для взвешивания отсчетов:

# 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 секунды - это значительная стоимость для более точных результатов.

обновление: 10 февраля 2016

Благодаря Северин и Энтони за предоставленную мне возможность заимствовать свои синапсы. Извините за задержку в последующем, у вас мало времени, чтобы опробовать оба ваших предложения.

Северин, вы правы в своих наблюдениях, что Revolution Analytics/MOR build быстрее для некоторых операций. Похоже, что это связано с библиотекой BLAS («Basic Linear Algebra Subprograms»), поставляемой с CRAN R. Это более точно, но медленнее. Таким образом, я оптимизировал BLAS на моей обработке с помощью собственной (но свободной от macs) Apple Accelerate vecLib, которая позволяет многопоточность (см. http://blog.quadrivio.com/2015/06/improved-r-performance-with-openblas.html). Это, казалось, сбривало некоторое время после операций, например. от 3 часов для свивы/свимэна чуть более 2 часов.

Энтони, которому не повезло с реплицируемым дизайном веса. type = "bootstrap" с репликами = 20 пробежал около 39 часов, прежде чем я ушел; type = «BRR» возвратила ошибку «Невозможно разделить с нечетными номерами блоков питания в страте», когда я установил параметры 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 выполнялись за одну ночь