[R] lmer() vs. lme() gave different variance component estimates
Peter Dalgaard
pdalgd at gmail.com
Sat Sep 18 10:35:45 CEST 2010
On 09/17/2010 10:50 PM, array chip 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).
>
For a nested design, the relation is quite straightforward: The residual
MS are the variances of sample means scaled to be comparable with the
residuals (so that in the absense of random components, all
MS are equal to within the F-ratio variability). So to get the id:eye
variance component, subtract the Within MS from the id:eye MS and divide
by the number of replicates (4 in this case since you have 640
observations on 160 eyes) (14.4 - 0.01875)/4 = 3.59, and similarly, the
id variance is the MS for id minus that for id:eye scaled by 8:
(42.482-14.4)/8 = 3.51.
I.e. it is reproducing the lmer results above, but of course not those
from your original post.
(Notice, by the way, that if you are only interested in the treatment
effect, you might as well have computed the average of all 8
measurements on each animal and computed a 1-way ANOVA).
--
Peter Dalgaard
Center for Statistics, Copenhagen Business School
Phone: (+45)38153501
Email: pd.mes at cbs.dk Priv: PDalgd at gmail.com
More information about the R-help
mailing list