[R-sig-ME] [R] lmer() vs. lme() gave different variance component estimates

[array chip]

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