1
$\begingroup$

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
2   Adult       M       X       210
2   Adult       M       X       211
2   Adult       M       X       215
2   Adult       M       X       220
3   Adult       F       Y       198
3   Adult       F       Y       190
3   Adult       F       Y       194
3   Adult       F       Y       186
4   Adult       F       Z       200
4   Adult       F       Z       250
4   Adult       F       Z       211
4   Adult       F       Z       230
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
  • X Adult vs. Y Adult vs. Z Adult

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.

$\endgroup$

1 Answer 1

2
$\begingroup$

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?

$\endgroup$
1
  • 2
    $\begingroup$ 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. $\endgroup$
    – Russ Lenth
    Apr 20, 2018 at 15:12

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.