Апостериорный тест на блеск

Question

Я анализирую свой биномиальный набор данных с помощью R с использованием обобщенной линейной смешанной модели (glmer, lme4-package). Я хотел сделать парные сравнения определенного фиксированного эффекта ("Звук") с помощью специального теста Тьюки (glht, multcomp-package).

Большинство из них работает нормально, но одна из моих фиксированных переменных эффектов ("SoundC") вообще не имеет дисперсии (96 раз "1" и ноль умножается на "0"), и кажется, что тест Тьюки не может справиться с этим. Все парные сравнения с этим "SoundC" дают p-значение 1.000, тогда как некоторые явно значимы.

В качестве проверки я изменил одно из 96 "1" на "0", и после этого я снова получил нормальные p-значения и существенные различия там, где я ожидал их, тогда как разница фактически стала меньше после моего ручного изменения.

У кого-нибудь есть решение? Если нет, то можно ли использовать результаты моего измененного набора данных и сообщать о моих изменениях вручную?

Воспроизводимый пример:

Response <- c(1,0,1,1,0,1,1,0,1,1,0,1,1,0,1,1,1,1,0,
              0,1,1,0,1,1,0,1,1,0,1,1,0,1,1,0,1,1,0,1,1,0,
              1,1,0,1,1,0,1,1,1,1,0,0,1,1,0,1,1,0,1,1,0,1,1,0,1)    
Data <- data.frame(Sound=rep(paste0('Sound',c('A','B','C')),22),
                   Response,
                   Individual=rep(rep(c('A','B'),2),rep(c(18,15),2)))


# Visual
boxplot(Response ~ Sound,Data)

# Mixed model
library (lme4)
model10 <- glmer(Response~Sound + (1|Individual), Data, family=binomial)

# Post-hoc analysis
library (multcomp)
summary(glht(model10, mcp(Sound="Tukey")))

Answer 1

Это граничит с вопросом CrossValidated; вы определенно видите полное разделение, где ваш ответ идеально разделен на 0 против 1 результатов. Это приводит к (1) бесконечным значениям параметров (они перечислены только как бесконечные из-за несовершенства вычислений) и (2) сумасшедшим / бесполезным значениям стандартных ошибок Вальда и соответствующим значениям $p$ (что вам и нужно видишь здесь). Обсуждение и решения приведены здесь, здесь и здесь, но я проиллюстрирую немного больше ниже.

Чтобы быть статистическим занудством на мгновение: вы действительно не должны пытаться установить случайный эффект только с 3 уровнями (см., Например, http://glmm.wikidot.com/faq)...

Ферт-скорректированная логистическая регрессия:

library(logistf)
L1 <- logistf(Response~Sound*Individual,data=Data,
        contrasts.arg=list(Sound="contr.treatment",
         Individual="contr.sum"))

                                 coef se(coef)            p
(Intercept)              3.218876e+00 1.501111 2.051613e-04 
SoundSoundB             -4.653960e+00 1.670282 1.736123e-05 
SoundSoundC             -1.753527e-15 2.122891 1.000000e+00 
IndividualB             -1.995100e+00 1.680103 1.516838e-01 
SoundSoundB:IndividualB  3.856625e-01 2.379919 8.657348e-01 
SoundSoundC:IndividualB  1.820747e+00 2.716770 4.824847e-01

Стандартные ошибки и p-значения теперь являются разумными (p-значение для сравнения A против C равно 1, потому что в буквальном смысле нет никакой разницы...)

Смешанная байесовская модель со слабыми априорами:

library(blme)
model20 <- bglmer(Response~Sound + (1|Individual), Data, family=binomial,
                  fixef.prior = normal(cov = diag(9,3)))

##              Estimate Std. Error    z value     Pr(>|z|)
## (Intercept)  1.711485   2.233667  0.7662221 4.435441e-01
## SoundSoundB -5.088002   1.248969 -4.0737620 4.625976e-05
## SoundSoundC  2.453988   1.701674  1.4421024 1.492735e-01

Спецификация diag(9,3) дисперсионно-ковариационной матрицы с фиксированным эффектом производит

$$ \left( \begin{array}{ccc} 9 & 0 & 0 \ 0 & 9 & 0 \ 0 & 0 & 9 \end{array} \right) $$

Другими словами, 3 указывает размерность матрицы (равную количеству параметров с фиксированным эффектом), а 9 определяет дисперсию - это соответствует стандартному отклонению 3 или диапазону 95% около $\pm. 6$, что является довольно большим / слабым / неинформативным для масштабированных ответов.

Это примерно соответствует (модель очень отличается)

library(multcomp)
summary(glht(model20, mcp(Sound="Tukey")))

##                     Estimate Std. Error z value Pr(>|z|)    
## SoundB - SoundA == 0   -5.088      1.249  -4.074 0.000124 ***
## SoundC - SoundA == 0    2.454      1.702   1.442 0.309216    
## SoundC - SoundB == 0    7.542      1.997   3.776 0.000397 ***

Как я уже сказал выше, я бы не рекомендовал смешанную модель в этом случае...