[R-sig-ME] lmer() vs. lme() gave different variance component estimates
Jarrod Hadfield
j.hadfield at ed.ac.uk
Sat Sep 18 12:09:04 CEST 2010
Both converge to the same answer on Fedora and Mac (Mac sessionInfo() below)
lmer(score~trt+(1|id/eye),dat, verbose=T, start=list(matrix(10),matrix(10)))
0: 429.50725: 10.0000 10.0000
1: 406.76151: 12.7962 11.2523
2: 403.69482: 13.3016 12.8538
3: 403.18712: 13.6918 13.5236
4: 403.15738: 13.8230 13.6823
5: 403.15696: 13.8432 13.6895
6: 403.15694: 13.8467 13.6854
7: 403.15694: 13.8475 13.6827
8: 403.15694: 13.8474 13.6826
lmer(score~trt+(1|id/eye),dat, verbose=T, start=list(matrix(.81),matrix(.57)))
1: 1304.7519: 2.42844 1.75398
2: 776.66857: 6.30805 0.00000
3: 657.56008: 7.73875 0.00000
4: 514.89384: 10.7783 0.000123881
5: 460.61372: 13.2310 0.00140974
6: 434.88979: 15.6779 0.00908817
7: 426.68762: 17.5388 0.0316599
8: 424.89288: 18.6843 0.0743135
9: 424.71437: 19.1269 0.132248
10: 424.70495: 19.2146 0.203365
11: 424.68492: 19.3175 0.420732
12: 424.52038: 19.7021 1.39710
13: 424.12511: 20.4042 3.37502
14: 422.97192: 20.4400 4.91187
15: 417.56194: 17.5366 5.92339
16: 406.81181: 13.4830 17.5342
17: 405.45194: 15.7947 12.4361
18: 404.19817: 15.0352 12.5415
19: 403.37771: 13.5806 13.0273
20: 403.37260: 13.7001 14.5562
21: 403.16170: 13.8571 13.8057
22: 403.15707: 13.8372 13.6724
23: 403.15694: 13.8483 13.6821
24: 403.15694: 13.8473 13.6827
25: 403.15694: 13.8474 13.6826
sessionInfo()
R version 2.10.1 (2009-12-14)
i386-apple-darwin9.8.0
locale:
[1] en_GB.UTF-8/en_GB.UTF-8/C/C/en_GB.UTF-8/en_GB.UTF-8
attached base packages:
[1] stats graphics grDevices utils datasets methods base
other attached packages:
[1] lme4_0.999375-35 Matrix_0.999375-43 lattice_0.18-3
loaded via a namespace (and not attached):
[1] grid_2.10.1 nlme_3.1-96 stats4_2.10.1
Quoting Peter Dalgaard <pdalgd at gmail.com>:
> On 09/18/2010 11:24 AM, Jarrod Hadfield wrote:
>> Fedora Core 13 seems fine:
>>
>> Linear mixed model fit by REML
>> Formula: score ~ trt + (1 | id/eye)
>> Data: dat
>> AIC BIC logLik deviance REMLdev
>> 425.2 474.2 -201.6 412.7 403.2
>> Random effects:
>> Groups Name Variance Std.Dev.
>> eye:id (Intercept) 3.59532 1.89613
>> id (Intercept) 3.51024 1.87356
>> Residual 0.01875 0.13693
>> Number of obs: 640, groups: eye:id, 160; id, 80
>
> Not in 32bit:
>
>> lmer(score~trt+(1|id/eye),dat)
> Linear mixed model fit by REML
> Formula: score ~ trt + (1 | id/eye)
> Data: dat
> AIC BIC logLik deviance REMLdev
> 446.7 495.8 -212.4 430.9 424.7
> Random effects:
> Groups Name Variance Std.Dev.
> eye:id (Intercept) 6.9208e+00 2.6307e+00
> id (Intercept) 4.0996e-12 2.0248e-06
> Residual 1.8750e-02 1.3693e-01
> Number of obs: 640, groups: eye:id, 160; id, 80
>
> Diddling the start values can change the result, though. Could you all try
>
> (a) setting verbose=T
> (b) try start=list(matrix(10),matrix(10))
> (c) same thing, but with .81 and .57
>
>
> In the (c) case, I get
>
>> l <- lmer(score~trt+(1|id/eye),dat,
> start=list(matrix(.81),matrix(.57)), verbose=T)
> 0: 2226.5319: 0.810000 0.570000
> 1: 1304.7507: 2.42844 1.75398
> 2: 776.66839: 6.30805 0.00000
> 3: 657.56008: 7.73875 0.00000
> 4: 514.89378: 10.7783 0.00000
> 5: 460.61371: 13.2310 0.00000
> 6: 434.88972: 15.6779 1.52026e-08
> 7: 426.68822: 17.5387 6.57868e-08
> 8: 424.89428: 18.6843 1.41900e-07
> 9: 424.71823: 19.1258 2.47916e-07
> 10: 424.71356: 19.2075 3.75766e-07
> 11: 424.71354: 19.2122 5.37851e-07
> 12: 424.71354: 19.2122 7.55161e-07
>
> which suggests an algorithm error where the parameters are getting stuck
> on the boundary. These values match the automatically generated starting
> values only to two digits, so I'm somewhat baffled that you apparently
> get something completely different in 64 bits!
>
> --
> Peter Dalgaard
> Center for Statistics, Copenhagen Business School
> Phone: (+45)38153501
> Email: pd.mes at cbs.dk Priv: PDalgd at gmail.com
>
>
--
The University of Edinburgh is a charitable body, registered in
Scotland, with registration number SC005336.
More information about the R-sig-mixed-models
mailing list