[R] Fitting mixed-effects models with lme with fixed error term variances

[Douglas Bates]

On 11/26/06, Gregor Gorjanc <gregor.gorjanc at bfro.uni-lj.si> wrote:
> Hello!
>
> Douglas Bates wrote:
> ...
> >> > Not quite.  There is now a capability in lmer to specify starting
> >> > estimates for the relative precision matrices which means that you can
> >> > use the estimates from one fit as starting estimates for another fit.
> >> > It looks like
> >> >
> >> > fm1 <- lmer(...)
> >> > fm2 <- lmer(y ~ x + (...), start = fm1 at Omega)
> >> >
> >> > I should say that in my experience this has not been as effective as I
> >> > had hoped it would be.  What I see in the optimization is that the
> >> > first few iterations reduce the deviance quite quickly but the
> >> > majority of the time is spent refining the estimates near the optimum.
> >> > So, for example, if it took 40 iterations to converge from the rather
> >> > crude starting estimates calculated within lmer it might instead take
> >> > 35 iterations if you give it good starting estimates.
> >>
> >> Yes, I also have the same experience with other programs, say VCE[1].
> >> What I was hopping for was just solutions from linear mixed model i.e.
> >> either via GLS or PLS representations and no optimization if values for
> >> variance components are passed as input - there should be a way to
> >> distinguish that from just passing starting values..
> >
> > The PLS representation in lmer is in terms of the relative
> > variance-covariance matrices of the random effects where "relative"
> > means relative to s^2.  (In fact, the Omega slot contains the relative
>
> What exactly is s^2?

Sorry, I confused a parameter and its estimate.  The
variance-covariance matrices are relative to the variance of the
per-observation noise which is assumed to be Gaussian with mean 0 and
variance-covariance matrix sigma^2 * I.  That is, they are relative to
sigma^2 which is estimated by s^2 the penalized residual sum of
squares divided by the number of observations.

> > precision matrices (inverses of the relative variance matrices) but
> > its the same idea.  All variance components are expressed relative to
> > s^2.)  If it is sufficient to fix these then you can easily do that.
> > Just follow the code in lmer until it gets to
> >
> >   .Call(mer_ECMEsteps, mer, cv$niterEM, cv$EMverbose)
> >   LMEoptimize(mer) <- cv
> >   return(new("lmer", mer,
> >
> > and replace the first two lines by
> >
> >   .Call("mer_factor", mer, PACKAGE = "lme4")
>
> Hmm. Sounds simple. I will try this. Do you find any value in adding
> this as an option in lmer or perhaps new? function, say lmeSol() or
> something similar?

I'll have to think about it.  It is the usual problem of determining
how to specify this option, documenting it, explaining it and
maintaining it along with the rest of the code.

> > A side-effect of performing the factorization is to evaluate the ML
> > and REML deviances and store the result in the deviance slot.
> >
> > The "follow the code in lmer" part will get a bit tricky because of
> > the number of functions that are hidden in the lme4 namespace but I'm
> > sure you know how to get around that.
>
> Thank you very much!
>
> Gregor
>
>
>
>