duplicated
For example,
I can assess whether this multilevel model is a singular fit or not by
using isSingular() function.
Likewise, is there any way that I can know whether this model is converged or failed to be converged?
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.
is there any good indicator that this model is converged or failed to converge? (other than noticing a warning message)
I am using the lme4 package and lmer() function.
For example, there is an example of failed convergence multilevel model
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)
This generates error Warning message: Model failed to converge with 1 negative eigenvalue: -2.3e+01
however, it seems that standard error seems to be estimated, although it is failed to converge.
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
(the above data and code is Not my model, I copied and pasted this data and code from one of stack overflow question.)
To sum up, my question is
(like, assessing singularity, isSingular() function gives clear indication)
the ultimate goal is for my simulation study, I will calculate the convergence rate.
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.
The model you showed has converged. You know this because of the message:
optimizer (nloptwrap) convergence code: 0 (OK)
If it had not converged you would see a warning like:
In checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
Model failed to converge: degenerate Hessian with 1 negative eigenvalues
However it has converged to a singular fit as shown in the next line:
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
I use the following helper function for that:
# 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)
}
which returns 1 if the model converged normally ie not to a singular fit, 0 if it converges to a singular fit and -1 if it fails to converge. Another approach is to promote the warnings to errors as per the comment by @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.
Moving on:
Why standard error still estimated while the model is failed to converge?
Well, as mentioned above, it has converged, and your advisor is wrong here:
My advisor said, if the model is failed to converge, the standard error will not be estimated.
If the model fails to converge it will output the estimates obtained on the last iteration before it gave up.