# How do I compare one variable across three groups and explore the impact of two fixed effects, and perform multiple comparisons between interactions?

I have a dataset in the form below:

ID  Age_Group   Gender  Class   Measurement
1   Teen        M       X       219
1   Teen        M       X       222
1   Teen        M       X       230
1   Teen        M       X       291
5   Teen        F       Z       219
5   Teen        F       Z       225
5   Teen        F       Z       211
5   Teen        F       Z       230
6   Teen        M       Y       199
6   Teen        M       Y       198
6   Teen        M       Y       200
6   Teen        M       Y       200


These are evidently repeated measurements on few subjects (ID).

My interest is the change in Measurement when looking at Age_Group between and within the Class categories.

In other words, I want to know if there are differences between:

• X Teen vs. X Adult
• Y Teen vs. Y Adult
• Z Teen vs. Z Adult

And then I also want to know if there are differences between:

• X Teen vs. Y Teen vs. Z Teen

And furthermore:

• X vs. Y. vs. Z (taking data from both Teen and Adult groups into consideration)

And lastly, in the above between-Class and within-Class comparisons, I also want to know if Gender plays a role in differences in measurements.

I have tried this in R:

test_data <- structure(list(ID = c(1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 3L, 3L, 3L, 3L, 4L, 4L, 4L, 4L, 5L, 5L, 5L, 5L, 6L, 6L, 6L, 6L), Age_Group = structure(c(2L, 2L, 2L, 2L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L), .Label = c("Adult", "Teen"), class = "factor"), Gender = structure(c(2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L ), .Label = c("F", "M"), class = "factor"), Class = structure(c(1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 2L, 2L, 2L, 2L), .Label = c("X", "Y", "Z" ), class = "factor"), Measurement = c(219L, 222L, 230L, 291L, 210L, 211L, 215L, 220L, 198L, 190L, 194L, 186L, 200L, 250L, 211L, 230L, 219L, 225L, 211L, 230L, 199L, 198L, 200L, 200L )), .Names = c("ID", "Age_Group", "Gender", "Class", "Measurement" ), class = "data.frame", row.names = c(NA, -24L))
test_data\$interactions <- with(test_data, interaction(Age_Group, Class, sep = " x "))
m1c <- lme4::lmer(Measurement ~ interactions + (1|ID), data = test_data, REML = FALSE)
ll <- multcomp::glht(m1c, linfct = mcp(interactions = "Tukey"))
summary(ll)


This approach performs more tests than necessary. Therefore, I am not sure if my p-values are biased in this way:

     Simultaneous Tests for General Linear Hypotheses

Multiple Comparisons of Means: Tukey Contrasts

Fit: lme4::lmer(formula = Measurement ~ interactions + (1 | ID), data = test_data,
REML = FALSE)

Linear Hypotheses:
Estimate Std. Error z value Pr(>|z|)
Teen x X - Adult x X == 0     26.50      10.47   2.531  0.11495
Adult x Y - Adult x X == 0   -22.00      10.47  -2.101  0.28642
Teen x Y - Adult x X == 0    -14.75      10.47  -1.409  0.72174
Adult x Z - Adult x X == 0     8.75      10.47   0.836  0.96091
Teen x Z - Adult x X == 0      7.25      10.47   0.692  0.98284
Adult x Y - Teen x X == 0    -48.50      10.47  -4.633  < 0.001 ***
Teen x Y - Teen x X == 0     -41.25      10.47  -3.940  0.00117 **
Adult x Z - Teen x X == 0    -17.75      10.47  -1.695  0.53478
Teen x Z - Teen x X == 0     -19.25      10.47  -1.839  0.44073
Teen x Y - Adult x Y == 0      7.25      10.47   0.692  0.98284
Adult x Z - Adult x Y == 0    30.75      10.47   2.937  0.03878 *
Teen x Z - Adult x Y == 0     29.25      10.47   2.794  0.05840 .
Adult x Z - Teen x Y == 0     23.50      10.47   2.245  0.21752
Teen x Z - Teen x Y == 0      22.00      10.47   2.101  0.28645
Teen x Z - Adult x Z == 0     -1.50      10.47  -0.143  0.99999
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
(Adjusted p values reported -- single-step method)


QUESTIONS

1. What is the correct statistical approach to explore this problem?

2. And additionally, is there a way to cover all the tests/comparisons I have mentioned in R efficiently?

It must be noted that in my original dataset I have 32 subjects in total. So this leaves me 5 subjects per each group of interactions of Age_Group and Class, and finally about 2-3 subjects per each group of interactions of Age_Group, Class, and Gender.

The multcomp package is difficult to use when there is more than one primary factor. I suggest using the emmeans package instead. First, re-fit the model with the factors of interest in there explicitly:

mod <- lme4::lmer(Measurement ~ Age_Group * Class + (1|ID), data = test_data,
REML = FALSE)


(why not REML = TRUE, though?)

Then obtain the estimated marginal means for each factor combination and do simple comparisons:

library("emmeans")
( emm <- emmeans(mod, ~ Age_Group * Class) )
pairs(emm, simple = "each")


This will give you pairwise comparisons of each factor at each level of the other factor.

Is there a reason you don't include Gender in the model?

• Hmmm. I tried this with your example, and there is some problem with the pbkrtest package in computing the adjusted covariances and the degrees of freedom. Suggest adding mode = "satter" to the emmeans() call to use the Satterthwaite d.f. instead. Apr 20, 2018 at 15:12