lmer2 slower?
I am forwarding (with permission) a conversation with Bill Auty who discovered that in one of his examples lmer2 is considerably slower than lmer, although they get to the same estimates. This can happen. As I mentioned in one reply I hope that forcing a supernodal Cholesky decomposition for models with crossed or partially crossed grouping factors for the random effects will help but I haven't had time to explore this yet. I'll post some other timing results in another message. ---------- Forwarded message ---------- From: Douglas Bates <bates at stat.wisc.edu> Date: Jan 26, 2007 3:38 PM Subject: Re: lmer2 slower? To: Bill Auty <bill at edmeasure.com>
On 1/26/07, Bill Auty <bill at edmeasure.com> wrote:
> class(xq2 at L)
[1] "dCHMsimpl" attr(,"package") [1] "Matrix"
Ah, good. That 's the answer that I wanted. Briefly, the CHOLMOD sparse matrix code allows for two types of Cholesky decompositions of sparse, symmetric matrices. These are called "simplicial" and "supernodal". In the lmer code I forced the supernodal decomposition. In the lmer2 code I allow the sparse matrix manipulation code to choose one or the other based on some criteria which I can "tune". I think that the default values that are currently being used are too conservative about switching to the supernodal decomposition so I may be able to gain back some of the loss of speed. Thanks again for checking this. May I send a copy of this reply to the R-SIG-mixed-models list?
>
I don't have the math chops to stay with you on the theory of this, but I'm glad to help by trying things on my data. Douglas Bates wrote:
Thanks. That's valuable information. Could you check something for me? If you still have the object fit by lmer2 could you send me the result of class(xq2 at L) The response should be either "dCHMsimpl" or "dCHMsuper" On 1/26/07, Bill Auty <bill at edmeasure.com> wrote:
I've attached the output from running lmer and lmer2 on the same model. The variables are: rit = scale score on a math test gr = grade (range 0 - 2) stu = student ID sch = school ID dist = district ID The students are partially crossed in school and district.
system.time(xq2<-lmer2(rit~gr+I(gr^2)+(gr|stu)+(gr|sch)+(gr|dist),e2math, + control = list(niterEM =0, gradient = FALSE),method="ML")) [1] 428.223 30.569 458.986 0.000 0.000
xq2
Linear mixed-effects model fit by maximum likelihood
Formula: rit ~ gr + I(gr^2) + (gr | stu) + (gr | sch) + (gr | dist)
Data: e2math
AIC BIC logLik MLdeviance REMLdeviance
977133 977251 -488554 977109 977117
Random effects:
Groups Name Variance Std.Dev. Corr
stu (Intercept) 78.60373 8.86587
gr 1.88712 1.37373 -0.227
sch (Intercept) 7.24479 2.69161
gr 0.75162 0.86696 -0.512
dist (Intercept) 4.30657 2.07523
gr 0.37241 0.61026 0.139
Residual 26.61534 5.15901
Number of obs: 136484, groups: stu, 80224; sch, 182; dist, 63
Fixed effects:
Estimate Std. Error t value
(Intercept) 209.65220 0.38078 550.6
gr 7.67870 0.13085 58.7
I(gr^2) -0.76376 0.01931 -39.6
Correlation of Fixed Effects:
(Intr) gr
gr -0.188
I(gr^2) 0.049 -0.411
system.time(q2<-lmer(rit~gr+I(gr^2)+(gr|stu)+(gr|sch)+(gr|dist),e2math, + control = list(niterEM =0, gradient = FALSE),method="ML")) [1] 241.323 24.070 265.412 0.000 0.000
q2
Linear mixed-effects model fit by maximum likelihood
Formula: rit ~ gr + I(gr^2) + (gr | stu) + (gr | sch) + (gr | dist)
Data: e2math
AIC BIC logLik MLdeviance REMLdeviance
977133 977251 -488554 977109 977117
Random effects:
Groups Name Variance Std.Dev. Corr
stu (Intercept) 78.60255 8.86581
gr 1.88767 1.37393 -0.227
sch (Intercept) 7.24380 2.69143
gr 0.75110 0.86666 -0.512
dist (Intercept) 4.30710 2.07536
gr 0.37265 0.61045 0.139
Residual 26.61352 5.15883
number of obs: 136484, groups: stu, 80224; sch, 182; dist, 63
Fixed effects:
Estimate Std. Error t value
(Intercept) 209.65218 0.38079 550.6
gr 7.67867 0.13086 58.7
I(gr^2) -0.76376 0.01931 -39.6
Correlation of Fixed Effects:
(Intr) gr
gr -0.188
I(gr^2) 0.049 -0.411