compare_lme4.html 1.05 MB
 Jonathon Love committed Sep 24, 2015 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60  Comparison of BayesFactor against other packages  Andreas Tille committed May 08, 2018 61 

Comparison of BayesFactor against other packages

This R markdown file runs a series of tests to ensure that the BayesFactor package is giving correct answers, and can gracefully handle probable input.

library(arm) library(lme4)

ANOVA

First we generate some data.

# 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)) ## ## Error: ID ## Df Sum Sq Mean Sq F value Pr(>F) ## Residuals 19 648 34.1 ## ## Error: ID:A ## Df Sum Sq Mean Sq F value Pr(>F) ## A 2 13.8 6.92 8.12 0.0012 ** ## Residuals 38 32.4 0.85 ## --- ## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 ## ## Error: ID:B ## Df Sum Sq Mean Sq F value Pr(>F) ## B 2 2.4 1.19 1.18 0.32 ## Residuals 38 38.4 1.01 ## ## Error: ID:C ## Df Sum Sq Mean Sq F value Pr(>F) ## C 2 5.5 2.767 2.95 0.064 . ## Residuals 38 35.6 0.937 ## --- ## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 ## ## Error: ID:A:B ## Df Sum Sq Mean Sq F value Pr(>F) ## A:B 4 1.6 0.402 0.41 0.8 ## Residuals 76 73.8 0.971 ## ## Error: ID:A:C ## Df Sum Sq Mean Sq F value Pr(>F) ## A:C 4 12.4 3.103 3.33 0.014 * ## Residuals 76 70.7 0.931 ## --- ## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 ## ## Error: ID:B:C ## Df Sum Sq Mean Sq F value Pr(>F) ## B:C 4 2.3 0.583 0.46 0.77 ## Residuals 76 96.6 1.271 ## ## Error: ID:A:B:C ## Df Sum Sq Mean Sq F value Pr(>F) ## A:B:C 8 12.6 1.58 1.4 0.2 ## Residuals 152 170.7 1.12 We can plot the data with standard errors: 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)) }
 Andreas Tille committed May 08, 2018 325 

 Jonathon Love committed Sep 24, 2015 326 327 328 329 330 331 332 333 334 335 336 337 338 339 340 341 342 343 

Bayes factor

Compute the Bayes factors, while testing the Laplace approximation

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") )
t.is
##    user  system elapsed  
Andreas Tille committed May 09, 2018  344                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                            ##   9.090   0.157   9.865 
Jonathon Love committed Sep 24, 2015  345     346     347     348     349     350                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                    
t.la
##    user  system elapsed  
Andreas Tille committed May 09, 2018  351                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                            ##   4.601   0.053   4.796 
Jonathon Love committed Sep 24, 2015  352     353     354     355     356     357     358     359                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                    
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)
 Andreas Tille committed May 08, 2018 360 

 Jonathon Love committed Sep 24, 2015 361 362 363 364 365 366 367 368 369 370 371 372 373 374 375 376 377 378 379 380 381 382 383 384 385 386 387 388 389 390 391 392 393 394 395 396 397 398 399 400 401 402 403 404 405 406 407 408 409 410 411 412 
bfs.is
## Bayes factor analysis ## -------------- ## [1] A + ID                                    : 9.02     ±0.93% ## [2] B + ID                                    : 0.059    ±1.76% ## [3] A + B + ID                                : 0.563    ±6.68% ## [4] A + B + A:B + ID                          : 0.00796  ±4.56% ## [5] C + ID                                    : 0.229    ±0.85% ## [6] A + C + ID                                : 2.29     ±5.06% ## [7] B + C + ID                                : 0.0133   ±1.08% ## [8] A + B + C + ID                            : 0.141    ±7.7% ## [9] A + B + A:B + C + ID                      : 0.00196  ±4.48% ## [10] A + C + A:C + ID                         : 2.49     ±13.1% ## [11] A + B + C + A:C + ID                     : 0.138    ±2.12% ## [12] A + B + A:B + C + A:C + ID               : 0.00206  ±2.86% ## [13] B + C + B:C + ID                         : 0.000251 ±1.44% ## [14] A + B + C + B:C + ID                     : 0.00253  ±2.14% ## [15] A + B + A:B + C + B:C + ID               : 3.67e-05 ±2.09% ## [16] A + B + C + A:C + B:C + ID               : 0.00264  ±1.66% ## [17] A + B + A:B + C + A:C + B:C + ID         : 4.11e-05 ±2.89% ## [18] A + B + A:B + C + A:C + B:C + A:B:C + ID : 8.47e-06 ±1.88% ##  ## Against denominator: ##   y ~ ID  ## --- ## Bayes factor type: BFlinearModel, JZS

Comparison to lmer and arm

We can use samples from the posterior distribution to compare BayesFactor with lmer and arm.

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)

Compare estimates of variance

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)
 Andreas Tille committed May 08, 2018 413 

 Jonathon Love committed Sep 24, 2015 414 415 416 417 418 419 420 421 422 

Compare estimates of participant effects:

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)
 Andreas Tille committed May 08, 2018 423 

 Jonathon Love committed Sep 24, 2015 424 425 426 427 428 429 430 431 432 433 434 435 436 437 438 439 440 441 442 443 444 445 446 447 448 449 450 

Compare estimates of fixed effects:

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)
 Andreas Tille committed May 08, 2018 451 

 Jonathon Love committed Sep 24, 2015 452 453 454 455 456 
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)
 Andreas Tille committed May 08, 2018 457 

 Jonathon Love committed Sep 24, 2015 458 459 460 461 462 463 464 
## 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")
 Andreas Tille committed May 08, 2018 465