"Model failed to converge" warning in lmer()

Question

With the following dataset, I wanted to see if the response (effect) changes with regard to sites, season, duration, and their interactions. Some online forums on statistics suggested me to go on with Linear Mixed-Effects Models, but the problem is that since replicates are randomised within each station, I have little chance to collect the sample from exactly the same spot in successive seasons (for example, repl-1 of s1 of post-monsoon may not be the same as that of monsoon). It is unlike the clinical trials (with within-subject design) where you measure the same subject repeatedly over seasons. However, considering sites and season as a random factor I ran the following commands and received a warning message:

Warning messages:
1: In checkConv(attr(opt, "derivs"), optpar,ctrl=controlpar,ctrl=controlcheckConv, 
: unable to evaluate scaled gradient
2: In checkConv(attr(opt, "derivs"), optpar,ctrl=controlpar,ctrl=controlcheckConv, 
: Model failed to converge: degenerate Hessian with 1 negative eigenvalues

Can anyone help me solve the issue? The codes are given below:

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


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

Answer 1

"Solving" the issue you experience in the sense of not receiving warnings about failed convergence is rather straightforward: you do not use the default BOBYQA optimiser but instead you opt to use the Nelder-Mead optimisation routine used by default in earlier 1.0.x previous versions. Or you install the package optimx so you can directly an L-BFGS-B routine or nlminb (same as lme4 versions prior to ver. 1). For example:

m1 = lmer(effect~duration+(1+duration|sites)+(1+duration|season), 
          data = dat1, REML = FALSE, 
          control = lmerControl(optimizer ="Nelder_Mead")
library(optimx)
m1 = lmer(effect~duration+(1+duration|sites)+(1+duration|season), 
          data = dat1, REML = FALSE, 
          control = lmerControl(
                           optimizer ='optimx', optCtrl=list(method='L-BFGS-B')))
m1 = lmer(effect~duration+(1+duration|sites)+(1+duration|season), 
          data = dat1, REML = FALSE, 
          control = lmerControl(
                           optimizer ='optimx', optCtrl=list(method='nlminb')))

all work fine (no warnings). The interesting questions are:

1. why you got these warnings to begin with and
2. why when you used REML = TRUE you got no warnings.

Succinctly, 1. you received those warnings because you defined duration both as a fixed effect as well as random slope for the factor sites as well as season. The model effectively ran-out of the degrees of freedom to estimate the correlations between the slopes and the intercepts you defined. If you used a marginally simpler model like:

m1 = lmer(effect~duration+ (1+duration|sites) + (0+duration|season) + (1|season),
          data=dat1, REML = FALSE)

you would experience no convergence issues. This model would effectively estimate uncorrelated random intercepts and random slopes for each season.

In addition, 2. when you defined REML = FALSE you used the Maximum Likelihood estimated instead of the Restricted Maximum Likelihood one. The REML estimates try to "factor out" the influence of the fixed effects $X$ before moving into finding the optimal random-effect variance structure (see the thread "What is "restricted maximum likelihood" and when should it be used?" for more detailed information on the matter). Computationally this procedure is essentially done by multiplying both parts of the original LME model equation $y = X\beta + Z\gamma + \epsilon$ by a matrix $K$ such that $KX = 0$, i.e. you change both the original $y$ to $Ky$ as well as the $Z$ to $KZ$. I strongly suspect that this effected the condition number of the design matrix $Z$ and as such help you out of the numerical hard-place you found yourself in the first place.

A final note is that I am not sure whether it makes sense to use season as a random effect to begin with. After all there are only so many seasons so you might as well treat them as fixed effects.

Answer 2

The question is statistical rather than technical. Actually, I used a random effect model instead of a fixed effect one.None of the factors, I think, should be treated as the random factor as we need at least 5 or 6 levels or replicates to treat a factor as random effect (see here What is the minimum recommended number of groups for a random effects factor?).

The above dataset contains only triplicate samples/site/season which is insufficient for a random effect model.In the data set the duration, 4-day and 7-day belong to two separate parallel experiments run under the same time. So spiting the data set by duration (4-day and 7-day) and performing a 2-way anova for each duration with season and sites as the factors would be sufficient to model the effect (response variable) here. The model should be the following:

lm(day_4_effect~sites*season, data=dat1)

lm(day_7_effect~sites*season, data=dat1)

Thanks to Bodo Winter (http://www.bodowinter.com/tutorial/bw_LME_tutorial2.pdf) and @usεr11852 who helped me solve the issue.

Answer 3

I've had success with the change of optimizer suggested here in larger samples, but this fix does not always work in relatively small samples in my experience. Matuschek et al., 2017 (Journal of Memory and Language) note that failed model convergence can also result from over specification of the random effects structure when the data cannot support maximal or near-maximal models. In this case, you will need to simplify your random effects structure, which Matuschek et al., 2017 offer some suggestions on. At least assuming I've interpreted their work correctly, mixed effects models are confusing.