[R-sig-ME] lmer() vs. lme() gave different variance component estimates
array chip
arrayprofile at yahoo.com
Fri Sep 17 20:39:17 CEST 2010
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
-------------- next part --------------
An embedded and charset-unspecified text was scrubbed...
Name: dat.txt
URL: <https://stat.ethz.ch/pipermail/r-sig-mixed-models/attachments/20100917/4662c616/attachment.txt>
More information about the R-sig-mixed-models
mailing list