Como se pode fazer um teste de hipótese do MCMC em um modelo de regressão de efeito misto com inclinações aleatórias?

12

A biblioteca languageR fornece um método (pvals.fnc) para realizar testes de significância do MCMC dos efeitos fixos em um modelo de regressão de efeito misto ajustado usando o lmer. No entanto, pvals.fnc comete um erro quando o modelo lmer inclui inclinações aleatórias.

Existe uma maneira de fazer um teste de hipótese do MCMC para esses modelos?
Se sim, como? (Para ser aceita, uma resposta deve ter um exemplo elaborado em R). Caso contrário, existe uma razão conceitual / computacional para que não haja como?

Esta pergunta pode estar relacionada a essa, mas eu não entendi o conteúdo o suficiente para ter certeza.

Edit 1 : Uma prova de conceito mostrando que pvals.fnc () ainda faz 'algo' com os modelos lme4, mas que não faz nada com os modelos de inclinação aleatória.

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)

Ele diz: Erro no pvals.fnc (primingHeid.lmer.rs): a amostragem MCMC ainda não foi implementada no lme4_0.999375 para modelos com parâmetros de correlação aleatória

Pergunta adicional: O pvals.fnc está executando como esperado para o modelo de interceptação aleatória? As saídas devem ser confiáveis?

russellpierce
fonte
3
(1) Estou surpreso que o pvals.fnc ainda funcione; Eu pensei que o mcmcsamp (no qual o pvals.fnc se baseia) não funcionou no lme4 por um bom tempo. Qual versão do lme4 você está usando? (2) Não há razão conceitual para que inclinações aleatórias quebrem o que está sendo feito para obter um teste de significância. (3) Combinar o teste de significância com o MCMC é um pouco incoerente, estatisticamente, embora eu compreenda o desejo de fazê-lo (obter confiança). intervalos é mais suportável) (4) a única relação entre este Q & 'MCMC' o outro é (nenhum ou seja, praticamente)
Ben Bolker
A versão do lme4 que eu uso depende do computador em que estou sentado. Esse console tem lme4_0.999375-32, mas raramente uso esse para análise. Há alguns meses, notei que o pvals.fnc () estava destruindo os resultados do lme4 após a análise - criei uma solução para isso na época, mas isso não parece mais ser um problema. Vou ter que fazer outra pergunta sobre o seu terceiro ponto no futuro próximo.
22410 russellpierce

Respostas:

4

Parece que sua mensagem de erro não é sobre variações de inclinação, é sobre efeitos aleatórios correlacionados. Você também pode ajustar o não correlacionado; isto é, um modelo de efeitos mistos com efeitos aleatórios independentes:

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

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

Patrick McCann
fonte
8

Aqui está (pelo menos a maioria) uma solução com MCMCglmm.

Primeiro, ajuste o modelo equivalente apenas de variação de interceptação com MCMCglmm:

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

Comparando ajustes entre MCMCglmme lmer, primeiro recuperando minha versão hackeada de 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)

Agora tente com inclinações aleatórias:

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

Isso fornece algum tipo de "valores-p do MCMC" ... você terá que explorar por si mesmo e ver se a coisa toda faz sentido ...

Ben Bolker
fonte
Muito obrigado Ben. Parece que isso me indicará a direção certa. Eu só preciso gastar algum tempo lendo a ajuda e os artigos relacionados do MCMCglmm para ver se consigo entender o que está acontecendo.
russellpierce
1
O modelo de declives aleatórios, neste caso, tem uma correlação entre o declive aleatório e a interceptação aleatória, ou, nesse contexto, essa idéia é absurda?
22813 russellpierce
Ajustei seu código aqui para facilitar a leitura, @Ben; se você não gostar, sinta-se à vontade para revertê-lo com minhas desculpas.
gung - Restabelece Monica