I actually think that if I understand the problem correctly, @JLC’s approach is worth considering (although @Solomon’s might in the end prove better). The important part here is the `(Condition|p|Subjs)`

term that appears in both formulas. This is somewhat convoluted syntax but it can be understood as:

- For both responses there is an unobserved “ability” of each subject to have correct answers in one task and have long reaction times in the other task (modelled as a varying effect)
- Those abilities are correlated between the tasks, so e.g. latent “ability” to respond correctly in the congruent condition is correlated with the “ability” to have long reaction times in the congruent condition.
- The amount and sign of correlation is estimated from data.

So here, the correlations between the subjects varying effects for the two tasks serve a similar purpose as the coefficients of the reaction times in @Solomon’s approach. IMHO the main conceptual difference is that when using the (average of) the reaction times to predict correct answers the predictors are on the response scale of the reaction time submodel and are directly influenced by the measurement noise (but not much since you are averaging). In the multivariate model, you look at correlations between variables on the latent linear scale. One would thus expect that inspecting the fitted correlations should give similar qualitative results as inspecting the fitted coefficients when using RT as a predictor.

A practical difference is also that correlations in my (limited) experience usually require more data to fit well than fixed effect coefficients.

Also, since the data from the two tasks don’t correspond directly to each other you can’t directly feed the data to `brms`

. But I think you could do something like (code not tested):

```
Merged2 <- rbind(Task.A %>% mutate(RT = NA, isA = TRUE),
Task.B %>% mutate(Correct = NA, isA = FALSE))
bf_correct <- bf(Correct | subset(isA) ~ Condition + (Condition|p|Subjs), family = bernoulli())
bf_rt <- bf(RT | subset(!isA) ~ Condition + (Condition|p|Subjs), family = "gamma")
fit <- brm(bf_correct + bf_rt + set_rescor(FALSE), data = Merged2)
```

Where the `subset`

term tells `brms`

to actually use only some rows for the relevant submodel. I am speculating here a bit, so feel free to ask for clarifications if you’d like to implement the model but encoutner difficulties.

**EDIT:** I just tested that the code above and after slight modifications it leads to a model that compiles and runs, although the simulated dataset is giving divergences (possibly because it is too small).

Best of luck with your model!