[R-sig-ME] lmer() vs. lme() gave different variance component estimates
Luca Borger
lborger at uoguelph.ca
Sat Sep 18 23:17:24 CEST 2010
> 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 case this is of any use, I get (Win-XP, 32-bit, sessionInfo() at the
end):
> mod_a <- lmer(score ~ trt+(1|id/eye), dat, verbose=T)
0: 2219.2409: 0.816497 0.577350
1: 1300.6006: 2.44030 1.76298
2: 774.29059: 6.33180 0.00000
3: 655.91611: 7.76290 0.00000
4: 514.09907: 10.8045 0.00000
5: 460.23606: 13.2559 0.00000
6: 434.74990: 15.6993 0.00000
7: 426.65081: 17.5537 0.00000
8: 424.88919: 18.6916 1.16103e-08
9: 424.71803: 19.1277 2.38735e-08
10: 424.71356: 19.2077 4.23682e-08
11: 424.71354: 19.2122 6.29032e-08
12: 424.71354: 19.2122 8.78525e-08
>
> mod_a
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) 1.4471e-16 1.2030e-08
Residual 1.8750e-02 1.3693e-01
Number of obs: 640, groups: eye:id, 160; id, 80
> mod_b <- lmer(score ~ trt+(1|id/eye), dat,
> start=list(matrix(10),matrix(10)), verbose=T)
0: 429.50725: 10.0000 10.0000
1: 406.77628: 12.7928 11.2494
2: 403.69604: 13.2997 12.8550
3: 403.18732: 13.6909 13.5238
4: 403.15739: 13.8229 13.6825
5: 403.15696: 13.8431 13.6896
6: 403.15694: 13.8467 13.6854
7: 403.15694: 13.8475 13.6828
8: 403.15694: 13.8474 13.6826
Warning message:
In is.na(x) : is.na() applied to non-(list or vector) of type 'NULL'
>
> mod_b
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 3.59532 1.89613
id 3.51024 1.87356
Residual 0.01875 0.13693
Number of obs: 640, groups: eye:id, 160; id, 80
> mod_c <- 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.56009: 7.73875 0.00000
4: 514.89375: 10.7783 0.00000
5: 460.61370: 13.2310 5.55866e-08
6: 434.88975: 15.6779 3.94635e-07
7: 426.68823: 17.5387 1.39786e-06
8: 424.89427: 18.6843 3.27134e-06
9: 424.71823: 19.1258 5.82021e-06
10: 424.71356: 19.2075 8.90695e-06
11: 424.71354: 19.2122 1.28100e-05
12: 424.71354: 19.2122 1.80876e-05
Warning message:
In is.na(x) : is.na() applied to non-(list or vector) of type 'NULL'
>
> mod_c
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 6.9208e+00 2.6307e+00
id 6.1342e-12 2.4767e-06
Residual 1.8750e-02 1.3693e-01
Number of obs: 640, groups: eye:id, 160; id, 80
> sessionInfo()
R version 2.11.1 (2010-05-31)
i386-pc-mingw32
locale:
[1] LC_COLLATE=English_United Kingdom.1252 LC_CTYPE=English_United
Kingdom.1252 LC_MONETARY=English_United Kingdom.1252
[4] LC_NUMERIC=C LC_TIME=English_United
Kingdom.1252
attached base packages:
[1] stats graphics grDevices utils datasets methods base
other attached packages:
[1] lme4_0.999375-35 Matrix_0.999375-44 lattice_0.19-11
loaded via a namespace (and not attached):
[1] grid_2.11.1 nlme_3.1-96 stats4_2.11.1 tools_2.11.1
>
HTH
Cheers,
Luca
----- Original Message -----
From: "Peter Dalgaard" <pdalgd at gmail.com>
To: "Jarrod Hadfield" <j.hadfield at ed.ac.uk>
Cc: <r-sig-mixed-models at r-project.org>
Sent: Saturday, September 18, 2010 5:58 AM
Subject: Re: [R-sig-ME] lmer() vs. lme() gave different variance component
estimates
> 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
>
> _______________________________________________
> 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