lmerTest::anova не показывает p-значения

Question

Я задаю новый вопрос, потому что dublicate ( anova () не отображает значение p при использовании с lmerTest) на самом деле не дает ответа: я столкнулся с той же проблемой, что lmerTest::anova не будет выводить степени свободы и p -значения для конкретной модели (которая гораздо менее сложна, чем та, что в посте, упомянутом выше):

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

Я заметил, что в модельном отчете содержится код сходимости 0 и было 2 предупреждения оптимизатора. Может ли быть так, что из-за этого lmerTest::anova не сообщает p-значения? Вывод самой модели выглядит нормально, только с очень высоким AIC, BIC и т. Д. (-10625). Кто-нибудь знает, что с этим делать? Я читал, что, возможно, следует выбрать другой оптимизатор, но это также решит проблему сходимости? Любая помощь приветствуется!

Изменить: Я загрузил свои данные на: http://www.sharecsv.com/s/1e303e9cef06d02af51ed5231385b759/TestData.csv Спасибо за вашу помощь!

Answer 1

Основная проблема заключается в том, что у вас единственная форма; оценочная матрица дисперсии-ковариации случайных эффектов находится на границе ее возможного пространства (эквивалентно, одному из внутренних параметров, которые lme4 используется для характеристики дисперсионно-ковариационной матрицы, которая должна быть неотрицательной, точно равна нулю). Это не проблема оптимизации как таковая; Обычно это означает, что ваша модель слишком сложна для ваших данных и должна быть упрощена (например, см. соответствующий раздел в FAQ по GLMM для получения дополнительной информации). В частности, хотя ваш экспериментальный дизайн в принципе позволяет вам соответствовать разным воздействиям picture_category а также ComparisonType надеяться на оценку матрицы дисперсии-ковариации 4х4 для случайного эффекта от 39 субъектов очень оптимистично. (Возможно, вы следуете совету Барра и его коллег "соблюдайте максимальные правила" ...)

То, что следует, может быть более техническим, чем вы на самом деле хотите, но я предоставляю его в качестве будущего справочника для устранения проблем такого рода проблем. Если вы хотите игнорировать тот факт, что ваша модель является единственной (что я бы не рекомендовал), и вы готовы исправить ошибку в текущей версии lmerTest [Я отправляю сопровождающему электронное письмо]), вы можете получить p-значения для этой задачи с помощью приближения Кенварда-Роджера: summary(m2, ddf="Kenward-Roger") работает, хотя и довольно медленно (163 секунды на моем ноутбуке).

---

Так как lmerTest перехватывает сообщения об ошибках, затрудняя понимание того, что произошло, я обычно пытаюсь устранить неполадки lmerTest проблемы, начиная с lme4,

Подходящая модель:

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 не является положительно определенным

По сути, мы не можем разложить по Холески матрицу Гессе (вторая производная), потому что она не является положительно определенной: как вы можете прочитать в Википедии, разложение Холецкого применимо к матрицам ПД.