I have a dataset with the following variables, which are measured between subjects (i.e. 1 measurement per individual):

  • ID: subject identifier
  • age: subject age
  • activity: overall activity level of that individual

Each individual completes 2 computer tasks (one easy, one hard) and reaction time (RT) was measured for each. So I have the following variables measured within subjects (repeated measures):

  • difficulty: task difficulty
  • RT: reaction time on the task

My data looks something like this:

df <- data.frame("ID" = c(1,1,2,2,3,3,4,4),
                 "age" = c(20,20,25,25,19,19,21,21),
                 "activity" = c(55,55,72,72,83,83,67,67),
                 "difficulty" = c(0,1,0,1,0,1,0,1),
                 "RT" = c(110,250,90,100,99,132,122,134))
df$difficulty <- factor(df$difficulty)

  ID age activity difficulty  RT
1  1  20       55          0 110
2  1  20       55          1 250
3  2  25       72          0  90
4  2  25       72          1 100
5  3  19       83          0  99
6  3  19       83          1 132
7  4  21       67          0 122
8  4  21       67          1 134

I'm expecting to see an interaction such that the slope between activity and RT will be different for each level of difficulty factor (while controlling for age). I have tried the following model in R:

lm(RT ~ activity*difficulty + age, data=df)

My concern here is the repeated nature of some of my variables. As you can see from my sample data there are twice as many rows as there should be. The values for my between subject variables are doubled and each participant was measured twice, which will affect my degrees of freedom, and thus my p values.

Is this a valid concern when testing regression interactions involving within subject variables?

Is there a more appropriate way to test this and what would it look like in R?

  • $\begingroup$ also wrth looking into is the plm-package for panel data $\endgroup$ Mar 8, 2019 at 0:36

1 Answer 1


The lm() model you proposed would make sense if each individual in your study provided a single reaction time (RT) value. However, each individual provides two RT values which are likely correlated for that individual. That correlation invalidates the assumption of independence of the RT values required by the lm() model, hence the results produced by that model.

What you need to use instead of the lm() function is the lmer() function from the lme4 package, which implements a linear mixed effects model. Here is an example of lmer syntax:



m <- lmer(RT ~ activity*difficulty + age + (1 + difficulty | ID), data = df)

where ID needs to be coded as a factor. The above model includes a random intercept and a random slope for difficulty, as well as a random grouping factor (i.e., individual). For a model with just a random intercept, you can use this syntax:

m <- lmer(RT ~ activity*difficulty + age + (1 | ID), data = df)

This tutorial will get you started with linear mixed effects models: http://www.bodowinter.com/tutorial/bw_LME_tutorial2.pdf.

As an aside, I have a feeling that RT might need to be transformed before being included in your model (e.g., via a log transformation).

  • $\begingroup$ Thats a great resource, thanks! One confusion remains...the random intercept model makes sense because we expect individuals to have different baseline RT, but the inclusion of difficulty in the first model is confusing. Thats testing whether the slope between difficulty and RT is different between individuals (I think). Why difficulty? If my main question is about activity and RT, wouldn't it make more sense to use (1+activity | ID) instead? Or is there something about the structure of my data that doesn't require the inclusion of activity in the random effects portion of the model? $\endgroup$
    – Simon
    Mar 8, 2019 at 5:34
  • $\begingroup$ also, when I try to run the first model on my data I get the following error and not sure how to interpret it: number of observations (=394) <= number of random effects (=394) for term (1 + difficulty | ID); the random-effects parameters and the residual variance (or scale parameter) are probably unidentifiable $\endgroup$
    – Simon
    Mar 8, 2019 at 5:59
  • $\begingroup$ @Simon: I did some digging and found this link (web.stanford.edu/~rag/stat209/lmer2wave.html) which states that the most recent version of lme4 objects to two-wave data (i.e., data where each individual produces two repeated values of the outcome variable). So you may have no option but to use the lme() function from the nlme package to fit your model or perhaps the mixed_model() from the GLMMadaptive package with the gaussian() distribution (cran.r-project.org/web/packages/GLMMadaptive/vignettes/…). $\endgroup$ Mar 9, 2019 at 17:55
  • $\begingroup$ How you define 'baseline' RT depends on what predictor variables are included in your model and what it means for those predictor variable to take the value zero. Based on your current model, 'baseline' RT for an individual means reaction time when age = 0, activity = 0 and difficulty = 0 (which is uninterpretable, unless you choose to center the continuous predictor variables). $\endgroup$ Mar 9, 2019 at 18:01
  • $\begingroup$ Any predictor measured at the lowest level of your data generating hierarchy can in principle have effects which vary across individuals. Since you measure each individual under two tasks/conditions, any predictor whose values change across the two conditions for an individual can be allowed its own random slope in the model (if you have reason to believe that the effect of that predictor can be heterogeneous across individuals AND have enough data to support its inclusion in the model). Both activity and difficulty seem to fit that description, hence my suggestion. $\endgroup$ Mar 9, 2019 at 18:06

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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