12

The library languageR provides a method (pvals.fnc) to do MCMC significance testing of the fixed effects in a mixed effect regression model fit using lmer. However, pvals.fnc gives an error when the lmer model includes random slopes.

Is there a way to do an MCMC hypothesis test of such models?
If so, how? (To be accepted an answer should have a worked example in R) If not, is there a conceptual/computation reason why there is no way?

This question might be related to this one but I didn't understand the content there well enough to be certain.

Edit 1: A proof of concept showing that pvals.fnc() still does 'something' with lme4 models, but that it doesn't do anything with random slope models.

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)

It says: Error in pvals.fnc(primingHeid.lmer.rs) : MCMC sampling is not yet implemented in lme4_0.999375 for models with random correlation parameters

Additional question: Is pvals.fnc performing as expected for random intercept model? Should the outputs be trusted?

Community
  • 1
russellpierce
  • 17,079
  • 16
  • 67
  • 98
  • 3
    (1) I'm surprised that pvals.fnc still works at all; I thought mcmcsamp (which pvals.fnc relies on) had been non-functional in lme4 for quite a while. What version of lme4 are you using? (2) There is no conceptual reason why having random slopes should break whatever one is doing to get a significance test (3) Combining significance testing with MCMC is a little bit incoherent, statistically, although I understand the urge to do so (getting confidence intervals is more supportable) (4) the only relationship between this Q & the other is 'MCMC' (i.e. none, practically) – Ben Bolker Dec 15 '10 at 16:42
  • The version of lme4 I use depends on the computer I am sitting at. This console has lme4_0.999375-32, but I seldom use this one for analysis. I noticed several months ago that pvals.fnc() was ripping apart the lme4 results after analysis - I built a work around for it at the time, but it doesn't seem to be an issue anymore. I'll have to ask another question on your 3rd point in the near future. – russellpierce Dec 16 '10 at 04:12

2 Answers2

8

Here's (at least most of) a solution with MCMCglmm.

First fit the equivalent intercept-variance-only model with MCMCglmm:

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

Comparing fits between MCMCglmm and lmer, first retrieving my hacked version of 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)

Now try it with random slopes:

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

This does give some sort of "MCMC p-values" ... you'll have to explore for yourself and see whether the whole thing makes sense ...

gung - Reinstate Monica
  • 132,789
  • 81
  • 357
  • 650
Ben Bolker
  • 34,308
  • 2
  • 93
  • 126
  • Thanks a lot Ben. This looks like it will point me in the right direction. I just need to spend some time reading through the help and related articles for MCMCglmm to see if I can wrap my head around what is happening. – russellpierce Dec 20 '10 at 08:26
  • 1
    Does the random slopes model in this case have a correlation between the random slope and random intercept, or in this framework is such an idea nonsensical? – russellpierce Feb 12 '13 at 15:34
  • I tweaked your code here to make it easier to read, @Ben; if you don't like it, feel free to roll it back w/ my apologies. – gung - Reinstate Monica Aug 24 '13 at 15:09
4

It looks like your error message isn't about varying slopes, it is about correlated random effects. You can fit the uncorrelated as well; that is, a mixed-effects model with independent random effects: 

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

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

gung - Reinstate Monica
  • 132,789
  • 81
  • 357
  • 650
Patrick McCann
  • 1,330
  • 7
  • 12