[R-sig-ME] lmer() vs. lme() gave different variance component estimates
Ista Zahn
izahn at psych.rochester.edu
Fri Sep 17 22:56:08 CEST 2010
I get the "PC" results on arch linux:
> 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) 3.0516e-17 5.5241e-09
Residual 1.8750e-02 1.3693e-01
Number of obs: 640, groups: eye:id, 160; id, 80
Fixed effects:
Estimate Std. Error t value
(Intercept) 2.5500 0.5885 4.333
trtB 0.9875 0.8322 1.187
trtC 0.8625 0.8322 1.036
trtControl 1.1500 0.8322 1.382
trtD 2.4125 0.8322 2.899
trtE 1.4375 0.8322 1.727
trtF 2.1750 0.8322 2.614
trtG -2.5500 0.8322 -3.064
> sessionInfo()
R version 2.11.1 (2010-05-31)
i686-pc-linux-gnu
locale:
[1] LC_CTYPE=en_US.utf8 LC_NUMERIC=C
[3] LC_TIME=en_US.utf8 LC_COLLATE=en_US.utf8
[5] LC_MONETARY=C LC_MESSAGES=en_US.utf8
[7] LC_PAPER=en_US.utf8 LC_NAME=C
[9] LC_ADDRESS=C LC_TELEPHONE=C
[11] LC_MEASUREMENT=en_US.utf8 LC_IDENTIFICATION=C
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
-Ista
On Sat, Sep 18, 2010 at 4:34 PM, Andrew Robinson
<A.Robinson at ms.unimelb.edu.au> wrote:
> So far the results split on platform lines: John's is pc, Daniel's and
> Reinhold's are Mac.
>
> I get the Mac result on FreeBSD.
>
>> 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
> 425.2 474.2 -201.6 412.7 403.2
> Random effects:
> Groups Name Variance Std.Dev.
> eye:id (Intercept) 3.59531 1.89613
> id (Intercept) 3.51025 1.87357
> Residual 0.01875 0.13693
> Number of obs: 640, groups: eye:id, 160; id, 80
>
> Fixed effects:
> Estimate Std. Error t value
> (Intercept) 2.5500 0.7287 3.499
> trtB 0.9875 1.0305 0.958
> trtC 0.8625 1.0305 0.837
> trtControl 1.1500 1.0305 1.116
> trtD 2.4125 1.0305 2.341
> trtE 1.4375 1.0305 1.395
> trtF 2.1750 1.0305 2.111
> trtG -2.5500 1.0305 -2.474
>
> Correlation of Fixed Effects:
> (Intr) trtB trtC trtCnt trtD trtE trtF
> trtB -0.707
> trtC -0.707 0.500
> trtControl -0.707 0.500 0.500
> trtD -0.707 0.500 0.500 0.500
> trtE -0.707 0.500 0.500 0.500 0.500
> trtF -0.707 0.500 0.500 0.500 0.500 0.500
> trtG -0.707 0.500 0.500 0.500 0.500 0.500 0.500
>> sessionInfo()
> R version 2.11.1 (2010-05-31)
> x86_64-unknown-freebsd8.0
>
> locale:
> [1] C
>
> 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.19-11
>
> loaded via a namespace (and not attached):
> [1] grid_2.11.1 nlme_3.1-96 stats4_2.11.1
>
>
>
> Cheers
>
> Andrew
>
>
> On Fri, Sep 17, 2010 at 04:57:59PM -0300, Daniel wrote:
>> I really don't know what is going wrong.
>>
>> sessionInfo()
>> R version 2.11.1 (2010-05-31)
>> x86_64-apple-darwin9.8.0
>>
>> locale:
>> [1] en_US.UTF-8/en_US.UTF-8/C/C/en_US.UTF-8/en_US.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-8
>>
>> loaded via a namespace (and not attached):
>> [1] grid_2.11.1 nlme_3.1-96 stats4_2.11.1
>>
>>
>> On Fri, Sep 17, 2010 at 4:33 PM, array chip <arrayprofile at yahoo.com> wrote:
>>
>> > Hi, Reinhold and Daniel,
>> >
>> > I just re-installed lme4 package, still got different results than yours.
>> >
>> > Here is my sessionInfo():
>> >
>> > R version 2.11.1 (2010-05-31)
>> > i386-pc-mingw32
>> >
>> > locale:
>> > [1] LC_COLLATE=English_United States.1252
>> > [2] LC_CTYPE=English_United States.1252
>> > [3] LC_MONETARY=English_United States.1252
>> > [4] LC_NUMERIC=C
>> > [5] LC_TIME=English_United States.1252
>> >
>> > attached base packages:
>> > [1] stats graphics grDevices datasets utils methods base
>> >
>> > other attached packages:
>> > [1] lme4_0.999375-35 Matrix_0.999375-43 lattice_0.18-8
>> >
>> > loaded via a namespace (and not attached):
>> > [1] grid_2.11.1 nlme_3.1-96 stats4_2.11.1
>> >
>> > John
>> >
>> >
>> >
>> >
>> >
>> >
>> >
>> >
>> >
>> > ----- Original Message ----
>> > From: Reinhold Kliegl <reinhold.kliegl at gmail.com>
>> > To: array chip <arrayprofile at yahoo.com>
>> > Cc: r-sig-mixed-models at r-project.org
>> > Sent: Fri, September 17, 2010 12:14:13 PM
>> > Subject: Re: [R-sig-ME] lmer() vs. lme() gave different variance component
>> > estimates
>> >
>> > Not on my computer. Perhaps you could provide sessionInfo()?
>> >
>> > > 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
>> > 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
>> >
>> > > sessionInfo()
>> > R version 2.11.1 Patched (2010-07-16 r52550)
>> > x86_64-apple-darwin9.8.0
>> >
>> > locale:
>> > [1] en_US.UTF-8/en_US.UTF-8/C/C/en_US.UTF-8/en_US.UTF-8
>> >
>> > attached base packages:
>> > [1] stats graphics grDevices utils datasets methods base
>> >
>> > other attached packages:
>> > [1] lme4_0.999375-34 Matrix_0.999375-43 lattice_0.18-8
>> >
>> > loaded via a namespace (and not attached):
>> > [1] grid_2.11.1 nlme_3.1-96 stats4_2.11.1 tools_2.11.1
>> >
>> > Reinhold Kliegl
>> >
>> >
>> > On Fri, Sep 17, 2010 at 8:39 PM, array chip <arrayprofile at yahoo.com>
>> > wrote:
>> > > Hi, I have a dataset of animals receiving some eye treatments. There are
>> > 8
>> > > treatments, each animal's right and left eye was measured with some
>> > scores
>> > > (ranging from 0 to 7) 4 times after treatment. So there are nesting
>> > groups
>> > eyes
>> > > within animal. Dataset attached
>> > >
>> > >> dat<-read.table("dat.txt",sep='\t',header=T,row.names=1)
>> > >> dat$id<-factor(dat$id)
>> > >> str(dat)
>> > > 'data.frame': 640 obs. of 5 variables:
>> > > $ score: int 0 2 0 7 4 7 0 2 0 7 ...
>> > > $ id : Factor w/ 80 levels "1","3","6","10",..: 7 48 66 54 18 26 38 52
>> > 39
>> > 63
>> > > ...
>> > > $ rep : int 1 1 1 1 1 1 1 1 1 1 ...
>> > > $ eye : Factor w/ 2 levels "L","R": 2 2 2 2 2 2 2 2 2 2 ...
>> > > $ trt : Factor w/ 8 levels "A","B","C","Control",..: 1 1 1 1 1 1 1 1 1
>> > 1 ...
>> > >
>> > > I fit a mixed model using both lmer() from lme4 package and lme() from
>> > nlme
>> > > package:
>> > >
>> > >> 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.630742315798
>> > > id (Intercept) 1.4471e-16 0.000000012030
>> > > Residual 1.8750e-02 0.136930641909
>> > > Number of obs: 640, groups: eye:id, 160; id, 80
>> > >
>> > >> summary(lme(score~trt, random=(~1|id/eye), dat))
>> > >
>> > > Linear mixed-effects model fit by REML
>> > > Data: dat
>> > > AIC BIC logLik
>> > > 425.1569 474.0947 -201.5785
>> > >
>> > > Random effects:
>> > > Formula: ~1 | id
>> > > (Intercept)
>> > > StdDev: 1.873576
>> > >
>> > > Formula: ~1 | eye %in% id
>> > > (Intercept) Residual
>> > > StdDev: 1.896126 0.1369306
>> > >
>> > > As you can see, the variance components estimates of random effects are
>> > quite
>> > > different between the 2 model fits. From the data, I know that the
>> > variance
>> > > component for "id" can't be near 0, which is what lmer() fit produced, so
>> > I
>> > > think the lme() fit is correct while lmer() fit is off. This can also be
>> > seen
>> > > from AIC, BIC etc. lme() fit has better values than lmer() fit.
>> > >
>> > >
>> > > I guess this might be due to lmer() didn't converge very well, is there
>> > anyway
>> > > to adjust to make lmer() converge better to get similar results as lme()?
>> > >
>> > > Thanks
>> > >
>> > > John
>> > >
>> > >
>> > >
>> > > _______________________________________________
>> > > R-sig-mixed-models at r-project.org mailing list
>> > > https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
>> > >
>> > >
>> >
>> >
>> >
>> >
>> >
>> >
>>
>>
>> --
>> Daniel Marcelino
>> Skype: dmsilv
>> http://bit.ly/pol4vc
>>
>> [[alternative HTML version deleted]]
>>
>> _______________________________________________
>> R-sig-mixed-models at r-project.org mailing list
>> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
>
> --
> Andrew Robinson
> Program Manager, ACERA
> Department of Mathematics and Statistics Tel: +61-3-8344-6410
> University of Melbourne, VIC 3010 Australia (prefer email)
> http://www.ms.unimelb.edu.au/~andrewpr Fax: +61-3-8344-4599
> http://www.acera.unimelb.edu.au/
>
> _______________________________________________
> R-sig-mixed-models at r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
>
--
Ista Zahn
Graduate student
University of Rochester
Department of Clinical and Social Psychology
http://yourpsyche.org
More information about the R-sig-mixed-models
mailing list