# Statistical test for two factors (each multiple levels) and repeated measures

I wondered what statistical test to use in the following case:

• one categorical factor/condition with three levels
• one ordinal factor/condition with three levels
• one numeric response variable
• a small population of individuals (N = 20)
• the individuals vary strongly with respect to the response variable
• each individual was sampled under all possible combinations of the conditions
• there could be interaction between the factors

Thus, the design looks as follows:

Now, I want to test

1. whether there are differences between treatments A, B, and C
2. whether there are trends in factor 2

I thought about multiple paired t-tests, but perhaps there is a better solution?

Here is one way you can do this, which will allow you to make comparisons between any of the factor levels from one fitted model.

First I'll simulate some data that matches your design:

library(dplyr)
library(lme4)
library(emmeans)

set.seed(1)

d <- expand.grid(
id = as.factor(1:20),
f1 = as.factor(c(0,1,2)),
f2 = as.factor(c(0,1,2))
)

design <- model.matrix(~ d$$f1*d$$f2 + d$id) coefs <- c(c(0,0.25,0.35,0.45,0.55), rnorm(19, 0, 0.5), c(0.1, 0.25, 0.30, 0.4)) d$y <- as.numeric(design %*% coefs + rnorm(180, 0, 1))

d <- as_tibble(d)


The data is in long format, i.e., there are 9 rows per individual (one for each combination of factor1 and factor2).

Here are the first 10 rows:

# A tibble: 180 × 4
id    f1    f2         y
<fct> <fct> <fct>  <dbl>
1 1     0     0      0.594
2 1     0     1      1.37
3 1     0     2      1.33
4 1     1     0      0.325
5 1     1     1     -1.19
6 1     1     2      1.72
7 1     2     0      0.294
8 1     2     1      0.894
9 1     2     2     -0.171
10 2     0     0     -0.791


Next fit a mixed effects model with a random intercept for individual:

m <- lmer(y ~ f1*f2 + (1 | id), data = d)


You can then estimate the mean value of the outcome for each combination of the factors:

m <- lmer(y ~ f1*f2 + (1 | id), data = d)

means <- emmeans(m, ~ c(f1, f2))



Which gives:

> means
f1 f2  emmean   SE  df lower.CL upper.CL
0  0  -0.1288 0.25 162  -0.6220    0.364
1  0   0.2108 0.25 162  -0.2824    0.704
2  0   0.1238 0.25 162  -0.3693    0.617
0  1   0.0759 0.25 162  -0.4173    0.569
1  1   0.8786 0.25 162   0.3854    1.372
2  1   1.0669 0.25 162   0.5737    1.560
0  2   0.4245 0.25 162  -0.0687    0.918
1  2   0.8619 0.25 162   0.3688    1.355
2  2   1.4285 0.25 162   0.9353    1.922


Then, you can make any kind of contrast you like. E.g., the following compares every combination of factors against all other combinations (showing only the first few lines):

pairs(means)

contrast          estimate    SE  df t.ratio p.value
f10 f20 - f11 f20  -0.3396 0.338 152  -1.003  0.9851
f10 f20 - f12 f20  -0.2526 0.338 152  -0.746  0.9980
f10 f20 - f10 f21  -0.2047 0.338 152  -0.605  0.9996
f10 f20 - f11 f21  -1.0074 0.338 152  -2.977  0.0795


Alternatively, you can compare the levels of factor1, averaging over the levels of factor2 (which is inadvisable when interactions are present):

> emmeans(m, ~f1) |> pairs()
NOTE: Results may be misleading due to involvement in interactions
contrast  estimate    SE  df t.ratio p.value
f10 - f11   -0.527 0.195 152  -2.695  0.0213
f10 - f12   -0.749 0.195 152  -3.834  0.0005
f11 - f12   -0.223 0.195 152  -1.139  0.4915

Results are averaged over the levels of: f2