(数据不是我的数据,而是来自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个级别实在太少了,无法很好地估计随机效应。