I've been putting a lot of work over the last few days into bring mixed effects models to bear on some behavioural data I've collected for my thesis, but it's occurred to me that I'm not 100% sure that this kind of model is actually appropriate for my data (I only came across them after starting the experiment).

In an experiment, 60 participants completed 28 trials of a reasoning task, consisting of 14 problems (call them "A"-"N"), with participants completing each in 2 conditions, x (conflict) and y (control).

Participant  Problem
          1  Ax/Ay    Bx/By    Cx/Cy    Dx/Dy  ...  Nx/Ny
          2  Ax/Ay    Bx/By    Cx/Cy    Dx/Dy  ...  Nx/Ny
          3  Ax/Ay    Bx/By    Cx/Cy    Dx/Dy  ...  Nx/Ny
          4  Ax/Ay    Bx/By    Cx/Cy    Dx/Dy  ...  Nx/Ny
          5  Ax/Ay    Bx/By    Cx/Cy    Dx/Dy  ...  Nx/Ny
        ...  ...      ...      ...      ...    ...  ...
         60  Ax/Ay    Bx/By    Cx/Cy    Dx/Dy  ...  Nx/Ny

I'm interested in the difference in a trial-by-trial variable (let's call it reaction time) between the conflict and control conditions, and would expect it to be higher for the conflict (y) trials.

Obviously, I would expect to find a within-subject correlation - some subjects are generally fast, some generally slow. I would also expect to find a within-problem correlation - some problems are answered faster than others, regardless of condition. To account for these, I include random intercepts for these two factors:

(1|participant) + (1|problem).

The difference between conflict and control conditions may or may not turn out to be the same for each subject, and for each problem. For this reason, I consider including random slopes as well:

(1 + condition|participant) + (1 + condition|problem).

Putting this together, I'm testing a model that looks either like:

null_model = lmer(reaction_time ~ condition(1|participant) + (1|problem), data=data)
condition_model = lmer(reaction_time ~ condition + (1|participant) + (1|problem), data=data)


null_model = lmer(reaction_time ~ (1 + condition|participant) + (1 + condition|problem), data=data)
condition_model = lmer(reaction_time ~ condition(1|participant) + (1|problem), data=data)


Please; have I horribly misunderstood how this is supposed to work?

Edit: The more traditional approach to analysing this data would be to average across problems within each participant, yielding two data points per participant: conflict condition mean and control condition mean, and then use a paired-samples t test. Reading this question, I thought that this approach should be largely the same as fitting

lmer(reaction_time ~ condition + (1|participant), data=data),

but trying this in R, it seems otherwise.

Edit #2: Bounty added.


This sounds like a two-level cross-classified model, where reaction time is predicted by condition, and variance in reaction time is structured within individuals, and independently within problems (i.e. the structure of your data is not "problems within persons" or "persons within problems"). You may also wish to include time in the model if problem fatigue/competence or problem ordering is a concern.

$y_{i,j_{1},j_{2}}$ is reaction time, indexed by $i$ (answers-within-persons-and-problems) which are within $j_{1}$ persons and $j_{2}$ problems.

$x$ is condition

$z_{1}$ is person id

$z_{2}$ is problem id

$y_{i,j_{1},j_{2}} = \beta_{0i,j_{1},j_{2}} + \beta_{x}x + \mu_{j_{1}} + \mu_{j_{2}} + \varepsilon_{i,j_{1},j_{2}}$,


$\beta_{0}$ is the intercept which varies randomly at levels $i$, $j_{1}$, and $j_{2}$,

$\beta_{x}$ is the fixed-effect of condition on $y$,

$\mu_{j_{1}}$ is the person-level random-effect, in this case contributing a person-level variance to the intercept.

$\mu_{j_{2}}$ is the problem-level random-effect, in case contributing a problem-level variance to the intercept.

$\varepsilon_{i,j_{1},j_{2}}$ is the answer-within-persons-and-places-level residual variance.

If you were also measuring person-level characteristics, or problem-level characteristics and wanted to add them to the model, you might then given them random-effects (i.e. random slopes) also. But the low-sample sizes make me think that's not the case in this example.

The low number of problems (14) means that you probably won't be performing much inference about problem-level variance, but it does mean you can account for problem-level variance while accounting for person-level variance.

A good multilevel modeling textbook (e.g. Goldstein, H. (2003). Multilevel statistical models. Oxford University Press, 3rd edition.) will give a good grounding in cross-classified models. I think the MLwiN manual givens good tutorials on them as well.

Update: I haven't yet dove into R for mixed effect models, but I think that the command for the cross-classified model you want is:

lmer(reaction_time ~ 1 + condition + (1|participant) + (1|problem), data)

  • $\begingroup$ Thanks, this seems a great answer, and the suggest about time is very welcome. Does this mean that I'm right in fitting this model as condition_model = lmer(reaction_time ~ condition + (1|participant) + (1|problem), data=data)? I'm also not 100% sure on when you should and shouldn't include random slopes (you've a typo there, by the way), but that may be a question for another time. $\endgroup$ – Eoin Apr 28 '14 at 15:00
  • $\begingroup$ I have corrected the typo, and provided my best guess as to the R syntax corresponding to the model I propose. You would include a random slope for a given level of data structure if you felt that the effect of a predictor ($\beta$) varied at that level (e.g. the effect of condition on reaction time differed between people—which is different that the person-level random intercept, which merely says that average reaction time varies between people). Advanced: Random slopes can also model heteroscedasticity of predictors at the residual level. $\endgroup$ – Alexis Apr 28 '14 at 15:44

As it turns out, my design is almost exactly that described by Baayen, Davidson and Bates (2008) ("Mixed-effects modeling with crossed random effects for subjects and items").

Data of this form is handled explicitly by R package lme4 (written by the third author), but not, apparently, by nlme.

My syntax in the question was (I think) correct: the random-intercept model is fitted by lmer(reaction_time ~ condition + (1|participant) + (1|problem), data=data).

Finally, I've found a useful overloaded version of the lmer function in package lmerTest which provides a number of methods for calculating p values for terms in mixed-models, and package MuMin provides an implementation of Nakagawa and Schielzeth's (2012) method of calculating $R^{2}$ for mixed-models (according to this blog).


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.