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?
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)