How do I interpret the results of lme4::rePCA towards diagnosing a singular mixed model?

Question

I'm having a difficult time diagnosing the reasons for singular fits in a business problem I'm working on.

The lme4::isSingular documentation recommends lme4::rePCA as a function to help diagnose singular fits. I can successfully run the function on my models, but I don't know how to interpret the results.

I created a reproducible example and included below. Some specific questions:

1. What is the meaning of Subject standard deviations? How do I know what components 1, 2, 3 of this vector correspond to in the model?

#> $Subject
#> Standard deviations (1, .., p=3):
#> [1] 9.8639459 0.2313664 0.0000000

2. What do the dimensions n, k correspond to in the PCA. Perhaps one of n, k corresponds to model terms and the other to principal components?

#> Rotation (n x k) = (3 x 3):
#>      [,1] [,2] [,3]
#> [1,]    1    0    0
#> [2,]    0    0   -1
#> [3,]    0   -1    0

3. How do I interpret those PCA rotation results to understand why this model is singular?
4. Are there any tools I can use for diagnosing singular lmer and glmer models?

Many thanks!

library(Matrix)
library(lme4)

model <- lmer(Reaction ~ 0 + Days + (1 | Subject) + (Days | Subject), data = sleepstudy)
#> boundary (singular) fit: see ?isSingular

print(model)
#> Linear mixed model fit by REML ['lmerMod']
#> Formula: Reaction ~ 0 + Days + (1 | Subject) + (Days | Subject)
#>    Data: sleepstudy
#> REML criterion at convergence: 1828.653
#> Random effects:
#>  Groups    Name        Std.Dev. Corr
#>  Subject   (Intercept) 252.446      
#>  Subject.1 (Intercept)   0.000      
#>            Days          5.921   NaN
#>  Residual               25.593      
#> Number of obs: 180, groups:  Subject, 18
#> Fixed Effects:
#>  Days  
#> 10.61  
#> optimizer (nloptwrap) convergence code: 0 (OK) ; 0 optimizer warnings; 1 lme4 warnings

print(isSingular(model))
#> [1] TRUE

print(rePCA(model))
#> $Subject
#> Standard deviations (1, .., p=3):
#> [1] 9.8639459 0.2313664 0.0000000
#> 
#> Rotation (n x k) = (3 x 3):
#>      [,1] [,2] [,3]
#> [1,]    1    0    0
#> [2,]    0    0   -1
#> [3,]    0   -1    0
#> 
#> attr(,"class")
#> [1] "prcomplist"

^{Created on 2022-04-22 by the reprex package (v2.0.1)}

Answer 1

I'll start by explaining what's wrong with this model, and then try to piece together what rePCA is doing.

If you look at the random effects component of the summary() you can see that there are two Subject/Intercept terms listed (one is estimated as zero, which also leads to the estimated correlation being NaN):

Random effects:
 Groups    Name        Variance Std.Dev. Corr
 Subject   (Intercept) 63729.17 252.446      
 Subject.1 (Intercept)     0.00   0.000      
           Days           35.06   5.921   NaN
 Residual                654.99  25.593      

The problem here is that the (Days|Subject) term automatically includes an intercept; that's why this type of model is usually written (1|Subject) + (0 + Day | Subject): see also here, here ...

OK, but what is rePCA doing? (In this case we don't need rePCA to figure it out, but in more complex situations it might be useful ...)

What order are the random effects terms coded in? You can tell this from the random effects report above (the first random effects term is the (1|Subject) term, which is a scalar random effect, corresponding to a single variance (a 1x1 matrix); the second is (Days|Subject), which is a vector-valued random effect (i.e., intercept and slope), corresponding to a 2x2 RE covariance matrix.

We could also get this information as follows:

getME(model, "cnms")
$Subject
[1] "(Intercept)"

$Subject
[1] "(Intercept)" "Days"       

rePCA combines all of the random-effect covariance matrices into a single block-diagonal matrix. I'm going to use some slightly tricky code to name the components (this should be built into rePCA but needs some work to make sure it handles more complex cases):

rr <- rePCA(model)
cc <- getME(model, "cnms")
nm <- unlist(cc)
dimnames(rr$Subject[[2]]) <- list(nm, nm)

Standard deviations (1, .., p=3):
  9.8639459   0.2313664   0.0000000 

Rotation (n x k) = (3 x 3):
            (Intercept) (Intercept) Days
(Intercept)           1           0    0
(Intercept)           0           0   -1
Days                  0          -1    0

The first thing this tells us is (again) that the model is singular (because the third singular value ['standard deviations'] is exactly zero). The second is that the problem arises in the second component (intercept 2 + Days).

I might dig around a little bit more with ?svd (singular value decomposition) and with the Bates et al. "Parsimonious mixed models" preprint (on ArXiv and available from the rePsychling github repo; most of the examples primarily use/interpret the standard deviations vector (to see how many components of the covariance matrix are exactly zero) and don't dig into the structure of the rotation matrix ...

---

Bates, Douglas, Reinhold Kliegl, Shravan Vasishth, and Harald Baayen. “Parsimonious Mixed Models.” ArXiv:1506.04967 [Stat], June 16, 2015. http://arxiv.org/abs/1506.04967.