[R-sig-ME] Constraining error variance to 0 in an lmer() model.
Rolf Turner
r.turner at auckland.ac.nz
Sun Mar 25 00:40:17 CET 2018
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
More information about the R-sig-mixed-models
mailing list