[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