inst/doc/compare_lme4.R

In BayesFactor: Computation of Bayes Factors for Common Designs

## ----echo=FALSE,message=FALSE,results='hide'----------------------------------
library(BayesFactor)

options(BFprogress = FALSE)
bfversion = BFInfo()
session = sessionInfo()[[1]]
rversion = paste(session$version.string," on ",session$platform,sep="")


options(markdown.HTML.stylesheet = 'extra/manual.css')
library(knitr)
opts_chunk$set(dpi = 200, out.width = "67%") 
options(digits=3)
require(graphics)
set.seed(2)

## ----message=FALSE,warning=FALSE----------------------------------------------
library(arm)
library(lme4)


## -----------------------------------------------------------------------------
# Number of participants
N <- 20
sig2 <- 1
sig2ID <- 1

# 3x3x3 design, with participant as random factor
effects <- expand.grid(A = c("A1","A2","A3"),
                       B = c("B1","B2","B3"),
                       C = c("C1","C2","C3"),
                       ID = paste("Sub",1:N,sep="")
)
Xdata <- model.matrix(~ A*B*C + ID, data=effects)
beta <- matrix(c(50,
          -.2,.2,
          0,0,
          .1,-.1,
          rnorm(N-1,0,sqrt(sig2ID)),
          0,0,0,0,
          -.1,.1,.1,-.1,
          0,0,0,0,
          0,0,0,0,0,0,0,0),
               ncol=1)
effects$y = rnorm(Xdata%*%beta,Xdata%*%beta,sqrt(sig2))

## -----------------------------------------------------------------------------
# Typical repeated measures ANOVA
summary(fullaov <- aov(y ~ A*B*C + Error(ID/(A*B*C)),data=effects))

## ----fig.width=10,fig.height=4------------------------------------------------
mns <- tapply(effects$y,list(effects$A,effects$B,effects$C),mean)
stderr = sqrt((sum(resid(fullaov[[3]])^2)/fullaov[[3]]$df.resid)/N)

par(mfrow=c(1,3),cex=1.1)
for(i in 1:3){
  matplot(mns[,,i],xaxt='n',typ='b',xlab="A",main=paste("C",i), 
          ylim=range(mns)+c(-1,1)*stderr,ylab="y")
  axis(1,at=1:3,lab=1:3)
  segments(1:3 + mns[,,i]*0,mns[,,i] + stderr,1:3 + mns[,,i]*0,mns[,,i] - stderr,col=rgb(0,0,0,.3))
}

## -----------------------------------------------------------------------------
t.is = system.time(bfs.is <- anovaBF(y ~ A*B*C + ID, data = effects, 
                                     whichRandom="ID")
)
t.la = system.time(bfs.la <- anovaBF(y ~ A*B*C + ID, data = effects, 
                                     whichRandom="ID",
                                     method = "laplace")
)

## ----fig.width=6,fig.height=6-------------------------------------------------
t.is
t.la

plot(log(extractBF(sort(bfs.is))$bf),log(extractBF(sort(bfs.la))$bf),
     xlab="Default Sampler",ylab="Laplace approximation",
     pch=21,bg=rgb(0,0,1,.2),col="black",asp=TRUE,cex=1.2)
abline(0,1)

bfs.is

## ----message=FALSE------------------------------------------------------------
chains <- lmBF(y ~ A + B + C + ID, data=effects, whichRandom = "ID", posterior=TRUE, iterations=10000)

lmerObj <- lmer(y ~ A + B + C + (1|ID), data=effects)
# Use arm function sim() to sample from posterior
chainsLmer = sim(lmerObj,n.sims=10000)

## -----------------------------------------------------------------------------
BF.sig2 <- chains[,colnames(chains)=="sig2"]
AG.sig2 <- (chainsLmer@sigma)^2
qqplot(log(BF.sig2),log(AG.sig2),pch=21,bg=rgb(0,0,1,.2),
       col=NULL,asp=TRUE,cex=1,xlab="BayesFactor samples",
       ylab="arm samples",main="Posterior samples of\nerror variance")
abline(0,1)

## -----------------------------------------------------------------------------
AG.raneff <- chainsLmer@ranef$ID[,,1]
BF.raneff <-  chains[,grep('ID-',colnames(chains),fixed='TRUE')]
plot(colMeans(BF.raneff),colMeans(AG.raneff),pch=21,bg=rgb(0,0,1,.2),col="black",asp=TRUE,cex=1.2,xlab="BayesFactor estimate",ylab="arm estimate",main="Random effect posterior means")
abline(0,1)


## ----tidy=FALSE---------------------------------------------------------------
AG.fixeff <- chainsLmer@fixef
BF.fixeff <-  chains[,1:10]

# Adjust AG results from reference cell to sum to 0
Z = c(1,  1/3,  1/3,  1/3,  1/3,  1/3,  1/3,
      0, -1/3, -1/3,    0,    0,    0,    0,
      0,  2/3, -1/3,    0,    0,    0,    0,
      0, -1/3,  2/3,    0,    0,    0,    0,
      0,     0,   0, -1/3, -1/3,    0,    0,
      0,     0,   0,  2/3, -1/3,    0,    0,
      0,     0,   0, -1/3,  2/3,    0,    0,
      0,     0,   0,    0,    0, -1/3, -1/3,
      0,     0,   0,    0,    0,  2/3, -1/3,
      0,     0,   0,    0,    0, -1/3,  2/3)
dim(Z) = c(7,10)
Z = t(Z)

AG.fixeff2 = t(Z%*%t(AG.fixeff))

## Our grand mean has heavier tails
qqplot(BF.fixeff[,1],AG.fixeff2[,1],pch=21,bg=rgb(0,0,1,.2),col=NULL,asp=TRUE,cex=1,xlab="BayesFactor estimate",ylab="arm estimate",main="Grand mean posterior samples")
abline(0,1)

plot(colMeans(BF.fixeff[,-1]),colMeans(AG.fixeff2[,-1]),pch=21,bg=rgb(0,0,1,.2),col="black",asp=TRUE,cex=1.2,xlab="BayesFactor estimate",ylab="arm estimate",main="Fixed effect posterior means")
abline(0,1)

## Compare posterior standard deviations
BFsd = apply(BF.fixeff[,-1],2,sd)
AGsd = apply(AG.fixeff2[,-1],2,sd)
plot(sort(AGsd/BFsd),pch=21,bg=rgb(0,0,1,.2),col="black",cex=1.2,ylab="Ratio of posterior standard deviations (arm/BF)",xlab="Fixed effect index")
## AG estimates are slightly larger, consistent with sig2 estimates
## probably due to prior


## ----message=FALSE,warning=FALSE----------------------------------------------
library(languageR)
library(xtable)

## -----------------------------------------------------------------------------
data(primingHeidPrevRT)

primingHeidPrevRT$lRTmin1 <- log(primingHeidPrevRT$RTmin1)

###Frequentist 

lr4 <- lmer(RT ~ Condition + (1|Word)+ (1|Subject) + lRTmin1 + RTtoPrime + ResponseToPrime + ResponseToPrime*RTtoPrime +BaseFrequency ,primingHeidPrevRT)
# Get rid rid of some outlying response times
INDOL <- which(scale(resid(lr4)) < 2.5)
primHeidOL <- primingHeidPrevRT[INDOL,]

## -----------------------------------------------------------------------------
# Center continuous variables
primHeidOL$BaseFrequency <- primHeidOL$BaseFrequency - mean(primHeidOL$BaseFrequency)
primHeidOL$lRTmin1 <- primHeidOL$lRTmin1 - mean(primHeidOL$lRTmin1)
primHeidOL$RTtoPrime <- primHeidOL$RTtoPrime - mean(primHeidOL$RTtoPrime)

## -----------------------------------------------------------------------------
# LMER
lr4b <- lmer(  RT ~ Condition + ResponseToPrime +  (1|Word)+ (1|Subject) + lRTmin1 + RTtoPrime + ResponseToPrime*RTtoPrime + BaseFrequency , primHeidOL)
# BayesFactor
B5out <- lmBF( RT ~ Condition + ResponseToPrime +     Word +    Subject  + lRTmin1 + RTtoPrime + ResponseToPrime*RTtoPrime + BaseFrequency  , primHeidOL , whichRandom = c("Word", "Subject"),  posterior = TRUE, iteration = 50000,columnFilter=c("Word","Subject"))

lmerEff <- fixef(lr4b)
bfEff <- colMeans(B5out[,1:10])

## ----results='asis'-----------------------------------------------------------
print(xtable(cbind("lmer fixed effects"=names(lmerEff))), type='html')

## ----tidy=FALSE---------------------------------------------------------------
# Adjust lmer results from reference cell to sum to 0
Z = c(1,   1/2, 1/2,    0,    0,    0,    0,
      0,  -1/2,   0,    0,    0,    0,    0,
      0,   1/2,   0,    0,    0,    0,    0,
      0,     0,-1/2,    0,    0,    0,    0,
      0,     0, 1/2,    0,    0,    0,    0,
      0,     0,   0,    1,    0,    0,    0,
      0,     0,   0,    0,    1,    0,  1/2,
      0,     0,   0,    0,    0,    1,    0,
      0,     0,   0,    0,    0,    0, -1/2,
      0,     0,   0,    0,    0,    0,  1/2)
dim(Z) = c(7,10)
Z = t(Z)

# Do reparameterization by pre-multimplying the parameter vector by Z
reparLmer <- Z %*% matrix(lmerEff,ncol=1)

# put results in data.frame for comparison
sideBySide <- data.frame(BayesFactor=bfEff,lmer=reparLmer)

## ----results='asis'-----------------------------------------------------------
print(xtable(sideBySide,digits=4), type='html')

## -----------------------------------------------------------------------------
# Notice Bayesian shrinkage
par(cex=1.5)
plot(sideBySide[-1,],pch=21,bg=rgb(0,0,1,.2),col="black",asp=TRUE,cex=1.2, main="fixed effects\n (excluding grand mean)")
abline(0,1, lty=2)