Как можно проверить гипотезу MCMC на регрессионной модели со смешанным эффектом со случайными наклонами?

12

Библиотека languageR предоставляет метод (pvals.fnc) для проверки значимости MCMC фиксированных эффектов в модели регрессии смешанного эффекта, подходящей с использованием lmer. Однако pvals.fnc выдает ошибку, когда модель lmer включает случайные наклоны.

Есть ли способ сделать проверку гипотезы MCMC таких моделей?
Если так, то как? (Чтобы быть принятым, ответ должен иметь проработанный пример в R). Если нет, есть ли концептуальная / вычислительная причина, почему нет пути?

Этот вопрос может быть связан с этим, но я не достаточно хорошо понял содержание там, чтобы быть уверенным.

Редактирование 1 : подтверждение концепции, показывающее, что pvals.fnc () все еще «что-то» делает с моделями lme4, но ничего не делает с моделями со случайным наклоном.

library(lme4)
library(languageR)
#the example from pvals.fnc
data(primingHeid) 
# remove extreme outliers
primingHeid = primingHeid[primingHeid$RT < 7.1,]
# fit mixed-effects model
primingHeid.lmer = lmer(RT ~ RTtoPrime * ResponseToPrime + Condition + (1|Subject) + (1|Word), data = primingHeid)
mcmc = pvals.fnc(primingHeid.lmer, nsim=10000, withMCMC=TRUE)
#Subjects are in both conditions...
table(primingHeid$Subject,primingHeid$Condition)
#So I can fit a model that has a random slope of condition by participant
primingHeid.lmer.rs = lmer(RT ~ RTtoPrime * ResponseToPrime + Condition + (1+Condition|Subject) + (1|Word), data = primingHeid)
#However pvals.fnc fails here...
mcmc.rs = pvals.fnc(primingHeid.lmer.rs)

Это говорит: Ошибка в pvals.fnc (primingHeid.lmer.rs): выборка MCMC еще не реализована в lme4_0.999375 для моделей со случайными параметрами корреляции

Дополнительный вопрос: работает ли pvals.fnc, как ожидается, для модели случайного перехвата? Следует ли доверять выходам?

russellpierce
источник
3
(1) Я удивлен, что pvals.fnc все еще работает; Я думал, что mcmcsamp (на который опирается pvals.fnc) долгое время не работал в lme4. Какую версию lme4 вы используете? (2) Нет концептуальной причины, по которой наличие случайных уклонов должно нарушать то, что вы делаете, чтобы получить тест на значимость. (3) Комбинация тестирования на значимость с MCMC немного статистически непоследовательна, хотя я понимаю желание сделать это (получить уверенность интервалы более приемлемы) (4) единственное отношение между этим вопросом и другим является «MCMC» (т.е. практически нет)
Бен Болкер
Используемая версия lme4 зависит от компьютера, на котором я сижу. Эта консоль имеет lme4_0.999375-32, но я редко использую ее для анализа. Несколько месяцев назад я заметил, что pvals.fnc () разбирает результаты lme4 после анализа - в то время я обходил его, но, похоже, это больше не проблема. Мне придется задать еще один вопрос по вашему третьему пункту в ближайшем будущем.
Расселпирс

Ответы:

4

Похоже, ваше сообщение об ошибке не об изменяющихся наклонах, а о связанных случайных эффектах. Вы можете соответствовать и некоррелированным; модель смешанных эффектов с независимыми случайными эффектами:

Linear mixed model fit by REML
Formula: Reaction ~ Days + (1 | Subject) + (0 + Days | Subject)
Data: sleepstudy

с http://www.stat.wisc.edu/~bates/IMPS2008/lme4D.pdf

Патрик Макканн
источник
8

Вот (по крайней мере, большинство) решение с MCMCglmm.

Сначала подберите эквивалентную модель только для перехвата MCMCglmm:

library(MCMCglmm)
primingHeid.MCMCglmm = MCMCglmm(fixed=RT ~ RTtoPrime * ResponseToPrime + Condition, 
                                random=~Subject+Word, data = primingHeid)

Сравнивая совпадения между MCMCglmmи lmer, сначала извлекая мою взломанную версию arm::coefplot:

source(url("http://www.math.mcmaster.ca/bolker/R/misc/coefplot_new.R"))
  ## combine estimates of fixed effects and variance components
pp  <- as.mcmc(with(primingHeid.MCMCglmm, cbind(Sol, VCV)))
  ## extract coefficient table
cc1 <- coeftab(primingHeid.MCMCglmm,ptype=c("fixef", "vcov"))
  ## strip fixed/vcov indicators to make names match with lmer output
rownames(cc1) <- gsub("(Sol|VCV).", "", rownames(cc1))
  ## fixed effects -- v. similar
coefplot(list(cc1[1:5,], primingHeid.lmer))
  ## variance components -- quite different.  Worth further exploration?
coefplot(list(cc1[6:8,], coeftab(primingHeid.lmer, ptype="vcov")),
         xlim=c(0,0.16), cex.pts=1.5)

Теперь попробуйте это со случайными наклонами:

primingHeid.rs.MCMCglmm = MCMCglmm(fixed=RT ~ RTtoPrime * ResponseToPrime + Condition,
                                   random=~Subject+Subject:Condition+Word, 
                                   data = primingHeid)        
summary(primingHeid.rs.MCMCglmm)

Это дает какие-то "значения p MCMC" ... вам придется изучить для себя и посмотреть, имеет ли смысл все это ...

Бен Болкер
источник
Большое спасибо, Бен. Похоже, это укажет мне правильное направление. Мне просто нужно потратить некоторое время на чтение справки и связанных статей для MCMCglmm, чтобы посмотреть, смогу ли я обернуться вокруг происходящего.
Расселпирс
1
Имеет ли модель случайных наклонов в этом случае корреляцию между случайным наклоном и случайным пересечением, или в этой структуре такая идея бессмысленна?
russellpierce
Я настроил ваш код здесь, чтобы его было легче читать, @Ben; если вам это не нравится, не стесняйтесь откатить его назад с моими извинениями.
gung - Восстановить Монику