21
$\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$
24
$\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")

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

$\endgroup$
4
  • 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$ 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$ Jun 18 '20 at 21:11
  • $\begingroup$ Thanks CrunchyTopping, I updated the links. The emmeans df can be specified with the emmeans argument, lmer.df = "satterthwaite". $\endgroup$ Nov 17 '20 at 21:50
12
$\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.

$\endgroup$
4
  • $\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$ Aug 6 '17 at 17:54
  • $\begingroup$ 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. $\endgroup$
    – kakarot
    Nov 4 '20 at 3:29
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)

$\endgroup$
2
  • 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$ Jun 8 '20 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 '20 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.