18
$\begingroup$

This is my data frame:

Group   <- c("G1","G1","G1","G1","G1","G1","G1","G1","G1","G1","G1","G1","G1","G1","G1","G2","G2","G2","G2","G2","G2","G2","G2","G2","G2","G2","G2","G2","G2","G2","G3","G3","G3","G3","G3","G3","G3","G3","G3","G3","G3","G3","G3","G3","G3")
Subject <- c("S1","S2","S3","S4","S5","S6","S7","S8","S9","S10","S11","S12","S13","S14","S15","S1","S2","S3","S4","S5","S6","S7","S8","S9","S10","S11","S12","S13","S14","S15","S1","S2","S3","S4","S5","S6","S7","S8","S9","S10","S11","S12","S13","S14","S15")
Value   <- c(9.832217741,13.62390117,13.19671612,14.68552076,9.26683366,11.67886655,14.65083473,12.20969772,11.58494621,13.58474896,12.49053635,10.28208078,12.21945867,12.58276212,15.42648969,9.466436017,11.46582655,10.78725485,10.66159358,10.86701127,12.97863424,12.85276916,8.672953949,10.44587257,13.62135205,13.64038394,12.45778874,8.655142642,10.65925259,13.18336949,11.96595556,13.5552118,11.8337142,14.01763101,11.37502161,14.14801305,13.21640866,9.141392359,11.65848845,14.20350364,14.1829714,11.26202565,11.98431285,13.77216009,11.57303893)

data <- data.frame(Group, Subject, Value)

Then I run a linear-mixed effects model to compare the 3 Groups' difference on "Value", where "Subject" is the random factor:

library(lme4)
library(lmerTest)
model <- lmer (Value~Group + (1|Subject), data = data)
summary(model)

The results are:

Fixed effects:
            Estimate Std. Error       df t value Pr(>|t|)    
(Intercept) 12.48771    0.42892 31.54000  29.114   <2e-16 ***
GroupG2     -1.12666    0.46702 28.00000  -2.412   0.0226 *  
GroupG3      0.03828    0.46702 28.00000   0.082   0.9353    

However, how to compare Group2 with Group3? What is the convention in academic article?

$\endgroup$
20
$\begingroup$

You could use emmeans::emmeans() or lmerTest::difflsmeans(), or multcomp::glht().

I prefer emmeans (previously lsmeans).

library(emmeans)
emmeans(model, list(pairwise ~ Group), adjust = "tukey")

Note difflsmeans cannot correct for multiple comparisons, and uses the Satterthwaite method for calculating degrees of freedom as default instead of the Kenward-Roger method used by emmeans.

library(lmerTest)
difflsmeans(model, test.effs = "Group")

The multcomp::glht() method is described in the other answer to this question, by Hack-R.

Also, you can get the ANOVA p-values by loading lmerTest and then using anova.

library(lmerTest)
anova(model)

Just to be clear, you intended for the Value to be assessed three times for each subject, right? It looks like Group is within-subjects, not between-subjects.

| cite | improve this answer | |
$\endgroup$
  • 1
    $\begingroup$ I just want to add to the response of Kayle Sawyer that the package lsmeans is being deprecated in favor of emmeans. $\endgroup$ – Downhiller Jun 14 '18 at 19:52
  • $\begingroup$ Note if you specify the library, you must use lmerTest::lmer(), not lme4::lmer() for anova() to show the p-values. $\endgroup$ – Kayle Sawyer Jun 28 '18 at 19:20
  • $\begingroup$ Link and uses the Satterthwaite ... instead of the Kenward-Roger has rotted; source for difflsmeans just says "based on the t-distribution using degrees of freedom based onSatterthwaites or Kenward-Roger methods" so it looks like you can specify one or the other but I can't find a similar for a ddf argument in ?emmeans $\endgroup$ – CrunchyTopping Jun 18 at 21:11
11
$\begingroup$

After you've fit your lmer model you can do ANOVA, MANOVA, and multiple comparison procedures on the model object, like this:

library(multcomp)
summary(glht(model, linfct = mcp(Group = "Tukey")), test = adjusted("holm"))
   Simultaneous Tests for General Linear Hypotheses

Multiple Comparisons of Means: Tukey Contrasts


Fit: lmer(formula = Value ~ Group + (1 | Subject), data = data)

Linear Hypotheses:
             Estimate Std. Error z value Pr(>|z|)  
G2 - G1 == 0 -1.12666    0.46702  -2.412   0.0378 *
G3 - G1 == 0  0.03828    0.46702   0.082   0.9347  
G3 - G2 == 0  1.16495    0.46702   2.494   0.0378 *
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
(Adjusted p values reported -- holm method)

As for the convention in academic papers, that's going to vary a lot by field, journal, and specific subject matter. So for that case just review related articles and see what they do.

| cite | improve this answer | |
$\endgroup$
  • $\begingroup$ Thank you. But which adjustment was actually used? Tukey or holm? Why both appears in the post-hoc test? $\endgroup$ – Ping Tang Sep 29 '16 at 2:41
  • $\begingroup$ @PingTang You're welcome. It's Bonferroni-Holm correction of all-pair multiple comparison. That's just one option, of course. You could also do summary(glht(model, linfct = mcp(Group = "Tukey"))). If you want to see the full academic / statistical descriptions of the various tests that can be performed check out the references in ?glht and multicomp more generally. I think Hsu 1996 would be the main one. $\endgroup$ – Hack-R Sep 29 '16 at 2:54
  • 4
    $\begingroup$ @PingTang , the mcp function, the Group = Tukey just means to compare all pairwise groups in the variable "Group". It doesn't mean a Tukey adjustment. $\endgroup$ – Sal Mangiafico Aug 6 '17 at 17:54
0
$\begingroup$

Why not just do a pairwise t.test, with either holm or bonferroni correction, between your groups, using the fitted values from the model, since you see that your group2 varies significantly in your linear model? You could then draw a comparison between all 3 groups from your data.

In which case, you could just write:

PT <- pairwise.t.test(fitted.Values,Group, method=bonferroni)

| cite | improve this answer | |
$\endgroup$
  • 1
    $\begingroup$ Welcome to the site. Was this intended as an answer to the OP's question, a comment requesting clarification from the OP or one of the answerers, or a new question of your own? Please only use the "Your Answer" field to provide answers to the original question. You will be able to comment anywhere when your reputation is >50. If you have a new question, click the gray ASK QUESTION at the top of the page & ask it there, then we can help you properly. Since you're new here, you may want to take our tour, which has information for new users. $\endgroup$ – gung - Reinstate Monica Jun 8 at 16:17
  • $\begingroup$ Hi, it was provided as a possible solution for a post-hoc test to this lmer example. However, it could also be interpreted as a question, since statistics is an on going discussion, and it's possible that a better solution exists than the one I proposed. Thanks for the clarification of where I should post questions $\endgroup$ – Andy Jun 9 at 9:36

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.