(the data is not my data, but from stack overflow website)
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 is default, so I did not put that.
one computer (which one is better one) give me this output
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
Because of this warning message, I thought this model failed to converge. However, my advisor said this is unclear because that eigenvalue is really close to 0.
A more confusing point is that, if I ran the same code on a different computer, the results are like this
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
In here, there is no "failed to converge" error message. So, I am really confused whether this is converged or failed to converge.
Added to that, in my previous question (How can I know whether the model is converged or failed to converge in lme4 without warning message in r?) @Robert Long gave me really helpful function to indicate whether certain model has converged or not
# 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)
}`
if I use this function it gives me 0, which means that it converges , but singular fit.
But again, due to "Model failed to converge with 1 negative eigenvalue: -2.3e+01" this warning message, I am so confused.
I need a function to indicate whether the certain model has converged or not. but I am not sure which element indicate whether the model converges or not (@optinfo$conv$lme4 was highly suspicious, but as you can see above, I am so confused)
(TMI: the ultimate goal is to calculate the convergence rate in my multilevel simulation study).
Because of this warning message, I thought this model failed to converge. However, my advisor said this is unclear because that eigenvalue is really close to 0.
No, the model has converged. I do not get the Model failed to converge with 1 negative eigenvalue message with your data.
It has converged to a singular fit, and this is because the model is overfitted. You have only 3 sites and 3 seasons and you are also fitting random slopes for duration over both grouping variables. Try this instead:
lmer(effect ~ duration + (1|sites) +(1|season), data = dat1)
However, 3 is really too few levels to get a good estimate of the random effects.