Brms - is my categorical predictor model specification/write up correct?

Interfaces brms
,
1

Hi,
I have completed my first analysis and write-up of some clinical trial data using brms - was hoping to get some experienced eyes to look at it. My model which has only categorical population-level predictors is specified as:

# Specify priors and run Stan model -------------------------------------------
prior1 <- c(prior(normal(0, 10), class = b),
            prior(normal(30, 5), class = b, coef = "Intercept")) 

model1 <- brm(formula = MADRS ~ 0+ Intercept + Group+Time+ Group:Time + (1|ID),
              family = gaussian,
              prior = prior1,
              data = dressMADRS,
              chains = 3, cores = 3,warmup = 500, iter = 10000,
              sample_prior = TRUE);

In my write-up I have defined the model as:

y \sim \operatorname{Normal}(\mu, \sigma_{t})
\mu = \beta_{0} + s_j + (\overrightarrow{\beta_{Time}} \times \overrightarrow{X_{Time}}) + (\overrightarrow{\beta_{Group}} \times \overrightarrow{X_{Group}}) + (\overrightarrow{\beta_{Int}} \times \overrightarrow{X_{Int}})
\overrightarrow{\beta_{Age}},\overrightarrow{\beta_{Time}},\overrightarrow{\beta_{Int}} \sim \operatorname{Normal}(0, 10)
\beta_{0} \sim \operatorname{Normal}(30, 5)
s,\sigma_{t},sd(s) \sim \operatorname{Half Student}(3,0,mad(y))

with the following text description:

Overarrows are used to indicate dummy coding of categorical fixed effect variables and sj
corresponds to the random efffects of participant j. Relatively non-committal priors are used. For the overall intercept a normal prior with mean of 30 with sd of 10 gives 95% coverage of the range 20-40. 20 is the minimum entry criteria for the study while a patient with score > 40 would typically be excluded. For fixed effects, normal priors are used with an sd of 10 which allows for a range of effect sizes of size comparable to the literature. sd priors are specified as default half student t priors as suggested by Gelman (2006) using the median absolute deviation (mad) of y as the scale factor (BĂĽrkner (2017)). For MCMC modelling Stan is acccessed through the brms interface (BĂĽrkner (2017)). For computing Bayes Factors, a ROPE of range |sd/10| is used (Kruschke and Liddell (2018)).

With thanks in advance!

1 Like
2

Hi,
sorry for not getting to you earlier. This looks mostly reasonable, I do have a few minor comments:

Shouldn’t those overarrows be only over the data (X)? The coefficents themselves are not coded in any way, right?

I think you are missing a definition of s_j in the model. I would put

s_j \sim Normal(0, sd(s))

I generally find it more readable to write the models with full indexing over observations to make it explicit, what changes and what is the same, e.g. something like:

y_i \sim Normal(\mu_i, \sigma_t) \\ \mu_i = \beta_{0} + s_{id(i)} + (\overrightarrow{\beta_{Time}} \times \overrightarrow{X_{Time,i}}) + ...

I guess you meant “sd of 5”?

A potentiall bigger issue is that I am not sure how Bayes factors work with interactions - if you test separately for the Group in ROPE and Group:Time in ROPE, you IMHO might get misleading results and unless the main effect and interaction are of interest separately you should take care to test them both at the same time, i.e. it is possible that you’ll get both the main effect and interaction consistent with no difference, but the data will not be consistent with no difference between Control at Time 2 and Treatment at Time 2…

Best of luck with your model!

3

Hi Martin,
Thanks for replying at all (!) and looking over it. All makes sense to me and thanks for picking up the mistakes. Just with regards to your comment on Bayes Factors. I was planning on mostly emphasising the 95% HDIs of the posterior beta parameter estimates but thought I should also include some Bayes Factors. I calculated them like this (using BayesTestR) which I think (?) is consistent with your comment and calculates BFs for all betas at the same time:

bf1 <- bayesfactor_parameters(model1, prior = NULL, null = c(-sd_y/10, sd_y/10))

Thanks,
Suresh

1 Like