模型在r (lme4)中收敛或不收敛

Question

(这些数据不是我的数据，而是来自堆栈溢出网站)

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)

REML=TRUE是默认的，所以我没有把它放进去。

一台计算机(哪一台更好)给我这个输出

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
Warning message:
Model failed to converge with 1 negative eigenvalue: -2.3e+01

因为这个警告信息，我认为这个模型不能收敛。然而，我的顾问说，这是不清楚的，因为特征值实际上接近0。

更令人困惑的一点是，如果我在另一台计算机上运行相同的代码，结果如下

Linear mixed model fit by REML ['lmerMod']
Formula: effect ~ duration + (1 + duration | sites) + (1 + duration |      season)
   Data: dat1
REML criterion at convergence: 968.9574
Random effects:
 Groups   Name        Std.Dev. Corr
 sites    (Intercept) 2834         
          duration7d  1285     1.00
 season   (Intercept)    0         
          duration7d  1084      NaN
 Residual             2301         
Number of obs: 54, groups:  sites, 3; season, 3
Fixed Effects:
(Intercept)   duration7d  
       4184         3266  
optimizer (nloptwrap) convergence code: 0 (OK) ; 0 optimizer warnings; 1 lme4 warnings 

在这里，没有“未能收敛”错误消息。所以，我真的很疑惑，究竟是融合了还是失败了。

此外，在我之前的问题(How can I know whether the model is converged or failed to converge in lme4 without warning message in r?)中，Robert Long给了我一个非常有用的函数来指示某个模型是否已经收敛

# 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)
}`

如果我使用这个函数，它给我0，这意味着它收敛，但奇异拟合。

但是，由于“模型未能收敛到一个负特征值：-2.3e+01”这个警告信息，我感到非常困惑。

我需要一个函数来表示某个模型是否已经收敛。但我不确定哪个元素表示模型是否收敛(@optinfo$conv$lme4 4非常可疑，但正如您在上面看到的，我非常困惑)。

(TMI:最终目标是在我的多级仿真研究中计算收敛速度)。

Answer 1

因为这个警告信息，我认为这个模型无法收敛。然而，我的顾问说，这是不清楚的，因为特征值实际上接近0。

不，模型已经收敛了。我没有从您的数据中获得Model failed to converge with 1 negative eigenvalue消息。

它已经收敛到一个奇异的拟合，这是因为模型被过度拟合。你只有3个地点和3个季节，你也是拟合随机斜率的持续时间超过两个分组变量。试一试：

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

然而，3的水平实在太少，无法很好地估计随机效应。