Using R studio to do Bayesian Multilevel Modeling

Dear Gang Chen,

I recently came across the paper titled *"Age of onset of obsessive-compulsive disorder differentially affects white matter microstructure", which utilized the Region-Based Analysis Program through Bayesian Multilevel Modeling to detect white matter alterations between different groups. I am interested in performing a similar analysis on my own data using the brms package in R.

For my analysis, I have two groups of participants labeled as "bl" and "hc." For each individual, I extracted FA measures for 72 tracts separately. My data is structured as follows:

my data

label Age Gender image_ID measurement value
hc 20 0 1 tract_1 0.45
hc 20 1 1 tract_2 0.40
...
bl 21 0 66 tract _1 0.40
bl 21 0 66 tract_2 0.42
...

My goal is to detect FA alterations in the "bl" group compared to the "hc" group, adjusting for age and gender as covariates. Additionally, I have included image_ID and measurement as random effects. Below is the model specification that I am currently using:

my prior

prior_1 <- c(
prior(normal(0, 0.02), class = "b"),
prior(normal(0.5, 0.005), class = "Intercept"),
prior(exponential(20), class = "sd"),
prior(exponential(30), class = "sigma"),
prior(lkj(2), class = "cor")
)

brms modeling

model_3 <- brms::brm(
value ~ label * Age * Gender + (1 + label | image_ID) + (1 | measurement),
data = data,
chains = 4,
warmup = 1000,
iter = 5000,
thin = 1,
prior = prior_1,
#control = list(max_treedepth = 15)
).

Since the article did not provide enough details, I would like to know should I included all variables in a single model to detect alterations across the tracts, or produced separate models for each tract? I would appreciate it if you could review my modeling approach and provide any feedback.

Thank you in advance for your help.

Best regards,
Shuangwei Chai

Hi Shuangwei,

Regarding the three variables of label, Age, and Gender, is the primary effect of interest the comparison between the bl and hc groups, with Age and Gender included as covariates?

If you have AFNI installed on your computer, you can directly use the RBA program, which streamlines both the modeling and effect extraction steps:

https://afni.nimh.nih.gov/pub/dist/doc/program_help/RBA.html

Gang Chen

Dear Gang Chen,

Thank you for your response. I am interested in detecting alterations in diffusion metrics (such as FA and MD) between the bl and hc groups. Specifically, I have extracted FA values from 72 white matter tracts for each participant. I would like to know whether I should include all 72 FA values in the model simultaneously or evaluate them separately.

Unfortunately, I have not installed AFNI on my computer, and I am hoping to replicate these analyses using R.

Thank you in advance for your assistance.

Best regards,
Shuangwei Chai

It would be easier to perform the Bayesian hierarchical modeling using the RBA program, provided you're open to installing AFNI.

Alternatively, if you prefer to analyze the data directly in R using the brms package, your model code looks fine:

model_3 <- brms::brm(
  value ~ label * Age * Gender + (1 + label | image_ID) + (1 | measurement),
  data = data,
  chains = 4,
  warmup = 1000,
  iter = 5000,
  thin = 1,
  prior = prior_1
  # control = list(max_treedepth = 15)
)

One subtle but important step is centering the age variable, which can improve the interpretability of the results. The more tedious part lies in extracting and assembling the effects of interest from the brms output.

Gang Chen