[R-sig-ME] lmer() vs. lme() gave different variance component estimates
array chip
arrayprofile at yahoo.com
Fri Sep 17 21:33:38 CEST 2010
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
>
>
More information about the R-sig-mixed-models
mailing list