lmerTest::anova没有显示p值

Question

我提出了一个新的问题，因为二元结构(anova()在与lmerTest一起使用时不显示p值。)并没有提供真正的答案:我遇到了与lmerTest::anova不会输出特定模型的自由度和p-值的问题(这比上面提到的模型要复杂得多)：

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

我注意到，模型报告包含收敛代码0，并且有两个优化器警告。难道正是因为这个，lmerTest::anova没有报告p值吗？模型本身的输出看起来正常，只是非常高的AIC，BIC等(-10625)。有人知道该怎么做吗？我已经读到，也许应该选择另一个优化器，但这是否也解决了收敛问题？任何帮助都是非常感谢的！

编辑：，我把我的数据上传到：http://www.sharecsv.com/s/1e303e9cef06d02af51ed5231385b759/TestData.csv，谢谢你的帮助！

Answer 1

基本问题是，您有一个奇异拟合；估计的随机效应方差协方差矩阵在其可行空间的边界上(等价地，lme4用来描述方差协方差矩阵的内部参数之一，必须是非负的，精确为零)。这本身并不是一个优化问题；它通常意味着您的模型对于您的数据来说太复杂了，应该简化(例如，请参阅GLMM常见问题中的相关部分以获得更多信息)。特别是，虽然你的实验设计在原则上允许你适应picture_category和ComparisonType效应中的主体变化，但希望为39个受试者的随机效应估计一个4x4方差协方差矩阵是非常乐观的。(你可能遵循Barr等人的“保持最大的建议”.)

下面的内容可能比你真正想要的更技术性，但我提供它作为将来解决这类问题的参考。如果您想忽略您的模型是单数的(我不推荐)，并且您愿意修复当前版本的lmerTest (我正在发送给维护人员的电子邮件)中的一个bug，那么您实际上可以通过Kenward-Roger近似获得这个问题的p-值：summary(m2, ddf="Kenward-Roger")工作，尽管速度很慢(在我的笔记本上有163秒)。

因为lmerTest拦截错误消息，因此很难看到发生了什么，所以我通常尝试从lme4开始来解决lmerTest问题。

Fit模型：

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

好吧，让我们看看发生了什么。summary(DirectionFit)，特别是VarCorr(DirectionFit)，并没有显示任何0的方差或+/-1的相关性，但拟合仍然是单数的：

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

(在实践中测试which(th<1e-10)可能更好)

那么这对lmerTest输出有什么影响呢？

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

当我们深入挖掘时，我们发现问题在calcApvar函数中，它有以下代码：

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)

如果我们尝试最后一个命令而没有沉默/捕获错误消息，我们将得到

chol.default(h)中的错误: 10阶前导小调不是正定的

基本上，我们不能对Hessian (二阶导数)矩阵进行Cholesky分解，因为它不是正定的:正如你可以在维基百科上看到的那样，Cholesky分解适用于PD矩阵。