19

I would like to share some of my thoughts when trying to improve the model fitting time of a linear mixed effects model in R using the lme4 package.

Dataset Size: The dataset consists, approximately, of 400.000 rows and 32 columns. Unfortunately, no information can be shared about the nature of the data.

Assumptions and Checks: It is assumed that the response variable comes from a Normal distribution. Prior to the model fitting process, variables were tested for collinearity and multicollinearity using correlation tables and the alias function provided in R.

Continuous variables were scaled in order to help convergence.

Model Structure: The model equation contains 31 fixed effects (including intercept) and 30 random effects (intercept is not included). Random effects are randomized for a specific factor variable that has 2700 levels. The covariance structure is Variance Components as it is assumed that there is independency between random effects.

Model equation example:

lmer(Response ~ 1 + Var1 + Var2 + ... + Var30 + (Var1-1| Group) + (Var2-1| Group) + ... + (Var30-1| Group), data=data, REML=TRUE)

Model was fitted successfully, however, it took about 3,1 hours to provide results. The same model in SAS took a few seconds. There is available literature on the web on how to reduce time by using the non-linear optimization algorithm nloptwrap and turnining off the time consuming derivative calculation that is performed after the optmization is finished calc.derivs = FALSE:

https://cran.r-project.org/web/packages/lme4/vignettes/lmerperf.html

Time was reduced by 78%.

Question: Is there any other alternative way to reduce the model fitting time by defining the lmer parameter inputs accordingly? There is so much difference between R and SAS in terms of model fitting time.

Any suggestion is appreciated.

mammask
  • 285
  • 1
  • 3
  • 10
  • My first move would be to subsample the dataset (rows) and fit that. – Roman Luštrik Aug 24 '15 at 09:03
  • 1
    Hi @Roman Luštrik. I need to obtain coefficient estimates using the entire dataset. – mammask Aug 24 '15 at 09:16
  • I think the algorithms used in `lmer` and SAS proc mixed are different, although you'd have to get @benbolker to confirm. You could also examine the [source code](https://github.com/lme4/lme4/tree/master/src) – alexwhitworth Aug 25 '15 at 20:30
  • 1
    a few seconds vs. 0.75 hours (after the optimization you've done) is really pretty surprising. Do the SAS and `lmer` results agree with each other, i.e. are you really fitting the same model? Are all of the `Varx` variables numeric? Can you provide a reproducible example (I know your data are confidential, but all you would really need to do is post simulation code that produces *approximately* the same structure as your data ...) – Ben Bolker Aug 29 '15 at 20:40
  • 1
    Hi Ben Bolker. Results are exactly the same and variables are also numeric. I will be back in a week so I will be glad to share approximately the same structure of the data. Thanks for your time. – mammask Aug 31 '15 at 13:19
  • Dear @BenBolker I am sending you an e-mail with approximately the same structure of my initial data. Thank you for your assistance. – mammask Sep 09 '15 at 15:03
  • this can be useful: https://cran.r-project.org/web/packages/lme4/vignettes/lmerperf.html – toto_tico Nov 15 '22 at 13:13

3 Answers3

13

lmer() determines the parameter estimates by optimizing the profiled log-likehood or profiled REML criterion with respect to the parameters in the covariance matrix of the random effects. In your example there will be 31 such parameters, corresponding to the standard deviations of the random effects from each of the 31 terms. Constrained optimizations of that size take time.

It is possible that SAS PROC MIXED has specific optimization methods or has more sophisticated ways of determining starting estimates. SAS being a closed-source system means we won't know what they do.

By the way, you can write the random effects as (1+Var1+Var2+...+Var30 || Group)

Axeman
  • 32,068
  • 8
  • 81
  • 94
Douglas Bates
  • 476
  • 3
  • 5
  • 1
    There is an aspect of the evaluation of the profiled log-likelihood in lmer that can be made faster for these types of models. When we wrote lme4 we solved for the fixed-effects coefficients and the modes of the random effects at each iteration. This involves operations on matrices with n rows and vectors of length n, where n is the number of observations. Recently I realized that this is not necessary. I have a [Julia](http://julialang.org) implementation of the modified algorithm but nothing in R yet. – Douglas Bates Sep 10 '15 at 17:48
  • 2
    I'm new to stackoverflow and apparently have missed out on some of the subtleties. I usually identify myself as Doug Bates, not user1864481 – Douglas Bates Sep 10 '15 at 17:50
  • 3
    Dear Douglas. E. Bates, your answer is a surprise. I managed to fit the same model in `Julia`. Results indicate that `Julia` is approximately 74 times faster than `R` when it comes to fit the specific model. This is a huge difference in performance. Also, results are aligned. Is there any plan to perform the same development in `R lme4` in the near future? Thank you for your assistance. Konstantinos Mammas – mammask Sep 15 '15 at 14:10
  • 4
    I have two [Julia](http://JuliaLang.org) packages for fitting linear mixed-effects models at present - an older algorithm in [MixedModels](https://github.com/dmbates/MixedModels.jl) and the new algorithm in [ReTerms](https://github.com/dmbates/ReTerms.jl). The plan is to fold ReTerms into MixedModels and document the new algorithm. I don't currently plan to implement the new algorithm in [lme4](https://github.com/lme4/lme4). Some of the key parts would, I think, be difficult to implement in R/Rcpp. – Douglas Bates Sep 16 '15 at 15:24
  • 2
    This is clear and also opens a completely new task for me; whether to use `Julia` instead of `R` in order to fit Generalized Linear Models or Linear Mixed Effects Models. Thank you very much for your valuable assistance. I think that this conversation will make a lot of people reconsider which statistical language to use in order to fit mixed models. – mammask Sep 17 '15 at 10:37
  • @mammask which Julia package did you finally use (the one that was 74 times faster than `R lme4`)? – elikesprogramming Apr 02 '17 at 17:36
  • 1
    sorry the late reply. The package I used is the following: https://github.com/dmbates/MixedModels.jl (implemented in Julia programming language). It is written from @DouglasBates so I hope it helps. – mammask Apr 14 '17 at 10:29
  • Any faster option? Something like speedglm or bigglm but with random effects? – skan May 08 '17 at 20:08
  • 1
    Hi @mammask , I am in a similar mess. Can you tell me what was the final package/approach that you used? LMER takes a lot of time but I dont know Julia. Shall I learn Julia to make my model faster or is there a way? – kawsleo Aug 05 '21 at 12:51
  • Hi @kawsleo, I think using Julia will help you speed up things a little bit. – mammask Aug 17 '21 at 21:55
6

We have implemented random intercepts regression assuming compound symmetry in the R package Rfast. The command is rint.reg. It is 30+ times faster than the corresponding lme4 function. I do not know if this helps, but just in case.

https://cran.r-project.org/web/packages/Rfast/index.html

Michail
  • 177
  • 1
  • 6
5

If you use glmer rather than lmer, there is a parameter nAGQ. I found that setting nAGQ=0 dramatically reduced the time it took to fit a fairly complex model (13 fixed effects, one random effect with varying intercept and slope, 300k rows). This basically tells glmer to use a less exact form of parameter estimation for GLMMs. See ?glmer for more details, or this post.

user2739472
  • 1,401
  • 17
  • 15