I want to run a mixed model (using lme4::lmer) on 60M observations of the following format; all predictor/dependent variables are categorical (factors) apart from the continuous dependent variable tc; patient is the grouping variable for a random intercept term. I have 64-bit R and 16Gb RAM and I'm working from a central server. RStudio is the most recent server version.
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
I'm getting a cannot allocate a vector of 25.5Gb error. I'm assigned 40Gb on the server and am using 25 so I guess that means I need another 10 or so. I don't think I can get any extra space assigned.
I don't know the first thing about parallel processing except that I'm using one of four cores at the moment. Can anyone suggest parallel code for this model, or perhaps a different fix?
As pointed out by Carl Witthoft, the standard parallelization tools in R work by making a copy of the data for every worker (I originally called this a shared memory model but I think that's not accurate), so running them all on a system with a single block of memory (as opposed, e.g., to a "simple network of workstations" (SNOW) where each process lives on a different machine) they will make things worse rather than better (their main purpose is to accelerate compute-bound jobs by using multiple processors).
In the short term, you might be able to save some memory by treating the categorical fixed-effect predictors (age, atc) as random effects but forcing their variances to be large. I don't know if this will be enough to save you or not; it will compress the fixed-effect model matrix a lot, but the model frame will still be stored/replicated with the model object ...
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)
Random effects are automatically sorted in order from largest
to smallest number of levels. Following the machinery described
in the ?modular help page:
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
You can ignore the reported variances for the categorical terms, and use
ranef() to recover their (unshrunk) estimates.
In the long term, the proper way to do this problem is probably to parallelize it with a distributed-memory model. In other words, you would want to parcel the data out in chunks to different servers; use the machinery described in ?modular to set up a likelihood function (actually a REML-criterion function) that gives the REML criterion for a subset of the data as a function of the parameters; then run a central optimizer that takes a set of parameters and evaluates the REML criterion by submitting the parameters to each server, retrieving the values from each server, and adding them. The only two problems I see with implementing this are (1) I don't actually know how to implement distributed-memory computation in R (based on this intro document it seems that the Rmpi/doMPI packages might be the right way to go); (2) in the default way that lmer is implemented, the fixed-effects parameters are profiled out rather than being explicitly part of the parameter vector.