lmer: problem in crossed random effect model with verydifferent variances

On Wed, Jun 17, 2009 at 4:02 PM, Douglas Bates<bates at stat.wisc.edu> wrote:

On Wed, Jun 17, 2009 at 2:48 PM, Doran, Harold<HDoran at air.org> wrote:

Michael

I'll take a quick stab at this, but there is really no way to know what
the issue is given that there is no real description of your data.
First, SAS and lmer use different algorithms for generating parameter
estimates, so it's no surprise that the world does not line up exactly
between the two. However, the estimates should be similar.

Lmer used REML by default. What did you use in SAS, ML or REML

The default is in SAS is also REML. So REML was used for both lmer and
PROC MIXED.

What are the values of the REML criterion (or the deviance, if you
used ML estimates) at convergence? ?If they are very close then it is
just a matter of the convergence criterion.

The REML criteria in fact were also close at convergence (SAS gives
448.45 while lmer says 448.5).

Consider that the
relative variance for the analyst in the SAS fit is less than 0.01.
When you estimate a variance relative to the residual variance as 0.01
it means "essentially zero"

You should read the result from lmer as an estimate of zero. ?For some
reason the optimization software doesn't like to converge on the
boundary and often ends up at very small but non-zero values.

Certainly in practice, the difference between 1500 and 10^-12 does not
really matter
when the error variance is 200,000. ?Still, 1500 seems to be more 'palatable'.
Anyway, it was just surprising to see the difference even though I
understand that lmer uses
difference algorithm and did not expect the answers to be exactly the same.

By the way, at log scale, SAS just gives 0 for analyst variance while
lmer says 10^-15.

I am attaching the code and data at the end of the email in case
anyone wants to try it out.

Thanks for the clarification.

Michael

R:
lmer (y ~ (1 | day) + (1 | analyst), data = tmp)

SAS:
proc mixed data = tmp;
? ? ? class day analyst;
? ? ? model y = / s;
? ? ? random day analyst;
run;

tmp:

? day analyst ? ?y
1 ? ?1 ? ? ? 1 5482
2 ? ?1 ? ? ? 1 3285
3 ? ?1 ? ? ? 1 4266
4 ? ?1 ? ? ? 1 3818
5 ? ?1 ? ? ? 1 4159
6 ? ?2 ? ? ? 1 3007
7 ? ?2 ? ? ? 1 3349
8 ? ?2 ? ? ? 1 3178
9 ? ?2 ? ? ? 1 3093
10 ? 2 ? ? ? 1 3242
11 ? 3 ? ? ? 1 2495
12 ? 3 ? ? ? 1 2687
13 ? 3 ? ? ? 1 2858
14 ? 3 ? ? ? 1 2090
15 ? 3 ? ? ? 1 2218
16 ? 1 ? ? ? 2 3882
17 ? 1 ? ? ? 2 4522
18 ? 1 ? ? ? 2 3647
19 ? 1 ? ? ? 2 3690
20 ? 1 ? ? ? 2 3754
21 ? 2 ? ? ? 2 3050
22 ? 2 ? ? ? 2 2858
23 ? 2 ? ? ? 2 3178
24 ? 2 ? ? ? 2 2901
25 ? 2 ? ? ? 2 2410
26 ? 3 ? ? ? 2 1855
27 ? 3 ? ? ? 2 3157
28 ? 3 ? ? ? 2 2111
29 ? 3 ? ? ? 2 2431
30 ? 3 ? ? ? 2 3114

str(dat)

summary(aov(y ~ day * analyst, dat))

summary(aov(y ~ day  + analyst, dat))