Hi there :) I have a vision for a model but I have doubts about its feasibility.
I ran two tasks:
in task A subjs answered whether the stimuli were words or not and I am interested in their accuracy (bernoulli link-function)
in task B subjs answered what’s the ink color of the stimuli and I am interested in their reaction times (lognormal link-function).
Both tasks had three conditions (Stroop’s congruent, neutral and incongruent) and all subjs did all the tasks.
Ultimately, I am interested in predicting the accuracy of task A with the reaction times of task B and with the conditions as another predictor.
I can create two models and extract the effects (or random effects) of task B and use these summarized observations as predictors, but it feels off. I’m worried that will make me lose a lot of information. Any suggestions/directions for which model I should use?
Rather than modeling the reaction times of B, you could just use those reaction times as un-modeled covariates for predicting A. This would require both A and B had the same three conditions.
Let’s walk this out a bit to make sure we’re on the same page. In the case of Merge.Tasks, you have two blocks of Stroop tasks, which are differentiated by whether they have Correct or RT as the DV. However, the Stroop tasks are related in that those DVs are from 36 trials within each of the 3 levels of Condition, yielding 36 * 3 = 108 trials for each subject. But since you have two kinds of Strop tasks, there are 108 * 2 = 216 total trials for each subject. Now I see your point that it would be odd to use the RT values from individual trials as predictors. The individual trials aren’t tightly linked between the DVs. I have two possible solutions:
You could just take the average RT within each of the three levels of Condition, for each subject. Thirty-six is a decent number of trials to compute a population estimate of a subject’s true reaction time. You could even account for uncertainty by computing the standard error for that mean and use the brms::me() function (see here).
Alternatively, you could use the 2-step model-based solution you alluded to in your initial post. To attenuate the loss of information, you could use the brms::me() function there, too. But this time you’d enter in the posterior mean and standard deviation of the random effects, rather than the sample statistics. Either way, you’d be relying on brms::me() to capture the uncertainty in your predictor.
Q2
As far as I can tell, both of the two solutions, above, would be robust to differences in trial numbers between the DVs. But as always, more trials are better than fewer.
Thank you @JLC! To be honest I am still quite a beginner and I never tried multivariate approaches. After looking at the more complex models of this vignette, I still can’t really see how this procedure will results in the data of bf_correct controlled/predicted (or something similar to controlled/predicted) by the data of bf_rt. I would be thankful for some elaboration and/or sources.
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).