False zero estimates as reported by Vincent Dorie

[dmbates]

Initial Comment:
I'm getting estimates of zero variance for the random effects when those values are demonstrably not the MLE. These appear to be artifacts of a) optimization on the scale of a standard deviation leaving the derivative of the log-likelihood with a trivial root of 0 and b) overzealous first step taken by the optimizer causing it to only explore near 0.

Specifically, for a univariate, balanced group model:

y_ij | theta_j ~ind N(mu + theta_j, sigma^2_eps), i = 1, ..., n,
theta_j ~iid N(0, sigma^2_the * sigma^2_eps), j = 1, ..., J.

The MLE for sigma_the^2 should be:
max( [(n - 1) / n] * R - 1 / n, 0),
where R is the ratio of the between group sum of squares to the within group:

R = S_b^2 / S_w^2,
S_b^2 = n * sum_j (y.bar_j - y.bar)^2,
S_w^2 = sum_ij (y_ij - y.bar_j)^2.

In this setting, code to reproduce the bug:

library(lme4);

sigma.eps <- 2;
sigma.the <- 0.75;
mu <- 2;

n <- 5;
J <- 10;
g <- gl(J, n);

set.seed(1);

theta <- rnorm(J, 0, sigma.eps * sigma.the);
y <- rnorm(n * J, mu + theta[g], sigma.eps);

lmerFit <- lmer(y ~ 1 + (1 | g), REML = FALSE, verbose=T);

y.bar <- mean(y);
y.bar.j <- sapply(1:J, function(j) mean(y[g == j]));
S.w <- sum((y - y.bar.j[g])^2);
S.b <- n * sum((y.bar.j - y.bar)^2);
R <- S.b / S.w;

sigma.the.hat <- sqrt(max((n - 1) * R / n - 1 / n, 0));

round(sigma.the.hat, 2); # 0.27
round(lmerFit@ST[[1]][1], 2); # 0

# refitting at analytic MLE to check that the formula is
# correct, or at least has a lower deviance
mleFit <- lmer(y ~ 1 + (1 | g), REML = FALSE, verbose=T,
               start = list(ST = list(as.matrix(sigma.the.hat))));

round(lmerFit@deviance[["ML"]] - mleFit@deviance[["ML"]], 2); #0.44
abs(sigma.the.hat - mleFit@ST[[1]][1]); # 0

Workaround:

It's not optimal, but drmnfb (used by S_nlminb_iterate) has a parameter LMAX0 (doc at: (http://www.netlib.org/port/dmng.f) that limits the size of the first step. Setting this to be below the minimum of the S components in the TSST' Cholesky decomposition seems to work, although is restrictive in higher dimensions. Walking back the first step if it hits the boundary is probably better.

[bbolker]

For what it's worth, this seems to work in the latest development version (0.999999911-0):

round(getME(lmerFit,"theta"),2)
g.(Intercept) 
         0.27