montaggio di un modello lineare misto su un set di dati molto grande

Question

Voglio eseguire un modello misto (utilizzando lme4::lmer) su 60M osservazioni del seguente formato; tutte le variabili predittore/dipendenti sono categoriali (fattori) a parte la variabile dipendente continua tc; patient è la variabile di raggruppamento per un termine di intercettazione casuale. Ho 64-bit R e 16 GB di RAM e sto lavorando da un server centrale. RStudio è la versione server più recente.

model <- lmer(tc~sex+age+lho+atc+(1|patient), 
       data=master,REML=TRUE) 

lho sex tc  age   atc patient 
18 M 16.61 45-54  H 628143 
7 F 10.52 12-15  G 2013855 
30 M 92.73 35-44  N 2657693 
19 M 24.92 70-74  G 2420965 
12 F 17.44 65-69  A 2833610 
31 F 7.03 75 and over A 1090322 
3 F 28.59 70-74  A 2718649 
29 F 4.09 75 and over C 384578 
16 F 67.22 65-69  R 1579355 
23 F 7.7  70-74  C 896374 

Ho riscontrato un errore cannot allocate a vector of 25.5Gb. Mi è stato assegnato 40 GB sul server e sto usando 25 quindi suppongo significhi che ho bisogno di altri 10 o giù di lì. Non penso di poter ottenere uno spazio extra assegnato.

Non conosco la prima cosa sull'elaborazione parallela, tranne che sto utilizzando uno dei quattro core al momento. Qualcuno può suggerire il codice parallelo per questo modello, o forse una correzione diversa?

Answer 1

Come rilevato da Carl Witthöft, gli strumenti di parallelizzazione standard in R utilizzano un modello di memoriacondivisa, in modo che si peggiorare le cose invece di migliorare (il loro scopo principale è quello di accelerare compute-bound lavori, usando multipla processori).

Nel breve termine, si potrebbe essere in grado di risparmiare un po 'di memoria trattando le categoriche predittori a effetti fissi (age, atc) come effetti casuali, ma costringendo loro varianze ad essere grandi. Non so se questo sarà sufficiente per salvarti o no; comprimerà la matrice a effetti fissi modello molto, ma la cornice modello verrà comunque memorizzato/replica con il modello oggetto ...

dd1 <- read.table(header=TRUE, 
text="lho sex tc  age   atc patient 
18 M 16.61 45-54  H 628143 
7 F 10.52 12-15  G 2013855 
30 M 92.73 35-44  N 2657693 
19 M 24.92 70-74  G 2420965 
12 F 17.44 65-69  A 2833610 
31 F 7.03 75_and_over A 1090322 
3 F 28.59 70-74  A 2718649 
29 F 4.09 75_and_over C 384578 
16 F 67.22 65-69  R 1579355 
23 F 7.7  70-74  C 896374") 
n <- 1e5 
set.seed(101) 
dd2 <- with(dd1, 
     data.frame(tc=rnorm(n,mean=mean(tc),sd=sd(tc)), 
       lho=round(runif(n,min=min(lho),max=max(lho))), 
       sex=sample(levels(sex),size=n,replace=TRUE), 
       age=sample(levels(age),size=n,replace=TRUE), 
       atc=sample(levels(atc),size=n,replace=TRUE), 
       patient=sample(1:1000,size=n,replace=TRUE))) 
library("lme4") 
m1 <- lmer(tc~sex+(1|lho)+(1|age)+(1|atc)+(1|patient), 
      data=dd2,REML=TRUE) 

effetti casuali sono ordinati automaticamente in ordine dal più grande al più piccolo numero di livelli. In seguito la macchina descritta nella pagina ?modular aiuto:

lmod <- lFormula(tc~sex+(1|lho)+(1|age)+(1|atc)+(1|patient), 
        data=dd2,REML=TRUE) 
names(lmod$reTrms$cnms) ## ordering 
devfun <- do.call(mkLmerDevfun, lmod) 
wrapfun <- function(tt,bigsd=1000) { 
    devfun(c(tt,rep(bigsd,3))) 
} 
wrapfun(1) 
opt <- optim(fn=wrapfun,par=1,method="Brent",lower=0,upper=1000) 
opt$fval <- opt$value ## rename/copy 
res <- mkMerMod(environment(devfun), opt, lmod$reTrms, fr=lmod$fr) 
res 

È possibile ignorare gli scostamenti rilevati per i termini categorici, e utilizzare ranef() per recuperare le loro stime (unshrunk).

A lungo termine, il modo corretto per risolvere questo problema è probabilmente quello di parallelizzarlo con un modello di memoria distribuita. In altre parole, vorrai suddividere i dati in blocchi su diversi server; utilizzare la macchina descritta in ?modular per impostare una funzione di verosimiglianza (in realtà una funzione del criterio REML) che fornisce il criterio REML per un sottoinsieme dei dati in funzione dei parametri; quindi eseguire un ottimizzatore centrale che accetta un set di parametri e valuta il criterio REML inviando i parametri a ciascun server, recuperando i valori da ciascun server e aggiungendoli. Gli unici due problemi che vedo con l'implementazione di questo sono (1) in realtà non so come implementare il calcolo della memoria distribuita in R (basato su this intro document sembra che i pacchetti Rmpi/doMPI potrebbero essere la strada giusta da percorrere); (2) nel modo predefinito in cui è implementato lmer, i parametri a effetti fissi vengono profilati anziché essere esplicitamente parte del vettore dei parametri.