I have read in several answers to questions on this site that the best way to choose the random structure for a mixed effects model is by using theoretical knowledge. On the other hand I have also read the advice in Barr et al (2013) to “keep it maximal”, that is, to fit the maximal random structure by including all fixed effects, including interactions, as random slopes. However, this seems to often lead to models which either won’t converge, or do converge, but with a warning of a “singular fit”. In the accepted answer to [this] (Is a singular fit with no correlations near +/- 1 or variances of zero, a false positive?) question it is stated that singular models should be simplified.

But how is this done when the reported correlations are not near +/-1 and there is no theoretical knowledge to help choose.

An example would be much appreciated.


A good approach to this kind of problem is outlined in Bates et al (2015).

But first a bit of background. Bates et al (2015) re-analysed several sets of experimental data where a maximal random structure was adopted. In particular they have re-analysed the dataset used by Barr et al (2013) that was used as an example of “keeping it maximal” and found that the model was severely overfitted. In Barr et al (2013) the authors fit a model with crossed random effects and random slopes for 8 fixed effects across both grouping factors. This means 8 variance components and 28 correlations between them, for /each/ grouping factor, that is a total of 72 parameters. Bearing in mind that the data had only 56 subjects who responded to 32 items, common sense should suggest that such a model would be severely overfitted. Bates, rather diplomatically assessed the idea that the data would support such a complex random structurel as "optimistic" ! However the model actually did converge without warnings, using lme4 in R, although as noted by Bates this was rather "unfortunate", as they went on to show that it was indeed overfitted, and they used principal components analysis to identify this. More recent versions of lme4 actually use very same PCA procedure explained below to determine whether the model has converged with a “singular fit” and produces a warning. Very often this is also accompanied by estimated correlations between the random effects of +1 or -1, and/or variance components estimated at zero, however when the random structure is complex (typically of dimension 3 or higher) then these "symptoms" can be absent.

In lme4, a Cholesky decomposition of the variance covariance (VCV) matrix is used during estimated. If the Cholesky factor (a lower triangular matrix) contains one or more columns of zero values, then it is rank deficient, which means there is no variability in one or more of the random effects. This is equivalent to having variance components with no variability. PCA is a dimensionality reduction procedure, and when applied to the estimated VCV matrix of random effects, will immediately indicate whether this matrix is of full rank. If we can reduce the dimensionality of the VCV matrix, that is, if the number of principal components that account for 100% of the variance is less than the number of columns in the VCV matrix, then we have prima facie evidence that the random effects structure is too complex to be supported by the data and can therefore be reduced.

Thus Bates suggests the following iterative procedure:

  1. Apply PCA to the VCV matrix to determine whether the model is overfitted (singular).
  2. Fit a “zero correlation parameter” (ZCP) which will identify random effects with zero, or very small, variance
  3. Remove these random effects from the model and fit a newly reduced model and check for any other near-zero random effects. Repeat as needed.
  4. Re-introduce correlations among the remaining random effects, and if a non-singular fit is obtained use a likelihood ratio test to compare this model with the previous one. If there is still a singular fit then go back to 2. At this point it is worth noting that lme4 now incorporates step 1 above during the fitting procedure and will produce a warning that the fit is singular. In models where the random structure is simple, such as random intercepts with a single random slope it is usually obvious where the problem lies and removing the random slope will usually cure the problem. It is important to note that this does not mean that there is no random slope in the population, only that the current data do not support it.

However, things can be a little confusing when lme4 reports that the fit is singular, but there are no correlations of +/- 1 or variance components of zero. But applying the above procedure can usually result in a more parsimonious model that is not singular. A worked example can demonstrate this:

This dataset has 3 variables to be considered as fixed effects: A, B and C, and one grouping factor group with 10 levels. The response variable is Y and there are 15 observations per group.

We begin by fitting the maximal model, as suggested by Barr et al (2013).

> library(lme4)

The data can be downloaded from: https://github.com/WRobertLong/Stackexchange/blob/master/data/singular.csv

Here they are loaded into R into the dataframe dt.

> m0 <- lmer(y ~ A * B * C + (A * B * C  | group), data = dt)
boundary (singular) fit: see ?isSingular

Note that this is a singular fit. However, if we inspect the VCV matrix we find no correlations near 1 or -1, nor any variance component near zeroL

> VarCorr(m0)
 Groups   Name        Variance Std.Dev. Corr                                     
 group    (Intercept) 3.710561 1.9263                                            
          A           4.054078 2.0135    0.01                                    
          B           7.092127 2.6631   -0.01 -0.03                              
          C           4.867372 2.2062   -0.05 -0.02 -0.22                        
          A:B         0.047535 0.2180   -0.05 -0.47 -0.83 -0.03                  
          A:C         0.049629 0.2228   -0.24 -0.51  0.47 -0.74  0.01            
          B:C         0.048732 0.2208   -0.17  0.08 -0.40 -0.77  0.50  0.44      
          A:B:C       0.000569 0.0239    0.24  0.43  0.37  0.65 -0.72 -0.63 -0.86
 Residual             3.905752 1.9763                                            
Number of obs: 150, groups:  group, 10

Now we apply PCA using the rePCA function in lme4:

> summary(rePCA(m0))
Importance of components:
                        [,1]  [,2]  [,3]  [,4]    [,5]     [,6]       [,7] [,8]
Standard deviation     1.406 1.069 1.014 0.968 0.02364 0.000853 0.00000322    0
Proportion of Variance 0.389 0.225 0.202 0.184 0.00011 0.000000 0.00000000    0
Cumulative Proportion  0.389 0.613 0.816 1.000 1.00000 1.000000 1.00000000    1

This shows that the VCV matrix has 8 columns, but is rank deficent, because the first 4 principal components explain 100% of the variance. Hence the singular fit, and this means it is over-fitted and we can remove parts of the random structure.

So next we fit a "Zero-correlation-parameter" model:

> m1 <- lmer(y ~ A * B * C + (A * B * C  || group), data = dt)
boundary (singular) fit: see ?isSingular

As we can see, this is also singular, however we can immediately see that several variance components are now very near zero:

> VarCorr(m1)
 Groups   Name        Variance     Std.Dev. 
 group    (Intercept) 3.2349037958 1.7985838
 group.1  A           0.9148149412 0.9564596
 group.2  B           0.4766785339 0.6904191
 group.3  C           1.0714133159 1.0350910
 group.4  A:B         0.0000000032 0.0000565
 group.5  A:C         0.0000000229 0.0001513
 group.6  B:C         0.0013923672 0.0373144
 group.7  A:B:C       0.0000000000 0.0000000
 Residual             4.4741626418 2.1152217

These happen to be all of the interaction terms. Moreover running PCA again, we find again that 4 components are superfluous:

> summary(rePCA(m1))
Importance of components:
                         [,1]   [,2]   [,3]    [,4]    [,5]      [,6]      [,7] [,8]
Standard deviation     0.8503 0.4894 0.4522 0.32641 0.01764 7.152e-05 2.672e-05    0
Proportion of Variance 0.5676 0.1880 0.1605 0.08364 0.00024 0.000e+00 0.000e+00    0
Cumulative Proportion  0.5676 0.7556 0.9161 0.99976 1.00000 1.000e+00 1.000e+00    1 

So now we remove the interactions from the random structure:

> m2 <- lmer(y ~ A * B *  C + (A + B + C  || group), data = dt)

The model now converges without warning, and PCA shows that the VCV is of full rank:

> summary(rePCA(m2))
Importance of components:
                         [,1]    [,2]    [,3]    [,4]
Standard deviation     1.5436 0.50663 0.45275 0.35898
Proportion of Variance 0.8014 0.08633 0.06894 0.04334
Cumulative Proportion  0.8014 0.88772 0.95666 1.00000

So we now re-introduce correlations:

m3 <- lmer(y ~ A * B * C + (A + B + C  | group), data = dt)
boundary (singular) fit: see ?isSingular

...and now the fit is singular again, meaning that at least one of the correlations are not needed. We could then proceed to further models with fewer correlations, but the previous PCA indicated that 4 components were not needed, so in this instance we will settle on the model with no interactions:

Random effects:
 Groups   Name        Variance Std.Dev.
 group    (Intercept) 10.697   3.271   
 group.1  A            0.920   0.959   
 group.2  B            0.579   0.761   
 group.3  C            1.152   1.073   
 Residual              4.489   2.119   

Fixed effects:
            Estimate Std. Error t value
(Intercept) -44.2911    30.3388   -1.46
A            12.9875     2.9378    4.42
B            13.6100     3.0910    4.40
C            13.3305     3.1316    4.26
A:B          -0.3998     0.2999   -1.33
A:C          -0.2964     0.2957   -1.00
B:C          -0.3023     0.3143   -0.96
A:B:C         0.0349     0.0302    1.16

We can also observe from the fixed effects estimates that the interaction terms have quite large standard errors, so in this instance we will also remove those, producing the final model:

> m4 <- lmer(y ~ A + B + C + (A + B + C  || group), data = dt)
> summary(m4)
Random effects:
 Groups   Name        Variance Std.Dev.
 group    (Intercept) 4.794    2.189   
 group.1  A           0.794    0.891   
 group.2  B           0.553    0.744   
 group.3  C           1.131    1.064   
 Residual             4.599    2.145   
Number of obs: 150, groups:  group, 10

Fixed effects:
            Estimate Std. Error t value
(Intercept)  -14.000      1.868    -7.5
A              9.512      0.301    31.6
B             10.082      0.255    39.5
C             10.815      0.351    30.8

I would also point out that I simulated this dataset with standard deviations of 2 for the residual error and random intercept, 1 for all the random slopes, no correlations between the slopes, -10 for the fixed intercept and 10 for each of the fixed effects, and no interactions.

So in this case, we have settled upon a model that has estimated all the parameters adequately.


Bates, D., Kliegl, R., Vasishth, S. and Baayen, H., 2015. Parsimonious mixed models. arXiv preprint arXiv:1506.04967. https://arxiv.org/pdf/1506.04967.pdf

Barr, D.J., Levy, R., Scheepers, C. and Tily, H.J., 2013. Random effects structure for confirmatory hypothesis testing: Keep it maximal. Journal of memory and language, 68(3), pp.255-278.

| cite | improve this answer | |
  • $\begingroup$ @rnorouzian I think so but I'm not near a computer to check. You can check pretty easily though, by just looking at the data from 1 group $\endgroup$ – Robert Long Sep 29 at 17:41
  • $\begingroup$ Thanks, Rob, I just did, they indeed are all level-1 predictors thanks. But why when we drop the correlation bet. intercepts & slopes, the small variances will reveal themselves? (I think it has to do with cov2cor() formula but am not sure). $\endgroup$ – rnorouzian Sep 29 at 20:31
  • $\begingroup$ @rnorouzian I'm also not sure. It could make for an interesting question! $\endgroup$ – Robert Long Sep 29 at 20:48
  • $\begingroup$ Rob, one more thought that keeps coming back to me is that the correlation among random-effects can become -+1 not only when some of the variances of random-effects are nearly zero or very small, but also when Covariance bet. two random-effects nearly equals the multiplication of the std. deviations of the two random-effects. This is because of $r_{x,y} = Cov(x,y)/sd(x)*sd(y)$. In such cases, we get a singular warning. BUT, that doesn't mean that variances in question are very small or nearly zero. Does that make sense? $\endgroup$ – rnorouzian Sep 30 at 23:35
  • $\begingroup$ I should add one more important comment. I can't cant exactly reproduce your results in this thread. In fact, for m3 is not singular at all! Also, all values for estimates are a bit different. packageVersion("lme4") [1] ‘1.1.23’ $\endgroup$ – rnorouzian Oct 1 at 0:36

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.