package comparison

## fitting methods
library(lme4)
library(blme)
library(glmmTMB)
library(brms)
## plotting/manipulation
library(broom)  ## install_github("bbolker/broom")
library(dotwhisker)
library(dplyr)
library(ggplot2)

Define a formula (which we’ll use repeatedly) and make a data frame that represents a fully crossed, randomized-block design with three factors for the fixed effects (3x3x2) and two random effects (id and item.new).

form <- contrast~c.con.tr*c.type.tr*c.diff.tr+(1|id)+(1|item.new)
mydata <- expand.grid(c.con.tr=factor(1:3),
                      c.type.tr=factor(1:2),
                      c.diff.tr=factor(1:2),
                      id=factor(1:10),
                      item.new=factor(1:10))

Simulate a Bernoulli response:

set.seed(101)
mydata$contrast <- simulate(form[-2],
                            newdata=mydata,
                            newparams=list(beta=rep(1,12),
                                           theta=rep(1,2)),
                            family=binomial,
                            weights=rep(1,nrow(mydata)))[[1]]

Fit the model with lme4, blme, glmmTMB, and brms (the first and third are un-penalized, the second and fourth set a zero-centered prior on the fixed effects):

riobglmer <- bglmer(form,
                    data=mydata, family=binomial (link='logit'),
                    fixef.prior= normal(cov = diag(9,12)))
rioglmer <- glmer(form,
                    data=mydata, family=binomial (link='logit'))
rioglmmTMB <- glmmTMB(form,
                    data=mydata, family=binomial (link='logit'))
prior <- get_prior(form,
                    data=mydata, family=bernoulli)
prior$prior[1] <- "normal(0,3)"
riobrms <- brm(contrast~c.con.tr*c.type.tr*c.diff.tr+(1|id)+(1|item.new),
               prior=prior,
               chains=3,
               data=mydata, family=bernoulli)
resList <- list(blme=riobglmer,lme4=rioglmer,glmmTMB=rioglmmTMB,
                brms=riobrms)
saveRDS(resList,file="bglmer_runs.rds")

Retrieve the stored data:

resList <- readRDS("bglmer_runs.rds")

Tweak labels a bit and put the results together:

Plot:

dwplot(resFrame)+theme_bw()+
    geom_vline(xintercept=1,lty=2)+
    scale_x_continuous(limits=c(-5,5),oob=scales::squish)+
    scale_colour_brewer(palette="Set1")

Conclusions:

• all approaches give generally similar answers
• even this fairly large data set isn’t enough to give accurate answers about interactions (unpenalized models give ridiculous answers)
• it’s a bit worrying that 2/12 parameters aren’t covered here (i.e., confidence intervals for c.con.tr3 and c.con.tr3:c.diff.tr2); this could be a weird single simulation (I don’t have the patience to redo this hundreds of times - brms in particular runs slowly)

Session info

sessionInfo()

## R Under development (unstable) (2017-05-26 r72742)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 14.04.5 LTS
## 
## Matrix products: default
## BLAS: /usr/local/lib/R/lib/libRblas.so
## LAPACK: /usr/local/lib/R/lib/libRlapack.so
## 
## locale:
##  [1] LC_CTYPE=en_CA.UTF8       LC_NUMERIC=C             
##  [3] LC_TIME=en_CA.UTF8        LC_COLLATE=en_CA.UTF8    
##  [5] LC_MONETARY=en_CA.UTF8    LC_MESSAGES=en_CA.UTF8   
##  [7] LC_PAPER=en_CA.UTF8       LC_NAME=C                
##  [9] LC_ADDRESS=C              LC_TELEPHONE=C           
## [11] LC_MEASUREMENT=en_CA.UTF8 LC_IDENTIFICATION=C      
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
##  [1] bindrcpp_0.1       dplyr_0.7.0        dotwhisker_0.2.6  
##  [4] gtable_0.2.0       gridExtra_2.2.1    broom_0.4.2       
##  [7] brms_1.7.0         ggplot2_2.2.1.9000 Rcpp_0.12.11      
## [10] glmmTMB_0.1.3      blme_1.0-4         lme4_1.1-14       
## [13] Matrix_1.2-10      rmarkdown_1.5     
## 
## loaded via a namespace (and not attached):
##  [1] rstan_2.15.1         TMB_1.7.10           reshape2_1.4.2      
##  [4] splines_3.5.0        lattice_0.20-35      colorspace_1.3-2    
##  [7] htmltools_0.3.6      stats4_3.5.0         loo_1.1.0           
## [10] yaml_2.1.14          rlang_0.1.1          nloptr_1.0.4        
## [13] foreign_0.8-68       glue_1.0.0           RColorBrewer_1.1-2  
## [16] bindr_0.1            matrixStats_0.52.2   plyr_1.8.4          
## [19] stringr_1.2.0        munsell_0.4.3        coda_0.19-1         
## [22] psych_1.7.5          evaluate_0.10        labeling_0.3        
## [25] inline_0.3.14        knitr_1.16           parallel_3.5.0      
## [28] bayesplot_1.2.0      rstantools_1.2.0     scales_0.4.1        
## [31] backports_1.1.0      StanHeaders_2.15.0-1 abind_1.4-5         
## [34] mnormt_1.5-5         digest_0.6.12        stringi_1.1.5       
## [37] grid_3.5.0           rprojroot_1.2        tools_3.5.0         
## [40] magrittr_1.5         lazyeval_0.2.0       tibble_1.3.3        
## [43] pkgconfig_2.0.1      tidyr_0.6.3          MASS_7.3-47         
## [46] assertthat_0.2.0     minqa_1.2.4          R6_2.2.1            
## [49] nlme_3.1-131         compiler_3.5.0