[R-sig-ME] Summary -- Re: Computational speed - MCMCglmm/lmer
David Atkins
datkins at u.washington.edu
Fri Jun 25 18:31:02 CEST 2010
First, thanks to Doug, Jarrod, Paul, and others who chimed in on this.
I will take a stab at a summary:
For the most part, it sounds like there are no simple solutions for
computational time. Specifically, beyond obvious differences in
hardware, there may be some minor gains from an optimized BLAS, but as
Doug notes, this is not a panacea.
Doug and Jarrod both commented on the fact that at its heart, there is a
single optimization process going on, so it's not clear how this could
be threaded or parallelized (and it would certainly not be trivial to try).
Finally, Doug has an option in the development version of (g)lmer that
allows "coarse" but faster fits. This might be quite useful for initial
model exploration.
Again, thanks all.
cheers, Dave
Dave Atkins, PhD
Research Associate Professor
Department of Psychiatry and Behavioral Science
University of Washington
datkins at u.washington.edu
Center for the Study of Health and Risk Behaviors (CSHRB)
1100 NE 45th Street, Suite 300
Seattle, WA 98105
206-616-3879
http://depts.washington.edu/cshrb/
(Mon-Wed)
Center for Healthcare Improvement, for Addictions, Mental Illness,
Medically Vulnerable Populations (CHAMMP)
325 9th Avenue, 2HH-15
Box 359911
Seattle, WA 98104
http://www.chammp.org
(Thurs)
Jarrod Hadfield wrote:
> Hi Dave,
>
> Sorry for the late response on this - I've been doing field work. I
> think there is little option for speeding up the time per iteration in
> MCMCglmm without considerable extra work. Most cpu are spent either
> solving the mixed model equations or generating random (normal)
> deviates. It may be possible to multi-thread the latter but how easy
> this would be and how much of an increase in speed you would get I'm not
> sure. For ZIP models at least, the long computing time is a result of
> poor mixing. The next version (as you know) will allow hurdle models to
> be fitted. These have much better mixing properties and in terms of
> model fit are very similar to ZIP models.
>
> Cheers,
>
> Jarrod
> On 19 Jun 2010, at 16:42, David Atkins wrote:
>
>>
>> Hi all--
>>
>> I use (g)lmer and MCMCglmm on a weekly basis, and I am wondering about
>> options for speeding up their computations. This is primarily an
>> issue with MCMCglmm, given the many necessary MCMC iterations to get
>> to convergence on some problems. But, even with glmer(), I have runs
>> that get into 20-30 minutes.
>>
>> First, let me be very clear that this is in no way a criticism of
>> Doug's and Jarrod's work (package developers for lme4 and MCMCglmm,
>> respectively). Their code has probably brought models/data into range
>> that would not have been possible.
>>
>> Second, I have included link to data and script below, along with some
>> timings on my computer: Mac Book Pro, 2.5GHz, with 4GB RAM. Here's
>> sessionInfo from my runs:
>>
>> > sessionInfo()
>> R version 2.11.1 (2010-05-31)
>> i386-apple-darwin9.8.0
>>
>> locale:
>> [1] en_US.UTF-8/en_US.UTF-8/C/C/en_US.UTF-8/en_US.UTF-8
>>
>> attached base packages:
>> [1] stats graphics grDevices utils datasets methods base
>>
>> other attached packages:
>> [1] MCMCglmm_2.05 corpcor_1.5.6 ape_2.5-3
>> [4] coda_0.13-5 tensorA_0.35 lme4_0.999375-34
>> [7] Matrix_0.999375-41 lattice_0.19-7
>>
>> loaded via a namespace (and not attached):
>> [1] gee_4.13-14 grid_2.11.1 nlme_3.1-96 stats4_2.11.1 tools_2.11.1
>>
>> Specific questions:
>>
>> 1. Would be curious to know timings on other people's set-ups. Jarrod
>> and I had an exchange one time where he was gracious enough to run a
>> zero-inflated model where I was concerned about convergence. He ran a
>> model with 1.3M iterations, which I think would take a number of days,
>> if not a week on my computer. This was part of what got me thinking
>> about this. Thus, my first interest is whether there is an "optimal"
>> hardware/OS configuration, or does it matter?
>>
>> Some things I see in R-help archives:
>>
>> 2. 32 vs. 64-bit: Seems like this is mostly an issue of data/model
>> size and whether you need to access more than 4GB of RAM. AFAICS,
>> 64-bit processors are not necessarily faster.
>>
>> 3. "Optimized" BLAS: There's a bit of discussion about optimized BLAS
>> (basis linear algebra... something). However, these discussions note
>> that there is no generally superior BLAS. Not sure whether specific
>> BLAS might be optimized for GLMM computations.
>>
>> 4. Parallel computing: With multi-core computers, looks like there are
>> some avenues for splitting intensive computations across processors.
>>
>> Finally, I'm asking here b/c I run into these issues with GLMM (and
>> zero-inflated mixed models), though most of the discussion I've seen
>> thus far about computation speed has been on R-help.
>>
>> The data below are self-reported drinks (alcohol) from college
>> students for up to the last 90 days. Distribution of counts is
>> zero-inflated. I run a Poisson GLMM with glmer, over-dispersed
>> Poisson GLMM with MCMCglmm, and then zero-inflated OD Poisson with
>> MCMCglmm and provide timings for my set-up.
>>
>> Okay, any and all thoughts welcomed.
>>
>> thanks, Dave
>>
>>
>> ### thinking through computational speed with lmer and MCMCglmm
>> #
>> ### read drinking data
>> drink.df <- read.table(file =
>> "http://depts.washington.edu/cshrb/newweb/stats%20documents/tlfb.txt",
>> header = TRUE, sep = "\t")
>> str(drink.df) # 57K rows
>>
>> ### id is id variable (shocking)
>> ### gender and weekday are 0/1 indicator variables
>> ### drinks has number of drinks consumed
>>
>> ### distribution of outcome
>> table(drink.df$drinks)
>> plot(table(drink.df$drinks), lwd=2) # zero-inflated
>>
>> ### how many people?
>> length(unique(drink.df$id)) # 990
>> sort(table(drink.df$id))
>>
>> ### NOTE: most people have max of 90, which they should
>> ### two folks with 180 and 435 (prob data errors)
>> ### long negative tail down from 90
>> #
>> ### NOTE: for this purpose, not worrying about the 180/435
>>
>> ### speed tests
>> #
>> ### fit random intercept and random slope for weekday
>> ### fixed effects for gender and weekday and interaction
>> #
>> ### Poisson GLMM with glmer()
>> library(lme4)
>>
>> system.time(
>> drk.glmer <- glmer(drinks ~ weekday*gender + (weekday | id),
>> data = drink.df, family = poisson,
>> verbose = TRUE)
>> )
>> summary(drk.glmer)
>>
>> ### timing
>> #
>> ### user system elapsed
>> ### 36.326 9.013 45.316
>>
>> ### over-dispersed Poisson GLMM with MCMCglmm()
>> library(MCMCglmm)
>>
>> prior <- list(R = list(V = 1, n = 1),
>> G = list(G1 = list(V = diag(2), n = 2)))
>> system.time(
>> drk.mcmc <- MCMCglmm(drinks ~ weekday*gender,
>> random = ~ us(1 + weekday):id,
>> data = drink.df, family = "poisson",
>> prior = prior)
>> )
>> summary(drk.mcmc) # NOTE: using summary.MCMCglmm in v 2.05 of package
>>
>> ### timing
>> #
>> ### user system elapsed
>> ### 1034.317 165.128 1203.536
>>
>> ### zero-inflated, over-dispersed Poisson GLMM with MCMCglmm()
>> #
>> ### NOTE: haven't run the following yet, other than a quick "toy run" to
>> ### sure that it is set up correctly.
>> ### NOTE: this only has random intercept in each portion of the model
>>
>> prior2 <- list(R = list(V = diag(2), n = 1, fix = 2),
>> G = list(G1 = list(V = 1, n = 1),
>> G2 = list(V = 1, n = 1)))
>> system.time(
>> drk.zimcmc <- MCMCglmm(drinks ~ -1 + trait*weekday*gender,
>> random = ~ idh(at.level(trait, 1)):id +
>> idh(at.level(trait, 2)):id,
>> rcov = ~ idh(trait):units,
>> data = drink.df, family = "zipoisson",
>> prior = prior2)
>> )
>> summary(drk.zimcmc)
>>
>> ### timing
>> #
>> ### user system elapsed
>> ### 2105.366 544.881 2640.030
>>
>> --
>> Dave Atkins, PhD
>> Research Associate Professor
>> Department of Psychiatry and Behavioral Science
>> University of Washington
>> datkins at u.washington.edu
>>
>> Center for the Study of Health and Risk Behaviors (CSHRB)
>> 1100 NE 45th Street, Suite 300
>> Seattle, WA 98105
>> 206-616-3879
>> http://depts.washington.edu/cshrb/
>> (Mon-Wed)
>>
>> Center for Healthcare Improvement, for Addictions, Mental Illness,
>> Medically Vulnerable Populations (CHAMMP)
>> 325 9th Avenue, 2HH-15
>> Box 359911
>> Seattle, WA 98104?
>> 206-897-4210
>> http://www.chammp.org
>> (Thurs)
>>
>> _______________________________________________
>> R-sig-mixed-models at r-project.org mailing list
>> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
>>
>
>
More information about the R-sig-mixed-models
mailing list