With Stan and frontend packages rstanarm
or brms
I can easily analyze data the Bayesian way as I did before with mixed-models such as lme
. While I have most of the book and articles by Kruschke-Gelman-Wagenmakers-etc on my desk, these don't tell me how to summarize results for a medical audience, torn between the Skylla of Bayesian's wrath and the Charybdis of medical reviewers ("we want significances, not that diffuse stuff").
An example: Gastric frequency (1/min) is measured in three groups; healthy controls are the reference. There are several measurements for each participant, so à la frequentist I used the following mixed-model lme
:
summary(lme(freq_min~ group, random = ~1|study_id, data = mo))
Slightly edited results:
Fixed effects: freq_min ~ group
Value Std.Error DF t-value p-value
(Intercept) 2.712 0.0804 70 33.7 0.0000
groupno_symptoms 0.353 0.1180 27 3.0 0.0058
groupwith_symptoms 0.195 0.1174 27 1.7 0.1086
For simplicity, I will use 2* std error as 95% CI.
In frequentist context, I would have summarized this as:
- In the control group the estimated frequency was 2.7/min (maybe add CI here, but I avoid this sometimes because of the confusion created by absolute and difference CI).
- In the no_symptoms group, the frequency was higher by 0.4/min, CI(0.11 to 0.59)/min, p = 0.006 than control.
- In the with_symptoms group, the frequency was higher by 0.2/min, CI(-0.04 to 0.4)/min, p = 0.11 than control.
This is about the maximum acceptable complexity for a medical publication, the reviewer will probably ask me to add "not significant" in the second case.
Here is the same with stan_lmer
and default priors.
freq_stan = stan_lmer(freq_min~ group + (1|study_id), data = mo)
contrast lower_CredI frequency upper_CredI
(Intercept) 2.58322 2.714 2.846
groupno_symptoms 0.15579 0.346 0.535
groupwith_symptoms -0.00382 0.188 0.384
where CredI are 90% credible intervals (see the rstanarm vignette why 90% is used as default.)
Questions:
- How to translate the above summary to the Bayesian world?
- To what extent is prior-discussion required? I am quite sure the paper will come back with the usual "subjective assumption" when I mention priors; or at least with "no technical discussion, please". But all Bayesian authorities request that interpretation is only valid in the context of priors.
- How can I deliver some "significance" surrogate in formulation, without betraying Bayesian concepts? Something like "credibly different" (uuuh...) or almost credibly different (buoha..., sounds like "at the brim of significance).
Jonah Gabry and Ben Goodrich (2016). rstanarm: Bayesian Applied Regression Modeling via Stan. R package version 2.9.0-3. https://CRAN.R-project.org/package=rstanarm
Stan Development Team (2015). Stan: A C++ Library for Probability and Sampling, Version 2.8.0. URL http://mc-stan.org/.
Paul-Christian Buerkner (2016). brms: Bayesian Regression Models using Stan. R package version 0.8.0. https://CRAN.R-project.org/package=brms
Pinheiro J, Bates D, DebRoy S, Sarkar D and R Core Team (2016). nlme: Linear and Nonlinear Mixed Effects Models. R package version 3.1-124, http://CRAN.R-project.org/package=nlme>.
mean(as.matrix(freq_stan)[,"groupwith_symptoms"] < 0)
. $\endgroup$group_nosymptoms
and then say the probability of it being negative is1 / draws
. But for the intercept, the chain is never going to wander into negative territory for these data, so I guess you could say the probability is less than1 / draws
. $\endgroup$