r - 在 r 中没有警告消息的情况下，我如何知道模型是否在 lme4 中收敛或未能收敛？

Question

重复

例如，我可以通过以下方式评估这个多级模型是否是单一拟合使用 isSingular() 函数。

同样，有什么办法可以知道这个模型是收敛还是收敛失败？

我的顾问说，如果模型无法收敛，就不会估计标准误差。然而，尽管下面未能收敛，但似乎估计了一个标准误差。

是否有任何好的指标表明该模型收敛或未能收敛？ (注意警告消息除外)

我正在使用 lme4 包和 lmer() 函数。

比如有收敛多级模型失败的例子

library(lme4)
read.table(textConnection("duration season  sites   effect
                          4d    mon s1  7305.91
                          4d    mon s2  856.297
                          4d    mon s3  649.93
                          4d    mon s1  10121.62
                          4d    mon s2  5137.85
                          4d    mon s3  3059.89
                          4d    mon s1  5384.3
                          4d    mon s2  5014.66
                          4d    mon s3  3378.15
                          4d    post    s1  6475.53
                          4d    post    s2  2923.15
                          4d    post    s3  554.05
                          4d    post    s1  7590.8
                          4d    post    s2  3888.01
                          4d    post    s3  600.07
                          4d    post    s1  6717.63
                          4d    post    s2  1542.93
                          4d    post    s3  1001.4
                          4d    pre s1  9290.84
                          4d    pre s2  2199.05
                          4d    pre s3  1149.99
                          4d    pre s1  5864.29
                          4d    pre s2  4847.92
                          4d    pre s3  4172.71
                          4d    pre s1  8419.88
                          4d    pre s2  685.18
                          4d    pre s3  4133.15
                          7d    mon s1  11129.86
                          7d    mon s2  1492.36
                          7d    mon s3  1375
                          7d    mon s1  10927.16
                          7d    mon s2  8131.14
                          7d    mon s3  9610.08
                          7d    mon s1  13732.55
                          7d    mon s2  13314.01
                          7d    mon s3  4075.65
                          7d    post    s1  11770.79
                          7d    post    s2  4254.88
                          7d    post    s3  753.2
                          7d    post    s1  11324.95
                          7d    post    s2  5133.76
                          7d    post    s3  2156.2
                          7d    post    s1  12103.76
                          7d    post    s2  3143.72
                          7d    post    s3  2603.23
                          7d    pre s1  13928.88
                          7d    pre s2  3208.28
                          7d    pre s3  8015.04
                          7d    pre s1  11851.47
                          7d    pre s2  6815.31
                          7d    pre s3  8478.77
                          7d    pre s1  13600.48
                          7d    pre s2  1219.46
                          7d    pre s3  6987.5
                          "),header=T)->dat1

lmer(effect ~ duration + (1+duration|sites) +(1+duration|season),
                                  data=dat1)

这会产生错误警告信息:模型未能收敛到 1 个负特征值:-2.3e+01

然而，尽管未能收敛，但似乎估计了标准误差。

Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
Formula: effect ~ duration + (1 + duration | sites) + (1 + duration |      season)
   Data: dat1

REML criterion at convergence: 969

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-2.0515 -0.6676  0.0075  0.5333  3.2161 

Random effects:
 Groups   Name        Variance Std.Dev. Corr
 sites    (Intercept) 8033602  2834         
          duration7d  1652488  1285     1.00
 season   (Intercept)       0     0         
          duration7d  1175980  1084      NaN
 Residual             5292365  2301         
Number of obs: 54, groups:  sites, 3; season, 3

Fixed effects:
            Estimate Std. Error       df t value Pr(>|t|)  
(Intercept) 4183.896   1695.252    2.008   2.468    0.132  
duration7d  3265.641   1155.357    3.270   2.827    0.060 .
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation of Fixed Effects:
           (Intr)
duration7d 0.520 
optimizer (nloptwrap) convergence code: 0 (OK)
boundary (singular) fit: see ?isSingular

(以上数据和代码不是我的模型，我从一个堆栈溢出问题中复制并粘贴了这些数据和代码。)

综上所述，我的问题是

1. 除了注意警告信息外，是否有明确的函数或方式通知该函数是否收敛或收敛失败

(比如，评估奇点，isSingular() 函数给出了明确的指示)

1. 为什么在模型无法收敛时仍然估计标准误差？

最终目的是为了我的模拟学习，我会计算收敛速度。

Answer 1

My advisor said, if the model is failed to converge, the standard error will not be estimated. However, although the below failed to converge, a standard error seems to be estimated.

您展示的模型已经收敛。你知道这一点是因为消息:

optimizer (nloptwrap) convergence code: 0 (OK)

如果它没有收敛，您会看到如下警告:

In checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv,  :
Model failed to converge: degenerate  Hessian with 1 negative eigenvalues

但是它已经收敛到一个奇异拟合，如下一行所示:

boundary (singular) fit: see ?isSingular

is there any clear function or way to notify whether this function is converged or failed converge, other than noticing warning message

我为此使用了以下辅助函数:

# helper function
# Has the model converged ?

hasConverged <- function (mm) {
  
  if ( !inherits(mm, "merMod")) stop("Error: must pass a lmerMod object")
  
  retval <- NULL
  
  if(is.null(unlist(mm@optinfo$conv$lme4))) {
    retval = 1
  }
  else {
    if (isSingular(mm)) {
      retval = 0
    } else {
      retval = -1
    }
  }
  return(retval)
}

如果模型正常收敛，即不收敛到奇异拟合，则返回 1；如果收敛到奇异拟合，则返回 0；如果收敛失败，则返回 -1。另一种方法是根据@SamR 的评论将警告提升为错误:

In general, if a warning is not enough, you can turn a warning into an error with options(warn=2), which means the operation will end so you should not get any standard errors or other output. Just remember to set warnings back to 1 afterwards.

继续:

Why standard error still estimated while the model is failed to converge?

嗯，如上所述，它已经收敛了，你的advisor在这里错了:

My advisor said, if the model is failed to converge, the standard error will not be estimated.

如果模型无法收敛，它将输出在放弃之前的最后一次迭代中获得的估计值。