正如在此答案中指出的那样，并且在其中一条评论中也提到过，问题似乎是局部最大值。为了更清楚地看到这一点，我在下面编写了一个简单的代码来计算该模型的负对数似然并使用optim(). 从不同的初始值开始会导致两种不同的解决方案：

# 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

我同意@DimitrisRizopoulos 的回答，还有几点要说。

• 首先我会说我很不高兴lmer没有找到最好的答案——尽管我怀疑这种情况可能仅限于小的、不寻常的（见下文）数据集。lme可能做得更好的原因之一是它适合对数标准偏差量表，这可能会使接近零的最小值“更广泛”。
• 您可以通过为缩放的标准偏差 ( ) 设置一个明确的、较低的起始值lmer来复制结果；基于下面的探索，或者任何较低的值都可以正常工作。对于它的价值，这将导致估计的随机效应方差为 0（以及“奇异拟合”消息，以及相当于完全忽略随机效应分量并使用...的答案）lmestart=...start=8lm()
• 在这种特殊情况下，使用“nloptwrap”优化器没有帮助；事实上，所有lmer可以使用的优化器，从默认起始值​​（（缩放标准偏差）= 1.0）开始，找到远离零的更高局部最小值。$\theta$
• 这是等效于lmer默认情况下用于查找起始值的方法的代码，当仅存在截距值随机效应时（请参见 此处）：

v0 <- with(multi,var(ave(y,gid)))  ## variance among group values  
v.e <- var(multi$y)-v0             ## residual var ~ total var - group variance
sqrt(v0/v.e)                       ## convert to scaled standard deviation

的起始值。$\theta=10.8$

• 我们可以系统地看到不同的起始值如何产生不同的结果：

m0 <- lmer(y~x+(1|(gid)), data=multi, REML=TRUE)
tvec2 <- seq(0,20,length=51)
ff <- function(t0) getME(update(m0,start=t0),"theta")
v <- sapply(tvec2,ff)
plot(tvec2,v)
abline(v=10.8,col="red")

• 我们还可以显式地可视化（负对数）似然面：

## helper function to capture fitting trajectory
cfun  <- function(...) {
    cc <- capture.output(x <- do.call(lmer,c(list(...),list(verbose=100))))
    gfun <- function(x,s) {
        as.numeric(gsub(s,"",grep(s,x,value=TRUE)))
    }
    it <- gfun(cc,"iteration: +")
    xval <- gfun(cc,"\tx = ")
    fval <- gfun(cc,"\tf\\(x\\) = +")
    attr(x,"optvals") <- data.frame(it,xval,fval)
    return(x)
}

c0 <- cfun(y~x+(1|(gid)), data=multi, REML=TRUE)
c1 <- cfun(y~x+(1|(gid)), data=multi, REML=FALSE)

f <- as.function(m0)
tvec <- seq(0,100,length=101)
dvec <- sapply(tvec,f)
m3 <- update(m0,REML=FALSE)
f2 <- as.function(m3)
dvec2 <- sapply(tvec,f2)
par(las=1,bty="l")
matplot(tvec,cbind(dvec,dvec2),type="l",
        ylab="deviance/REMLcrit",
        xlab="scaled standard dev")
with(attr(c0,"optvals"),text(xval,fval,it))
with(attr(c1,"optvals"),text(xval,fval,it,col=2))
legend("bottomright",c("REML","ML"),
       col=1:2,lty=1:2)

数字显示尝试的值的顺序。我们可以看到，只是ML 曲线的形状略有不同，它使优化器倾向于边界拟合而不是内部拟合。

• 这些数据是人为的吗？下面的左图按组显示数据；右图显示减去组均值的值。每组内的 5 个值之间几乎没有变化......

• 如果我们模拟具有相同属性的数据（从估计的系数开始），但变化实际上是高斯的，我们根本不会得到相同类型的多峰曲面：

multi_sim <- transform(multi,y=simulate(m0,seed=101)[[1]])
f3 <- as.function(update(m0,data=multi_sim))
dvec3 <- sapply(tvec,f3)
plot(tvec,dvec3,type="l")