12
$\begingroup$

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?

$\endgroup$
  • 3
    $\begingroup$ (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) $\endgroup$ – Ben Bolker Dec 15 '10 at 16:42
  • $\begingroup$ 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. $\endgroup$ – russellpierce Dec 16 '10 at 4:12
4
$\begingroup$

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

$\endgroup$
8
$\begingroup$

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 ...

$\endgroup$
  • $\begingroup$ 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. $\endgroup$ – russellpierce Dec 20 '10 at 8:26
  • 1
    $\begingroup$ 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? $\endgroup$ – russellpierce Feb 12 '13 at 15:34
  • $\begingroup$ 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. $\endgroup$ – gung Aug 24 '13 at 15:09

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.