6
$\begingroup$

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)
$\endgroup$
  • $\begingroup$ I've used MCMCglmm for Bayesian linear mixed models. I haven't used the binomial family very much but it does have that as an option. $\endgroup$ – TPM Nov 15 '18 at 21:00
  • $\begingroup$ I use rstan. Check out maths.bath.ac.uk/~jjf23/stan for helpful material on the linear model. Extending it to the logistic model is relatively trivial. You should also keep the Stan manual handy mc-stan.org/users/documentation $\endgroup$ – Heteroskedastic Jim Nov 16 '18 at 3:46
7
$\begingroup$

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

| cite | improve this answer | |
$\endgroup$
  • 1
    $\begingroup$ good answer. Also rstanarm as an alternative to brms $\endgroup$ – Ben Bolker Nov 17 '18 at 13:42

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.