[R-sig-ME] How does lmer obtain ML estimates that are not stationary points?

[Ben Bolker]

Asaf Weinstein <asafw.at.wharton at ...> writes:

> 
> Hi,
> 
> I am looking at a two-way random-effects ANOVA layout as a particular case
> of the general mixed-model,
> 
> y|b ~ N(X beta + Zb, sigsq I)
> b ~ N(0,sigsq Gamma_theta )    [Gamma diagonal],
> 
> I am trying to compute ML estimates for theta under a KNOWN sigsq (i.e.,
> error variance is known). I derived my own ML estimates since lmer()
> estimates sigsq rather than assuming it is known.
> As long as the variance components (or theta's) are "big" (far from zero),
> the output of my algorithm is consistent with the output of lmer (when
> plugging in the estimated sigsq into my functions); but I run into problems
> when the ML estimates for at least one of the theta's is close zero.
> 
> Before i roll up my sleeves, I would like to know how lmer handles the case
> in which the value of theta which nullifies the derivative is negative (or
> does not exist at all), i.e., the ML estimate cannot be obtained by
> searching for a root of the derivative of the profile likelihood.
> 
> (If my point was not clear, I refer to the case similar to what happens in
> the balanced two-way ANOVA when the row (or column) sum of squares minus
> the estimated error component is negative..).

  lmer fits models on a constrained space where the variances are not
allowed to be negative.  So it would give the best fit on the boundary
of the feasible space (although I would be very slightly suspicious of
the results in this case; it is easy to misconverge / run into
optimization difficulties when the results are on the boundary,
although I don't have any concrete examples of where lmer gets this 
wrong).

  With the development version of R, you can get the deviance function
and evaluate it yourself over a range of parameter values to see what
it does.