lmerTest::anova not showing p-values

Question

I am asking a new question because the dublicate (anova() does not display p-value when used with lmerTest) is not really providing an answer: I ran into the same problem that lmerTest::anova will not output degrees of freedom and p-values for a specific model (that is much less complicated than the one in the post mentioned above):

DirectionFit <- lmer(Similarity ~  picture_category * ComparisonType + (1 + picture_category + ComparisonType|Subject), data = DirectionData, REML=FALSE)

I noticed, that the model report includes the convergence code 0 and that there have been 2 optimizer warnings. Could it be that it is because of this that lmerTest::anova is not reporting p-values? The output of the model itself looks normal, just with very high AIC, BIC etc (-10625). Does anyone know what to do about this? I have read that maybe another optimizer should be chosen, but would that also solve the convergence problem? Any help is appreciated!

Edit: I uploaded my data on: http://www.sharecsv.com/s/1e303e9cef06d02af51ed5231385b759/TestData.csv Thank you for your help!

Answer 1

The basic problem is that you have a singular fit; the estimated random-effects variance-covariance matrix is on the boundary of its feasible space (equivalently, one of the internal parameters that lme4 uses to characterize the variance-covariance matrix, which must be non-negative, is exactly zero). This is not an optimization problem per se; it generally means your model is too complex for your data, and should be simplified (e.g. see the relevant section in the GLMM FAQ for more information). In particular, although your experimental design in principle allows you to fit among-subject variation in effects of picture_category and ComparisonType, hoping to estimate a 4x4 variance-covariance matrix for the random effect from 39 subjects is very optimistic. (You may be following Barr et al's "keep it maximal advice" ...)

What follows may be more technical than you really want, but I'm providing it as future reference for trouble-shooting these kinds of problems. If you want to ignore the fact that your model is singular (which I wouldn't recommend), and you are willing to fix a bug in the current version of lmerTest [I'm sending the maintainer an e-mail]), you can actually get p-values for this problem via Kenward-Roger approximation: summary(m2, ddf="Kenward-Roger") works, although it's quite slow (163 seconds on my laptop).

---

Because lmerTest intercepts error messages, making it hard to see what has happened, I generally try to trouble-shoot lmerTest problems by starting with lme4.

Fit model:

DirectionData <- read.csv("TestData.csv")
library(lme4)
DirectionFit <- lmer(Similarity ~  picture_category * ComparisonType +
           (1 + picture_category + ComparisonType|Subject),
                     data = DirectionData, REML=FALSE)
## Warning messages:
## 1: In checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv,  :
##   unable to evaluate scaled gradient
## 2: In checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv,  :
##   Model failed to converge: degenerate  Hessian with 1 negative eigenvalues

OK, Let's look at what happened. summary(DirectionFit), and particularly VarCorr(DirectionFit), doesn't show us any 0 variances or +/-1 correlations, but the fit is nonetheless singular:

th <- getME(DirectionFit,"theta")
which(th==0)
## Subject.picture_categoryland 
##                           8 

(it might be better in practice to test for which(th<1e-10))

So how does this affect the lmerTest output?

library(lmerTest)
m2 <- as(DirectionFit,"merModLmerTest")
trace("summary",sig="merModLmerTest",browser)  ## debug
summary(m2)

As we dig down, we find that the trouble is in the calcApvar function, which has this code:

dd <- devfunTheta(rho$model)  ## set up deviance function
h <- myhess(dd, c(rho$thopt, sigma = rho$sigma)) ## compute Hessian
ch <- try(chol(h), silent = TRUE)

If we try the last command without silencing/catching the error message we get

Error in chol.default(h) : the leading minor of order 10 is not positive definite

Basically, we can't Cholesky-decompose the Hessian (second-derivative) matrix because it's not positive definite: as you can read on Wikipedia, the Cholesky decomposition applies to PD matrices.