Help with setting up a random effects model

I am slowly becoming familiar with mixed models. I am facing, however, some difficulties in defining a model for a new data set I am working on.

The dataset contains the reaction times RT for 50 subjects when they performed 3 tasks (the three tasks are A, B, and C). The task order was randomized between subjects. That is, subject 1 performed A, B, C; subject 2 performed B, A, C, and so forth. Each subject did the task only once, that is I have 3 RT measurements for each subject.

The data set could look like:

RT Subject Task Order
1.5 1 A I
1.1 1 B II
1.8 1 C III
1.0 2 A II
1.2 2 B I
1.5 2 C III

I would like to model the distribution of reaction times with a mixed model. I would like to know if there is any difference in reaction time for the task type, and the overall variability between subjects. One model could be:

RT ~ task + (task | subject)

I believe, however, that the order of the task was presented may have an effect on the reaction time. Therefore, by including the categorical variable order, I hope I would get some hints if there is a learning effect. The model would like:

RT ~ task * order + (task * order | subject)

Would this model be correct? Is this the way one would define a model to take into account the learning effect?

Is this question related to Stan or one of its interfaces?

It is a “general” question really. Once I understand how to set up a model, I would implement it in STAN or brms. I am sorry if this kind of question should not be ask in this forum.

If you are trying to get this model run in brms or rstanarm that it qualifies as a valid question I guess.

The model in general looks sensible, but you will likely not be able to fit it as you have just 3 RTs per person (if I understood you correctly), but 3 * 3 = 9 varying effects per person.

In fact, I tried to fit it in brms with this formula and with the lognormal family. The model was very slow to run and got many divergences and high r_hat values. What would be your suggestion? To remove some random effects (e.g., remove the interactions)?

You seem to have too few data to reliably fit more than 2 varying effects, right now you have 9.

Ok, that’s good information anyway, thank you.

You counted 3 \cdot 3=9 varying effects per person because both task and order have 3 levels? Thus, my design matrix X for the fixed effects would have 9 columns:


Then, the Z matrix for the random effects would have a 9 coefficients for each subject as well (i.e., 9 \cdot 50 = 450 total columns).

Given the measurements I have, you suggest that I can afford only 2 random effects per subject, which are too few even for the simplest model with only the task type as variable, because I have 3 coefficients to estimate.


Do you have any suggestions on how to get around this issue?

Run more trials per person, trial and order. Otherwise, you have to use
RT ~ task * order + (1 | subject)

1 Like

Thank you, I really appreciate your insights.