3
$\begingroup$

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:

Design

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?

$\endgroup$

1 Answer 1

3
$\begingroup$

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 
$\endgroup$

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.