Adjustment for multiple comparison in bayesian multivariate regression model (using brms)

Question

I am investigating age and timepoint effects on different (correlated) EEG parameters in a repeated measurements structure. I chose to use the brms R package to fit a bayesian multivariate model with random intercepts for each subject:

[EEG dvs] ~ age *  timepoint + (1|id)

I am particularly interested on which of the 6 outcomes there is an effect of age. Since I don't have any well informed prior, I only want to interpret the 95% credibility interval and not calculate any Bayes factors.

How can I correct for the multiple comparisons (6 EEG parameters) in this case? I don't see how I can build a hierarchical model like the one proposed in Gelman's paper (http://www.stat.columbia.edu/~gelman/research/published/multiple2f.pdf) for the different EEG parameters.

Answer 1

If I correctly understand the goal here, it sounds as though you may need to restructure your model and data slightly.

First the model:

EEG ~ 1 + age*timepoint + (1|id) + (1 + age*timepoint|node)

where node is a label that indicates which sensor position the EEG value maps onto.

Now the data you might feed this model then needs to be in "stacked" or sometimes referred to as "long" format. Each row would represent the EEG value for individual $i$ derived from the $j$th sensor.

In this scenario, I imagine you might be primarily interested in the credibility interval around the random slopes specified by the (... + age*timepoint|node) portion of the model.

For a more structural equation-y way to use brms and achieve (again what I understand) the goal to be, you could create a series of equations as follows:

form_EEG_1 <- bf(
    EEG_1 ~ 1 + age*timepoint + (1|id)
) + gaussian() #assuming this is the appropriate family for EEG values

form_EEG_2 <- bf(
    EEG_2 ~ 1 + age*timepoint + (1|id)
) + gaussian()

.
.
.

form_EEG_6 <- bf(
    EEG_6 ~ 1 + age*timepoint + (1|id)
) + gaussian()

and then in the brm function you would add the formulae together as follows:

fit_EEG <- brm(formula = form_EEG_1 + form_EEG_2 + ... form_EEG_6 + set_rescor(TRUE), 
               [other brm control arguments]
               )

The downside with this second approach to your model is that you would not be able to include residual correlations among your DVs if you cannot assume a multivariate normal distribution for your EEG values. So depending on what is a reasonable prior for these scores, this may or may not be an option for you.

Note these approaches are not equivalent, and they handle the dependency in your data in slightly different ways so be sure to select an approach that best maps on to your analytic goals and your data properties.