Based on pilot data, I have several significant links between different types of y values and different types of x values. To reduce the need for multiple comparison correction, I was wondering about integrating this in a model hierarchically, but I did not manage. Instead, the solution so far I have found (if it is correct…) is only using model comparison.

```
# x values
x_vars = data.frame(x1 = rnorm(50,mean=0,sd=1),
x2 = rnorm(50,mean=0,sd=1),
x3 = rnorm(50,mean=0,sd=1))
# regression weights that will link x_vars to y
x_regs=data.frame(x1.y1 = rnorm(50,mean=1,sd=1),
x2.y1 = rnorm(50,mean=1,sd=1),
x3.y1 = 0,
x1.y2 = rnorm(50,mean=-1,sd=1),
x2.y2 =0,
x3.y2 = rnorm(50,mean=-1,sd=1),
x1.y3 = 0,
x2.y3 = rnorm(50,mean=2,sd=1),
x3.y3 = rnorm(50,mean=2,sd=1))
y = data.frame(y1 = x_vars$x1* x_regs$x1.y1 + x_vars$x2* x_regs$x2.y1 ,
y2 = x_vars$x1* x_regs$x1.y2 + 0+ x_vars$x3* x_regs$x3.y2,
y3 = 0 + x_vars$x2* x_regs$x2.y3 + x_vars$x3* x_regs$x3.y3)
# bind y and x together to give to brms
df = cbind(y,x_vars) %>% dplyr::mutate_all(scale)
fit.selective = brm(bf(y1~x1+x2)+ bf(y2~x1+x3)+bf(y3~x2+x3)+set_rescor(TRUE),
data=df,prior=prior(normal(0,1), class = b),save_pars = save_pars(all = TRUE))
fit.null= brm(bf(mvbind(y1,y2,y3) ~ 1)+set_rescor(TRUE),data=df2, save_pars = save_pars(all = TRUE))
bayes_factor(fit.selective,fit.null)
```

And then I would say if this Bayes_factor suggests that indeed including the x-values in the regressions is better, I follow-up report regression weights for individual regressors. (This is just a toy example, in my real data, I have more nested groupings, I would test each with another bayes_factor analysis if the ‘hierarchically next up’ test is significant.)

Is this valid? And is there a way to include knowledge from the pilot data not just about which y-x links to test, but also about the direction? I.e. if suddenly the sign of x1-y1 was the other way round, that would be bad and not align with my hypothesis. E.g. could I use priors centred on the pilot data mean for each regression weight?

Is there an alternative model that could be done where instead of Bayes factor, with a single regressor that expresses the overall significance? I was thinking of somehow maybe nesting x1.y1, x2.y2 together, but I didn’t know how to implement this correctly (I tried row binding the ys and including a regressor for type of y and then e.g. x1|y.type, but those models did not converge)

I’d be very grateful for any pointers

Jacquie