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