Skip to main content
tidied up code for better readability
Source Link
gung - Reinstate Monica
  • 147.5k
  • 89
  • 406
  • 717

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

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

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

Source Link
Ben Bolker
  • 47.3k
  • 3
  • 130
  • 182

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