Dependence between binary responses in multivariate brms

Modeling
1

I have data in which participants chose between 4 possible responses on each trial. The 4 responses are best understood, however, as two binary choices, i.e., “w” vs “b” and “2” vs “1”. What I want to get from the model is how making one of these choices depends on the other, i.e., how choosing “w” depends on choosing “2”, etc. for each of the 4 combinations. I think it makes no sense to treat the responses as multinomial (4 values) nor as ordinal. Example code treating this as a multivariate model without any dependence between responses is displayed below. The two response choices are “ncRsp” (“w” vs “b”) and “nsRsp” (“2” vs “1”). There are three continuous predictors, scVDuration, scTDuration, and scSpeed. Any advice about how to test for dependence between responses and how that might depend on the values of the predictors would be very welcome.

exp2l_0_14sep26 ← brm(mvbind(ncRsp, nsRsp) ~ scVDuration +
scTDuration +
scSpeed +
(1 + scVDuration + scTDuration + scSpeed | q | subject),
family = bernoulli(link = “logit”),
data = exp2l,
prior = priors,
warmup = 1000, iter = 5000, chains = 4, cores = 8, threads = threading(2),
save_pars = save_pars(all = TRUE),
control = list(adapt_delta = 0.99,
max_treedepth = 15),
seed = 101457,
file = “exp2l_0_14sep26”,
file_refit = “on_change”,
backend = “cmdstanr”)

2

I’ve done a little digging since putting up this post. Gemini suggested that I could assess the extent to which responses are correlated with one another by setting up the model as shown below:

formula_ncRsp ← bf(ncRsp ~ scVDuration +
scTDuration +
scSpeed +
(1 + scVDuration +
scTDuration +
scSpeed | q | participantFact))

formula_nsRsp ← bf(nsRsp ~ scVDuration +
scTDuration +
scSpeed +
(1 + scVDuration +
scTDuration +
scSpeed | q | participantFact))

exp2l_0probit_18sep26 ← brm(formula_ncRsp +
formula_nsRsp +
set_rescor(FALSE),
family = bernoulli(link = “probit”),
data = exp2l,
prior = priors,
init = “0”,
warmup = 1000, iter = 5000, chains = 4, cores = 8, threads = threading(2),
save_pars = save_pars(all = TRUE),
control = list(adapt_delta = 0.999,
max_treedepth = 15),
seed = 101457,
file = “exp2l_0probit_18sep26”,
file_refit = “on_change”,
backend = “cmdstanr”)

and then looking at the correlations between the two responses with respect to the intercept and the population-level effects.

I have a question and a problem.

The question: Is this the right way to proceed? Related: Can’t I just look at those correlations in the output of the model I first posted? By the way, both this approach and that one produce correlations for which the 95% CIs include 0.

The problem: This alternative analysis throws divergences in the chains on 25% of transitions.

3

Can you say more about the structure of the data? Do you have more than one response for ncRsp and nsRsp from each participant? Do all your predictors vary within participant? What would your model formula look like if you were fitting a model for ncRsp only?

4

I have multiple responses for each response variable from each participant (10 for each possible combination of the predictors), and all the predictors vary within participant. The model formula would therefore look the same for each response variable, namely:

ncRsp ~ scVDuration + scTDuration + scSpeed + (1 + scVDuration + scTDuration + scSpeed | participantFact

5

The only difference between your first model and the second one is that in the second you have set_rescor =F, while the first defaults to T. It is surprising to me that not estimating residual correlations correlations makes sampling more difficult; I would expect the opposite.

My understanding is that interpreting residual correlations without taking into account correlations among participants would be misleading but I don’t see why the reverse would be problematic. That is, if correlations among participants is what you are interested in, whether you also estimate residual correlations or not should not matter for your among participants correlations estimates.

6

Thanks!

set_rescor = TRUE doesn’t work for binary responses.

I need to check Gemini’s suggestion in other sources before proceeding.