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”)
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.