Zwiększanie szybkości zapytań przy użyciu złożonego projektu ankiety w R

Question

Mam duży zestaw danych (ponad 20 milionów obs.), Który analizuję z pakietem survey, a przeprowadzenie wyszukiwania prostych zapytań zajmuje mi wiele wieków. Próbowałem znaleźć sposób na przyspieszenie mojego kodu, ale chciałbym wiedzieć, czy istnieją lepsze sposoby na zwiększenie wydajności.

W moim teście, ja porównać szybkość trzech poleceniach svyby/svytotal:

1. proste polecenie svyby/svytotal
2. obliczeń równoległych z foreachdopar użyciu 7 rdzeni
3. Skompilowany wersją opcji 2

Spoiler: Opcja 3 jest ponad dwukrotnie szybsza niż pierwsza opcja, ALE nie jest odpowiednia dla dużych zestawów danych, ponieważ opiera się na obliczeniach równoległych, które szybko osiągają limity pamięci w przypadku dużych zbiorów danych. Ja też mam ten problem, mimo 16GB pamięci RAM. Jest kilka solutions to this memory limitation, ale żadna z nich nie ma zastosowania do obiektów projektu ankiety.

Jakieś pomysły na to, jak przyspieszyć i nie zawieść z powodu ograniczeń pamięci?

My kod z powtarzalnym np

# Load Packages 
library(survey) 
library(data.table) 
library(compiler) 
library(foreach) 
library(doParallel) 
options(digits=3) 

# Load Data 
data(api) 

# Convert data to data.table format (mostly to increase speed of the process) 
apiclus1 <- as.data.table(apiclus1) 

# Multiplicate data observations by 1000 
apiclus1 <- apiclus1[rep(seq_len(nrow(apiclus1)), 1000), ] 

# create a count variable 
apiclus1[, Vcount := 1] 

# create survey design 
dclus1 <- svydesign(id=~dnum, weights=~pw, data=apiclus1, fpc=~fpc) 

1) prosty kod

t1 <- Sys.time() 
table1 <- svyby(~Vcount, 
       ~stype+dnum+cname, 
       design = dclus1, 
       svytotal) 
T1 <- Sys.time() - t1 

2) Obliczenia równoległe z foreach dopar stosując 7 rdzeni

# in this option, I create a list with different subsets of the survey design 
# that will be passed to different CPU cores to work at the same time 

subdesign <- function(i){ subset(dclus1, dnum==i)} 
groups <- unique(apiclus1$dnum) 
list_subsets <- lapply(groups[], subdesign) # apply function and get all  subsets in a list 
i <- NULL 

# Start Parallel 
registerDoParallel(cores=7) 

t2 <- Sys.time() 
table2 <- foreach (i = list_subsets, .combine= rbind, .packages="survey")  %dopar% { 
    options(survey.lonely.psu = "remove") 
    svyby(~Vcount, 
     ~stype+dnum+cname, 
     design = i, 
     svytotal)} 
T2 <- Sys.time() - t2 

3. Skompilowana wersja opcji 2

# make a function of the previous query 
query2 <- function (list_subsets) { foreach (i = list_subsets, .combine=  rbind, .packages="survey") %dopar% { 
    svyby(~Vcount, 
     ~stype+dnum+cname, 
     design = i, 
     svytotal)}} 

# Compile the function to increase speed 
query3 <- cmpfun(query2) 

t3 <- Sys.time() 
table3 <- query3 (list_subsets) 
T3 <- Sys.time() - t3 

Wyniki

>T1: 1.9 secs 
>T2: 1.13 secs 
>T3 0.58 secs 

barplot(c(T1, T2, T3), 
     names.arg = c("1) simple table", "2) parallel", "3) compiled parallel"), 
     ylab="Seconds") 

Answer 1

dzięki do układania na to pytanie tak dobrze. wydajna praca z dużymi zestawami danych badania w R prawdopodobnie wymaga podstawowej składni SQL (która jest o wiele łatwiejsza do nauczenia niż R). MonetDB jest jedyną opcją dużych zbiorów danych zgodną z pakietem survey, więc eksplorowanie innych pakietów o wysokiej wydajności (prawdopodobnie) nie przyniesie rezultatów. na ogół podczas eksploracji ogromnego zbioru danych piszę bezpośrednio w zapytaniach SQL, zamiast korzystać z pakietu ankiet, ponieważ standardowe obliczenia błędów są obciążające obliczeniowo (a wariancje nie są tak użyteczne podczas interaktywnego eksplorowania danych). zauważ, jak końcowy znacznik czasu SQL usuwa wszystkie inne opcje. obliczyć szybki średnią ważoną, użyć czegoś podobnego "SELECT by_column , SUM(your_column * the_weight)/SUM(the_weight) FROM yourdata GROUP BY by_column"

kiedy potrzebują standardowych błędów interaktywnie, linearyzacji (svydesign) jest często bardziej obliczeniowo niż replikacji (svrepdesign), ale czasami tworzenia wzorów replikacji (co mam zrobić z poniższym jk1w_dclus1) wymaga większej znajomości metod ankiet niż niektórzy użytkownicy są w stanie z nimi wygodnie pracować.

# Load Packages 
library(MonetDB.R) 
library(MonetDBLite) 
library(DBI) # suggested in comments and needed on OSX 
library(survey) 

# Load Data 
data(api) 

# Multiplicate data observations by 10000 
apiclus1 <- apiclus1[rep(seq_len(nrow(apiclus1)), 10000), ] 

# create a count variable 
apiclus1$vcount <- 1 

# create survey design 
dclus1 <- svydesign(id=~dnum, weights=~pw, data=apiclus1, fpc=~fpc) 


dbfolder <- tempdir() 

db <- dbConnect(MonetDBLite() , dbfolder) 
dbWriteTable(db , 'apiclus1' , apiclus1) 


db_dclus1 <- 
    svydesign(
     weight = ~pw , 
     id = ~dnum , 
     data = "apiclus1" , 
     dbtype = "MonetDBLite" , 
     dbname = dbfolder , 
     fpc = ~fpc 
    ) 

# you provided a design without strata, 
# so type="JK1" matches that most closely. 
# but see survey:::as.svrepdesign for other linearization-to-replication options 
jk1w <- jk1weights(psu = apiclus1$dnum , fpc = apiclus1$fpc) 

# after the replicate-weights have been constructed, 
# here's the `svrepdesign` call.. 
jk1w_dclus1 <- 
    svrepdesign(
     weight = ~pw , 
     type = "JK1" , 
     repweights = jk1w$repweights , 
     combined.weights = FALSE , 
     scale = jk1w$scale , 
     rscales = jk1w$rscales , 
     data = 'apiclus1' , 
     dbtype = "MonetDBLite" , 
     dbname = dbfolder 
    ) 

# slow 
system.time(res1 <- svyby(~vcount,~stype+dnum+cname,design = dclus1,svytotal)) 
# > system.time(res1 <- svyby(~vcount,~stype+dnum+cname,design = dclus1,svytotal)) 
    # user system elapsed 
    # 17.40 2.86 20.27 


# faster 
system.time(res2 <- svyby(~vcount,~stype+dnum+cname,design = db_dclus1,svytotal)) 
# > system.time(res2 <- svyby(~vcount,~stype+dnum+cname,design = db_dclus1,svytotal)) 
    # user system elapsed 
    # 13.00 1.20 14.18 


# fastest 
system.time(res3 <- svyby(~vcount,~stype+dnum+cname,design = jk1w_dclus1,svytotal)) 
# > system.time(res3 <- svyby(~vcount,~stype+dnum+cname,design = jk1w_dclus1,svytotal)) 
    # user system elapsed 
    # 10.75 1.19 11.96 

# same standard errors across the board 
all.equal(SE(res1) , SE(res2)) 
all.equal(SE(res2) , SE(res3)) 
# NOTE: the replicate-weighted design will be slightly different 
# for certain designs. however this technique is defensible 
# and gets used in 
# https://github.com/ajdamico/asdfree/tree/master/Censo%20Demografico 


# at the point you do not care about standard errors, 
# learn some sql: 
system.time(res4 <- dbGetQuery(db , "SELECT stype , dnum , cname , SUM(pw) FROM apiclus1 GROUP BY stype , dnum , cname")) 
# because this is near-instantaneous, no matter how much data you have. 

# same numbers as res1: 
all.equal(as.numeric(sort(coef(res1))) , sort(res4$L1)) 
# > system.time(res4 <- dbGetQuery(db , "SELECT stype , dnum , cname , SUM(pw) FROM apiclus1 GROUP BY stype , dnum , cname")) 
    # user system elapsed 
    # 0.15 0.20 0.23