Solved – Completely different results from lme() and lmer()

As it was noted in this answer, and also mentioned in one of the comments, the problem seems to be a local maximum. To see this more clearly I have written below a simple code to calculate the negative log-likelihood of this model and do the optimization using optim(). Starting with different initial values leads to the two different solutions:

# data multi <- structure(list(x = c(4.9, 4.84, 4.91, 5, 4.95, 3.94, 3.88, 3.95,                                4.04, 3.99, 2.97, 2.92, 2.99, 3.08, 3.03, 2.01, 1.96, 2.03, 2.12,                                2.07, 1.05, 1, 1.07, 1.16, 1.11),                          y = c(3.2, 3.21, 3.256, 3.25,                                3.256, 3.386, 3.396, 3.442, 3.436, 3.442, 3.572, 3.582, 3.628,                                3.622, 3.628, 3.758, 3.768, 3.814, 3.808, 3.814, 3.944, 3.954,                                4, 3.994, 4),                          pid = 1:25,                          gid = c(1L, 1L, 1L, 1L, 1L, 2L, 2L,                                  2L, 2L, 2L, 3L, 3L, 3L, 3L, 3L, 4L, 4L, 4L, 4L, 4L, 5L, 5L, 5L,                                  5L, 5L)), class = "data.frame", row.names = c(NA, -25L))  # function to calculate the negative log-likelihood of the random intercepts model library("mvtnorm") logLik <- function (thetas, y, X, id) {     ncX <- ncol(X)     betas <- thetas[seq_len(ncX)]     sigma_b <- exp(thetas[ncX + 1])     sigma <- exp(thetas[ncX + 2])     eta <- c(X %*% betas)     unq_id <- unique(id)     n <- length(unq_id)     lL <- numeric(n)     for (i in seq_len(n)) {         id_i <- id == unq_id[i]         n_i <- sum(id_i)         V_i <- matrix(sigma_b^2, n_i, n_i)         diag(V_i) <- diag(V_i) + sigma^2         lL[i] <- dmvnorm(y[id_i], mean = eta[id_i], sigma = V_i, log = TRUE)     }     - sum(lL, na.rm = TRUE) }  # optimization using as initial values 0 for the fixed effects,  # and 1 for the variance components  opt <- optim(rep(0, 4), logLik, method = "BFGS",               y = multi$y, X = cbind(1, multi$x), id = multi$gid)  opt$par[1:2] # fixed effects #> [1] 2.855872 0.250341 exp(opt$par[3]) # sd random intercepts #> [1] 0.6029724 exp(opt$par[4]) # sd error terms #> [1] 0.01997889  # optimization using as initial values 4 & -0.2 for the fixed effects,  # and 0.0003 and 0.034 for the variance components  opt2 <- optim(c(4, -0.2, -8, -3.4), logLik, method = "BFGS",                y = multi$y, X = cbind(1, multi$x), id = multi$gid)  opt2$par[1:2] # fixed effects #> [1]  4.1846965 -0.1928397 exp(opt2$par[3]) # sd random intercepts #> [1] 0.000270746 exp(opt2$par[4]) # sd error terms #> [1] 0.03239167