[R-sig-ME] Constraining error variance to 0 in an lmer() model.

[Rolf Turner]

In a longish thread that he initiated ("What is the lmer/nlme equivalent 
of the REPEATED subcommand in SPSS's MIXED procedure?") Maarten Jung
answered a question that has long been bothering me.

Using an example data set provided by Ben Pelzer he showed how to fit
a model in which in effect the error variance is 0.  Many years ago I 
was struggling to analyse a data set which in its simplest form had
exactly the same structure as Ben Pelzer's data, essentially repeated 
measures for each of a number of subjects, the number of repetitions
("occasions" in Ben Pelzer's data) being the same for each subject 
("person").

I asked for help with this problem, but got no really satisfactory answer.

The model that I wanted to fit was (adapting Maarten Jung's notation 
slightly to suit my own tastes):

     (*) lmer(test ~ 0+person + (0+person | occasion),data=Dat)

The reason that I wanted to fit this model was to try to verify that I 
understood what was going on (in this minimal instance) before 
proceeding to deal with more complicated models).  The point is that the
foregoing model can be "fitted" without recourse to lmer() simply by
thinking of the data as being multivariate (normal).  I.e. each "person"
has an associated tri-variate vector of observations, assumed to be 
independent from person to person.  If we form the data matrix, say "X" 
(with 20 rows and 3 columns) and calculate the column means and
var(X)/20 we should get the same results as would be produced by lmer(),
or so I thought (given that I understood correctly what lmer() was doing.)

The problem that I encountered (back in the dim dark past) was that the 
code in (*) falls over.  This is understandable since there is 
redundancy in the model, as Maarten Jung explained.  He also explained
how to get around the problem by adding a control argument:

     lmer(test ~ 0+person + (0+person | occasion),data=Dat,
                 control=lmerControl(check.nobs.vs.nRE = "ignore"))

does not fall over.  And it indeed gives the same results as does the 
multivariate analysis, to at least 5 decimal places in all instances:

> lmer results:
> fixef(fit): > occasion1 occasion2 occasion3
>     29.80     30.05     33.30
>
> vcov(fit):
> 3 x 3 Matrix of class "dpoMatrix"
>           occasion1 occasion2 occasion3
> occasion1 1.7821054  1.150526 0.5768422
> occasion2 1.1505265  2.249868 3.0623683
> occasion3 0.5768422  3.062368 5.8531577
> 
> Multivariate normal results:
> column means: >
> occasion1 occasion2 occasion3 
>     29.80     30.05     33.30
>
> var(X)/20:
>           occasion1 occasion2 occasion3
> occasion1 1.7821053  1.150526 0.5768421
> occasion2 1.1505263  2.249868 3.0623684
> occasion3 0.5768421  3.062368 5.8531579

Unfortunately lmer() still gives an off-putting warning message:

> Warning message:
> In checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv,  :
>   Model is nearly unidentifiable: large eigenvalue ratio
>  - Rescale variables?

I have no idea what is actually going on here, but it seems to me that 
this should not be happening.  The setting is dead simple, and the 
answer that lmer() obtains is "exactly" correct, as we can check in this 
instance by other means.

What we have here is a "boundary" or "edge" or "corner" case, and I have 
heard it asserted by someone knowledgeable and respected (can't remember 
who, it was a long while ago) that in testing your software it is 
crucial to see how it behaves with such "boundary" cases.  In the 
current instance lmer() is handling this boundary case alright, but 
perhaps worrying the user unnecessarily.

Can anything be done?  I would judge by the reactions to my previous 
enquiries on this matter that the answer is "no", but I thought I'd 
throw the question out to the cognoscenti again, prompted by being shown 
by Maarten Jung how to at least get an answer out of lmer().

cheers,

Rolf Turner

-- 
Technical Editor ANZJS
Department of Statistics
University of Auckland
Phone: +64-9-373-7599 ext. 88276