22

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?

Karolis Koncevičius
  • 4,282
  • 7
  • 30
  • 47
Ping Tang
  • 425
  • 1
  • 5
  • 7

3 Answers3

24

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")

The next option is difflsmeans. 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 default by emmeans, so it might be best to explicitly specify the method you prefer.

library(lmerTest)
difflsmeans(model, test.effs = "Group", ddf="Kenward-Roger")

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)
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."

Kayle Sawyer
  • 666
  • 6
  • 11
  • 1
    I just want to add to the response of Kayle Sawyer that the package *lsmeans* is being deprecated in favor of *emmeans*. – Downhiller Jun 14 '18 at 19:52
  • Note if you specify the library, you must use lmerTest::lmer(), not lme4::lmer() for anova() to show the p-values. – Kayle Sawyer Jun 28 '18 at 19:20
  • 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 ` – CrunchyTopping Jun 18 '20 at 21:11
  • Thanks CrunchyTopping, I updated the links. The emmeans df can be specified with the `emmeans` argument, `lmer.df = "satterthwaite"`. – Kayle Sawyer Nov 17 '20 at 21:50
12

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.

Hack-R
  • 848
  • 9
  • 24
  • Thank you. But which adjustment was actually used? Tukey or holm? Why both appears in the post-hoc test? – Ping Tang Sep 29 '16 at 02:41
  • @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. – Hack-R Sep 29 '16 at 02:54
  • 4
    @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. – Sal Mangiafico Aug 06 '17 at 17:54
  • 2
    Note that for `lmer()` models, the default pvalues from `glht()` and `emmeans()` will be different. This is because `emmeans()` uses the K-R estimate of degrees of freedom, while `glht()` defaults to a normal approximation (z-score). See `?glht`. You can enter the degrees of freedom directly via `glht(...,df=28)` to make the p-values agree. – kakarot Nov 04 '20 at 03:29
0

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)

Andy
  • 116
  • 1
  • 1
    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. – gung - Reinstate Monica Jun 08 '20 at 16:17
  • 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 – Andy Jun 09 '20 at 09:36