This is a repost from StackOverflow, I was advised to post here instead, thank you in advance!

I've seen many lmers , but not so many of them of longitudinal data (which is my case right now). The thing is, I have scores (outcome) collected at 2021 and at 2022 (predictor: Year) and in two languages (Spanish and English, predictor: Language). Therefore, my data look like this:


SUBJ01  200   2021  SPANISH
SUBJ01  250   2022  ENGLISH
SUBJ01  230   2021  SPANISH
SUBJ01  300   2022  ENGLISH

So each subject has 4 scores, 2 per year (2021/2022) and 2 per language (Span/Eng). I'm fitting a lmer but I'm wondering if I should include random slopes as well or stay with random intercepts.

So, basically, I'm thinking about:

lmer(SCORE ~ YEAR * LANGUAGE + (1|SUBJECT), data = data)



I'm wondering if I should use the random slopes as well or not. Any thoughts would be much appreciated, thanks in advance :)


I've ran two more models,

modYearRandomSlope <- lmer(SCORE ~ LANGUAGE * YEAR + (YEAR|SUBJECT), data = data)


modLangRandomSlope <- lmer(SCORE ~ LANGUAGE * YEAR + (LANGUAGE|SUBJECT), data = data)
error: boundary (singular) fit: see help('isSingular')

why do I get an error fitting random slopes for language but not for year? Ps: If I only fit (YEAR|ID) as random slopes, do I still account for the variability within language as well?

I've plotted these two graphs:



(ps: take 2020 as 2021 and 2021 as 2022) and



If I run an anova between these models and a model with just the subject random intercept, it selects the model with no random slopes, any guesses why?

ModOnlyRandomIntercept <- lmer(SCORE ~ LANGUAGE * YEAR + (1|ID)

enter image description here

I'm trying to model the best fit to account for all these within subjects variability, maybe it's clear now, thanks in advance

  • 1
    $\begingroup$ This isn't really longitudinal data as there are only two years: time is either time1 or time2. Same with languages. So neither YEAR nor LANGUAGE make a good candidate to be treated as "random effects". I wouldn't consider the second model. $\endgroup$
    – dipetkov
    Commented Aug 25, 2022 at 7:00
  • $\begingroup$ The first model has fixed YEAR and LANGUAGE effects with an interaction and random intercepts for the students. You can think about adding random slopes as well. For example (YEAR|SUBJECT) adds a random slope for the effect of time, so that the performance of each student can change from one year to the next. $\endgroup$
    – dipetkov
    Commented Aug 25, 2022 at 7:05
  • $\begingroup$ Do you have only one data point per (YEAR,LANGUAGE) pair for each student? $\endgroup$
    – dipetkov
    Commented Aug 25, 2022 at 7:06
  • $\begingroup$ @dipetkov , each student has 4 data points: 2 scores for LANGUAGE for each YEAR (i.e, 1 score for English (2021), 1 score for Spanish (2021), 1 score for English (2022) and 1 score for Spanish (2022). Is it clearer now? I really appreaciate your help, thank you $\endgroup$ Commented Aug 25, 2022 at 11:10
  • $\begingroup$ @dipetkov , if you could, would you take a look at my edit, please? $\endgroup$ Commented Aug 25, 2022 at 13:39

1 Answer 1


You want to analyze test score data with a mixed effects model. Each student was tested four times: in two languages (Spanish and English) and in two consecutive years (2020 and 2021). Since there are only two time points, the most complex form of the time effect is a line. Conceptually, it is easier to think of the problem not as longitudinal but in terms of the four unique combinations of two two-level factors.

Mathematically, you can let both (or either) the language effect and the year effect vary by student.

I introduce notation to make this more precise. Let $E(Y_{jk})$ be the mean score for language $j$ in year $k$. In the population (of students), the mean scores are modeled by the fixed effects (an intercept, two main effects and an interaction) and represented by SCORE ~ YEAR * LANGUAGE in the model formula.

For student $i$, let $Y_{j[i],k[i]}$ denote their score for language $j$ in year $k$. Here are the random components that seem appropriate for your test score data. The random effects (intercept, slopes) are differences from the population mean to account for additional variability from student to student.

  • (1 | SUBJECT): random intercept $Y_{j[i],k[i]} = E(Y_{jk}) + a_i + \text{error}$
  • (YEAR | SUBJECT): random intercept and random year slope $Y_{j[i],k[i]} = E(Y_{jk}) + a_i + b_j\left(\text{YEAR}\right) + \text{error}$
  • (LANGUAGE | SUBJECT): random intercept and random language slope $Y_{j[i],k[i]} = E(Y_{jk}) + a_i + b_k\left(\text{LANGUAGE}\right) + \text{error}$
  • (YEAR + LANGUAGE | SUBJECT): random intercept, with random slopes for both year and language $Y_{j[i],k[i]} = E(Y_{jk}) + a_i + b_j\left(\text{YEAR}\right) + b_k\left(\text{LANGUAGE}\right) + \text{error}$

Aside: You don't have enough observations per student to include a random interaction with (YEAR * LANGUAGE | SUBJECT).

So there are four options for the random component in the mixed effects model. How do you choose?

If you can fit all four models successfully — in particular, without getting a boundary (singular) fit error — you can use AIC or the likelihood ratio test to select the best model. Based on the anova() table you report at the end, you will choose the random intercept (simplest) model. The danger of this approach is that you pick the "best" model for the sample, not the "best" model for the population. That is, you overfit to the data.

Alternatively, you can skip model selection altogether and choose the "maximal" model, the one with two random slopes, as this is the model you believed a priori to be appropriate. (I infer this from the fact that you wrote this question). These random effects have meaningful and intuitive interpretation: the student-specific year slope means that the improvement from one year to the next varies by student; similarly, the student-specific language slope means that the language aptitude varies by student as well.

I prefer the second approach. Unfortunately, it may not be possible to put in practice because each student was tested only four times. It's challenging to estimate how scores vary by student when there are two observations per year/language.

You ask specifically how to formulate a mixed effects model. There are alternative approaches to analyze structured data like yours (structured = observations are not independent but correlated within clusters such a student taking a language test multiple times). For example, generalized least squares (GLS) might work nicely.

  • $\begingroup$ thank you very much for taking your time to help me, I REALLY appreciate it. I was thinking, in order to account for overfitting, maybe I could run the complete model (random intercept + 2 random slopes), which will give me a boundary error for lang and then stick with the model with the random intercept + random slope for year ? I could argue that the best model would be the complete one, but since there seems to be little variation within language, I decided to stick with the one which still accounts (even partially) to overfitting rather than the one with just the intercept? $\endgroup$ Commented Aug 26, 2022 at 11:55
  • $\begingroup$ Sounds good. I think to some extent you are overthinking it. I tried some simulations last night (though I wasn't sure how to generate data that looks like yours) and when there are four observations per student, the standard errors on the fixed effects are almost the same in all four models. This means that the inference (statistical conclusions) about the fixed effects would be about the same. $\endgroup$
    – dipetkov
    Commented Aug 26, 2022 at 11:57
  • $\begingroup$ I was wondering another thing (if it's better, I may open another question for that). Why I don't get an error for this: aov(SCORE ~ LANGUAGE * YEAR + Error(SUBJECT/LANGUAGE * YEAR), data = data), but for the model (lmer(SCORE ~ (YEAR * LANGUAGE | SUBJECT), data = data)) I do get one (i.e, the boundary issue). $\endgroup$ Commented Aug 26, 2022 at 12:04
  • $\begingroup$ Yes, better ask a new question. I never use ANOVA myself. $\endgroup$
    – dipetkov
    Commented Aug 26, 2022 at 12:05
  • $\begingroup$ All right, I'll do that. Many thanks once more!!! 😁😁😁 $\endgroup$ Commented Aug 26, 2022 at 12: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.