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:
- Apply PCA to the VCV matrix to determine whether the model is overfitted (singular).
- Fit a “zero correlation parameter” (ZCP) which will identify random effects with zero, or very small, variance
- 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.
- 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))
$`group`
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))
$`group`
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))
$`group`
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.
References:
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.
