lmer() vs. lme() gave different variance component estimates
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 <pdalgd at gmail.com>
To: array chip <arrayprofile at yahoo.com>
Cc: r-help at 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.mes at cbs.dk Priv: PDalgd at gmail.com