Bayesian GLM – What Would Be a Bayesian Equivalent for Mixed-Effects Logistic Regression Model

bayesiangeneralized linear modellme4-nlme

I have been using "glmer" in R to model a binary outcome for approximately 500 persons, in two groups, each measured at three points in time. My questions of interest are a) whether the change over time differs between groups and b) whether there is an overall higher response in one group.

Here is a slightly simplified version of my current model. 

m <- glmer(Shop ~ Time + Group + Time:Group + (1 | subj), 
           data = Shopping, family = binomial, 
           control = glmerControl(optimizer = "bobyqa"), nAGQ = 10)

I've gotten a result that I find plausible from this mixed-effects logistic model but would like to see if a Bayesian approach might suggest any different conclusion. To be honest, mostly I want to learn about Bayesian alternatives to the sorts of "Frequentist" models I've used for years.

What R procedure should I delve into here? I know there's one called "blme" in an R package of the same name. Is that a good one for a newcomer to Bayesian modelling to apply to my example problem?

EDIT

The brms suggestion was very apt. I loaded the brms, rstan and loo packages and was able to compare the loo and kfold types of AIC-like statistics to the fit statistics given by PROC GLIMMIX (SAS is my usual working tool and is where this model was originally run).

I also extended it to compare simpler (no random intercept) and over-fitted straw man alternative models, just to familiarize myself with how variant models are present in GLIMMIX and in the Stan world. 

require(brms)
require(rstan)
require(loo)

options(mc.cores = parallel::detectCores())
rstan_options(auto_write = TRUE)

m_stan <- brm(FMshop ~ Time*Group + ( 1 | subj) + 
                    age_BL, 
           data = Shopping, family = bernoulli)

m_stan_nomixed <- brm(FMshop ~ Time*Group + 
                               age_BL, 
           data = Shopping, family = bernoulli)

m_stan_rantime <- brm(FMshop ~ Time*Group + ( 1 + Time | subj ) +
                               age_BL, 
           data = Shopping, family = bernoulli)

m_stan_kfold<-kfold(m_stan,K=10)
m_stan_kfold_nomixed<-kfold(m_stan_nomixed,K=10)
compare(m_stan_kfold,m_stan_kfold_nomixed)
m_stan_kfold_rantime<-kfold(m_stan_rantime,K=10)
compare(m_stan_kfold,m_stan_kfold_rantime)

Best Answer

Stan is the state-of-the-art in Bayesian model fitting. It has an official R interface through rstan. With rstan you would need to learn how to write your models in the Stan language. Alternatively, Stan also provides the rstanarm package (hat-tip to @ben-bolker for pointing out the omission), through which you can write your models in the familiar lme4-style syntax. An equally user-friendly interface for Stan is the R package brms which is in addition very flexible to handle models that should satisfy basic and moderately advanced users.

For example, in your case the syntax would be exactly the same:

m <- brm(Shop ~ Time + Group + Time:Group + (1 | subj), 
       data = Shopping, family = binomial)

or more concretely (same would work with glmer as well)

m <- brm(Shop ~ Time*Group + (1 | subj), 
       data = Shopping, family = binomial)

This model in brms will assume reasonable defaults for the prior distributions but you are encouraged to select your own.

The syntax for basic models such as the one you give as an example is going to be the same between rstanarm and brms. The advantage of using rstanarmto fit these basic models is that it comes with pre-compiled Stan code so it is going to run faster than brms that needs to compile its Stan code for every model. To name a few distinguishing features, brms shines due to its extended support for different distributions (e.g. "zero-inflated beta", "von Mises", categorical), its extended syntax to cover cases where the user needs to model e.g. predictor or outcome (as in meta-analyses) measurement error, and its ability to fit distributional regressions, non-linear models, or mixture models. For a more extensive comparison of R packages for Bayesian analysis have a look at Bürkner 2018.

Since you are a newcomer to Bayesian models, I would also highly encourage you to read the book "Statistical Rethinking" which also comes with its own R package, rethinking that is also an excellent choice, although not as remarkably user-friendly and flexible as brms. There's even a version of the book adapted for brms.

**References** [Paul-Christian Bürkner, The R Journal (2018) 10:1, pages 395-411.][1]
Related Question