Repeated Measures in R Using Linear Mixed‑Effects Models

By /

Analyzing repeated‑measures data is a common task in psychology, education, basic science, health sciences, and many other fields where the same experimental units are measured across multiple time points or conditions. Traditionally, researchers have relied on repeated‑measures ANOVA to analyze these designs. While repeated‑measures ANOVA is familiar and easy to implement, it comes with several restrictive assumptions—most notably sphericity, equal spacing of measurements, and the requirement that all participants provide complete data. Violations of these assumptions are common in real‑world studies and may lead to biased results or loss of statistical power.

Linear mixed‑effects models (LMMs) offer a more flexible and robust alternative. Instead of forcing all experimental units to follow the same trajectory, mixed models allow each experimental unit to have their own random intercepts and, when appropriate, random slopes. This framework naturally accounts for the correlation among repeated observations and does so without requiring sphericity. Mixed models also handle unbalanced data, missing time points, and unequal spacing between measurements—situations where repeated‑measures ANOVA either fails or requires ad‑hoc fixes.

In this post, I will demonstrate how to perform one-way repeated‑measures analysis in R using linear mixed effects models, focusing on practical implementation with the lme4 package.

Repeated Measures using Linear Mixed Models in R

Research scenario: A school psychologist wants to examine whether participating in a school‑based yoga program helps reduce students’ math anxiety over time. Twenty‑five students identified with elevated math anxiety complete a standardized anxiety questionnaire at the start of the program, again at three months, and once more at six months. Because these repeated measurements come from the same students, the observations are correlated and may vary in their individual trajectories. The data for this example can be downloaded as a long format CSV file from here or directly by the R code in Listing 1.

For this example, our research question is: Does students’ math anxiety decrease over time, and can a linear mixed‑effects model capture both the overall trend and the differences across three measurements?

To formalize the linear mixed model, let yij represent the math‑anxiety score for student i at time point j. A basic random‑intercept linear mixed‑effects model can be written as:

yij=β0+β1Timeij+b0i+εij,b0i∼𝒩(0,σb2),εij∼𝒩(0,σ2).\begin{aligned} y_{ij} &= \beta_0 + \beta_1 \,\text{Time}_{ij} + b_{0i} + \varepsilon_{ij}, \\ b_{0i} &\sim \mathcal{N}(0, \sigma_b^2), \\ \varepsilon_{ij} &\sim \mathcal{N}(0, \sigma^2). \end{aligned}

where:

β0 is the population‑level intercept (average baseline anxiety), β1 is the fixed effect of time (average change in anxiety across the sample), b0i∼N(0,σb2) is the random intercept for student i, capturing individual differences in baseline anxiety, and εij∼N(0,σ2) is the residual error term. We run this model using the R code in Listing 1.

Listing 1: R code to run repeated measures using linear mixed modeling.
# Load libraries
library(lme4)
library(lmerTest)   # Gives p-values for fixed effects
library(car) 
library(emmeans) # Posthoc tests

# Download the data
url <- "https://statisticswithr.com/datasets/df_math_anxiety_long.csv"
df_math_anxiety_long <- read.csv(url)

# Build model
anxietyModel <- lmer(anxiety_score ~ time + (1 | student_id), data = df_math_anxiety_long)
summary(anxietyModel)


# Posthoc
emmeans(anxietyModel, pairwise ~ time, adjust = "bonferroni")
Linear Mixed Model Results
# Output 1
# summary(anxietyModel)
# Output 2
Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
Formula: anxiety_score ~ time + (1 | student_id)
   Data: df_math_anxiety_long

REML criterion at convergence: 530.5

Scaled residuals: 
     Min       1Q   Median       3Q      Max 
-2.66254 -0.48781 -0.01677  0.58582  2.21391 

Random effects:
 Groups     Name        Variance Std.Dev.
 student_id (Intercept) 27.10    5.206   
 Residual               61.22    7.824   
Number of obs: 75, groups:  student_id, 25

Fixed effects:
                Estimate Std. Error      df t value Pr(>|t|)    
(Intercept)       75.960      1.880  60.589  40.413  < 2e-16 ***
timeanxiety_t02  -20.520      2.213  48.000  -9.272 2.83e-12 ***
timeanxiety_t03  -30.960      2.213  48.000 -13.990  < 2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation of Fixed Effects:
            (Intr) tmn_02
timnxty_t02 -0.589       
timnxty_t03 -0.589  0.500

# Output 2
# Post hoc
Post hoc Tests Results
# Output 2
# Posthoc test

$emmeans
 time        emmean   SE   df lower.CL upper.CL
 anxiety_t01   76.0 1.88 60.6     72.2     79.7
 anxiety_t02   55.4 1.88 60.6     51.7     59.2
 anxiety_t03   45.0 1.88 60.6     41.2     48.8

Degrees-of-freedom method: kenward-roger 
Confidence level used: 0.95 

$contrasts
 contrast                  estimate   SE df t.ratio p.value
 anxiety_t01 - anxiety_t02     20.5 2.21 48   9.272 <0.0001
 anxiety_t01 - anxiety_t03     31.0 2.21 48  13.990 <0.0001
 anxiety_t02 - anxiety_t03     10.4 2.21 48   4.718 <0.0001

Degrees-of-freedom method: kenward-roger 
P value adjustment: bonferroni method for 3 tests

The linear mixed‑effects model revealed a strong effect of time on students’ math‑anxiety scores. At baseline, students reported high anxiety levels (intercept = 75.96, SE = 1.88, t(60.59) = 40.41, p < .001). Anxiety decreased substantially at each subsequent measurement occasion, with scores dropping by 20.52 points at three months (SE = 2.21, t(48) = –9.27, p = 2.83 × 10⁻¹²) and by 30.96 points at six months (SE = 2.21, t(48) = –13.99, p < .001). The random‑intercept variance (27.10) indicated meaningful individual differences in baseline anxiety, while the residual variance (61.22) reflected within‑student variability across time. Overall, the model provides strong evidence that math anxiety declined over the six‑month period.

In addition, estimated marginal means showed a clear downward trend across the three time points: 76.0 at baseline, 55.4 at three months, and 45.0 at six months (all SEs = 1.88). Bonferroni‑adjusted pairwise comparisons confirmed that each time point differed significantly from the others. Anxiety decreased by 20.5 points from baseline to three months (t(48) = 9.27, p < .0001), by 31.0 points from baseline to six months (t(48) = 13.99, p < .0001), and by 10.4 points from three to six months (t(48) = 4.72, p < .0001). These post‑hoc results reinforce the conclusion that students experienced steady and meaningful reductions in math anxiety over time.

References

Fox, J. (2016). Applied regression analysis and generalized linear models (3rd ed.). SAGE Publications.

Hothorn, T., & Everitt, B. S. (2014). A handbook of statistical analyses using R (3rd ed.). Chapman and Hall/CRC. https://doi.org/10.1201/b17081

Quinn, G. P., & Keough, M. J. (2023). Experimental Design and Data Analysis for Biologists (2nd ed.). Cambridge: Cambridge University Press.

Wang, Z., & Goonewardene, L. A. (2004). The use of mixed models in the analysis of animal experiments with repeated measures data. Canadian Journal of Animal Science, 84(1), 1–11. https://doi.org/10.4141/A03-123

Scroll to Top