fitting a lmer with longitudinal data for 2 categorical predictors

Question

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:

SUBJECT SCORE YEAR LANGUAGE

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)

or

lmer(SCORE ~ YEAR * LANGUAGE + (1|SUBJECT)+ (1|YEAR:SUBJECT) + (1|LANGUAGE:SUBJECT)

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

EDIT:

I've ran two more models,

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

and

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:

SCORE against YEAR by SUBJECT

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

SCORE against LANGUAGE by SUBJECT

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)

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

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.