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



(数据不是我的数据,而是来自stack overflow网站)

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 

在这里,没有"未能收敛"。错误消息。所以,我真的很困惑这是收敛的还是不收敛的

补充说,在我之前的问题中(我如何知道模型在lme4中是否收敛或未能收敛,而在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,这意味着它是收敛的,但是奇异拟合。

但同样,由于"模型未能收敛于1个负特征值:-2.3e+01"这个警告信息,我很困惑。

我需要一个函数来表示某个模型是否收敛。但我不确定哪个元素表明模型是否收敛(@optinfo$conv$lme4是高度可疑的,但正如你在上面看到的,我很困惑)

(TMI:最终目标是在我的多级模拟研究中计算收敛率)。

由于这个警告消息,我认为这个模型没有收敛。然而,我的导师说这是不清楚的,因为特征值确实接近于0。

不,模型已经收敛。我没有收到Model failed to converge with 1 negative eigenvalue消息与您的数据。

它已经收敛到一个奇异拟合,这是因为模型是过拟合的。你只有3个地点和3个季节,而且你还在两个分组变量上拟合持续时间的随机斜率。试试这个:

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

然而,3个级别实在太少了,无法很好地估计随机效应。

相关内容

  • 没有找到相关文章

最新更新