I confirm John's problems with lmer. I'm using R 2.11.1. on Windows XP.
R> m4 <- lmer(score~trt+(1|id/eye),dat)
R> m4
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
One other note--the variance components from asreml agree with *lme*:
ma <- asreml(score~trt , random=~id + id:eye, data=dat)
R> sqrt(ma$gammas * ma$sigma2)
id!id.var id:eye!id.var R!variance
1.8735643 1.8961309 0.1369306
Kevin
On Fri, Sep 17, 2010 at 3:50 PM, array chip <arrayprof...@yahoo.com> wrote:
>
>
> Thank you Peter. Actually 3 people from mixed model mailing list tried my
> code
> using lmer(). They got the same results as what I got from lme4(). So they
> couldn't replicate my lmer() results:
>
> 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
>
> The only difference they can think of is they are using Mac and FreeBSD
> while
> mine is PC. I can't imagine this can be the reason. I re-install lme4
> package,
> but still got weird results with lmer().
>
> Per your suggestion, here is the results for aov()
>
> summary(aov(score~trt+Error(id/eye), data=dat))
>
> Error: id
> Df Sum Sq Mean Sq F value Pr(>F)
> trt 7 1353.6 193.378 4.552 0.0002991 ***
> Residuals 72 3058.7 42.482
> ---
> Signif. codes: 0 *** 0.001 ** 0.01 * 0.05 . 0.1 1
>
> Error: id:eye
> Df Sum Sq Mean Sq F value Pr(>F)
> Residuals 80 1152 14.4
>
> Error: Within
> Df Sum Sq Mean Sq F value Pr(>F)
> Residuals 480 9 0.01875
>
> As can be seen, thr within subject variance estimate (0.01875) is the same
> between aov, lmer and lme. But I am not sure how to relate results between
> aov
> and lmer/lme for the other 2 variance components (id and eye%in%id).
>
> Thanks
>
> John
>
>
>
>
>
>
> ----- Original Message ----
> From: Peter Dalgaard <pda...@gmail.com>
> To: array chip <arrayprof...@yahoo.com>
> Cc: r-help@r-project.org
> Sent: Fri, September 17, 2010 1:05:27 PM
> Subject: Re: [R] lmer() vs. lme() gave different variance component
> estimates
>
> On 09/17/2010 09:14 PM, array chip wrote:
> > Hi, I asked this on mixed model mailing list, but that list is not very
> active,
> >
> > so I'd like to try the general R mailing list. Sorry if anyone receives
> the
> > double post.
> >
> >
> > 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()?
>
> That's your guess... I'd be more careful about jumping to conclusions.
> If this is a balanced data set, perhaps you could supply the result of
>
> summary(aov(score~trt+Error(id/eye), data=dat))
>
> The correct estimates should be computable from the ANOVA table.
>
>
> --
> Peter Dalgaard
> Center for Statistics, Copenhagen Business School
> Phone: (+45)38153501
> Email: pd....@cbs.dk Priv: pda...@gmail.com
>
>
>
>
>
> ______________________________________________
> R-help@r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide
> http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.
>
--
Kevin Wright
[[alternative HTML version deleted]]
______________________________________________
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
